""" 教育用再現コード: 2021年 統計データ分析コンペティション 審査員奨励賞（大学生・一般の部） ================================================================================= 論文タイトル：地方からの人口流出：転出行動の決定要因パネル分析 受賞区分 ：審査員奨励賞 ［大学生・一般の部］ 【分析概要】 データ：SSDSE-B-2026.csv（都道府県別パネルデータ, 2012〜2023年度） 対象 ：全47都道府県 × 12年（2012〜2023） 本研究は「地方からの人口流出（転出行動）」の決定要因を 都道府県パネルデータにより分析する。 1. 転出率・純移動率の時系列推移（地域別） 2. 転出率の都道府県別ランキング（2021年断面） 3. PanelOLS 固定効果モデルによる転出率の決定要因推定 4. 求人倍率 vs 純移動率 散布図（2021年, 都道府県ラベル） 【変数定義】 転出率 = 転出者数（日本人移動者）/ 総人口 × 1000 純移動率 = (転入者数 - 転出者数) / 総人口 × 1000 求人倍率 = 月間有効求人数（一般）/ 月間有効求職者数（一般） 住宅地価_log = log(標準価格（平均価格）（住宅地）) 高齢化率 = 65歳以上人口 / 総人口 × 100 大学学生数_log = log(大学学生数) 【推定手法】 linearmodels PanelOLS: entity_effects=True, cov_type='clustered' Hausman検定（固定効果 vs ランダム効果）も実施 【使用データ】 SSDSE-B-2026.csv（統計数理研究所 教育用標準データセット） https://www.ism.ac.jp/editsection/ssdse/ 【図の出力】 html/figures/2021_U5_1_fig1_timeseries.png ... 転出率の時系列推移（地域別平均） html/figures/2021_U5_1_fig2_ranking.png ... 転出率都道府県別ランキング（2021年） html/figures/2021_U5_1_fig3_fe_coef.png ... PanelOLS FE 係数プロット html/figures/2021_U5_1_fig4_scatter.png ... 求人倍率 vs 純移動率 散布図（2021年） ================================================================================= """ # ============================================================ # 【データの準備】実行前に以下のデータファイルを用意してください # # 必要ファイル: # ・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/2021_U5_1_shorei.py # ============================================================ import os import warnings import numpy as np import pandas as pd import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt import matplotlib.patches as mpatches import statsmodels.api as sm from scipy import stats from linearmodels.panel import PanelOLS, RandomEffects warnings.filterwarnings('ignore') # ── パス設定 ────────────────────────────────────────────────────────────────── FIG_DIR = 'html/figures' DATA_B = 'data/raw/SSDSE-B-2026.csv' 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, }) # ── 地域区分マップ ───────────────────────────────────────────────────────────── REGION_MAP = { '北海道': '北海道・東北', '青森県': '北海道・東北', '岩手県': '北海道・東北', '宮城県': '北海道・東北', '秋田県': '北海道・東北', '山形県': '北海道・東北', '福島県': '北海道・東北', '茨城県': '関東', '栃木県': '関東', '群馬県': '関東', '埼玉県': '関東', '千葉県': '関東', '東京都': '関東', '神奈川県': '関東', '新潟県': '中部', '富山県': '中部', '石川県': '中部', '福井県': '中部', '山梨県': '中部', '長野県': '中部', '岐阜県': '中部', '静岡県': '中部', '愛知県': '中部', '三重県': '近畿', '滋賀県': '近畿', '京都府': '近畿', '大阪府': '近畿', '兵庫県': '近畿', '奈良県': '近畿', '和歌山県': '近畿', '鳥取県': '中国・四国', '島根県': '中国・四国', '岡山県': '中国・四国', '広島県': '中国・四国', '山口県': '中国・四国', '徳島県': '中国・四国', '香川県': '中国・四国', '愛媛県': '中国・四国', '高知県': '中国・四国', '福岡県': '九州・沖縄', '佐賀県': '九州・沖縄', '長崎県': '九州・沖縄', '熊本県': '九州・沖縄', '大分県': '九州・沖縄', '宮崎県': '九州・沖縄', '鹿児島県': '九州・沖縄', '沖縄県': '九州・沖縄', } REGION_COLORS = { '北海道・東北': '#4e9af1', '関東': '#e05c5c', '中部': '#f0a500', '近畿': '#5cb85c', '中国・四国': '#9b59b6', '九州・沖縄': '#f39c12', } # ═══════════════════════════════════════════════════════════════════════════════ # ■ Step1. 実データ読み込み（SSDSE-B-2026） # ═══════════════════════════════════════════════════════════════════════════════ print("=" * 65) print("■ SSDSE-B-2026 実データ読み込み") print("=" * 65) df_raw = pd.read_csv(DATA_B, encoding='cp932', header=1) # 都道府県レベルのみ（地域コード = R + 5桁数字） df = df_raw[df_raw['地域コード'].str.match(r'^R\d{5}$', na=False)].copy() df['年度'] = df['年度'].astype(int) df = df.sort_values(['都道府県', '年度']).reset_index(drop=True) print(f" 読み込み完了: {df['都道府県'].nunique()} 都道府県 × {df['年度'].nunique()} 年度") print(f" 年度範囲: {df['年度'].min()} 〜 {df['年度'].max()}") # ── 変数計算 ───────────────────────────────────────────────────────────────── pop = pd.to_numeric(df['総人口'], errors='coerce') migout = pd.to_numeric(df['転出者数（日本人移動者）'], errors='coerce') migin = pd.to_numeric(df['転入者数（日本人移動者）'], errors='coerce') eff_req = pd.to_numeric(df['月間有効求人数（一般）'], errors='coerce') eff_seek = pd.to_numeric(df['月間有効求職者数（一般）'], errors='coerce') land_res = pd.to_numeric(df['標準価格（平均価格）（住宅地）'], errors='coerce') pop65 = pd.to_numeric(df['65歳以上人口'], errors='coerce') univ_st = pd.to_numeric(df['大学学生数'], errors='coerce') df['転出率'] = migout / pop * 1000 df['純移動率'] = (migin - migout) / pop * 1000 df['求人倍率'] = eff_req / eff_seek.replace(0, np.nan) df['住宅地価_log'] = np.log(land_res.replace(0, np.nan)) df['高齢化率'] = pop65 / pop * 100 df['大学学生数_log'] = np.log(univ_st.replace(0, np.nan)) df['地域'] = df['都道府県'].map(REGION_MAP) # 都道府県省略名（2文字） df['都道府県2'] = df['都道府県'].str[:2] print(f"\n 変数計算後 欠損確認:") for v in ['転出率', '純移動率', '求人倍率', '住宅地価_log', '高齢化率', '大学学生数_log']: n_miss = df[v].isna().sum() print(f" {v:<16}: 欠損 {n_miss} 件") # ── 記述統計 ────────────────────────────────────────────────────────────────── print("\n 2021年断面の記述統計（転出率）:") d21 = df[df['年度'] == 2021].copy() print(f" 平均={d21['転出率'].mean():.2f}, SD={d21['転出率'].std():.2f}, " f"最小={d21['転出率'].min():.2f}, 最大={d21['転出率'].max():.2f}") # ═══════════════════════════════════════════════════════════════════════════════ # ■ Step2. PanelOLS 固定効果モデル # 被説明変数: 転出率 # 説明変数 : 求人倍率, 住宅地価_log, 高齢化率, 大学学生数_log # 固定効果 : 都道府県固定効果 (entity_effects=True) # クラスター標準誤差: 都道府県クラスター # ═══════════════════════════════════════════════════════════════════════════════ print("\n" + "=" * 65) print("■ PanelOLS 固定効果モデル推定") print("=" * 65) # パネルデータ用マルチインデックス設定 EXOG_VARS = ['求人倍率', '住宅地価_log', '高齢化率', '大学学生数_log'] ENDOG_VAR = '転出率' df_panel = df[['都道府県', '年度', ENDOG_VAR] + EXOG_VARS].dropna().copy() df_panel = df_panel.set_index(['都道府県', '年度']) y_fe = df_panel[ENDOG_VAR] X_fe = sm.add_constant(df_panel[EXOG_VARS]) X_fe.columns = ['const'] + EXOG_VARS # PanelOLS 固定効果 model_fe = PanelOLS(y_fe, X_fe, entity_effects=True) res_fe = model_fe.fit(cov_type='clustered', cluster_entity=True) print("\n 【固定効果モデル（FE）推定結果】") print(f" 観測数: {res_fe.nobs}, 都道府県数: {df_panel.index.get_level_values(0).nunique()}") print(f" Within R²: {res_fe.rsquared_within:.4f}") print() param_names = [v for v in EXOG_VARS] print(f" {'変数':<18} {'係数':>10} {'SE':>10} {'t':>8} {'p値':>8}") print(f" {'-'*56}") for vname in EXOG_VARS: coef = res_fe.params[vname] se = res_fe.std_errors[vname] tval = res_fe.tstats[vname] pval = res_fe.pvalues[vname] sig = '***' if pval < 0.01 else ('**' if pval < 0.05 else ('*' if pval < 0.1 else '')) print(f" {vname:<18} {coef:>10.4f} {se:>10.4f} {tval:>8.3f} {pval:>8.4f} {sig}") # ── Hausman 検定（FE vs RE）────────────────────────────────────────────────── print("\n" + "=" * 65) print("■ Hausman 検定（固定効果 vs ランダム効果の選択）") print("=" * 65) model_re = RandomEffects(y_fe, X_fe) res_re = model_re.fit(cov_type='robust') # Hausman 統計量の手計算 b_fe = res_fe.params[EXOG_VARS].values b_re = res_re.params[EXOG_VARS].values # FEの分散共分散行列のEXOG部分のみ V_fe = res_fe.cov.loc[EXOG_VARS, EXOG_VARS].values V_re = res_re.cov.loc[EXOG_VARS, EXOG_VARS].values diff = b_fe - b_re V_diff = V_fe - V_re try: # 正定値でない場合の対処: 対角要素のみで近似 V_diff_diag = np.diag(np.abs(np.diag(V_diff))) H_stat = float(diff @ np.linalg.pinv(V_diff_diag) @ diff) H_df = len(EXOG_VARS) H_pval = 1 - stats.chi2.cdf(H_stat, df=H_df) except Exception: H_stat, H_df, H_pval = np.nan, len(EXOG_VARS), np.nan print(f"\n Hausman 統計量 H = {H_stat:.4f}") print(f" 自由度 df = {H_df}") print(f" p値 = {H_pval:.4f}") if not np.isnan(H_pval): if H_pval < 0.05: print(" → p < 0.05: 固定効果モデル（FE）を支持") else: print(" → p >= 0.05: ランダム効果モデル（RE）も有力") # ═══════════════════════════════════════════════════════════════════════════════ # ■ Step3. 図の生成（4枚） # ═══════════════════════════════════════════════════════════════════════════════ # ──────────────────────────────────────────────────────────────────────────── # 図1：転出率の時系列推移（地域別平均, 2012-2023） # ──────────────────────────────────────────────────────────────────────────── print("\n図1: 転出率の時系列推移を作成中...") ts = (df.groupby(['年度', '地域'])['転出率'] .mean() .reset_index() .pivot(index='年度', columns='地域', values='転出率')) fig1, ax1 = plt.subplots(figsize=(11, 5.5)) for region, color in REGION_COLORS.items(): if region in ts.columns: ax1.plot(ts.index, ts[region], marker='o', markersize=4, linewidth=2, color=color, label=region) # 最終年ラベル ax1.annotate( region, xy=(ts.index[-1], ts[region].iloc[-1]), xytext=(3, 0), textcoords='offset points', fontsize=8, color=color, va='center' ) ax1.axvline(2021, color='gray', linestyle='--', linewidth=1.0, alpha=0.7, label='2021年（論文対象）') ax1.set_xlabel('年度', fontsize=12) ax1.set_ylabel('転出率（人口千対）', fontsize=12) ax1.set_title('転出率の時系列推移（地域別平均）\n（実データ：SSDSE-B-2026, 47都道府県）', fontsize=13, fontweight='bold') ax1.set_xticks(range(df['年度'].min(), df['年度'].max() + 1, 2)) ax1.legend(fontsize=8, loc='upper left', ncol=2) ax1.grid(True, alpha=0.25) plt.tight_layout() fig1.savefig(os.path.join(FIG_DIR, '2021_U5_1_fig1_timeseries.png'), bbox_inches='tight', dpi=150) plt.close(fig1) print(" → 2021_U5_1_fig1_timeseries.png 保存完了") # ──────────────────────────────────────────────────────────────────────────── # 図2：転出率の都道府県別ランキング（2021年） # ──────────────────────────────────────────────────────────────────────────── print("図2: 転出率都道府県別ランキングを作成中...") d21_sorted = d21.sort_values('転出率', ascending=True).copy() colors_bar = [REGION_COLORS.get(REGION_MAP.get(p, ''), '#aaa') for p in d21_sorted['都道府県']] fig2, ax2 = plt.subplots(figsize=(9, 11)) bars = ax2.barh( d21_sorted['都道府県'], d21_sorted['転出率'], color=colors_bar, alpha=0.85, edgecolor='white', height=0.75 ) national_mean = d21['転出率'].mean() ax2.axvline(national_mean, color='#333', linestyle='--', linewidth=1.3, label=f'全国平均 {national_mean:.1f}‰') ax2.set_xlabel('転出率（人口千対, ‰）', fontsize=11) ax2.set_title('都道府県別 転出率ランキング（2021年）\n（実データ：SSDSE-B-2026）', fontsize=12, fontweight='bold') # 凡例（地域別色） patches = [mpatches.Patch(color=c, label=r) for r, c in REGION_COLORS.items()] ax2.legend(handles=patches, fontsize=8, loc='lower right') ax2.grid(axis='x', alpha=0.25) # 値ラベル for bar, val in zip(bars, d21_sorted['転出率']): ax2.text(val + 0.05, bar.get_y() + bar.get_height() / 2, f'{val:.1f}', va='center', fontsize=7) plt.tight_layout() fig2.savefig(os.path.join(FIG_DIR, '2021_U5_1_fig2_ranking.png'), bbox_inches='tight', dpi=150) plt.close(fig2) print(" → 2021_U5_1_fig2_ranking.png 保存完了") # ──────────────────────────────────────────────────────────────────────────── # 図3：PanelOLS FE 係数プロット（転出率の決定要因） # ──────────────────────────────────────────────────────────────────────────── print("図3: FE 係数プロットを作成中...") coef_df = pd.DataFrame({ '変数': EXOG_VARS, '係数': [res_fe.params[v] for v in EXOG_VARS], 'SE': [res_fe.std_errors[v] for v in EXOG_VARS], 'p値': [res_fe.pvalues[v] for v in EXOG_VARS], }).sort_values('係数', ascending=True).reset_index(drop=True) coef_df['CI_lo'] = coef_df['係数'] - 1.96 * coef_df['SE'] coef_df['CI_hi'] = coef_df['係数'] + 1.96 * coef_df['SE'] coef_df['sig'] = coef_df['p値'] < 0.05 colors_fe = ['#C62828' if sig else '#90A4AE' for sig in coef_df['sig']] fig3, ax3 = plt.subplots(figsize=(9, 5)) ax3.barh(coef_df['変数'], coef_df['係数'], color=colors_fe, alpha=0.85, xerr=1.96 * coef_df['SE'], capsize=5, error_kw={'elinewidth': 1.5, 'ecolor': '#555'}, edgecolor='white', height=0.55) ax3.axvline(0, color='black', linewidth=1.2) for _, row in coef_df.iterrows(): label = f"{row['係数']:.4f}" if row['p値'] < 0.01: label += '***' elif row['p値'] < 0.05: label += '**' elif row['p値'] < 0.1: label += '*' offset = 0.0005 if row['係数'] >= 0 else -0.0005 ha = 'left' if row['係数'] >= 0 else 'right' ax3.text(row['係数'] + offset, row['変数'], label, va='center', fontsize=9, fontweight='bold', ha=ha) ax3.set_xlabel('回帰係数（転出率, 単位：人口千対）', fontsize=11) ax3.set_title('PanelOLS 固定効果モデル（FE）推定結果\n被説明変数：転出率 *** p<0.01, ** p<0.05, * p<0.1', fontsize=12, fontweight='bold') ax3.grid(axis='x', alpha=0.25) sig_patch = mpatches.Patch(color='#C62828', alpha=0.85, label='p<0.05（有意）') ns_patch = mpatches.Patch(color='#90A4AE', alpha=0.85, label='p≥0.05（非有意）') ax3.legend(handles=[sig_patch, ns_patch], fontsize=9, loc='lower right') plt.tight_layout() fig3.savefig(os.path.join(FIG_DIR, '2021_U5_1_fig3_fe_coef.png'), bbox_inches='tight', dpi=150) plt.close(fig3) print(" → 2021_U5_1_fig3_fe_coef.png 保存完了") # ──────────────────────────────────────────────────────────────────────────── # 図4：求人倍率 vs 純移動率 散布図（2021年, 都道府県ラベル） # ──────────────────────────────────────────────────────────────────────────── print("図4: 求人倍率 vs 純移動率 散布図を作成中...") d21_sc = d21[['都道府県', '都道府県2', '地域', '求人倍率', '純移動率']].dropna().copy() colors_sc = [REGION_COLORS.get(r, '#aaa') for r in d21_sc['地域']] fig4, ax4 = plt.subplots(figsize=(11, 7)) ax4.scatter(d21_sc['求人倍率'], d21_sc['純移動率'], c=colors_sc, s=70, alpha=0.85, zorder=3, edgecolors='white', linewidths=0.5) for _, row in d21_sc.iterrows(): ax4.annotate( row['都道府県2'], xy=(row['求人倍率'], row['純移動率']), xytext=(3, 3), textcoords='offset points', fontsize=7.5, color='#333', zorder=4 ) # 回帰直線 x_vals = d21_sc['求人倍率'].values y_vals = d21_sc['純移動率'].values coef_sc = np.polyfit(x_vals, y_vals, 1) x_range = np.linspace(x_vals.min(), x_vals.max(), 200) r_val, p_val = stats.pearsonr(x_vals, y_vals) ax4.plot(x_range, np.polyval(coef_sc, x_range), color='black', linewidth=1.5, linestyle='--', label=f'回帰直線 (r={r_val:.3f}, p={p_val:.3f})', zorder=2) ax4.axhline(0, color='#aaa', linewidth=1.0, linestyle='-') ax4.axvline(d21_sc['求人倍率'].mean(), color='#aaa', linewidth=1.0, linestyle=':', label=f'求人倍率 全国平均 ({d21_sc["求人倍率"].mean():.2f})') ax4.set_xlabel('求人倍率（月間有効求人数÷月間有効求職者数）', fontsize=12) ax4.set_ylabel('純移動率（人口千対, ‰）\n（転入率 − 転出率）', fontsize=12) ax4.set_title('求人倍率 vs 純移動率（2021年, 都道府県別）\n（実データ：SSDSE-B-2026）', fontsize=13, fontweight='bold') patches = [mpatches.Patch(color=c, label=r) for r, c in REGION_COLORS.items() if r in d21_sc['地域'].values] ax4.legend(handles=patches + [ plt.Line2D([0], [0], color='black', linestyle='--', label=f'回帰直線 (r={r_val:.3f})') ], fontsize=8, loc='upper left') ax4.grid(True, alpha=0.2) plt.tight_layout() fig4.savefig(os.path.join(FIG_DIR, '2021_U5_1_fig4_scatter.png'), bbox_inches='tight', dpi=150) plt.close(fig4) print(" → 2021_U5_1_fig4_scatter.png 保存完了") # ═══════════════════════════════════════════════════════════════════════════════ # ■ 最終サマリー # ═══════════════════════════════════════════════════════════════════════════════ print("\n" + "=" * 65) print("✓ 全図の生成完了（4枚）") print("=" * 65) print("\n【最終サマリー】") print(f" データ: SSDSE-B-2026 実公的データ") print(f" 47都道府県 × {df['年度'].nunique()} 年度 = {len(df_panel)} 有効観測（パネル）") print(f"\n 転出率（2021年）: 全国平均 {d21['転出率'].mean():.2f}‰") print(f" 最大: {d21.loc[d21['転出率'].idxmax(), '都道府県']} {d21['転出率'].max():.2f}‰") print(f" 最小: {d21.loc[d21['転出率'].idxmin(), '都道府県']} {d21['転出率'].min():.2f}‰") print(f"\n FE モデル Within R²: {res_fe.rsquared_within:.4f}") print(f"\n FE 推定結果（主要変数）:") for v in EXOG_VARS: coef = res_fe.params[v] pval = res_fe.pvalues[v] sig = '***' if pval < 0.01 else ('**' if pval < 0.05 else ('*' if pval < 0.1 else 'ns')) print(f" {v:<20}: {coef:>+.4f} ({sig})") print(f"\n 求人倍率 vs 純移動率 相関係数: r={r_val:.3f}, p={p_val:.4f}") print(f"\n Hausman 検定: H={H_stat:.4f}, p={H_pval:.4f}")