地方創生政策が都道府県間人口移動に与えた影響

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()