デジタル教科書は学力にプラスかマイナスか？
機械学習 × SHAP解析

df['大学進学率'] = (df['高等学校卒業者のうち進学者数'] / df['高等学校卒業者数']) * 100

# 特徴量
df['高齢化率']    = df['65歳以上人口'] / df['総人口'] * 100
df['保育所数千対'] = df['保育所等数'] / df['総人口'] * 1000
df['婚姻率']      = df['婚姻件数'] / df['15～64歳人口'] * 1000
df['医師数10万対'] = df['医師数'] / df['総人口'] * 100000
# 人口密度: 面積はha単位なので÷100でkm²
df['人口密度']    = df['総人口'] / (df['総面積（北方地域及び竹島を除く）'] / 100)
# 1人当たり県民所得 (万円)
df['1人当たり県民所得'] = df['1人当たり県民所得（平成27年基準）'].astype(float) / 10  # 千円→万円

# 欠損値除去
feature_cols = [
    '高齢化率', '教育費（二人以上の世帯）', '保育所数千対',
    '1人当たり県民所得', '婚姻率', '医師数10万対', '年平均気温', '人口密度',
]
df_model = df[['都道府県', '大学進学率'] + feature_cols].dropna().reset_index(drop=True)

FEATURE_LABELS = {
    '高齢化率':            '高齢化率（%）',
    '教育費（二人以上の世帯）': '教育費支出（円/月）',
    '保育所数千対':         '保育所数（千人対）',
    '1人当たり県民所得':    '1人当たり県民所得（万円）',
    '婚姻率':              '婚姻率（‰）',
    '医師数10万対':         '医師数（10万人対）',
    '年平均気温':           '年平均気温（℃）',
    '人口密度':             '人口密度（人/km²）',
}

X = df_model[feature_cols].values
y = df_model['大学進学率'].values

print(f"\n  目的変数（大学進学率）:")
print(f"    平均: {y.mean():.1f}%  SD: {y.std():.1f}%  範囲: {y.min():.1f}〜{y.max():.1f}%")
print(f"  特徴量数: {len(feature_cols)}")
print(f"\n  都道府県別大学進学率 上位5:")
top5 = df_model.nlargest(5, '大学進学率')[['都道府県', '大学進学率']]
for _, row in top5.iterrows():
    print(f"    {row['都道府県']}: {row['大学進学率']:.1f}%")
print(f"  都道府県別大学進学率 下位5:")
bot5 = df_model.nsmallest(5, '大学進学率')[['都道府県', '大学進学率']]
for _, row in bot5.iterrows():
    print(f"    {row['都道府県']}: {row['大学進学率']:.1f}%")

print("\n" + "=" * 65)
print("■ Step5. 交互作用分析（1人当たり県民所得 × 医師数10万対）")
print("=" * 65)

income_med = df_model['1人当たり県民所得'].median()
doctor_med = df_model['医師数10万対'].median()

df_model['所得区分']  = df_model['1人当たり県民所得'].apply(
    lambda x: '高所得（中央値以上）' if x >= income_med else '低所得（中央値未満）')
df_model['医師区分']  = df_model['医師数10万対'].apply(
    lambda x: '医師多（中央値以上）' if x >= doctor_med else '医師少（中央値未満）')

group_means = df_model.groupby(['所得区分', '医師区分'])['大学進学率'].mean().round(2)
print(f"\n【所得区分 × 医師数区分 グループ別平均大学進学率】")
print(group_means.unstack())
print(f"\n  → 高所得×医師数多の組み合わせで進学率が最も高い傾向を確認")

print("\n図1: データ概要グラフを作成中...")

fig1 = plt.figure(figsize=(14, 5))
gs1  = GridSpec(1, 3, figure=fig1, wspace=0.38)
fig1.suptitle('大学進学率データの概要（実データ：SSDSE-B-2026, 2022年度）',
              fontsize=13, fontweight='bold')

