"""
教育用再現コード: 2022年 統計データ分析コンペティション 審査員奨励賞 [高校生の部]
==========================================================================
論文タイトル：高齢化と医療費：都道府県別パネルデータによる要因分析
手法：パネルOLS固定効果、相関分析、時系列、散布図

【分析概要】
  データ：SSDSE-B-2026.csv（都道府県×年度パネルデータ）
          47都道府県 × 12年（2012〜2023年）= 564観測
  目的変数：保健医療費（二人以上の世帯）[円/月] — 医療費の家計負担代理
  説明変数：
    ・高齢化率（65歳以上人口 / 総人口 × 100）[%]
    ・一般病院数 [施設]
    ・一般診療所数 [施設]
    ・消費支出_log（log変換）
    ・保育所等数

  Fig1. 保健医療費の時系列推移（地域別平均, 2012〜2023年）
  Fig2. 高齢化率 vs 保健医療費 散布図（2022年, 都道府県ラベル付き, 地域色分け）
  Fig3. パネルOLS固定効果 係数プロット
  Fig4. 一般病院数 vs 保健医療費 散布図（2022年）

【データ出典】
  SSDSE-B-2026.csv: 社会・人口統計体系（都道府県データ）
  統計センターより公表の実データ（2012〜2023年）

【データサイエンス学習ポイント】
  1. 高齢化率の定義と計算（超高齢社会の統計的定義）
  2. パネルデータの構造（個体×時点）と固定効果モデルの直感的説明
  3. 医療費の代理変数としての保健医療費（家計調査と医療費統計の違い）
  4. 人口高齢化の政策課題（統計から見える社会保障の持続可能性）
==========================================================================
"""

# ============================================================
# 【データの準備】実行前に以下のデータファイルを用意してください
#
#   必要ファイル:
#     ・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_H5_6_shorei.py
# ============================================================


import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
from matplotlib.patches import Patch
import warnings
warnings.filterwarnings('ignore')

import statsmodels.api as sm
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)

# ================================================================
# ■ 実データ読み込み（SSDSE-B-2026: 都道府県パネル）
# ================================================================
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)

print(f"パネルデータ形状: {df_b.shape}  "
      f"({df_b['都道府県'].nunique()} 都道府県 × {df_b['年度'].nunique()} 年度)")
print(f"年度範囲: {df_b['年度'].min()} 〜 {df_b['年度'].max()}")

# ================================================================
# ■ 変数生成
# ================================================================
df_b['高齢化率']     = df_b['65歳以上人口'] / df_b['総人口'] * 100
df_b['消費支出_log'] = np.log(df_b['消費支出（二人以上の世帯）'].clip(lower=1))
df_b['保健医療費']   = pd.to_numeric(df_b['保健医療費（二人以上の世帯）'], errors='coerce')
df_b['一般病院数']   = pd.to_numeric(df_b['一般病院数'], errors='coerce')
df_b['一般診療所数'] = pd.to_numeric(df_b['一般診療所数'], errors='coerce')
df_b['保育所等数']   = pd.to_numeric(df_b['保育所等数'], errors='coerce')

# 保健医療費の対数（分析用）
df_b['保健医療費_log'] = np.log(df_b['保健医療費'].clip(lower=1))

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

# ================================================================
# ■ 記述統計
# ================================================================
print("\n" + "=" * 60)
print("■ 記述統計（全期間 2012〜2023年）")
print("=" * 60)
desc_cols = ['高齢化率', '保健医療費', '一般病院数', '一般診療所数', '保育所等数', '消費支出_log']
print(df_b[desc_cols].describe().round(2))

# ================================================================
# ■ 図1: 保健医療費の時系列推移（地域別平均）
# ================================================================
print("\n[図1] 保健医療費 時系列プロット作成中...")

ts_region = df_b.groupby(['年度', '地域'])['保健医療費'].mean().reset_index()
ts_all    = df_b.groupby('年度')['保健医療費'].mean().reset_index()
ts_all.columns = ['年度', '全国平均']

fig1, ax1 = plt.subplots(figsize=(10, 5.5))
for region in region_order:
    rdata = ts_region[ts_region['地域'] == region]
    ax1.plot(rdata['年度'], rdata['保健医療費'],
             'o-', color=region_colors[region], linewidth=1.8,
             markersize=5, label=region, alpha=0.85)

