""" 教育用再現コード: 2023年 統計データ分析コンペティション 審査員奨励賞 [高校生の部] ================================================================= 論文タイトル：都道府県ごとの学力の差 著者：（高校生の部 審査員奨励賞） 【分析概要】 データ：SSDSE-B-2026.csv（都道府県データ） 2022年度 47都道府県 目的変数：大学進学率（E4602/E4601 × 100） 高校卒業後大学進学者数 ÷ 高校卒業者数 Step1. 都道府県間の大学進学率の可視化（横棒グラフ） Step2. 各説明変数との相関係数（ピアソン） Step3. ステップワイズ変数選択（バックワードAIC）による重回帰 Step4. 散布図（教育費割合 vs 大学進学率） 【説明変数】 教育費割合 : L322108/L3221 × 100（家計の教育費が消費支出に占める割合） 消費支出 : L3221（所得水準の代理指標） 小学教員/児童比 : E2401/E2501（小さいほど手厚い指導） 中学教員/生徒比 : E3401/E3501 住宅地価格 : C5401（富裕度の指標） 高齢化率 : A1303/A1101 × 100 宿泊者per capita: G7101/A1101（都市・観光地の代理指標） 年平均気温 : B4101 【データ出典】 SSDSE-B-2026.csv: 社会・人口統計体系（都道府県データ） 統計センターより公表の実データ 【データサイエンス学習ポイント】 1. 都道府県間の格差の記述統計・可視化 2. ピアソン相関係数と無相関検定 3. ステップワイズ変数選択（バックワードAIC） 4. 標準化係数（βeta coefficient）による効果量比較 ================================================================= """ # ============================================================ # 【データの準備】実行前に以下のデータファイルを用意してください # # 必要ファイル: # ・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/2023_H5_2_shorei.py # ============================================================ import pandas as pd import numpy as np import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt import statsmodels.api as sm from scipy import stats import warnings warnings.filterwarnings('ignore') # ────────────────────────────────────────────────────────────── # 共通設定 # ────────────────────────────────────────────────────────────── plt.rcParams['font.family'] = 'Hiragino Sans' plt.rcParams['axes.unicode_minus'] = False plt.rcParams['figure.dpi'] = 150 import os FIG_DIR = 'html/figures' DATA_DIR = 'data/raw' os.makedirs(FIG_DIR, exist_ok=True) # ================================================================ # ■ 実データ読み込み（SSDSE-B-2026: 都道府県データ） # ================================================================ df_raw = pd.read_csv( os.path.join(DATA_DIR, 'SSDSE-B-2026.csv'), encoding='cp932', header=None ) # 行0がコード行・行1がラベル行・行2以降がデータ df_raw.columns = df_raw.iloc[0] df_all = df_raw.iloc[2:].copy().reset_index(drop=True) # 都道府県コード（R\d{5}形式）に絞る mask_pref = df_all['Code'].str.match(r'^R\d{5}$', na=False) df_pref = df_all[mask_pref].copy() # 2022年度のみ抽出 df_pref['年度_str'] = df_pref['SSDSE-B-2026'].astype(str).str.strip() df = df_pref[df_pref['年度_str'] == '2022'].copy().reset_index(drop=True) print("=" * 60) print(f"■ 対象都道府県数: N={len(df)}（2022年度）") print("=" * 60) # ── 必要列を数値変換 ───────────────────────────────────────────── NUM_COLS = [ 'E4602', 'E4601', # 大学進学者数・高校卒業者数 'L322108', 'L3221', # 教育費・消費支出 'E2401', 'E2501', # 小学教員数・児童数 'E3401', 'E3501', # 中学教員数・生徒数 'C5401', # 住宅地標準価格 'A1303', 'A1101', # 65歳以上人口・総人口 'G7101', # 延べ宿泊者数 'B4101', # 年平均気温 ] for c in NUM_COLS: df[c] = pd.to_numeric(df[c], errors='coerce') # ── 目的変数 ───────────────────────────────────────────────────── df['大学進学率'] = df['E4602'] / df['E4601'] * 100 # ── 説明変数 ───────────────────────────────────────────────────── df['教育費割合'] = df['L322108'] / df['L3221'] * 100 df['消費支出'] = df['L3221'] df['小学教員児童比'] = df['E2401'] / df['E2501'] df['中学教員生徒比'] = df['E3401'] / df['E3501'] df['住宅地価格'] = df['C5401'] df['高齢化率'] = df['A1303'] / df['A1101'] * 100 df['宿泊者per capita'] = df['G7101'] / df['A1101'] df['年平均気温'] = df['B4101'] VAR_NAMES = [ '教育費割合', '消費支出', '小学教員児童比', '中学教員生徒比', '住宅地価格', '高齢化率', '宿泊者per capita', '年平均気温' ] VAR_LABELS = { '教育費割合': '教育費割合（%）', '消費支出': '消費支出（円/月）', '小学教員児童比': '小学教員/児童比', '中学教員生徒比': '中学教員/生徒比', '住宅地価格': '住宅地価格（千円/m²）', '高齢化率': '高齢化率（%）', '宿泊者per capita': '宿泊者per capita', '年平均気温': '年平均気温（℃）', } # 都道府県名 pref_names = df['Prefecture'].values # 有効データ確認 y = df['大学進学率'].values X_df = df[VAR_NAMES].copy() valid_mask = ( np.isfinite(y) & X_df.apply(lambda col: np.isfinite(col)).all(axis=1).values ) df_valid = df[valid_mask].reset_index(drop=True) y = df_valid['大学進学率'].values X_df = df_valid[VAR_NAMES].copy() pref_names = df_valid['Prefecture'].values N = len(y) print(f"有効都道府県数（欠損除外後）: N={N}") print() print(df_valid[['Prefecture', '大学進学率'] + VAR_NAMES].describe().round(3).to_string()) # ================================================================ # ■ Step1. 相関係数（ピアソン） # ================================================================ print("\n" + "=" * 60) print("■ Step1. ピアソン相関係数（目的変数：大学進学率）") print("=" * 60) print(f"\n {'変数':<20} {'r':>8} {'p値':>10} {'有意':>6}") print(" " + "-" * 48) corr_vals = {} pval_vals = {} for name in VAR_NAMES: r, p = stats.pearsonr(X_df[name].values, y) corr_vals[name] = r pval_vals[name] = p sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else 'n.s.' print(f" {name:<20} {r:>8.4f} {p:>10.4f} {sig:>6}") # ================================================================ # ■ Step2. バックワードAIC変数選択 # ================================================================ def backward_aic(X_df, y): """バックワードAIC変数選択。AICが改善しなくなるまで変数を除外する。""" cols = list(X_df.columns) step = 0 print("\n" + "=" * 60) print("■ Step2. バックワードAIC変数選択") print("=" * 60) while len(cols) > 1: full_model = sm.OLS(y, sm.add_constant(X_df[cols])).fit() full_aic = full_model.aic aic_drop = {} for c in cols: sub = [x for x in cols if x != c] m = sm.OLS(y, sm.add_constant(X_df[sub])).fit() aic_drop[c] = m.aic best_drop = min(aic_drop, key=aic_drop.get) if aic_drop[best_drop] < full_aic: step += 1 print(f" Step{step}: '{best_drop}' を除外" f"（AIC: {full_aic:.3f} → {aic_drop[best_drop]:.3f}）") cols.remove(best_drop) else: break return cols selected_cols = backward_aic(X_df, y) print(f"\n 最終選択変数: {selected_cols}") # ── 最終モデルの推定 ────────────────────────────────────────────── X_final = sm.add_constant(X_df[selected_cols]) final_model = sm.OLS(y, X_final).fit() print(f"\n R² = {final_model.rsquared:.4f}") print(f" Adj. R² = {final_model.rsquared_adj:.4f}") print(f" AIC = {final_model.aic:.3f}") print() print(final_model.summary2()) # ── 標準化係数（βeta） ────────────────────────────────────────── y_std = (y - y.mean()) / y.std() X_sel_std = X_df[selected_cols].copy() for c in selected_cols: X_sel_std[c] = (X_sel_std[c] - X_sel_std[c].mean()) / X_sel_std[c].std() std_model = sm.OLS(y_std, sm.add_constant(X_sel_std)).fit() beta_coefs = std_model.params[1:].values # 定数項除く beta_ses = std_model.bse[1:].values beta_pvals = std_model.pvalues[1:].values print("\n標準化係数（β）:") for name, b, p in zip(selected_cols, beta_coefs, beta_pvals): sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else 'n.s.' print(f" {name:<20}: β={b:.4f} p={p:.4f} {sig}") # ================================================================ # ■ 図の生成（4枚） # ================================================================ # ──────────────────────────────────────────────────────────────── # 図1: 47都道府県の大学進学率 横棒グラフ（降順、上位10強調） # ──────────────────────────────────────────────────────────────── fig1, ax1 = plt.subplots(figsize=(10, 12)) sort_idx = np.argsort(y) y_sorted = y[sort_idx] names_sorted = pref_names[sort_idx] top10_threshold = np.sort(y)[-10] colors_bar = ['#C62828' if v >= top10_threshold else '#1565C0' for v in y_sorted] alphas_bar = [0.95 if v >= top10_threshold else 0.65 for v in y_sorted] bars = ax1.barh( np.arange(N), y_sorted, color=colors_bar, alpha=0.75, edgecolor='white', height=0.7 ) for bar, alpha in zip(bars, alphas_bar): bar.set_alpha(alpha) ax1.set_yticks(np.arange(N)) ax1.set_yticklabels(names_sorted, fontsize=9) ax1.set_xlabel('大学進学率（%）', fontsize=12) ax1.set_title('都道府県別 大学進学率（2022年度）\n（赤：上位10都道府県）', fontsize=13, fontweight='bold', pad=12) ax1.axvline(y.mean(), color='#FF8F00', linestyle='--', linewidth=1.5, label=f'全国平均 {y.mean():.1f}%') ax1.grid(axis='x', alpha=0.3) ax1.legend(fontsize=10) for i, (v, name) in enumerate(zip(y_sorted, names_sorted)): ax1.text(v + 0.3, i, f'{v:.1f}%', va='center', fontsize=7.5, color='#C62828' if v >= top10_threshold else '#555') plt.tight_layout() fig1.savefig(os.path.join(FIG_DIR, '2023_H5_2_fig1_bar.png'), bbox_inches='tight', dpi=150) plt.close(fig1) print("\n図1保存: 2023_H5_2_fig1_bar.png") # ──────────────────────────────────────────────────────────────── # 図2: 相関係数棒グラフ（全説明変数 vs 大学進学率） # ──────────────────────────────────────────────────────────────── fig2, ax2 = plt.subplots(figsize=(9, 5)) corr_list = [corr_vals[n] for n in VAR_NAMES] pval_list = [pval_vals[n] for n in VAR_NAMES] colors2 = ['#C62828' if p < 0.05 else '#9E9E9E' for p in pval_list] y_pos2 = np.arange(len(VAR_NAMES)) ax2.barh(y_pos2, corr_list, color=colors2, alpha=0.75, edgecolor='white', height=0.6) ax2.axvline(0, color='gray', linestyle='--', linewidth=1.0) ax2.set_yticks(y_pos2) ax2.set_yticklabels(VAR_NAMES, fontsize=10) ax2.set_xlabel('ピアソン相関係数', fontsize=11) ax2.set_title('各説明変数と大学進学率の相関係数\n（無相関検定: 赤=p<0.05）', fontsize=12, fontweight='bold') ax2.invert_yaxis() ax2.grid(axis='x', alpha=0.3) from matplotlib.patches import Patch ax2.legend( handles=[ Patch(color='#C62828', alpha=0.75, label='有意（p<0.05）'), Patch(color='#9E9E9E', alpha=0.75, label='非有意'), ], fontsize=9 ) for i, (r, p) in enumerate(zip(corr_list, pval_list)): sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else '' offset = 0.01 if r >= 0 else -0.01 ha = 'left' if r >= 0 else 'right' ax2.text(r + offset, i, f' {r:.3f}{sig}', va='center', ha=ha, fontsize=8.5) plt.tight_layout() fig2.savefig(os.path.join(FIG_DIR, '2023_H5_2_fig2_corr.png'), bbox_inches='tight', dpi=150) plt.close(fig2) print("図2保存: 2023_H5_2_fig2_corr.png") # ──────────────────────────────────────────────────────────────── # 図3: 最終OLSモデルの標準化係数 # ──────────────────────────────────────────────────────────────── fig3, ax3 = plt.subplots(figsize=(9, max(4, len(selected_cols) * 0.9))) colors3 = ['#C62828' if p < 0.01 else '#FF8F00' if p < 0.05 else '#9E9E9E' for p in beta_pvals] y_pos3 = np.arange(len(selected_cols)) ax3.barh(y_pos3, beta_coefs, color=colors3, alpha=0.75, edgecolor='white', height=0.6) ax3.errorbar(beta_coefs, y_pos3, xerr=1.96 * beta_ses, fmt='none', color='#333', capsize=4, linewidth=1.2) ax3.axvline(0, color='gray', linestyle='--', linewidth=1.0) ax3.set_yticks(y_pos3) ax3.set_yticklabels(selected_cols, fontsize=10) ax3.set_xlabel('標準化係数 β（±1.96SE）', fontsize=11) ax3.set_title( f'バックワードAIC選択後の重回帰 — 標準化係数\n' f'R²={final_model.rsquared:.3f} Adj.R²={final_model.rsquared_adj:.3f} N={N}都道府県', fontsize=12, fontweight='bold' ) ax3.invert_yaxis() ax3.grid(axis='x', alpha=0.3) ax3.legend( handles=[ Patch(color='#C62828', alpha=0.75, label='p<0.01'), Patch(color='#FF8F00', alpha=0.75, label='p<0.05'), Patch(color='#9E9E9E', alpha=0.75, label='n.s.'), ], fontsize=9 ) for i, (b, p) in enumerate(zip(beta_coefs, beta_pvals)): sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else '' offset = 0.01 if b >= 0 else -0.01 ha = 'left' if b >= 0 else 'right' ax3.text(b + offset, i, f' {b:.3f}{sig}', va='center', ha=ha, fontsize=8.5) plt.tight_layout() fig3.savefig(os.path.join(FIG_DIR, '2023_H5_2_fig3_beta.png'), bbox_inches='tight', dpi=150) plt.close(fig3) print("図3保存: 2023_H5_2_fig3_beta.png") # ──────────────────────────────────────────────────────────────── # 図4: 散布図（教育費割合 vs 大学進学率）+ 回帰直線 + 都道府県ラベル # ──────────────────────────────────────────────────────────────── fig4, ax4 = plt.subplots(figsize=(11, 7)) x4 = X_df['教育費割合'].values r4, p4 = stats.pearsonr(x4, y) ax4.scatter(x4, y, color='#1565C0', s=50, alpha=0.7, zorder=3, edgecolors='white', linewidth=0.5) # 回帰直線 z4 = np.polyfit(x4, y, 1) xs4 = np.linspace(x4.min() - 0.05, x4.max() + 0.05, 200) ax4.plot(xs4, np.poly1d(z4)(xs4), color='#C62828', linewidth=1.8, alpha=0.85, label=f'回帰直線 r={r4:.3f}（p={p4:.4f}）') # 都道府県ラベル for xi, yi, name in zip(x4, y, pref_names): # 短縮ラベル label = name.replace('都', '').replace('道', '').replace('府', '').replace('県', '') ax4.text(xi + 0.01, yi + 0.15, label, fontsize=7, color='#333', ha='left', va='bottom', alpha=0.85) ax4.set_xlabel('教育費割合（消費支出に占める割合 %）', fontsize=12) ax4.set_ylabel('大学進学率（%）', fontsize=12) ax4.set_title('教育費割合と大学進学率の関係（2022年度・47都道府県）', fontsize=13, fontweight='bold') ax4.legend(fontsize=10) ax4.grid(True, alpha=0.3) plt.tight_layout() fig4.savefig(os.path.join(FIG_DIR, '2023_H5_2_fig4_scatter.png'), bbox_inches='tight', dpi=150) plt.close(fig4) print("図4保存: 2023_H5_2_fig4_scatter.png") print("\n全図の生成完了（4枚）") print(" fig1_bar.png : 都道府県別大学進学率（横棒グラフ）") print(" fig2_corr.png : 各変数との相関係数棒グラフ") print(" fig3_beta.png : 最終OLSモデル標準化係数") print(" fig4_scatter.png: 教育費割合 vs 大学進学率 散布図")