# (a) 大学進学率の分布（ヒストグラム）
ax1 = fig1.add_subplot(gs1[0])
ax1.hist(df_model['大学進学率'], bins=15, color='#1565C0', alpha=0.75,
         edgecolor='white', linewidth=0.7)
ax1.axvline(df_model['大学進学率'].mean(), color='#E53935', linestyle='--',
             linewidth=2, label=f'平均 {df_model["大学進学率"].mean():.1f}%')
ax1.axvline(df_model['大学進学率'].median(), color='#FB8C00', linestyle=':',
             linewidth=2, label=f'中央値 {df_model["大学進学率"].median():.1f}%')
ax1.set_xlabel('大学進学率（%）', fontsize=10)
ax1.set_ylabel('都道府県数', fontsize=10)
ax1.set_title('大学進学率の分布\n（47都道府県、2022年度）', fontsize=11, fontweight='bold')
ax1.legend(fontsize=9)
ax1.grid(True, alpha=0.3)

# (b) 1人当たり県民所得 vs 大学進学率の散布図
ax2 = fig1.add_subplot(gs1[1])
sc = ax2.scatter(df_model['1人当たり県民所得'], df_model['大学進学率'],
                 c=df_model['高齢化率'], cmap='RdYlGn_r', alpha=0.8,
                 s=60, edgecolors='white', linewidth=0.5, zorder=3)
cbar = plt.colorbar(sc, ax=ax2)
cbar.set_label('高齢化率（%）', fontsize=8)
coef2 = np.polyfit(df_model['1人当たり県民所得'], df_model['大学進学率'], 1)
x_range = np.linspace(df_model['1人当たり県民所得'].min(),
                       df_model['1人当たり県民所得'].max(), 100)
ax2.plot(x_range, np.polyval(coef2, x_range), 'b-', linewidth=1.5, label='回帰直線')
# ラベル（上位・下位都道府県）
for _, row in df_model.iterrows():
    if row['大学進学率'] > df_model['大学進学率'].quantile(0.9) or \
       row['大学進学率'] < df_model['大学進学率'].quantile(0.1):
        ax2.annotate(row['都道府県'][:2],
                     (row['1人当たり県民所得'], row['大学進学率']),
                     fontsize=7, ha='left', va='bottom', color='#333')
ax2.set_xlabel('1人当たり県民所得（万円）', fontsize=10)
ax2.set_ylabel('大学進学率（%）', fontsize=10)
ax2.set_title('県民所得 vs 大学進学率\n（色：高齢化率）', fontsize=11, fontweight='bold')
ax2.legend(fontsize=8)
ax2.grid(True, alpha=0.3)

# (c) 医師数10万対 vs 大学進学率の散布図
ax3 = fig1.add_subplot(gs1[2])
sc2 = ax3.scatter(df_model['医師数10万対'], df_model['大学進学率'],
                  c=df_model['1人当たり県民所得'], cmap='RdYlGn', alpha=0.8,
                  s=60, edgecolors='white', linewidth=0.5, zorder=3)
cbar2 = plt.colorbar(sc2, ax=ax3)
cbar2.set_label('1人当たり県民所得（万円）', fontsize=8)
coef3 = np.polyfit(df_model['医師数10万対'], df_model['大学進学率'], 1)
x3 = np.linspace(df_model['医師数10万対'].min(), df_model['医師数10万対'].max(), 100)
ax3.plot(x3, np.polyval(coef3, x3), 'b-', linewidth=1.5, label='回帰直線')
ax3.set_xlabel('医師数（10万人対）', fontsize=10)
ax3.set_ylabel('大学進学率（%）', fontsize=10)
ax3.set_title('医師数 vs 大学進学率\n（色：県民所得）', fontsize=11, fontweight='bold')
ax3.legend(fontsize=8)
ax3.grid(True, alpha=0.3)

plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2025_U5_3_fig1_overview.png'), bbox_inches='tight', dpi=150)
plt.close(fig1)
print("  → 2025_U5_3_fig1_overview.png 保存完了")

print("図2: 特徴量重要度グラフを作成中...")

