コロナ禍が人口移動に与えた影響
転入・転出データのパネル分析

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)

# ── 地域マップ ─────────────────────────────────────────────────
region_map = {
    '北海道': '北海道・東北', '青森県': '北海道・東北', '岩手県': '北海道・東北',
    '宮城県': '北海道・東北', '秋田県': '北海道・東北', '山形県': '北海道・東北',
    '福島県': '北海道・東北', '茨城県': '関東', '栃木県': '関東', '群馬県': '関東',
    '埼玉県': '関東', '千葉県': '関東', '東京都': '関東', '神奈川県': '関東',
    '新潟県': '中部', '富山県': '中部', '石川県': '中部', '福井県': '中部',
    '山梨県': '中部', '長野県': '中部', '岐阜県': '中部', '静岡県': '中部', '愛知県': '中部',
    '三重県': '近畿', '滋賀県': '近畿', '京都府': '近畿', '大阪府': '近畿',
    '兵庫県': '近畿', '奈良県': '近畿', '和歌山県': '近畿',
    '鳥取県': '中国・四国', '島根県': '中国・四国', '岡山県': '中国・四国',
    '広島県': '中国・四国', '山口県': '中国・四国', '徳島県': '中国・四国',
    '香川県': '中国・四国', '愛媛県': '中国・四国', '高知県': '中国・四国',
    '福岡県': '九州・沖縄', '佐賀県': '九州・沖縄', '長崎県': '九州・沖縄',
    '熊本県': '九州・沖縄', '大分県': '九州・沖縄', '宮崎県': '九州・沖縄',
    '鹿児島県': '九州・沖縄', '沖縄県': '九州・沖縄'
}
region_colors = {
    '北海道・東北': '#4e9af1', '関東': '#e05c5c', '中部': '#f0a500',
    '近畿': '#5cb85c', '中国・四国': '#9b59b6', '九州・沖縄': '#f39c12'
}

# ── 変数作成 ──────────────────────────────────────────────────
df_b['転入率']   = df_b['転入者数（日本人移動者）'] / df_b['総人口'] * 1000
df_b['転出率']   = df_b['転出者数（日本人移動者）'] / df_b['総人口'] * 1000
df_b['純移動率'] = (df_b['転入者数（日本人移動者）'] - df_b['転出者数（日本人移動者）']) / df_b['総人口'] * 1000
df_b['COVID']    = df_b['年度'].isin([2020, 2021]).astype(int)
df_b['地域']     = df_b['都道府県'].map(region_map)

# ── 都市型・地方型の分類（2019年純移動率の上位10・下位10）──────
pivot_2019 = df_b[df_b['年度'] == 2019].set_index('都道府県')['純移動率']
top10_prefs  = pivot_2019.nlargest(10).index.tolist()
bot10_prefs  = pivot_2019.nsmallest(10).index.tolist()
df_b['都市型'] = df_b['都道府県'].apply(
    lambda x: '都市型（上位10）' if x in top10_prefs
    else ('地方型（下位10）' if x in bot10_prefs else 'その他')
)

# ── 統計値の計算（HTMLへの埋め込み用）────────────────────────
# 地域別平均純移動率の時系列
region_ts = df_b.groupby(['年度', '地域'])['純移動率'].mean().unstack()

# 都道府県別2022年純移動率
nmr_2022 = df_b[df_b['年度'] == 2022].set_index('都道府県')['純移動率'].sort_values(ascending=False)

# 2019→2020の変化量
nmr_2019 = df_b[df_b['年度'] == 2019].set_index('都道府県')['純移動率']
nmr_2020 = df_b[df_b['年度'] == 2020].set_index('都道府県')['純移動率']
delta_2019_2020 = (nmr_2020 - nmr_2019).dropna()

# 差の差分析（DiD）
# 都市型vs地方型 × COVID前後
did_df = df_b[df_b['都市型'] != 'その他'].copy()
did_df['post'] = (did_df['年度'] >= 2020).astype(int)
did_df['treat'] = (did_df['都市型'] == '都市型（上位10）').astype(int)
did_df['did'] = did_df['treat'] * did_df['post']
X_did = sm.add_constant(did_df[['treat', 'post', 'did']])
ols_did = sm.OLS(did_df['純移動率'], X_did).fit()

