地方移住と地域活性化
転入増加要因の分析

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

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

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

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