ax1.plot(ts_all['年度'], ts_all['全国平均'],
         'k--', linewidth=2.5, markersize=0, label='全国平均', zorder=10)

ax1.set_xlabel('年度', fontsize=11)
ax1.set_ylabel('保健医療費（円/月、二人以上世帯）', fontsize=11)
ax1.set_title('保健医療費の時系列推移（地域別平均、2012〜2023年）\n'
              'データ：SSDSE-B-2026 実データ（47都道府県）',
              fontsize=12, fontweight='bold')
ax1.legend(fontsize=9, loc='upper left', framealpha=0.9)
ax1.grid(axis='y', alpha=0.35)
ax1.xaxis.set_major_locator(mticker.MultipleLocator(2))
ax1.yaxis.set_major_formatter(mticker.FuncFormatter(lambda x, _: f'{int(x):,}'))

# コロナ注釈
y_range = ax1.get_ylim()
try:
    covid_val = ts_all.loc[ts_all['年度'] == 2020, '全国平均'].values[0]
    ax1.annotate('コロナ禍\n(2020)', xy=(2020, covid_val),
                 xytext=(2021, covid_val - 800),
                 fontsize=8.5, color='#555',
                 arrowprops=dict(arrowstyle='->', color='#555', lw=1))
except Exception:
    pass

plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2022_H5_6_fig1_ts.png'), dpi=150, bbox_inches='tight')
plt.close(fig1)
print("  -> 保存: 2022_H5_6_fig1_ts.png")

# ================================================================
# ■ 図2: 高齢化率 vs 保健医療費 散布図（2022年、都道府県ラベル付き）
# ================================================================
print("\n[図2] 高齢化率 vs 保健医療費 散布図（2022年）作成中...")

df_2022 = df_b[df_b['年度'] == 2022].dropna(subset=['高齢化率', '保健医療費']).copy()

fig2, ax2 = plt.subplots(figsize=(12, 8))

for region in region_order:
    rdf = df_2022[df_2022['地域'] == region]
    ax2.scatter(rdf['高齢化率'], rdf['保健医療費'],
                color=region_colors[region], s=60, alpha=0.85,
                label=region, zorder=4)
    for _, row in rdf.iterrows():
        # 都道府県名を短縮（県・府・都・道を削除）
        short = row['都道府県'].replace('都', '').replace('道', '').replace('府', '').replace('県', '')
        ax2.annotate(short,
                     xy=(row['高齢化率'], row['保健医療費']),
                     xytext=(2, 2), textcoords='offset points',
                     fontsize=6.5, alpha=0.8)

# 回帰直線
x_fit = df_2022['高齢化率'].values
y_fit = df_2022['保健医療費'].values
mask_fit = np.isfinite(x_fit) & np.isfinite(y_fit)
z = np.polyfit(x_fit[mask_fit], y_fit[mask_fit], 1)
xs = np.linspace(x_fit[mask_fit].min(), x_fit[mask_fit].max(), 100)
ax2.plot(xs, np.poly1d(z)(xs), 'k--', linewidth=1.5, alpha=0.7, label='回帰直線', zorder=5)

r, p = stats.pearsonr(x_fit[mask_fit], y_fit[mask_fit])
ax2.set_xlabel('高齢化率（65歳以上人口比率 %）', fontsize=12)
ax2.set_ylabel('保健医療費（円/月、二人以上世帯）', fontsize=12)
ax2.set_title(f'高齢化率 vs 保健医療費（2022年、都道府県別）\n'
              f'相関係数 r = {r:.3f}（p = {p:.4f}）',
              fontsize=13, fontweight='bold')
ax2.legend(fontsize=9, loc='upper left', framealpha=0.9)
ax2.grid(True, alpha=0.3)
ax2.yaxis.set_major_formatter(mticker.FuncFormatter(lambda x, _: f'{int(x):,}'))

plt.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2022_H5_6_fig2_scatter.png'), dpi=150, bbox_inches='tight')
plt.close(fig2)
print(f"  -> 保存: 2022_H5_6_fig2_scatter.png  (r={r:.3f}, p={p:.4f})")

# ================================================================
# ■ パネルOLS 固定効果モデル
# ================================================================
print("\n[パネル分析] 固定効果モデル推定中...")