print("=== 差の差分析 結果 ===")
print(ols_did.summary())

# 重回帰分析（パネル固定効果なし簡易版）
reg_df = df_b.copy()
pref_dummies = pd.get_dummies(reg_df['都道府県'], drop_first=True, dtype=float)
X_reg = sm.add_constant(pd.concat([reg_df[['COVID', '転入率']], pref_dummies], axis=1))
ols_reg = sm.OLS(reg_df['純移動率'], X_reg).fit()

print("\n=== 重回帰分析（主要係数） ===")
print(f"COVID係数: {ols_reg.params['COVID']:.4f}, p={ols_reg.pvalues['COVID']:.4f}")
print(f"転入率係数: {ols_reg.params['転入率']:.4f}, p={ols_reg.pvalues['転入率']:.4f}")

# 東京都の純移動率推移
tokyo_ts = df_b[df_b['都道府県'] == '東京都'].set_index('年度')['純移動率']
print(f"\n=== 東京都の純移動率 ===\n{tokyo_ts.to_string()}")

# 2020年の全国平均純移動率と2019年の差
avg_2020 = df_b[df_b['年度'] == 2020]['純移動率'].mean()
avg_2019 = df_b[df_b['年度'] == 2019]['純移動率'].mean()
print(f"\n全国平均純移動率 2019: {avg_2019:.4f}, 2020: {avg_2020:.4f}, 差: {avg_2020-avg_2019:.4f}")

# DiD係数
did_coef = ols_did.params['did']
did_pval = ols_did.pvalues['did']
treat_coef = ols_did.params['treat']
post_coef  = ols_did.params['post']
print(f"\nDiD係数: {did_coef:.4f}, p={did_pval:.4f}")
print(f"treat係数: {treat_coef:.4f}, post係数: {post_coef:.4f}")

# 2019基準の変化トップ・ワースト
print(f"\n=== 2019→2020 純移動率変化量 Top5 ===")
print(delta_2019_2020.nlargest(5).to_string())
print(f"\n=== 2019→2020 純移動率変化量 Bottom5 ===")
print(delta_2019_2020.nsmallest(5).to_string())

# ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━

fig, ax = plt.subplots(figsize=(10, 5.5))

region_order = ['関東', '近畿', '中部', '北海道・東北', '中国・四国', '九州・沖縄']
line_styles   = ['-', '--', '-.', ':', '-', '--']
markers       = ['o', 's', '^', 'D', 'v', 'P']

for i, region in enumerate(region_order):
    if region in region_ts.columns:
        ax.plot(region_ts.index, region_ts[region],
                label=region, color=region_colors[region],
                linestyle=line_styles[i], marker=markers[i],
                markersize=5, linewidth=1.8)

# COVID期シェーディング
ax.axvspan(2019.5, 2021.5, color='#ffd700', alpha=0.18, label='COVID-19流行期（2020–2021年）')
ax.axhline(0, color='gray', linewidth=0.8, linestyle='--')

ax.set_xlabel('年度', fontsize=12)
ax.set_ylabel('純移動率（‰：人口1000人あたり）', fontsize=12)
ax.set_title('図1 地域ブロック別 純移動率の推移（2012〜2023年）', fontsize=13, fontweight='bold', pad=12)
ax.legend(loc='upper right', fontsize=10, framealpha=0.9)
ax.set_xticks(range(2012, 2024))
ax.set_xticklabels([str(y) for y in range(2012, 2024)], rotation=45, fontsize=9)
ax.grid(axis='y', alpha=0.3)
plt.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2020_H1_fig1.png'))
plt.close(fig)
print("fig1 saved")

# ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━

nmr_2022_sorted = nmr_2022.sort_values(ascending=True)
colors_fig2 = [region_colors[region_map.get(pref, '関東')]
               for pref in nmr_2022_sorted.index]

