医療アクセスと健康アウトカム：
操作変数法によるパネル分析

import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')

import statsmodels.api as sm
from linearmodels.iv import IV2SLS
from scipy import stats

plt.rcParams['font.family'] = 'Hiragino Sans'
plt.rcParams['axes.unicode_minus'] = False
plt.rcParams['figure.dpi'] = 150

import os
FIG_DIR = 'html/figures'
DATA_B  = 'data/raw/SSDSE-B-2026.csv'
os.makedirs(FIG_DIR, exist_ok=True)

df_b = pd.read_csv(DATA_B, encoding='cp932', header=1)
df_b = df_b[df_b['地域コード'].str.match(r'^R\d{5}$', na=False)].copy()
df_b['年度'] = df_b['年度'].astype(int)
df_b = df_b.sort_values(['都道府県', '年度']).reset_index(drop=True)

# 変数生成
df_b['保健医療費'] = df_b['保健医療費（二人以上の世帯）']
df_b['医療費_log'] = np.log(df_b['保健医療費'].clip(lower=1))
df_b['診療所数_log'] = np.log(df_b['一般診療所数'].clip(lower=1))
df_b['降水量_log'] = np.log(df_b['降水量（年間）'].clip(lower=1))
df_b['高齢化率'] = df_b['65歳以上人口'] / df_b['総人口'] * 100
df_b['消費支出_log'] = np.log(df_b['消費支出（二人以上の世帯）'].clip(lower=1))
df_b['人口_log'] = np.log(df_b['総人口'].clip(lower=1))

region_map = {
    '北海道': '北海道・東北', '青森県': '北海道・東北', '岩手県': '北海道・東北',
    '宮城県': '北海道・東北', '秋田県': '北海道・東北', '山形県': '北海道・東北',
    '福島県': '北海道・東北',
    '茨城県': '関東', '栃木県': '関東', '群馬県': '関東', '埼玉県': '関東',
    '千葉県': '関東', '東京都': '関東', '神奈川県': '関東',
    '新潟県': '中部', '富山県': '中部', '石川県': '中部', '福井県': '中部',
    '山梨県': '中部', '長野県': '中部', '岐阜県': '中部', '静岡県': '中部', '愛知県': '中部',
    '三重県': '近畿', '滋賀県': '近畿', '京都府': '近畿', '大阪府': '近畿',
    '兵庫県': '近畿', '奈良県': '近畿', '和歌山県': '近畿',
    '鳥取県': '中国・四国', '島根県': '中国・四国', '岡山県': '中国・四国',
    '広島県': '中国・四国', '山口県': '中国・四国', '徳島県': '中国・四国',
    '香川県': '中国・四国', '愛媛県': '中国・四国', '高知県': '中国・四国',
    '福岡県': '九州・沖縄', '佐賀県': '九州・沖縄', '長崎県': '九州・沖縄',
    '熊本県': '九州・沖縄', '大分県': '九州・沖縄', '宮崎県': '九州・沖縄',
    '鹿児島県': '九州・沖縄', '沖縄県': '九州・沖縄',
}
df_b['地域'] = df_b['都道府県'].map(region_map)
region_colors = {
    '北海道・東北': '#4e9af1', '関東': '#e05c5c', '中部': '#f0a500',
    '近畿': '#5cb85c', '中国・四国': '#9b59b6', '九州・沖縄': '#f39c12',
}

fig, ax = plt.subplots(figsize=(10, 5))
yearly = df_b.groupby(['年度', '地域'])['保健医療費'].mean().reset_index()
for reg, grp in yearly.groupby('地域'):
    ax.plot(grp['年度'], grp['保健医療費'] / 1000, marker='o', markersize=4,
            label=reg, color=region_colors.get(reg, 'gray'))
ax.set_xlabel('年度', fontsize=12)
ax.set_ylabel('保健医療費（千円）', fontsize=12)
ax.set_title('地域別 保健医療費の推移（2012〜2023年）', fontsize=14, fontweight='bold')
ax.legend(fontsize=9)
ax.grid(alpha=0.3)
plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2022_U5_12_fig1_ts.png'), dpi=150, bbox_inches='tight')
plt.close()
print("Fig1 saved")

df_2022 = df_b[df_b['年度'] == 2022].copy()
df_sc = df_2022[['都道府県', '一般診療所数', '保健医療費', '地域']].dropna()
fig, ax = plt.subplots(figsize=(9, 6))
for reg, grp in df_sc.groupby('地域'):
    ax.scatter(grp['一般診療所数'], grp['保健医療費'] / 1000,
               color=region_colors.get(reg, 'gray'), label=reg, s=60, alpha=0.8)
for _, row in df_sc.iterrows():
    ax.annotate(row['都道府県'][:2], (row['一般診療所数'], row['保健医療費'] / 1000),
                fontsize=6, alpha=0.6, xytext=(2, 2), textcoords='offset points')
slope, intercept, r, p, _ = stats.linregress(df_sc['一般診療所数'], df_sc['保健医療費'] / 1000)
xr = np.linspace(df_sc['一般診療所数'].min(), df_sc['一般診療所数'].max(), 100)
ax.plot(xr, slope * xr + intercept, 'k--', linewidth=1.5,
        label=f'回帰直線 (r={r:.2f}, p={p:.3f})')