EXOG_VARS = ['高齢化率', '一般病院数', '一般診療所数', '消費支出_log', '保育所等数']
ENDOG_VAR  = '保健医療費'

df_panel = df_b[['都道府県', '年度', ENDOG_VAR] + EXOG_VARS].dropna().copy()
print(f"  分析対象: {len(df_panel)} 観測 "
      f"({df_panel['都道府県'].nunique()} 都道府県 × {df_panel['年度'].nunique()} 年度)")

fe_success = False
fe_params  = None
fe_bse     = None
fe_pvalues = None
fe_rsquared = None

# PanelOLS 固定効果モデル
try:
    from linearmodels.panel import PanelOLS
    df_pi = df_panel.set_index(['都道府県', '年度'])
    fe_model = PanelOLS(
        df_pi[ENDOG_VAR],
        sm.add_constant(df_pi[EXOG_VARS]),
        entity_effects=True,
        time_effects=False,
    )
    fe_res = fe_model.fit(cov_type='clustered', cluster_entity=True)
    print("\n[固定効果モデル推定結果（クラスター標準誤差）]")
    print(fe_res.summary.tables[1])

    # パラメータ抽出（const 除く）
    params_s   = fe_res.params.drop('const', errors='ignore')
    bse_s      = fe_res.std_errors.drop('const', errors='ignore')
    pvalues_s  = fe_res.pvalues.drop('const', errors='ignore')

    fe_params   = params_s.values
    fe_bse      = bse_s.values
    fe_pvalues  = pvalues_s.values
    fe_var_names = params_s.index.tolist()
    fe_rsquared  = fe_res.rsquared_between if hasattr(fe_res, 'rsquared_between') else float('nan')
    fe_success   = True
    print(f"  固定効果モデル成功（R²_within ≒ {fe_res.rsquared:.4f}）")

except ImportError:
    print("  linearmodels 未インストール -> クロスセクションOLSにフォールバック")

# フォールバック: 2022年クロスセクションOLS
if not fe_success:
    df_ols = df_b[df_b['年度'] == 2022][['保健医療費'] + EXOG_VARS].dropna()
    X_ols  = sm.add_constant(df_ols[EXOG_VARS])
    ols_res = sm.OLS(df_ols['保健医療費'], X_ols).fit()
    print("\n[クロスセクションOLS（2022年）]")
    print(ols_res.summary2())
    fe_params    = ols_res.params[EXOG_VARS].values
    fe_bse       = ols_res.bse[EXOG_VARS].values
    fe_pvalues   = ols_res.pvalues[EXOG_VARS].values
    fe_var_names = EXOG_VARS
    fe_rsquared  = ols_res.rsquared

# ================================================================
# ■ 図3: パネルOLS固定効果 係数プロット
# ================================================================
print("\n[図3] 係数プロット作成中...")

# 係数の色（p値によって色分け）
colors3 = ['#C62828' if p < 0.01 else '#FF8F00' if p < 0.05 else '#9E9E9E'
           for p in fe_pvalues]
y_pos3 = np.arange(len(fe_var_names))

fig3, ax3 = plt.subplots(figsize=(9, max(4.5, len(fe_var_names) * 0.9)))

ax3.barh(y_pos3, fe_params, color=colors3, alpha=0.78,
         edgecolor='white', height=0.6)
ax3.errorbar(fe_params, y_pos3, xerr=1.96 * fe_bse,
             fmt='none', color='#333', capsize=4, linewidth=1.2)
ax3.axvline(0, color='gray', linestyle='--', linewidth=1.2)
ax3.set_yticks(y_pos3)
ax3.set_yticklabels(fe_var_names, fontsize=11)
ax3.set_xlabel('回帰係数（±1.96 SE）', fontsize=11)
model_label = 'パネルOLS 固定効果（クラスター SE）' if fe_success else 'OLS（2022年クロスセクション）'
ax3.set_title(f'保健医療費の要因分析 — {model_label}\n'
              f'R² = {fe_rsquared:.4f}  |  N = {len(df_panel)} 観測',
              fontsize=12, fontweight='bold')