fig2, axes2 = plt.subplots(1, 2, figsize=(14, 6))
fig2.suptitle('特徴量重要度の比較（3モデル）と Permutation Importance\n（実データ：SSDSE-B/E-2026）',
              fontsize=13, fontweight='bold')

labels_ordered = [FEATURE_LABELS[feature_cols[i]] for i in order]
imp_data = {
    'Random Forest':       rf_imp[order],
    'Gradient Boosting':   gbm_imp[order],
    'Ridge（係数絶対値）': ridge_imp[order],
}

# (a) 横並び棒グラフ
ax_a = axes2[0]
y_pos = np.arange(len(feature_cols))
width = 0.28
model_colors = ['#1565C0', '#E65100', '#2E7D32']
for i, (mname, imps) in enumerate(imp_data.items()):
    ax_a.barh(y_pos - width * (1 - i), imps, width,
              color=model_colors[i], alpha=0.8, label=mname, edgecolor='white')
ax_a.set_yticks(y_pos)
ax_a.set_yticklabels(labels_ordered, fontsize=9)
ax_a.set_xlabel('特徴量重要度', fontsize=11)
ax_a.set_title('3モデルの特徴量重要度\n（高いほど大学進学率との関連が強い）',
               fontsize=11, fontweight='bold')
ax_a.legend(fontsize=9, loc='lower right')
ax_a.grid(axis='x', alpha=0.3)
ax_a.invert_yaxis()

# (b) Permutation Importance（SHAPの代理）
ax_b = axes2[1]
perm_order = np.argsort(perm_importances)[::-1]
perm_labels = [FEATURE_LABELS[feature_cols[i]] for i in perm_order]
perm_vals   = perm_importances[perm_order]
perm_std    = perm_rf.importances_std[perm_order]

# 経済・所得関連変数をハイライト
econ_features = ['1人当たり県民所得', '医師数10万対', '教育費（二人以上の世帯）']
bar_colors_perm = ['#E53935' if feature_cols[i] in econ_features else '#1565C0'
                   for i in perm_order]
ax_b.barh(np.arange(len(perm_vals)), perm_vals, xerr=perm_std, capsize=3,
          color=bar_colors_perm, alpha=0.8, edgecolor='white',
          error_kw={'elinewidth': 1.2})
ax_b.set_yticks(np.arange(len(perm_vals)))
ax_b.set_yticklabels(perm_labels, fontsize=9)
ax_b.set_xlabel('Permutation Importance（±1SD）', fontsize=11)
ax_b.set_title('Permutation Importance\n(SHAP値の代理指標, 赤=所得・医療関連変数)',
               fontsize=11, fontweight='bold')
ax_b.axvline(0, color='gray', linewidth=0.8, linestyle='--')
ax_b.grid(axis='x', alpha=0.3)
ax_b.invert_yaxis()

plt.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2025_U5_3_fig2_importance.png'), bbox_inches='tight', dpi=150)
plt.close(fig2)
print("  → 2025_U5_3_fig2_importance.png 保存完了")

print("図3: SHAPサマリープロット（代理）を作成中...")

fig3, axes3 = plt.subplots(1, 2, figsize=(14, 6))
fig3.suptitle('SHAP値による特徴量効果の解釈（Random Forest）\n（実データ：SSDSE-B/E-2026）',
              fontsize=13, fontweight='bold')

pred_vals = rf.predict(X)

# パネル(a): 1人当たり県民所得 × RF予測値 Dependence Plot
ax_a = axes3[0]
income_idx = feature_cols.index('1人当たり県民所得')
feat_income = df_model['1人当たり県民所得'].values
color_arr   = df_model['医師数10万対'].values
sc = ax_a.scatter(feat_income, pred_vals,
                  c=color_arr, cmap='RdYlGn', alpha=0.8,
                  s=60, edgecolors='white', linewidth=0.5, zorder=3)