ax.set_xlabel('一般診療所数（件）', fontsize=12)
ax.set_ylabel('保健医療費（千円）', fontsize=12)
ax.set_title('一般診療所数 vs 保健医療費（2022年, N=47）', fontsize=13, fontweight='bold')
ax.legend(fontsize=8, ncol=2)
ax.grid(alpha=0.3)
plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2022_U5_12_fig2_scatter.png'), dpi=150, bbox_inches='tight')
plt.close()
print("Fig2 saved")

iv_vars = ['高齢化率', '消費支出_log', '人口_log']
df_iv = df_b[['医療費_log', '診療所数_log', '降水量_log'] + iv_vars].dropna()

# OLS
X_ols = sm.add_constant(df_iv[['診療所数_log'] + iv_vars])
res_ols = sm.OLS(df_iv['医療費_log'], X_ols).fit()

# 2SLS using linearmodels IV2SLS
try:
    exog = sm.add_constant(df_iv[iv_vars])
    endog = df_iv[['診療所数_log']]
    instruments = df_iv[['降水量_log']]
    mod_iv = IV2SLS(df_iv['医療費_log'], exog, endog, instruments)
    res_iv = mod_iv.fit(cov_type='unadjusted')

    # Compare 診療所数_log coefficient
    ols_coef = res_ols.params['診療所数_log']
    ols_se = res_ols.bse['診療所数_log']
    iv_coef = res_iv.params['診療所数_log']
    iv_se = res_iv.std_errors['診療所数_log']

    fig, ax = plt.subplots(figsize=(8, 5))
    methods = ['OLS（内生性無視）', '2SLS（IV推定）']
    coefs_c = [ols_coef, iv_coef]
    ses_c = [ols_se, iv_se]
    colors_c = ['#888888', '#e05c5c']
    ax.barh(methods, coefs_c, xerr=[1.96 * s for s in ses_c], color=colors_c, alpha=0.8,
            error_kw={'elinewidth': 2, 'capsize': 6})
    ax.axvline(0, color='black', linewidth=0.8)
    ax.set_xlabel('診療所数（log）の係数', fontsize=12)
    ax.set_title('OLS vs 2SLS：診療所数が保健医療費に与える効果\n（内生性バイアスの補正比較）', fontsize=12, fontweight='bold')
    ax.grid(axis='x', alpha=0.3)
    for i, (c, s) in enumerate(zip(coefs_c, ses_c)):
        ax.text(c + 0.02, i, f'{c:.3f}±{1.96*s:.3f}', va='center', fontsize=9)
    plt.tight_layout()
    plt.savefig(os.path.join(FIG_DIR, '2022_U5_12_fig3_iv_compare.png'), dpi=150, bbox_inches='tight')
    plt.close()
    print(f"Fig3 saved: OLS={ols_coef:.3f}, 2SLS={iv_coef:.3f}")
    print(res_iv.summary)
except Exception as e:
    print(f"2SLS error: {e}")
    # Fallback: just OLS coefficients
    coefs = res_ols.params.drop('const')
    ses = res_ols.bse.drop('const')
    pvals = res_ols.pvalues.drop('const')
    fig, ax = plt.subplots(figsize=(8, 5))
    colors_c = ['#e05c5c' if p < 0.05 else '#888888' for p in pvals]
    ax.barh(range(len(coefs)), coefs, xerr=1.96 * ses, color=colors_c, alpha=0.8,
            error_kw={'elinewidth': 1.5, 'capsize': 4})
    ax.set_yticks(range(len(coefs)))
    ax.set_yticklabels(coefs.index, fontsize=10)
    ax.axvline(0, color='black', linewidth=0.8)
    ax.set_xlabel('OLS係数', fontsize=12)
    ax.set_title('保健医療費の決定要因（OLS）', fontsize=12, fontweight='bold')
    ax.grid(axis='x', alpha=0.3)
    plt.tight_layout()
    plt.savefig(os.path.join(FIG_DIR, '2022_U5_12_fig3_iv_compare.png'), dpi=150, bbox_inches='tight')
    plt.close()

df_iv_sc = df_b[df_b['年度'] == 2022][['都道府県', '降水量（年間）', '一般診療所数', '地域']].dropna()
fig, ax = plt.subplots(figsize=(9, 6))
for reg, grp in df_iv_sc.groupby('地域'):
    ax.scatter(grp['降水量（年間）'], grp['一般診療所数'],
               color=region_colors.get(reg, 'gray'), label=reg, s=60, alpha=0.8)
for _, row in df_iv_sc.iterrows():
    ax.annotate(row['都道府県'][:2], (row['降水量（年間）'], row['一般診療所数']),
                fontsize=6, alpha=0.6, xytext=(2, 2), textcoords='offset points')
x_v = df_iv_sc['降水量（年間）']
y_v = df_iv_sc['一般診療所数']
slope, intercept, r, p, _ = stats.linregress(x_v, y_v)
xr = np.linspace(x_v.min(), x_v.max(), 100)
ax.plot(xr, slope * xr + intercept, 'k--', linewidth=1.5,
        label=f'相関: r={r:.2f} (p={p:.3f})')
ax.set_xlabel('降水量・年間（mm）', fontsize=12)
ax.set_ylabel('一般診療所数（件）', fontsize=12)
ax.set_title('操作変数の妥当性確認：降水量 vs 一般診療所数（2022年）\n（第1段階の関連性）', fontsize=12, fontweight='bold')
ax.legend(fontsize=8, ncol=2)
ax.grid(alpha=0.3)
plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2022_U5_12_fig4_iv_valid.png'), dpi=150, bbox_inches='tight')
plt.close()
print("Fig4 saved")
print("All done!")