ax3.invert_yaxis()
ax3.grid(axis='x', alpha=0.3)
ax3.legend(handles=[Patch(color='#C62828', alpha=0.78, label='p<0.01'),
                    Patch(color='#FF8F00', alpha=0.78, label='p<0.05'),
                    Patch(color='#9E9E9E', alpha=0.78, label='n.s.')],
           fontsize=9, loc='lower right')

# 係数値と有意性ラベルを棒の外側に表示
for i, (coef, pval) in enumerate(zip(fe_params, fe_pvalues)):
    sig = '***' if pval < 0.001 else '**' if pval < 0.01 else '*' if pval < 0.05 else ''
    offset = abs(fe_bse[i]) * 2.2 if fe_bse[i] > 0 else 5
    ha = 'left' if coef >= 0 else 'right'
    x_label = coef + offset if coef >= 0 else coef - offset
    ax3.text(x_label, i, f'{coef:.1f}{sig}', va='center', ha=ha, fontsize=8.5)

plt.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2022_H5_6_fig3_coef.png'), dpi=150, bbox_inches='tight')
plt.close(fig3)
print("  -> 保存: 2022_H5_6_fig3_coef.png")

# ================================================================
# ■ 図4: 一般病院数 vs 保健医療費 散布図（2022年）
# ================================================================
print("\n[図4] 一般病院数 vs 保健医療費 散布図（2022年）作成中...")

df_2022b = df_b[df_b['年度'] == 2022].dropna(subset=['一般病院数', '保健医療費']).copy()

fig4, ax4 = plt.subplots(figsize=(10, 6.5))

for region in region_order:
    rdf = df_2022b[df_2022b['地域'] == region]
    ax4.scatter(rdf['一般病院数'], rdf['保健医療費'],
                color=region_colors[region], s=65, alpha=0.85,
                label=region, zorder=4)
    for _, row in rdf.iterrows():
        short = row['都道府県'].replace('都', '').replace('道', '').replace('府', '').replace('県', '')
        ax4.annotate(short,
                     xy=(row['一般病院数'], row['保健医療費']),
                     xytext=(2, 2), textcoords='offset points',
                     fontsize=6.5, alpha=0.8)

# 回帰直線
x_h = df_2022b['一般病院数'].values
y_h = df_2022b['保健医療費'].values
mask_h = np.isfinite(x_h) & np.isfinite(y_h)
z_h = np.polyfit(x_h[mask_h], y_h[mask_h], 1)
xs_h = np.linspace(x_h[mask_h].min(), x_h[mask_h].max(), 100)
ax4.plot(xs_h, np.poly1d(z_h)(xs_h), 'k--', linewidth=1.5, alpha=0.7, label='回帰直線')

r_h, p_h = stats.pearsonr(x_h[mask_h], y_h[mask_h])
ax4.set_xlabel('一般病院数（施設）', fontsize=12)
ax4.set_ylabel('保健医療費（円/月、二人以上世帯）', fontsize=12)
ax4.set_title(f'一般病院数 vs 保健医療費（2022年、都道府県別）\n'
              f'相関係数 r = {r_h:.3f}（p = {p_h:.4f}）',
              fontsize=13, fontweight='bold')
ax4.legend(fontsize=9, loc='upper right', framealpha=0.9)
ax4.grid(True, alpha=0.3)
ax4.yaxis.set_major_formatter(mticker.FuncFormatter(lambda x, _: f'{int(x):,}'))

plt.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2022_H5_6_fig4_hospital.png'), dpi=150, bbox_inches='tight')
plt.close(fig4)
print(f"  -> 保存: 2022_H5_6_fig4_hospital.png  (r={r_h:.3f}, p={p_h:.4f})")

# ================================================================
# ■ 完了サマリー
# ================================================================
print("\n" + "=" * 60)
print("■ 全図の生成完了（4枚）")
print("=" * 60)
print("  fig1: 2022_H5_6_fig1_ts.png       — 保健医療費 時系列（地域別平均）")
print("  fig2: 2022_H5_6_fig2_scatter.png  — 高齢化率 vs 保健医療費（2022年）")
print("  fig3: 2022_H5_6_fig3_coef.png     — パネルOLS 固定効果 係数プロット")
print("  fig4: 2022_H5_6_fig4_hospital.png — 一般病院数 vs 保健医療費（2022年）")
print()
print(f"  保存先: {os.path.abspath(FIG_DIR)}")