"""
教育用再現コード: 2022年 統計データ分析コンペティション 審査員奨励賞 [大学生・一般の部]
=================================================================
論文タイトル：医療アクセスと健康アウトカム：操作変数法によるパネル分析
受賞：審査員奨励賞（大学生・一般の部）

【分析概要】
  データ：SSDSE-B-2026.csv（都道府県別パネルデータ, 2012〜2023年度）
  対象：全47都道府県 × 12年（2012〜2023）

  目的変数：保健医療費（二人以上の世帯）→ 健康コスト・医療利用の代理変数
  内生変数：一般診療所数（医療アクセスの指標）
  操作変数：降水量（年間）→ 外出・受診行動への影響（気象条件）

  分析の流れ
  1. 時系列：保健医療費の推移（地域別）
  2. 散布図：医療機関数 vs 保健医療費（OLS偏相関）
  3. 2SLS：第1段階（IV → 診療所数）+ 第2段階（診療所数 → 医療費）
  4. OLS vs 2SLS 係数比較グラフ

【被説明変数】
  保健医療費（二人以上の世帯）(円)

【内生変数（第2段階）】
  一般診療所数

【操作変数（第1段階）】
  降水量（年間）— 気象条件が受診行動・医療施設立地に影響（外生性の仮定）

【その他の制御変数】
  高齢化率、消費支出_log、総人口

【推定手法】
  statsmodels IV2SLS（2段階最小二乗法）+ OLS比較

【データ出典】
  SSDSE-B-2026.csv: 社会・人口統計体系（都道府県データ）
=================================================================
"""

# ============================================================
# 【データの準備】実行前に以下のデータファイルを用意してください
#
#   必要ファイル:
#     ・SSDSE-B-2026.csv
#       → data/raw/SSDSE-B-2026.csv に配置
#
#   ダウンロード先:
#     https://www.nstac.go.jp/use/literacy/ssdse/
#     （SSDSE-B（社会・人口統計体系 都道府県データ） の CSV をダウンロード）
#
#   フォルダ配置（プロジェクトルートからの相対パス）:
#     code/                ← このスクリプトの場所
#     data/raw/            ← CSV ファイルをここに配置
#     html/figures/        ← 図の出力先（自動生成）
#
#   実行方法（ファイルを一切編集せず実行可能）:
#     python3 code/2022_U5_12_shorei.py
# ============================================================


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',
}

# Fig1: 保健医療費の時系列推移（地域別）
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")

# Fig2: 散布図 診療所数 vs 保健医療費
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")

# Fig3: 2SLS vs OLS 係数比較
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()

# Fig4: 操作変数の妥当性（降水量 vs 診療所数）
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!")