パネルデータを用いた進学と就職時による
人口流出の要因分析

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)}")

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

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

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

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/ に保存しました。")

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)

PRED_COLS = [
    '大学進学率',
    '有効求人倍率',
    '消費支出（二人以上の世帯）',
    '標準価格（平均価格）（住宅地）',
    '高齢化率',
    '合計特殊出生率',
]
PRED_LABELS = [
    '大学進学率\n（%）',
    '有効求人倍率',
    '消費支出\n（円/月）',
    '住宅地\n標準価格（円）',
    '高齢化率\n（%）',
    'TFR',
]
SHORT_LABELS = ['大学進学率', '有効求人倍率', '消費支出', '住宅地価格', '高齢化率', 'TFR']

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回帰 目的変数②: 純移動率（%）")

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の推定値が乖離 → 固定効果モデルを採用")