""" 教育用再現コード: 2022年 統計データ分析コンペティション 審査員奨励賞 [高校生の部] ================================================================= 論文タイトル：地方移住と地域活性化：転入増加要因の分析 受賞：審査員奨励賞（高校生の部） 【分析概要】 データ：SSDSE-B-2026.csv（都道府県別パネルデータ 2012-2023） 目的変数：転入率 = 転入者数（日本人移動者）/ 総人口 × 1000 説明変数：求人倍率、住宅地価_log、保育所密度、高齢化率、消費支出_log Fig1. 純移動率の時系列推移（地域別、東京圏の一人勝ちとCOVID後変化） Fig2. 求人倍率 vs 転入率 散布図（2022年、都道府県ラベル付き） Fig3. PanelOLS固定効果係数プロット（転入率の決定要因） Fig4. COVID前後（2019 vs 2022）の純移動率変化ランキング（上位・下位10都道府県） 【変数定義】 転入率 = 転入者数（日本人移動者）/ 総人口 × 1000 純移動率 = (転入者数 - 転出者数) / 総人口 × 1000 求人倍率 = 月間有効求人数（一般）/ 月間有効求職者数（一般） 住宅地価_log = log(標準価格（平均価格）（住宅地）) 保育所密度 = 保育所等数 / 総人口 × 10000 高齢化率 = 65歳以上人口 / 総人口 × 100 消費支出_log = log(消費支出（二人以上の世帯）) 【注意】 本スクリプトは実データ（SSDSE-B-2026）のみを使用。合成データは一切使用しない。 """ # ============================================================ # 【データの準備】実行前に以下のデータファイルを用意してください # # 必要ファイル: # ・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/2022_H5_13_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.font_manager as fm warnings.filterwarnings('ignore') # ── パス設定 ───────────────────────────────────────────────────────────────── FIG_DIR = 'html/figures' DATA_B = 'data/raw/SSDSE-B-2026.csv' os.makedirs(FIG_DIR, exist_ok=True) # ── フォント設定 ────────────────────────────────────────────────────────────── def get_japanese_font(): candidates = ['Hiragino Sans', 'Hiragino Kaku Gothic ProN', 'AppleGothic', 'Noto Sans CJK JP', 'IPAexGothic'] available = [f.name for f in fm.fontManager.ttflist] for c in candidates: if c in available: return c return 'sans-serif' FONT = get_japanese_font() plt.rcParams['font.family'] = FONT plt.rcParams['axes.unicode_minus'] = False # ── データ読み込み ──────────────────────────────────────────────────────────── df = pd.read_csv(DATA_B, encoding='cp932', header=0, skiprows=[1]) df.rename(columns={'SSDSE-B-2026': 'year'}, inplace=True) # 都道府県のみ（県コード R + 5桁 = 6文字、末尾000） df = df[df['Code'].str.startswith('R') & df['Code'].str.endswith('000')].copy() df['year'] = df['year'].astype(int) # ── 変数構築 ────────────────────────────────────────────────────────────────── df['転入率'] = df['A5101'] / df['A1101'] * 1000 df['純移動率'] = (df['A5101'] - df['A5102']) / df['A1101'] * 1000 df['求人倍率'] = df['F3103'] / df['F3102'] # 月間有効求人 / 求職 df['住宅地価_log'] = np.log(df['C5401'].replace(0, np.nan)) df['保育所密度'] = df['J2503'] / df['A1101'] * 10000 df['高齢化率'] = df['A1303'] / df['A1101'] * 100 df['消費支出_log'] = np.log(df['L3221'].replace(0, np.nan)) # ── 地域グループ定義 ────────────────────────────────────────────────────────── tokyo_area = ['東京都', '神奈川県', '埼玉県', '千葉県'] kinki_area = ['大阪府', '京都府', '兵庫県', '奈良県', '滋賀県', '和歌山県'] tohoku = ['青森県', '岩手県', '宮城県', '秋田県', '山形県', '福島県'] kyushu = ['福岡県', '佐賀県', '長崎県', '熊本県', '大分県', '宮崎県', '鹿児島県', '沖縄県'] def region_label(pref): if pref in tokyo_area: return '東京圏' elif pref in kinki_area: return '近畿圏' elif pref in tohoku: return '東北' elif pref in kyushu: return '九州・沖縄' else: return 'その他' df['region'] = df['Prefecture'].apply(region_label) REGION_COLORS = { '東京圏': '#E53935', '近畿圏': '#1E88E5', '東北': '#43A047', '九州・沖縄': '#FB8C00', 'その他': '#90A4AE', } print(f"データ行数: {len(df)}, 年度: {sorted(df['year'].unique())}") print(f"都道府県数（2022年）: {len(df[df['year']==2022])}") # ═══════════════════════════════════════════════════════════════════════════════ # Fig1: 純移動率の時系列推移（地域別） # ═══════════════════════════════════════════════════════════════════════════════ fig1, ax1 = plt.subplots(figsize=(10, 6), dpi=150) region_ts = df.groupby(['year', 'region'])['純移動率'].mean().reset_index() years_all = sorted(df['year'].unique()) for region, color in REGION_COLORS.items(): sub = region_ts[region_ts['region'] == region].sort_values('year') lw = 2.8 if region == '東京圏' else 1.8 ls = '-' if region in ['東京圏', '近畿圏', '東北'] else '--' ax1.plot(sub['year'], sub['純移動率'], color=color, lw=lw, ls=ls, marker='o', markersize=4, label=region) # COVID前後の境界線 ax1.axvline(2019.5, color='gray', lw=1.2, ls=':', alpha=0.7) ax1.axvline(2020.5, color='gray', lw=1.2, ls=':', alpha=0.7) ax1.axhline(0, color='black', lw=0.8, alpha=0.5) ax1.annotate('COVID-19\n感染拡大', xy=(2020, ax1.get_ylim()[1] if ax1.get_ylim()[1] > 0 else 2), xytext=(2020.3, 3.5), fontsize=9, color='gray', ha='left') ax1.set_title('純移動率の時系列推移（地域別平均）\n東京圏の一人勝ちとCOVID後の変化', fontsize=13, fontweight='bold', pad=12) ax1.set_xlabel('年度', fontsize=11) ax1.set_ylabel('純移動率（‰）', fontsize=11) ax1.legend(loc='upper left', fontsize=10, framealpha=0.9) ax1.grid(axis='y', alpha=0.3) ax1.set_xticks(years_all) ax1.tick_params(axis='x', rotation=45) plt.tight_layout() fig1_path = os.path.join(FIG_DIR, '2022_H5_13_fig1_timeseries.png') fig1.savefig(fig1_path, dpi=150, bbox_inches='tight') plt.close(fig1) print(f"Fig1 saved: {fig1_path}") # ═══════════════════════════════════════════════════════════════════════════════ # Fig2: 求人倍率 vs 転入率 散布図（2022年） # ═══════════════════════════════════════════════════════════════════════════════ df22 = df[df['year'] == 2022].dropna(subset=['求人倍率', '転入率']).copy() fig2, ax2 = plt.subplots(figsize=(10, 7), dpi=150) for region, color in REGION_COLORS.items(): sub = df22[df22['region'] == region] ax2.scatter(sub['求人倍率'], sub['転入率'], color=color, s=60, alpha=0.8, label=region, zorder=3) # 都道府県ラベル（省略なし — ただし重複を避けるため主要県のみ表示） label_prefs = ['東京都', '大阪府', '愛知県', '沖縄県', '秋田県', '青森県', '福岡県', '神奈川県', '埼玉県', '京都府', '宮城県', '北海道'] for _, row in df22.iterrows(): pref = row['Prefecture'] if pref in label_prefs: ax2.annotate(pref.replace('都','').replace('道','').replace('府','').replace('県',''), xy=(row['求人倍率'], row['転入率']), xytext=(4, 4), textcoords='offset points', fontsize=8, color='#333', zorder=4) # 回帰直線 from numpy.polynomial.polynomial import polyfit as npfit mask = df22[['求人倍率', '転入率']].notna().all(axis=1) x_v = df22.loc[mask, '求人倍率'].values y_v = df22.loc[mask, '転入率'].values coef = np.polyfit(x_v, y_v, 1) xline = np.linspace(x_v.min(), x_v.max(), 100) ax2.plot(xline, np.polyval(coef, xline), 'k--', lw=1.5, alpha=0.6, label='回帰直線') from scipy import stats as spstats r, p = spstats.pearsonr(x_v, y_v) ax2.text(0.97, 0.05, f'r = {r:.3f} (p = {p:.3f})', transform=ax2.transAxes, ha='right', va='bottom', fontsize=10, bbox=dict(boxstyle='round', facecolor='#EEE', alpha=0.8)) ax2.set_title('求人倍率 vs 転入率（2022年・47都道府県）', fontsize=13, fontweight='bold', pad=12) ax2.set_xlabel('求人倍率（月間有効求人 / 求職者数）', fontsize=11) ax2.set_ylabel('転入率（転入者数 / 総人口 × 1000）', fontsize=11) ax2.legend(loc='upper left', fontsize=9, framealpha=0.9) ax2.grid(alpha=0.3) plt.tight_layout() fig2_path = os.path.join(FIG_DIR, '2022_H5_13_fig2_scatter.png') fig2.savefig(fig2_path, dpi=150, bbox_inches='tight') plt.close(fig2) print(f"Fig2 saved: {fig2_path}") # ═══════════════════════════════════════════════════════════════════════════════ # Fig3: PanelOLS固定効果係数プロット（転入率の決定要因） # ═══════════════════════════════════════════════════════════════════════════════ # パネルデータ構築 panel_vars = ['転入率', '求人倍率', '住宅地価_log', '保育所密度', '高齢化率', '消費支出_log'] df_panel = df[['year', 'Prefecture'] + panel_vars].dropna().copy() var_labels = { '求人倍率': '求人倍率', '住宅地価_log': '住宅地価（log）', '保育所密度': '保育所密度\n（人口1万対）', '高齢化率': '高齢化率', '消費支出_log': '消費支出（log）', } exog_vars = list(var_labels.keys()) try: from linearmodels.panel import PanelOLS import statsmodels.api as sm df_panel = df_panel.set_index(['Prefecture', 'year']) endog = df_panel['転入率'] exog = sm.add_constant(df_panel[exog_vars]) model = PanelOLS(endog, exog, entity_effects=True, time_effects=False, drop_absorbed=True) result = model.fit(cov_type='clustered', cluster_entity=True) print(result.summary) coefs = result.params.drop('const', errors='ignore') se = result.std_errors.drop('const', errors='ignore') pvals = result.pvalues.drop('const', errors='ignore') ci_low = coefs - 1.96 * se ci_hi = coefs + 1.96 * se method_label = 'PanelOLS 固定効果（クラスター標準誤差）' except Exception as e: print(f"PanelOLS failed ({e}), falling back to OLS") import statsmodels.api as sm df_ols = df[panel_vars].dropna().copy() X = sm.add_constant(df_ols[exog_vars]) res_ols = sm.OLS(df_ols['転入率'], X).fit() coefs = res_ols.params.drop('const', errors='ignore') se = res_ols.bse.drop('const', errors='ignore') pvals = res_ols.pvalues.drop('const', errors='ignore') ci_low = coefs - 1.96 * se ci_hi = coefs + 1.96 * se method_label = 'OLS（固定効果なし）' # 標準化係数でなく生係数を表示（スケール差があるため正規化して可視化） # 視覚的比較のため z-score 標準化 df_std = df[panel_vars].dropna().copy() for v in exog_vars: df_std[v] = (df_std[v] - df_std[v].mean()) / df_std[v].std() import statsmodels.api as sm_std X_std = sm_std.add_constant(df_std[exog_vars]) res_std = sm_std.OLS(df_std['転入率'], X_std).fit() coefs_std = res_std.params.drop('const') se_std = res_std.bse.drop('const') pvals_std = res_std.pvalues.drop('const') ci_low_std = coefs_std - 1.96 * se_std ci_hi_std = coefs_std + 1.96 * se_std fig3, ax3 = plt.subplots(figsize=(8, 5), dpi=150) y_pos = np.arange(len(exog_vars)) colors_bar = ['#E53935' if p < 0.05 else '#90A4AE' for p in pvals_std] ax3.barh(y_pos, coefs_std.values, xerr=se_std.values, color=colors_bar, alpha=0.85, height=0.6, capsize=4) ax3.axvline(0, color='black', lw=1.0) ax3.set_yticks(y_pos) ax3.set_yticklabels([var_labels[v] for v in exog_vars], fontsize=10) ax3.set_xlabel('標準化回帰係数（β）\n赤：p < 0.05', fontsize=10) ax3.set_title(f'転入率の決定要因\n（{method_label}）', fontsize=12, fontweight='bold', pad=10) for i, (c, p, cl, ch) in enumerate(zip(coefs_std.values, pvals_std.values, ci_low_std.values, ci_hi_std.values)): sig = '***' if p < 0.001 else ('**' if p < 0.01 else ('*' if p < 0.05 else '')) label = f'β={c:.3f}{sig}' x_off = 0.01 if c >= 0 else -0.01 ha = 'left' if c >= 0 else 'right' ax3.text(c + x_off, i, label, va='center', ha=ha, fontsize=8.5) ax3.grid(axis='x', alpha=0.3) ax3.invert_yaxis() plt.tight_layout() fig3_path = os.path.join(FIG_DIR, '2022_H5_13_fig3_panel.png') fig3.savefig(fig3_path, dpi=150, bbox_inches='tight') plt.close(fig3) print(f"Fig3 saved: {fig3_path}") # ═══════════════════════════════════════════════════════════════════════════════ # Fig4: COVID前後（2019 vs 2022）の純移動率変化ランキング # ═══════════════════════════════════════════════════════════════════════════════ df19 = df[df['year'] == 2019][['Prefecture', '純移動率']].set_index('Prefecture') df22b = df[df['year'] == 2022][['Prefecture', '純移動率']].set_index('Prefecture') change = (df22b['純移動率'] - df19['純移動率']).dropna().sort_values(ascending=False) change.name = '純移動率変化（2022−2019）' # 上位10・下位10 top10 = change.head(10) bottom10 = change.tail(10) plot_data = pd.concat([top10, bottom10]).sort_values(ascending=True) fig4, ax4 = plt.subplots(figsize=(9, 8), dpi=150) bar_colors = ['#E53935' if v > 0 else '#1E88E5' for v in plot_data.values] bars = ax4.barh(plot_data.index, plot_data.values, color=bar_colors, alpha=0.85, height=0.7) ax4.axvline(0, color='black', lw=1.0) for bar, val in zip(bars, plot_data.values): x_off = 0.02 if val >= 0 else -0.02 ha = 'left' if val >= 0 else 'right' ax4.text(val + x_off, bar.get_y() + bar.get_height()/2, f'{val:+.2f}‰', va='center', ha=ha, fontsize=8.5) ax4.set_xlabel('純移動率の変化（2022年 − 2019年） [‰]', fontsize=10) ax4.set_title('COVID前後の純移動率変化\n上位10都道府県（赤）・下位10都道府県（青）', fontsize=12, fontweight='bold', pad=10) ax4.grid(axis='x', alpha=0.3) # 注釈 ax4.text(0.98, 0.98, '赤：流入増加\n青：流入減少 / 流出増加', transform=ax4.transAxes, ha='right', va='top', fontsize=9, bbox=dict(boxstyle='round', facecolor='#FFF9C4', alpha=0.9)) plt.tight_layout() fig4_path = os.path.join(FIG_DIR, '2022_H5_13_fig4_covid.png') fig4.savefig(fig4_path, dpi=150, bbox_inches='tight') plt.close(fig4) print(f"Fig4 saved: {fig4_path}") print("\n全図の生成が完了しました。") print(f" Fig1: {os.path.basename(fig1_path)}") print(f" Fig2: {os.path.basename(fig2_path)}") print(f" Fig3: {os.path.basename(fig3_path)}") print(f" Fig4: {os.path.basename(fig4_path)}")