""" 2018_U1_daijin.py 地方創生政策が都道府県間人口移動に与えた影響：差分の差分法による政策評価 2018年度 統計データ分析コンペティション 総務大臣賞（大学生・一般の部） 使用データ: 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/2018_U1_daijin.py # ============================================================ import os import numpy as np import pandas as pd import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt import statsmodels.api as sm from scipy import stats plt.rcParams['font.family'] = 'Hiragino Sans' plt.rcParams['axes.unicode_minus'] = False plt.rcParams['figure.dpi'] = 150 FIG_DIR = 'html/figures' DATA_B = 'data/raw/SSDSE-B-2026.csv' os.makedirs(FIG_DIR, exist_ok=True) # ── データ読み込み ────────────────────────────────────────────────────────── df_b = pd.read_csv(DATA_B, encoding='cp932', header=1) df_b = df_b[df_b['地域コード'].str.match(r'^R\d{5}', na=False)].copy() df_b['年度'] = df_b['年度'].astype(int) print("列名一覧:") print(df_b.columns.tolist()) print(f"\n年度範囲: {df_b['年度'].min()}〜{df_b['年度'].max()}") print(f"都道府県数: {df_b['都道府県'].nunique()}") # ── 都道府県名の正規化（県・都・道・府 を除去）──────────────────────────── def normalize_pref(name): for suffix in ['県', '都', '道', '府']: if name.endswith(suffix): return name[:-1] return name df_b['都道府県_短'] = df_b['都道府県'].apply(normalize_pref) # ── 地域マップ ────────────────────────────────────────────────────────────── region_map = { '北海道': '北海道・東北', '青森': '北海道・東北', '岩手': '北海道・東北', '宮城': '北海道・東北', '秋田': '北海道・東北', '山形': '北海道・東北', '福島': '北海道・東北', '茨城': '関東', '栃木': '関東', '群馬': '関東', '埼玉': '関東', '千葉': '関東', '東京': '関東', '神奈川': '関東', '新潟': '中部', '富山': '中部', '石川': '中部', '福井': '中部', '山梨': '中部', '長野': '中部', '岐阜': '中部', '静岡': '中部', '愛知': '中部', '三重': '近畿', '滋賀': '近畿', '京都': '近畿', '大阪': '近畿', '兵庫': '近畿', '奈良': '近畿', '和歌山': '近畿', '鳥取': '中国・四国', '島根': '中国・四国', '岡山': '中国・四国', '広島': '中国・四国', '山口': '中国・四国', '徳島': '中国・四国', '香川': '中国・四国', '愛媛': '中国・四国', '高知': '中国・四国', '福岡': '九州・沖縄', '佐賀': '九州・沖縄', '長崎': '九州・沖縄', '熊本': '九州・沖縄', '大分': '九州・沖縄', '宮崎': '九州・沖縄', '鹿児島': '九州・沖縄', '沖縄': '九州・沖縄' } region_colors = { '北海道・東北': '#4e9af1', '関東': '#e05c5c', '中部': '#f0a500', '近畿': '#5cb85c', '中国・四国': '#9b59b6', '九州・沖縄': '#f39c12' } df_b['地域'] = df_b['都道府県_短'].map(region_map) # ── アウトカム変数：転入超過率（転入者数 - 転出者数）/ 総人口 × 1000 ─────── df_b['転入超過数'] = df_b['転入者数（日本人移動者）'] - df_b['転出者数（日本人移動者）'] df_b['転入超過率'] = df_b['転入超過数'] / df_b['総人口'] * 1000 # 千人あたり # 高齢化率 df_b['高齢化率'] = df_b['65歳以上人口'] / df_b['総人口'] # 求人率（月間有効求人数 / 総人口 × 1万） df_b['求人率'] = df_b['月間有効求人数（一般）'] / df_b['総人口'] * 10000 # ── DiD設定 ──────────────────────────────────────────────────────────────── # 分析期間: 2012-2018（政策前: 2012-2014、政策後: 2015-2018） df_panel = df_b[df_b['年度'].between(2012, 2018)].copy() # 治療群の定義：2014年の総人口下位1/3（過疎地） # 対照群：総人口上位1/3（大都市圏） df_2014 = df_panel[df_panel['年度'] == 2014][['都道府県', '総人口']].copy() q33 = df_2014['総人口'].quantile(0.33) q67 = df_2014['総人口'].quantile(0.67) treat_prefs = df_2014[df_2014['総人口'] <= q33]['都道府県'].tolist() control_prefs = df_2014[df_2014['総人口'] >= q67]['都道府県'].tolist() print(f"\n治療群（小規模県, N={len(treat_prefs)}）: {[normalize_pref(p) for p in treat_prefs]}") print(f"対照群（大規模県, N={len(control_prefs)}）: {[normalize_pref(p) for p in control_prefs]}") # DiD用データ（治療群 + 対照群のみ） df_did = df_panel[df_panel['都道府県'].isin(treat_prefs + control_prefs)].copy() df_did['treat'] = df_did['都道府県'].isin(treat_prefs).astype(int) df_did['post'] = (df_did['年度'] >= 2015).astype(int) df_did['DiD'] = df_did['treat'] * df_did['post'] # ── 図1：転入超過率の時系列（治療群vs対照群） ──────────────────────────── fig, ax = plt.subplots(figsize=(10, 5.5)) treat_mean = df_did[df_did['treat']==1].groupby('年度')['転入超過率'].mean() ctrl_mean = df_did[df_did['treat']==0].groupby('年度')['転入超過率'].mean() ax.plot(treat_mean.index, treat_mean.values, color='#e05c5c', marker='o', linewidth=2.2, markersize=7, label='治療群（小規模県・過疎地）') ax.plot(ctrl_mean.index, ctrl_mean.values, color='#1565C0', marker='s', linewidth=2.2, markersize=7, label='対照群（大規模県・大都市圏）') ax.axvline(x=2014.5, color='#2e7d32', linewidth=2, linestyle='--', alpha=0.8) ax.text(2014.6, ax.get_ylim()[0] + (ax.get_ylim()[1]-ax.get_ylim()[0])*0.05, '地方創生政策\n開始（2015年）', color='#2e7d32', fontsize=10, va='bottom') ax.fill_between(treat_mean.index, df_did[df_did['treat']==1].groupby('年度')['転入超過率'].mean() - df_did[df_did['treat']==1].groupby('年度')['転入超過率'].std(), df_did[df_did['treat']==1].groupby('年度')['転入超過率'].mean() + df_did[df_did['treat']==1].groupby('年度')['転入超過率'].std(), alpha=0.12, color='#e05c5c') ax.set_xlabel('年度', fontsize=12) ax.set_ylabel('転入超過率（‰, 千人あたり）', fontsize=12) ax.set_title('図1：転入超過率の推移（治療群 vs 対照群）\n2015年の地方創生政策前後の比較', fontsize=13, fontweight='bold') ax.legend(fontsize=11, loc='upper right') ax.grid(axis='y', alpha=0.3) ax.set_xticks(range(2012, 2019)) plt.tight_layout() plt.savefig(os.path.join(FIG_DIR, '2018_U1_fig1.png'), dpi=150, bbox_inches='tight') plt.close() print("fig1 saved") # ── 図2：平行トレンド確認（政策前2012-2014） ───────────────────────────── fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5)) # 左：pre-period トレンド df_pre = df_did[df_did['年度'] <= 2014].copy() treat_pre = df_pre[df_pre['treat']==1].groupby('年度')['転入超過率'].mean() ctrl_pre = df_pre[df_pre['treat']==0].groupby('年度')['転入超過率'].mean() ax1.plot(treat_pre.index, treat_pre.values, color='#e05c5c', marker='o', linewidth=2.2, markersize=8, label='治療群') ax1.plot(ctrl_pre.index, ctrl_pre.values, color='#1565C0', marker='s', linewidth=2.2, markersize=8, label='対照群') # トレンド線 for grp, col, lab in [(1, '#e05c5c', '治療群'), (0, '#1565C0', '対照群')]: y_vals = df_pre[df_pre['treat']==grp].groupby('年度')['転入超過率'].mean().values x_vals = np.array([2012, 2013, 2014]) slope, intercept, r, p_val, se = stats.linregress(x_vals, y_vals) x_line = np.array([2011.8, 2014.2]) ax1.plot(x_line, intercept + slope * x_line, linestyle=':', color=col, alpha=0.6, linewidth=1.5) ax1.set_title('（A）政策前トレンド（2012-2014年）\n平行トレンド仮定の確認', fontsize=12, fontweight='bold') ax1.set_xlabel('年度', fontsize=11) ax1.set_ylabel('転入超過率（‰）', fontsize=11) ax1.legend(fontsize=10) ax1.grid(axis='y', alpha=0.3) ax1.set_xticks([2012, 2013, 2014]) # 右：DiDグラフ（2×2） categories = ['政策前（2012-2014）', '政策後（2015-2018）'] treat_vals = [ df_did[(df_did['treat']==1) & (df_did['post']==0)]['転入超過率'].mean(), df_did[(df_did['treat']==1) & (df_did['post']==1)]['転入超過率'].mean() ] ctrl_vals = [ df_did[(df_did['treat']==0) & (df_did['post']==0)]['転入超過率'].mean(), df_did[(df_did['treat']==0) & (df_did['post']==1)]['転入超過率'].mean() ] x = np.array([0, 1]) ax2.plot(x, treat_vals, color='#e05c5c', marker='o', linewidth=2.5, markersize=10, label='治療群（過疎地）') ax2.plot(x, ctrl_vals, color='#1565C0', marker='s', linewidth=2.5, markersize=10, label='対照群（大都市圏）') # 反事実（counterfactual） ctrl_change = ctrl_vals[1] - ctrl_vals[0] cf_val = treat_vals[0] + ctrl_change ax2.plot([0, 1], [treat_vals[0], cf_val], color='#e05c5c', linewidth=1.5, linestyle='--', alpha=0.5, label='治療群の反事実（推定）') # DiD矢印 ax2.annotate('', xy=(1, treat_vals[1]), xytext=(1, cf_val), arrowprops=dict(arrowstyle='<->', color='#2e7d32', lw=2)) did_val = treat_vals[1] - cf_val ax2.text(1.05, (treat_vals[1] + cf_val) / 2, f'DiD効果\n{did_val:+.2f}‰', color='#2e7d32', fontsize=10, va='center') ax2.set_xticks([0, 1]) ax2.set_xticklabels(['政策前', '政策後'], fontsize=11) ax2.set_ylabel('転入超過率（‰）', fontsize=11) ax2.set_title('（B）DiDの直感的理解（2×2表）\n反事実との差がDiD推定量', fontsize=12, fontweight='bold') ax2.legend(fontsize=9, loc='lower left') ax2.grid(axis='y', alpha=0.3) plt.suptitle('図2：平行トレンド確認とDiDの概念図', fontsize=13, fontweight='bold', y=1.01) plt.tight_layout() plt.savefig(os.path.join(FIG_DIR, '2018_U1_fig2.png'), dpi=150, bbox_inches='tight') plt.close() print("fig2 saved") # ── パネル固定効果モデル ──────────────────────────────────────────────────── y_col = '転入超過率' x_cols = ['DiD', 'treat', 'post', '高齢化率', '求人率'] df_fe = df_did.dropna(subset=[y_col] + x_cols).copy() # 固定効果：都道府県ダミー try: from linearmodels.panel import PanelOLS df_panel_idx = df_fe.set_index(['都道府県', '年度']) mod = PanelOLS(df_panel_idx[y_col], df_panel_idx[x_cols], entity_effects=True) res = mod.fit(cov_type='clustered', cluster_entity=True) print("\n=== Fixed Effects Model (linearmodels) ===") print(res.summary) use_lm = True except Exception as e: print(f"linearmodels not available: {e} → fallback OLS with dummies") use_lm = False # フォールバック: OLS with prefecture dummies dummies = pd.get_dummies(df_fe['都道府県'], drop_first=True) X_fe = sm.add_constant(pd.concat([df_fe[x_cols], dummies], axis=1).astype(float)) res_ols = sm.OLS(df_fe[y_col].astype(float), X_fe).fit(cov_type='HC3') print("\n=== FE via OLS+Dummies ===") print(res_ols.summary().tables[1]) # 係数と95%CI（DiD, treat, post, 高齢化率, 求人率） if use_lm: coef = res.params[x_cols] ci_low = res.conf_int()['lower'][x_cols] ci_high = res.conf_int()['upper'][x_cols] pvals = res.pvalues[x_cols] else: coef = res_ols.params[x_cols] ci_low = res_ols.conf_int()[0][x_cols] ci_high = res_ols.conf_int()[1][x_cols] pvals = res_ols.pvalues[x_cols] # ── 図3：固定効果モデル係数プロット ────────────────────────────────────── var_labels = { 'DiD': 'DiD効果\n（地方創生×過疎地）', 'treat': '治療群ダミー\n（過疎地）', 'post': '政策後ダミー\n（2015年〜）', '高齢化率': '高齢化率', '求人率': '求人率\n（月間有効求人数/人口）' } fig, ax = plt.subplots(figsize=(9, 5.5)) colors = [] for var in x_cols: p = pvals[var] if p < 0.05: colors.append('#C62828' if coef[var] < 0 else '#1565C0') else: colors.append('#9E9E9E') y_pos = np.arange(len(x_cols)) ax.barh(y_pos, coef[x_cols].values, xerr=[coef[x_cols].values - ci_low[x_cols].values, ci_high[x_cols].values - coef[x_cols].values], color=colors, alpha=0.85, height=0.55, error_kw=dict(elinewidth=1.8, capsize=5, ecolor='#333')) ax.axvline(0, color='black', linewidth=1.2, linestyle='-') ax.set_yticks(y_pos) ax.set_yticklabels([var_labels[v] for v in x_cols], fontsize=11) ax.set_xlabel('係数（転入超過率への影響, ‰）', fontsize=11) ax.set_title('図3：固定効果モデルの係数（95%信頼区間付き）\n' + '赤：有意な負効果　青：有意な正効果　灰：非有意', fontsize=12, fontweight='bold') # p値注釈 for i, var in enumerate(x_cols): p = pvals[var] star = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else '' if star: x_annot = ci_high[var] + abs(ci_high[var] - ci_low[var]) * 0.05 ax.text(x_annot, i, star, va='center', fontsize=12, color='#333', fontweight='bold') # DiD係数を強調 did_idx = x_cols.index('DiD') ax.get_children()[did_idx].set_edgecolor('#FFD700') ax.get_children()[did_idx].set_linewidth(2.5) ax.grid(axis='x', alpha=0.3) plt.tight_layout() plt.savefig(os.path.join(FIG_DIR, '2018_U1_fig3.png'), dpi=150, bbox_inches='tight') plt.close() print("fig3 saved") # ── 図4：地域別 転入超過率の箱ひげ図（最新年） ─────────────────────────── latest_year = df_b['年度'].max() df_latest = df_b[df_b['年度'] == latest_year].copy() df_latest['地域'] = df_latest['都道府県_短'].map(region_map) region_order = ['北海道・東北', '関東', '中部', '近畿', '中国・四国', '九州・沖縄'] fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 5.5)) # 左：転入超過率の箱ひげ図 data_by_region = [ df_latest[df_latest['地域'] == r]['転入超過率'].dropna().values for r in region_order ] bp = ax1.boxplot(data_by_region, patch_artist=True, vert=True, medianprops=dict(color='black', linewidth=2)) for patch, region in zip(bp['boxes'], region_order): patch.set_facecolor(region_colors[region]) patch.set_alpha(0.75) ax1.set_xticklabels(region_order, fontsize=8.5, rotation=20, ha='right') ax1.axhline(0, color='#999', linewidth=1, linestyle='--') ax1.set_ylabel('転入超過率（‰）', fontsize=11) ax1.set_title(f'（A）地域別 転入超過率（{latest_year}年）', fontsize=11, fontweight='bold') ax1.grid(axis='y', alpha=0.3) # 右：総人口の箱ひげ図 data_pop = [ df_latest[df_latest['地域'] == r]['総人口'].dropna().values / 10000 for r in region_order ] bp2 = ax2.boxplot(data_pop, patch_artist=True, vert=True, medianprops=dict(color='black', linewidth=2)) for patch, region in zip(bp2['boxes'], region_order): patch.set_facecolor(region_colors[region]) patch.set_alpha(0.75) ax2.set_xticklabels(region_order, fontsize=8.5, rotation=20, ha='right') ax2.set_ylabel('総人口（万人）', fontsize=11) ax2.set_title(f'（B）地域別 総人口（{latest_year}年）', fontsize=11, fontweight='bold') ax2.grid(axis='y', alpha=0.3) plt.suptitle(f'図4：地域別（6地域）人口・転入移動指標の箱ひげ図（{latest_year}年）', fontsize=12, fontweight='bold', y=1.01) # 凡例 handles = [plt.Rectangle((0,0),1,1, facecolor=region_colors[r], alpha=0.75) for r in region_order] fig.legend(handles, region_order, loc='lower center', ncol=3, bbox_to_anchor=(0.5, -0.12), fontsize=9) plt.tight_layout() plt.savefig(os.path.join(FIG_DIR, '2018_U1_fig4.png'), dpi=150, bbox_inches='tight') plt.close() print("fig4 saved") # ── Hausman検定（簡易）──────────────────────────────────────────────────── print("\n=== Hausman検定（FE vs RE）===") try: from linearmodels.panel import RandomEffects df_panel_idx2 = df_fe.set_index(['都道府県', '年度']) mod_re = RandomEffects(df_panel_idx2[y_col], df_panel_idx2[x_cols]) res_re = mod_re.fit(cov_type='robust') # Hausman統計量 b_fe = res.params[x_cols].values b_re = res_re.params[x_cols].values V_fe = res.cov.values[:len(x_cols), :len(x_cols)] V_re = res_re.cov.values[:len(x_cols), :len(x_cols)] diff = b_fe - b_re V_diff = V_fe - V_re try: hausman_stat = diff @ np.linalg.inv(V_diff) @ diff df_h = len(x_cols) p_hausman = 1 - stats.chi2.cdf(hausman_stat, df_h) print(f"H統計量 = {hausman_stat:.4f}, 自由度 = {df_h}, p値 = {p_hausman:.4f}") if p_hausman < 0.05: print("→ H0棄却: 固定効果モデル（FE）を採用") else: print("→ H0非棄却: 変量効果モデル（RE）が適切") except np.linalg.LinAlgError: print("Hausman: V_diff が特異行列のため簡易版で計算") hausman_stat = float(np.dot(diff, diff) / (np.trace(V_fe) - np.trace(V_re) + 1e-8)) print(f"簡易 H = {hausman_stat:.4f}") except Exception as e_h: print(f"Hausman検定スキップ: {e_h}") # ── 記述統計 ───────────────────────────────────────────────────────────── print("\n=== 記述統計（転入超過率, 2012-2018） ===") desc = df_did.groupby(['treat', 'post'])['転入超過率'].agg(['mean', 'std', 'count']).round(3) desc.index = desc.index.map(lambda x: ('治療群' if x[0]==1 else '対照群', '政策後' if x[1]==1 else '政策前')) print(desc) print("\nDONE: 2018_U1_daijin")