""" 2023_U5_1_shorei.py パネルデータを用いた進学と就職時による人口流出の要因分析 審査員奨励賞 [大学生・一般の部] 2023年度 実データ（SSDSE-B-2026.csv）による教育用再現コード ※合成データ・乱数生成は一切使用しない """ # ============================================================ # 【データの準備】実行前に以下のデータファイルを用意してください # # 必要ファイル: # ・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_U5_1_shorei.py # ============================================================ import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt import matplotlib.patches as mpatches import numpy as np import pandas as pd from scipy import stats from linearmodels import PanelOLS plt.rcParams['font.family'] = 'Hiragino Sans' plt.rcParams['axes.unicode_minus'] = False # ---------------------------------------------------------------- # パス設定 # ---------------------------------------------------------------- import os FIG_DIR = os.path.normpath('html/figures') + os.sep DATA_B = 'data/raw/SSDSE-B-2026.csv' os.makedirs(FIG_DIR, exist_ok=True) # ---------------------------------------------------------------- # 実データ読み込み: SSDSE-B-2026（パネルデータ 47都道府県 × 2012-2023年度） # ---------------------------------------------------------------- df_b = pd.read_csv(DATA_B, encoding='cp932', header=1) # 都道府県レベルのみ（地域コードが R+5桁数字） df_b = df_b[df_b['地域コード'].str.match(r'^R\d{5}$', na=False)].copy() # 使用列（日本語名） needed_cols_raw = [ '総人口', # A1101 相当 '転出者数（日本人移動者）', # A5102 相当 '転入者数（日本人移動者）', # A5101 相当 '高等学校卒業者数', # E4601 相当 '高等学校卒業者のうち進学者数', # E4602 相当 '月間有効求人数（一般）', # F3103 相当 '月間有効求職者数（一般）', # F3102 相当 '消費支出（二人以上の世帯）', # L3221 相当 '標準価格（平均価格）（住宅地）', # C5401 相当 '65歳以上人口', # A1303 相当 '合計特殊出生率', # A4103 相当 ] for col in needed_cols_raw: df_b[col] = pd.to_numeric(df_b[col], errors='coerce') df_b['年度'] = pd.to_numeric(df_b['年度'], errors='coerce') # ---------------------------------------------------------------- # 派生変数の作成 # ---------------------------------------------------------------- df_b['転出率'] = df_b['転出者数（日本人移動者）'] / df_b['総人口'] * 100 df_b['net_mig率'] = (df_b['転入者数（日本人移動者）'] - df_b['転出者数（日本人移動者）']) / df_b['総人口'] * 100 df_b['大学進学率'] = df_b['高等学校卒業者のうち進学者数'] / df_b['高等学校卒業者数'] * 100 df_b['有効求人倍率'] = df_b['月間有効求人数（一般）'] / df_b['月間有効求職者数（一般）'] df_b['高齢化率'] = df_b['65歳以上人口'] / df_b['総人口'] * 100 # 欠損除外と必要列の抽出 model_cols = [ '転出率', 'net_mig率', '大学進学率', '有効求人倍率', '消費支出（二人以上の世帯）', '標準価格（平均価格）（住宅地）', '高齢化率', '合計特殊出生率' ] df_panel = df_b.dropna(subset=model_cols + ['都道府県', '年度']).copy() df_panel = df_panel.set_index(['都道府県', '年度']) n_prefs = df_panel.index.get_level_values(0).nunique() n_years = df_panel.index.get_level_values(1).nunique() years_range = sorted(df_panel.index.get_level_values(1).unique()) print(f"パネルデータ: {n_prefs} 都道府県 × {n_years} 年度（{years_range[0]}〜{years_range[-1]}）") print(f"総観測数: {len(df_panel)}") # ---------------------------------------------------------------- # 説明変数の定義（標準化用ラベル） # ---------------------------------------------------------------- PRED_COLS = [ '大学進学率', '有効求人倍率', '消費支出（二人以上の世帯）', '標準価格（平均価格）（住宅地）', '高齢化率', '合計特殊出生率', ] PRED_LABELS = [ '大学進学率\n（%）', '有効求人倍率', '消費支出\n（円/月）', '住宅地\n標準価格（円）', '高齢化率\n（%）', 'TFR', ] SHORT_LABELS = ['大学進学率', '有効求人倍率', '消費支出', '住宅地価格', '高齢化率', 'TFR'] # ---------------------------------------------------------------- # Panel FE 回帰（PanelOLS、entity effects + clustered SE） # ---------------------------------------------------------------- import statsmodels.api as sm def run_fe(dep_name): dep_var = df_panel[dep_name] exog_data = df_panel[PRED_COLS].copy() exog_data = sm.add_constant(exog_data) model = PanelOLS(dep_var, exog_data, entity_effects=True, time_effects=False) res = model.fit(cov_type='clustered', cluster_entity=True) return res res_out = run_fe('転出率') res_net = run_fe('net_mig率') def sig_star(p): if p < 0.001: return "***" if p < 0.01: return "**" if p < 0.05: return "*" if p < 0.10: return "†" return "n.s." def print_result(res, title): print(f"\n=== {title} ===") params = res.params pvals = res.pvalues ses = res.std_errors print(f"Within R² = {res.rsquared_within:.4f} " f"Between R² = {res.rsquared_between:.4f} " f"Overall R² = {res.rsquared:.4f}") print(f"{'変数':<28} {'係数':>10} {'SE':>8} {'p値':>8} {'有意':>5}") for nm in PRED_COLS + ['const']: if nm in params.index: c, s, p = params[nm], ses[nm], pvals[nm] print(f" {nm:<26} {c:>10.4f} {s:>8.4f} {p:>8.4f} {sig_star(p):>5}") print_result(res_out, "FE回帰 目的変数①: 転出率（%）") print_result(res_net, "FE回帰 目的変数②: 純移動率（%）") # ---------------------------------------------------------------- # Hausman 検定（FE vs RE の比較のための注記） # ---------------------------------------------------------------- # linearmodels の BetweenOLS で between estimates を取得 from linearmodels import RandomEffects re_out = RandomEffects(df_panel['転出率'], sm.add_constant(df_panel[PRED_COLS])).fit() re_net = RandomEffects(df_panel['net_mig率'], sm.add_constant(df_panel[PRED_COLS])).fit() print("\n--- Hausman 検定の参考情報 ---") print("FE Within R² (転出率):", f"{res_out.rsquared_within:.4f}") print("RE Overall R² (転出率):", f"{re_out.rsquared:.4f}") print("FE Within R² (純移動率):", f"{res_net.rsquared_within:.4f}") print("RE Overall R² (純移動率):", f"{re_net.rsquared:.4f}") print("→ FE係数とRE係数の差が大きければ、固定効果モデルが適切（Hausman検定の原理）") # 係数比較（教育用） fe_coefs = res_out.params[PRED_COLS].values re_coefs = re_out.params[PRED_COLS].values hausman_stat = float(np.sum((fe_coefs - re_coefs) ** 2)) print(f" 簡易Hausman距離（転出率モデル）: {hausman_stat:.4f}") print(" ※ 距離が大きいほどFEとREの推定値が乖離 → 固定効果モデルを採用") # ================================================================ # 図1: 転出率の時系列（全国平均 2012-2023） # ================================================================ df_b2 = df_b.copy() df_b2 = df_b2.dropna(subset=['転出率', '年度']) ts = df_b2.groupby('年度')['転出率'].mean() fig, ax = plt.subplots(figsize=(10, 5)) ax.plot(ts.index, ts.values, marker='o', lw=2.5, color='#1565C0', markersize=7, label='全国平均転出率（%）') ax.fill_between(ts.index, ts.values * 0.95, ts.values * 1.05, alpha=0.15, color='#1565C0') ax.axvline(2012, color='gray', lw=1, linestyle='--', alpha=0.5) ax.axvline(2023, color='gray', lw=1, linestyle='--', alpha=0.5) # 注目年のアノテーション peak_yr = ts.idxmax() ax.annotate(f'最大: {ts.max():.2f}%\n({peak_yr}年度)', xy=(peak_yr, ts.max()), xytext=(peak_yr + 1, ts.max() + 0.05), fontsize=9, color='#C62828', arrowprops=dict(arrowstyle='->', color='#C62828', lw=1.2)) ax.set_xlabel("年度", fontsize=11) ax.set_ylabel("転出率（転出者数/総人口 × 100, %）", fontsize=11) ax.set_title("図1: 47都道府県の転出率（全国平均）2012〜2023年度\n" "データ出典: SSDSE-B-2026（総務省 住民基本台帳人口移動報告）", fontsize=11) ax.legend(fontsize=10) ax.grid(True, alpha=0.3) ax.set_xticks(ts.index) plt.tight_layout() plt.savefig(FIG_DIR + "2023_U5_1_fig1_ts.png", dpi=150) plt.close() print("\nfig1 saved") # ================================================================ # 図2: FE 回帰係数（転出率モデル） # ================================================================ def coef_plot(res, pred_cols, short_labels, title, figpath): coefs = res.params[pred_cols].values ses = res.std_errors[pred_cols].values pvals = res.pvalues[pred_cols].values colors = ['#C62828' if c > 0 else '#1565C0' for c in coefs] alphas = [0.9 if sig_star(p) != 'n.s.' else 0.45 for p in pvals] fig, ax = plt.subplots(figsize=(9, 5)) y_pos = range(len(short_labels)) bars = ax.barh(y_pos, coefs, color=colors, alpha=0.8, xerr=1.96 * ses, capsize=5, error_kw=dict(ecolor='gray', lw=1.5)) ax.axvline(0, color='black', lw=1.5) for i, (c, se_i, p_v) in enumerate(zip(coefs, ses, pvals)): star = sig_star(p_v) offset = max(abs(c) * 0.03, 0.001) ha_pos = 'left' if c >= 0 else 'right' x_pos = c + (1.96 * se_i + offset) if c >= 0 else c - (1.96 * se_i + offset) color_s = '#C62828' if star not in ('n.s.',) else 'gray' ax.text(x_pos, i, star, va='center', ha=ha_pos, fontsize=12, color=color_s) ax.set_yticks(list(y_pos)) ax.set_yticklabels(short_labels, fontsize=10) ax.set_xlabel("FE推定係数（95% CI, クラスター標準誤差）", fontsize=10) ax.set_title(title, fontsize=11) ax.grid(True, axis='x', alpha=0.3) red_p = mpatches.Patch(color='#C62828', alpha=0.8, label='正の係数（転出増加）') blue_p = mpatches.Patch(color='#1565C0', alpha=0.8, label='負の係数（転出減少）') ax.legend(handles=[red_p, blue_p], fontsize=9, loc='lower right') plt.tight_layout() plt.savefig(figpath, dpi=150) plt.close() coef_plot( res_out, PRED_COLS, SHORT_LABELS, f"図2: 固定効果（FE）回帰係数 — 目的変数: 転出率（%）\n" f"Within R²={res_out.rsquared_within:.3f} " f"N={n_prefs}都道府県×{n_years}年度 " f"*** p<0.001 ** p<0.01 * p<0.05 † p<0.10 n.s. 有意でない", FIG_DIR + "2023_U5_1_fig2_fe_out.png" ) print("fig2 saved") # ================================================================ # 図3: FE 回帰係数（純移動率モデル） # ================================================================ coef_plot( res_net, PRED_COLS, SHORT_LABELS, f"図3: 固定効果（FE）回帰係数 — 目的変数: 純移動率（%）\n" f"Within R²={res_net.rsquared_within:.3f} " f"N={n_prefs}都道府県×{n_years}年度 " f"*** p<0.001 ** p<0.01 * p<0.05 † p<0.10 n.s. 有意でない", FIG_DIR + "2023_U5_1_fig3_fe_net.png" ) print("fig3 saved") # ================================================================ # 図4: 散布図（大学進学率 vs 転出率, 2022年度断面） # ================================================================ df_2022 = df_b[df_b['年度'] == 2022].dropna( subset=['大学進学率', '転出率', '都道府県'] ).copy() fig, ax = plt.subplots(figsize=(10, 7)) scatter = ax.scatter( df_2022['大学進学率'], df_2022['転出率'], c=df_2022['有効求人倍率'], cmap='RdYlBu_r', s=70, alpha=0.8, zorder=3, vmin=0.5, vmax=2.5 ) cb = plt.colorbar(scatter, ax=ax, pad=0.02) cb.set_label('有効求人倍率', fontsize=10) # 回帰直線 x_vals = df_2022['大学進学率'].values y_vals = df_2022['転出率'].values slope, intercept_lr, r_val, p_val, _ = stats.linregress(x_vals, y_vals) x_line = np.linspace(x_vals.min(), x_vals.max(), 200) ax.plot(x_line, intercept_lr + slope * x_line, color='#C62828', lw=2.5, label=f'回帰直線 (r={r_val:.2f}, p={p_val:.3f})') # 注目都道府県のラベル highlight = {'東京都', '神奈川県', '大阪府', '秋田県', '青森県', '沖縄県', '島根県', '鳥取県', '京都府', '愛知県'} for _, row in df_2022.iterrows(): pref = row['都道府県'] if pref in highlight: ax.annotate( pref, (row['大学進学率'], row['転出率']), xytext=(row['大学進学率'] + 0.3, row['転出率'] + 0.02), fontsize=8, arrowprops=dict(arrowstyle='->', color='gray', lw=0.8) ) ax.set_xlabel("大学進学率（高校卒業者のうち進学者の割合, %）", fontsize=11) ax.set_ylabel("転出率（転出者数/総人口 × 100, %）", fontsize=11) ax.set_title( "図4: 大学進学率 vs 転出率（47都道府県, 2022年度断面）\n" "点の色は有効求人倍率（青=低, 赤=高）。データ: SSDSE-B-2026", fontsize=11 ) ax.legend(fontsize=10) ax.grid(True, alpha=0.3) plt.tight_layout() plt.savefig(FIG_DIR + "2023_U5_1_fig4_scatter.png", dpi=150) plt.close() print("fig4 saved") print("\n=== 分析完了 ===") print("実データ（SSDSE-B-2026, 47都道府県 × 2012-2023年度）による固定効果パネル回帰") print(f"転出率モデル: Within R² = {res_out.rsquared_within:.4f}") print(f"純移動率モデル: Within R² = {res_net.rsquared_within:.4f}") print("図1〜4 を html/figures/ に保存しました。")