高齢化と医療費
都道府県別パネルデータによる要因分析

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))

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))

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")

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})")

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")

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)

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

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)}")