""" 教育用再現コード: 2025年 統計データ分析コンペティション 審査員奨励賞（大学生） ================================================================================= 論文タイトル：デジタル教科書は学力にプラスかマイナスか？機械学習×SHAP解析 受賞区分 ：審査員奨励賞 ［大学生・一般の部］ 【分析概要】 データ：SSDSE-B-2026 + SSDSE-E-2026（47都道府県、実公的統計データ） 目的変数：大学進学率 = 高等学校卒業者のうち進学者数 / 高等学校卒業者数 × 100 ※実際の教育成果を示す公式統計量 Step1. 実データ読み込み・特徴量構築 Step2. 特徴量エンジニアリング Step3. 3モデル（RF / GradientBoosting / Ridge）の学習・評価 Step4. 特徴量重要度 + Permutation Importance（SHAP代理） Step5. 交互作用分析（1人当たり県民所得 × 医師数10万対） 【使用モデル】 RandomForestRegressor（sklearn） GradientBoostingRegressor（sklearn, XGBoost代替） Ridge regression（sklearn） 【使用データ】 SSDSE-B-2026.csv + SSDSE-E-2026.csv （統計数理研究所 教育用標準データセット） https://www.ism.ac.jp/editsection/ssdse/ 【データサイエンス学習ポイント】 1. アンサンブル学習（バギング=RF / ブースティング=GBM） 2. Permutation Importance による機械学習モデルの解釈（SHAP代理） 3. 特徴量重要度の比較（不純度減少 vs Permutation） 4. 交互作用効果の可視化（Dependence Plot + 等高線） 5. 実公的統計データの結合・前処理 【出力図】 html/figures/2025_U5_3_fig1_overview.png ... データ概要と目的変数の分布 html/figures/2025_U5_3_fig2_importance.png ... 特徴量重要度（3モデル比較） html/figures/2025_U5_3_fig3_shap.png ... SHAP サマリープロット（RF） html/figures/2025_U5_3_fig4_interact.png ... 所得×医師数の交互作用 ================================================================================= """ # ============================================================ # 【データの準備】実行前に以下のデータファイルを用意してください # # 必要ファイル: # ・SSDSE-B-2026.csv # → data/raw/SSDSE-B-2026.csv に配置 # ・SSDSE-E-2026.csv # → data/raw/SSDSE-E-2026.csv に配置 # # ダウンロード先: # https://www.nstac.go.jp/use/literacy/ssdse/ # （SSDSE-B（社会・人口統計体系 都道府県データ） の CSV をダウンロード） # （SSDSE-E（社会・人口統計体系 都道府県の指標2） の CSV をダウンロード） # # フォルダ配置（プロジェクトルートからの相対パス）: # code/ ← このスクリプトの場所 # data/raw/ ← CSV ファイルをここに配置 # html/figures/ ← 図の出力先（自動生成） # # 実行方法（ファイルを一切編集せず実行可能）: # python3 code/2025_U5_3_shorei.py # ============================================================ import os import numpy as np import pandas as pd import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt import matplotlib.patches as mpatches from matplotlib.gridspec import GridSpec from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor from sklearn.linear_model import Ridge from sklearn.inspection import permutation_importance from sklearn.model_selection import cross_val_score, KFold from sklearn.preprocessing import StandardScaler from sklearn.metrics import r2_score, mean_squared_error import warnings warnings.filterwarnings('ignore') # ── パス設定 ───────────────────────────────────────────────────────────────── FIG_DIR = 'html/figures' DATA_DIR = 'data/raw' os.makedirs(FIG_DIR, exist_ok=True) plt.rcParams.update({ 'font.family': 'Hiragino Sans', 'axes.unicode_minus': False, 'figure.dpi': 150, 'axes.spines.top': False, 'axes.spines.right': False, }) # ═══════════════════════════════════════════════════════════════════════════════ # ■ Step1. 実データ読み込み # ═══════════════════════════════════════════════════════════════════════════════ print("=" * 65) print("■ 実データ読み込み（SSDSE-B-2026 + SSDSE-E-2026）") print("=" * 65) # ── SSDSE-B-2026（2022年度、47都道府県） ────────────────────────────────── df_b = pd.read_csv(os.path.join(DATA_DIR, 'SSDSE-B-2026.csv'), encoding='cp932', header=1) mask_b = (df_b['年度'] == 2022) & df_b['地域コード'].str.match(r'^R\d{5}$', na=False) df_b22 = df_b[mask_b].copy().reset_index(drop=True) num_cols_b = [ '総人口', '65歳以上人口', '15～64歳人口', '保育所等数', '年平均気温', '高等学校卒業者のうち進学者数', '高等学校卒業者数', '婚姻件数', '教育費（二人以上の世帯）', ] for c in num_cols_b: df_b22[c] = pd.to_numeric(df_b22[c], errors='coerce') print(f" SSDSE-B 2022年: {len(df_b22)} 都道府県") # ── SSDSE-E-2026（47都道府県） ───────────────────────────────────────────── df_e_raw = pd.read_csv(os.path.join(DATA_DIR, 'SSDSE-E-2026.csv'), encoding='cp932') df_e = df_e_raw.iloc[2:].copy() df_e.columns = df_e_raw.iloc[1].values df_e = df_e[df_e['都道府県'] != '全国'].reset_index(drop=True) num_cols_e = [ '医師数', '1人当たり県民所得（平成27年基準）', '総面積（北方地域及び竹島を除く）', '総人口', ] for c in num_cols_e: df_e[c] = pd.to_numeric(df_e[c], errors='coerce') print(f" SSDSE-E: {len(df_e)} 都道府県") # ── データ結合（都道府県名でマージ） ────────────────────────────────────── df = df_b22[['都道府県'] + num_cols_b].merge( df_e[['都道府県', '医師数', '1人当たり県民所得（平成27年基準）', '総面積（北方地域及び竹島を除く）']], on='都道府県', how='inner' ) print(f" 結合後: {len(df)} 都道府県") # ═══════════════════════════════════════════════════════════════════════════════ # ■ Step2. 特徴量エンジニアリング # ═══════════════════════════════════════════════════════════════════════════════ # 目的変数：大学進学率（%） 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}%") # ═══════════════════════════════════════════════════════════════════════════════ # ■ Step3. 3モデルの学習・評価 # ═══════════════════════════════════════════════════════════════════════════════ print("\n" + "=" * 65) print("■ Step3. 3モデルの学習・LOO-CV 評価（少標本のため Leave-One-Out）") print("=" * 65) from sklearn.model_selection import LeaveOneOut def loo_metrics(model, X_in, y_in): """LOO予測をループで収集し R² と RMSE を返す（nan を回避）""" loo = LeaveOneOut() preds = np.empty(len(y_in)) for train_idx, test_idx in loo.split(X_in): model.fit(X_in[train_idx], y_in[train_idx]) preds[test_idx] = model.predict(X_in[test_idx]) loo_r2 = r2_score(y_in, preds) loo_rmse = np.sqrt(mean_squared_error(y_in, preds)) return loo_r2, loo_rmse # ── Random Forest ───────────────────────────────────────────────────────────── rf = RandomForestRegressor(n_estimators=300, max_depth=4, min_samples_leaf=3, random_state=42, n_jobs=-1) rf.fit(X, y) rf_loo_r2, rf_loo_rmse = loo_metrics( RandomForestRegressor(n_estimators=300, max_depth=4, min_samples_leaf=3, random_state=42, n_jobs=-1), X, y) # ── Gradient Boosting（XGBoost代替） ────────────────────────────────────────── gbm = GradientBoostingRegressor(n_estimators=200, max_depth=3, learning_rate=0.05, subsample=0.8, random_state=42) gbm.fit(X, y) gbm_loo_r2, gbm_loo_rmse = loo_metrics( GradientBoostingRegressor(n_estimators=200, max_depth=3, learning_rate=0.05, subsample=0.8, random_state=42), X, y) # ── Ridge Regression ────────────────────────────────────────────────────────── scaler = StandardScaler() X_sc = scaler.fit_transform(X) ridge = Ridge(alpha=1.0) ridge.fit(X_sc, y) ridge_loo_r2, ridge_loo_rmse = loo_metrics( Ridge(alpha=1.0), X_sc, y) print(f"\n {'モデル':<30} {'LOO R²':>8} {'LOO RMSE':>10} {'Train R²':>10}") print(" " + "-" * 62) print(f" {'Random Forest':<30} {rf_loo_r2:>8.4f} {rf_loo_rmse:>10.4f} {r2_score(y, rf.predict(X)):>10.4f}") print(f" {'Gradient Boosting (XGB代替)':<30} {gbm_loo_r2:>8.4f} {gbm_loo_rmse:>10.4f} {r2_score(y, gbm.predict(X)):>10.4f}") print(f" {'Ridge Regression':<30} {ridge_loo_r2:>8.4f} {ridge_loo_rmse:>10.4f} {r2_score(y, ridge.predict(X_sc)):>10.4f}") # ═══════════════════════════════════════════════════════════════════════════════ # ■ Step4. 特徴量重要度 + Permutation Importance（SHAP代理） # ═══════════════════════════════════════════════════════════════════════════════ print("\n" + "=" * 65) print("■ Step4. 特徴量重要度分析（Permutation Importance）") print("=" * 65) rf_imp = rf.feature_importances_ gbm_imp = gbm.feature_importances_ # Ridge: 係数の絶対値を重要度の代理とする（標準化済み） ridge_imp = np.abs(ridge.coef_) / np.abs(ridge.coef_).sum() # Permutation Importance（SHAP値の代理） perm_rf = permutation_importance(rf, X, y, n_repeats=30, random_state=42, n_jobs=-1) perm_importances = perm_rf.importances_mean order = np.argsort(rf_imp)[::-1] print(f"\n {'特徴量':<24} {'RF重要度':>10} {'GBM重要度':>10} {'Ridge係数':>10} {'Perm.Imp.':>10}") print(" " + "-" * 70) for idx in order: fn = feature_cols[idx] print(f" {FEATURE_LABELS[fn]:<24} {rf_imp[idx]:>10.4f} {gbm_imp[idx]:>10.4f}" f" {ridge_imp[idx]:>10.4f} {perm_importances[idx]:>10.4f}") # ═══════════════════════════════════════════════════════════════════════════════ # ■ Step5. 交互作用分析：1人当たり県民所得 × 医師数10万対 # ═══════════════════════════════════════════════════════════════════════════════ 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 → 高所得×医師数多の組み合わせで進学率が最も高い傾向を確認") # ═══════════════════════════════════════════════════════════════════════════════ # ■ 図の生成（4枚） # ═══════════════════════════════════════════════════════════════════════════════ # ──────────────────────────────────────────────────────────────────────────── # 図1：データ概要と目的変数の分布 # ──────────────────────────────────────────────────────────────────────────── 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 保存完了") # ──────────────────────────────────────────────────────────────────────────── # 図2：特徴量重要度（3モデル比較） # ──────────────────────────────────────────────────────────────────────────── 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 保存完了") # ──────────────────────────────────────────────────────────────────────────── # 図3：SHAP サマリープロット（代理可視化 + Dependence Plot） # ──────────────────────────────────────────────────────────────────────────── 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 保存完了") # ──────────────────────────────────────────────────────────────────────────── # 図4：1人当たり県民所得 × 医師数10万対 交互作用分析 # ──────────────────────────────────────────────────────────────────────────── 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 保存完了") # ═══════════════════════════════════════════════════════════════════════════════ # ■ 最終サマリー # ═══════════════════════════════════════════════════════════════════════════════ print("\n" + "=" * 65) print("✓ 全図の生成完了（4枚）") print("=" * 65) print(f"\n データ: SSDSE-B-2026 + SSDSE-E-2026（実公的統計データ）") print(f" 総観測数: {len(df_model)} 都道府県（2022年度）") print(f" 目的変数: 大学進学率 平均 {y.mean():.1f}%（SD: {y.std():.1f}%）") print(f"\n 【モデル評価（LOO-CV R²）】") print(f" Random Forest: {rf_loo_r2:.4f}") print(f" Gradient Boosting: {gbm_loo_r2:.4f}") print(f" Ridge Regression: {ridge_loo_r2:.4f}") print(f"\n 【上位特徴量（Permutation Importance, RF）】") perm_order_all = np.argsort(perm_importances)[::-1] for i in perm_order_all[:5]: print(f" {FEATURE_LABELS[feature_cols[i]]:<24} {perm_importances[i]:.4f}") print(f"\n 【主要知見】") print(f" - 1人当たり県民所得と医師数10万対が大学進学率の主要な予測因子") print(f" - 高所得×医師数多の都道府県で大学進学率が最も高い傾向") print(f" - 高齢化率は大学進学率と負の相関（人口構造の影響）")