cbar = plt.colorbar(sc, ax=ax_a, shrink=0.85)
cbar.set_label('医師数（10万人対）', fontsize=9)
# 都道府県ラベル
for i, (xi, yi) in enumerate(zip(feat_income, pred_vals)):
    if yi > np.percentile(pred_vals, 85) or yi < np.percentile(pred_vals, 15):
        ax_a.annotate(df_model.loc[i, '都道府県'][:2], (xi, yi),
                      fontsize=7, ha='left', va='bottom', color='#333')
coef_inc = np.polyfit(feat_income, pred_vals, 2)
x_inc    = np.linspace(feat_income.min(), feat_income.max(), 200)
ax_a.plot(x_inc, np.polyval(coef_inc, x_inc), 'b-', linewidth=2, label='二次フィット')
ax_a.set_xlabel(FEATURE_LABELS['1人当たり県民所得'], fontsize=11)
ax_a.set_ylabel('RF予測値（大学進学率 %）', fontsize=11)
ax_a.set_title('Dependence Plot\n（県民所得 × RF予測値, 色：医師数）',
               fontsize=11, fontweight='bold')
ax_a.legend(fontsize=9)
ax_a.grid(True, alpha=0.2)

# パネル(b): 医師数10万対（高 vs 低）で予測値分布を比較
ax_b = axes3[1]
doc_idx = feature_cols.index('医師数10万対')
doc_med = np.median(X[:, doc_idx])
high_mask = X[:, doc_idx] >= doc_med
low_mask  = ~high_mask

y_high = pred_vals[high_mask]
y_low  = pred_vals[low_mask]

# バイオリンプロット風にKDE近似で表示（stripplot風）
jitter_scale = 0.08
jp = np.full(high_mask.sum(), 1.0)
jq = np.full(low_mask.sum(),  0.0)
ax_b.scatter(jp + (np.arange(high_mask.sum()) % 5 - 2) * jitter_scale * 0.4,
             y_high, color='#E53935', alpha=0.6, s=40,
             label=f'医師数: 高（≥中央値 {doc_med:.0f}人）')
ax_b.scatter(jq + (np.arange(low_mask.sum()) % 5 - 2) * jitter_scale * 0.4,
             y_low, color='#1565C0', alpha=0.6, s=40,
             label=f'医師数: 低（<中央値 {doc_med:.0f}人）')
ax_b.axhline(y_high.mean(), color='#E53935', linestyle='--', linewidth=2,
             label=f'高群平均: {y_high.mean():.1f}%')
ax_b.axhline(y_low.mean(),  color='#1565C0',  linestyle='--', linewidth=2,
             label=f'低群平均: {y_low.mean():.1f}%')
ax_b.set_xticks([0, 1])
ax_b.set_xticklabels(['医師数: 低', '医師数: 高'], fontsize=11)
ax_b.set_ylabel('RF予測値（大学進学率 %）', fontsize=11)
ax_b.set_title('医師数（高 vs 低）の予測値分布\n→ 医師数が多い都道府県で進学率予測が高い',
               fontsize=11, fontweight='bold')
ax_b.legend(fontsize=8.5, loc='upper left')
ax_b.grid(axis='y', alpha=0.3)

plt.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2025_U5_3_fig3_shap.png'), bbox_inches='tight', dpi=150)
plt.close(fig3)
print("  → 2025_U5_3_fig3_shap.png 保存完了")

print("図4: 所得×医師数 交互作用図を作成中...")

fig4, axes4 = plt.subplots(1, 2, figsize=(14, 6))
fig4.suptitle('交互作用分析：1人当たり県民所得 × 医師数（10万人対）\n（実データ：SSDSE-B/E-2026）',
              fontsize=13, fontweight='bold')

# (a) 線形交互作用プロット（グループ別平均大学進学率）
ax4a = axes4[0]
income_labels3 = ['低所得\n（低1/3）', '中所得\n（中1/3）', '高所得\n（高1/3）']
income_q33 = df_model['1人当たり県民所得'].quantile(0.33)
income_q67 = df_model['1人当たり県民所得'].quantile(0.67)