fig, ax = plt.subplots(figsize=(9, 11))
bars = ax.barh(range(len(nmr_2022_sorted)), nmr_2022_sorted.values,
               color=colors_fig2, edgecolor='white', linewidth=0.4)
ax.set_yticks(range(len(nmr_2022_sorted)))
ax.set_yticklabels(nmr_2022_sorted.index, fontsize=9.5)
ax.axvline(0, color='black', linewidth=0.8)
ax.set_xlabel('純移動率（‰：人口1000人あたり）', fontsize=11)
ax.set_title('図2 都道府県別 純移動率ランキング（2022年）\n（+は転入超過、−は転出超過）',
             fontsize=12, fontweight='bold', pad=10)

# 凡例
from matplotlib.patches import Patch
legend_elements = [Patch(facecolor=region_colors[r], label=r) for r in region_order]
ax.legend(handles=legend_elements, loc='lower right', fontsize=9, title='地域ブロック')

ax.grid(axis='x', alpha=0.3)
plt.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2020_H1_fig2.png'))
plt.close(fig)
print("fig2 saved")

# ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━

common_prefs = nmr_2019.index.intersection(nmr_2020.index)
x_vals = nmr_2019[common_prefs]
y_vals = nmr_2020[common_prefs]
delta_vals = (y_vals - x_vals)

fig, ax = plt.subplots(figsize=(10, 8))

for pref in common_prefs:
    region = region_map.get(pref, '関東')
    color  = region_colors.get(region, '#888888')
    ax.scatter(x_vals[pref], y_vals[pref], color=color, s=55, zorder=3, edgecolors='white', linewidth=0.5)
    offset_x, offset_y = 0.03, 0.02
    ax.annotate(pref.replace('県', '').replace('都', '').replace('府', '').replace('道', ''),
                (x_vals[pref], y_vals[pref]),
                xytext=(x_vals[pref] + offset_x, y_vals[pref] + offset_y),
                fontsize=6.5, color='#333333')

# 対角線（変化なし線）
lims = [min(x_vals.min(), y_vals.min()) - 0.5,
        max(x_vals.max(), y_vals.max()) + 0.5]
ax.plot(lims, lims, 'k--', linewidth=0.9, alpha=0.5, label='変化なし（y=x）')

ax.set_xlabel('2019年 純移動率（‰）', fontsize=12)
ax.set_ylabel('2020年 純移動率（‰）', fontsize=12)
ax.set_title('図3 COVID前後の純移動率変化（2019年→2020年）\n（対角線より下：転出超過が悪化、上：改善）',
             fontsize=12, fontweight='bold', pad=10)

legend_elements = [Patch(facecolor=region_colors[r], label=r) for r in region_order]
legend_elements.append(plt.Line2D([0], [0], color='k', linestyle='--', label='変化なし（y=x）'))
ax.legend(handles=legend_elements, loc='upper left', fontsize=9)

ax.axhline(0, color='gray', linewidth=0.6, linestyle=':')
ax.axvline(0, color='gray', linewidth=0.6, linestyle=':')
ax.grid(alpha=0.2)
plt.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2020_H1_fig3.png'))
plt.close(fig)
print("fig3 saved")

# ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━

fig, axes = plt.subplots(1, 2, figsize=(11, 5.5), sharey=False)

box_data = {}
for city_type in ['都市型（上位10）', '地方型（下位10）']:
    sub = did_df[did_df['都市型'] == city_type]
    box_data[city_type] = {
        'pre':  sub[sub['年度'] < 2020]['純移動率'].values,
        'post': sub[sub['年度'] >= 2020]['純移動率'].values
    }

colors_box = {'pre': '#90CAF9', 'post': '#EF9A9A'}
labels_box  = {'pre': 'COVID前\n(2012〜2019)', 'post': 'COVID後\n(2020〜2023)'}