def income_group(x):
    if x < income_q33:
        return '低所得\n（低1/3）'
    elif x < income_q67:
        return '中所得\n（中1/3）'
    else:
        return '高所得\n（高1/3）'

df_plot = df_model.copy()
df_plot['所得3区分'] = df_plot['1人当たり県民所得'].apply(income_group)
df_plot['医師グループ'] = df_plot['医師数10万対'].apply(
    lambda x: '医師多（中央値以上）' if x >= doctor_med else '医師少（中央値未満）')

for grp, col, ls in [('医師多（中央値以上）', '#E53935', '-'), ('医師少（中央値未満）', '#1565C0', '--')]:
    d = df_plot[df_plot['医師グループ'] == grp].groupby('所得3区分', observed=True)['大学進学率'].mean()
    vals = [d.get(lab, np.nan) for lab in income_labels3]
    ax4a.plot(income_labels3, vals, 'o-', color=col, linestyle=ls,
              linewidth=2.5, markersize=9, markerfacecolor='white',
              markeredgewidth=2.5, label=grp)

ax4a.set_xlabel('1人当たり県民所得（区分）', fontsize=11)
ax4a.set_ylabel('平均大学進学率（%）', fontsize=11)
ax4a.set_title('交互作用プロット\n（医師多 vs 少）', fontsize=11, fontweight='bold')
ax4a.legend(fontsize=10)
ax4a.grid(True, alpha=0.3)

# (b) 2D等高線プロット（RF予測面）
ax4b = axes4[1]
income_grid = np.linspace(df_model['1人当たり県民所得'].min(),
                           df_model['1人当たり県民所得'].max(), 50)
doctor_grid = np.linspace(df_model['医師数10万対'].min(),
                           df_model['医師数10万対'].max(), 50)
INC, DOC = np.meshgrid(income_grid, doctor_grid)

# 平均的な観測値を固定し、income と doctor だけ変化させる
X_grid = np.tile(X.mean(axis=0), (50 * 50, 1))
inc_idx = feature_cols.index('1人当たり県民所得')
doc_idx = feature_cols.index('医師数10万対')
X_grid[:, inc_idx] = INC.ravel()
X_grid[:, doc_idx] = DOC.ravel()

pred_grid = rf.predict(X_grid).reshape(50, 50)

cs = ax4b.contourf(INC, DOC, pred_grid, levels=20, cmap='RdYlGn')
cbar4 = plt.colorbar(cs, ax=ax4b)
cbar4.set_label('RF予測 大学進学率（%）', fontsize=9)

# 実データ点をオーバーレイ
ax4b.scatter(df_model['1人当たり県民所得'], df_model['医師数10万対'],
             c=df_model['大学進学率'], cmap='RdYlGn', s=40, alpha=0.7,
             edgecolors='black', linewidth=0.5, zorder=5)
# 都道府県ラベル（上位・下位）
for _, row in df_model.iterrows():
    if row['大学進学率'] > df_model['大学進学率'].quantile(0.85) or \
       row['大学進学率'] < df_model['大学進学率'].quantile(0.15):
        ax4b.annotate(row['都道府県'][:2],
                      (row['1人当たり県民所得'], row['医師数10万対']),
                      fontsize=7, color='black', ha='left', va='bottom')

ax4b.set_xlabel('1人当たり県民所得（万円）', fontsize=11)
ax4b.set_ylabel('医師数（10万人対）', fontsize=11)
ax4b.set_title('RF予測面（等高線）\n県民所得 × 医師数の交互作用',
               fontsize=11, fontweight='bold')
ax4b.text(0.05, 0.95, '→ 右上（高所得×医師数多）\nで大学進学率の予測が最大',
          transform=ax4b.transAxes, fontsize=9.5, va='top',
          bbox=dict(boxstyle='round', facecolor='white', alpha=0.85, edgecolor='#aaa'))

plt.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2025_U5_3_fig4_interact.png'), bbox_inches='tight', dpi=150)
plt.close(fig4)
print("  → 2025_U5_3_fig4_interact.png 保存完了")