for ax, city_type in zip(axes, ['都市型（上位10）', '地方型（下位10）']):
    bp = ax.boxplot(
        [box_data[city_type]['pre'], box_data[city_type]['post']],
        patch_artist=True,
        widths=0.45,
        medianprops=dict(color='black', linewidth=2),
        whiskerprops=dict(linewidth=1.2),
        capprops=dict(linewidth=1.2),
        flierprops=dict(marker='o', markerfacecolor='#666', markersize=4, alpha=0.5)
    )
    for patch, col in zip(bp['boxes'], [colors_box['pre'], colors_box['post']]):
        patch.set_facecolor(col)

    # t検定
    t_stat, t_pval = stats.ttest_ind(
        box_data[city_type]['pre'],
        box_data[city_type]['post']
    )
    sig_str = f"p = {t_pval:.3f}" + (" *" if t_pval < 0.05 else " (n.s.)")

    ax.set_xticks([1, 2])
    ax.set_xticklabels([labels_box['pre'], labels_box['post']], fontsize=10)
    ax.set_ylabel('純移動率（‰）', fontsize=11)
    ax.set_title(f'{city_type}\n（{sig_str}）', fontsize=11, fontweight='bold')
    ax.axhline(0, color='gray', linewidth=0.7, linestyle='--')
    ax.grid(axis='y', alpha=0.3)

fig.suptitle('図4 都市型・地方型 COVID前後の純移動率分布比較（箱ひげ図）',
             fontsize=12, fontweight='bold', y=1.01)
plt.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2020_H1_fig4.png'), bbox_inches='tight')
plt.close(fig)
print("fig4 saved")

# ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
# 統計サマリー出力
# ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━
print("\n=== 統計サマリー（HTML埋め込み用） ===")

tokyo_2019 = df_b[(df_b['都道府県'] == '東京都') & (df_b['年度'] == 2019)]['純移動率'].values[0]
tokyo_2020 = df_b[(df_b['都道府県'] == '東京都') & (df_b['年度'] == 2020)]['純移動率'].values[0]
tokyo_2021 = df_b[(df_b['都道府県'] == '東京都') & (df_b['年度'] == 2021)]['純移動率'].values[0]
tokyo_2022 = df_b[(df_b['都道府県'] == '東京都') & (df_b['年度'] == 2022)]['純移動率'].values[0]

print(f"東京都 純移動率: 2019={tokyo_2019:.3f}‰, 2020={tokyo_2020:.3f}‰, 2021={tokyo_2021:.3f}‰, 2022={tokyo_2022:.3f}‰")
print(f"東京 2019→2020 変化: {tokyo_2020-tokyo_2019:.3f}‰")

osaka_2020 = df_b[(df_b['都道府県'] == '大阪府') & (df_b['年度'] == 2020)]['純移動率'].values[0]
osaka_2019 = df_b[(df_b['都道府県'] == '大阪府') & (df_b['年度'] == 2019)]['純移動率'].values[0]
print(f"大阪府 純移動率: 2019={osaka_2019:.3f}‰, 2020={osaka_2020:.3f}‰, 変化={osaka_2020-osaka_2019:.3f}‰")

# DiD
print(f"\nDiD分析:")
print(f"  都市型平均 COVID前: {did_df[(did_df['都市型']=='都市型（上位10）')&(did_df['post']==0)]['純移動率'].mean():.3f}‰")
print(f"  都市型平均 COVID後: {did_df[(did_df['都市型']=='都市型（上位10）')&(did_df['post']==1)]['純移動率'].mean():.3f}‰")
print(f"  地方型平均 COVID前: {did_df[(did_df['都市型']=='地方型（下位10）')&(did_df['post']==0)]['純移動率'].mean():.3f}‰")
print(f"  地方型平均 COVID後: {did_df[(did_df['都市型']=='地方型（下位10）')&(did_df['post']==1)]['純移動率'].mean():.3f}‰")
print(f"  DiD係数={did_coef:.4f}, p={did_pval:.4f}")

# 2022年ランキング
print(f"\n2022年 純移動率 Top5:")
print(nmr_2022.head(5).to_string())
print(f"2022年 純移動率 Bottom5:")
print(nmr_2022.tail(5).to_string())

print("\n全図の生成完了。")