人口減少対策の効果検証
移住促進・子育て支援政策の地域分析

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'
}

# ── 変数の準備 ──────────────────────────────────────────────────────────────
# 2012年・2022年断面
df12 = df_b[df_b['年度'] == 2012].set_index('都道府県').copy()
df22 = df_b[df_b['年度'] == 2022].set_index('都道府県').copy()

# 共通47都道府県
pref_47 = df22.index.intersection(df12.index)
df22 = df22.loc[pref_47].copy()
df12 = df12.loc[pref_47].copy()

# 人口増加率（2012→2022）
df22['人口増加率'] = (df22['総人口'] - df12['総人口']) / df12['総人口'] * 100

# 人口千人あたり指標
df22['保育所密度']  = df22['保育所等数'] / df22['総人口'] * 1000
df22['転入率']     = df22['転入者数（日本人移動者）'] / df22['総人口'] * 100
df22['転出率']     = df22['転出者数（日本人移動者）'] / df22['総人口'] * 100
df22['純移動率']   = df22['転入率'] - df22['転出率']
df22['住居費割合'] = df22['住居費（二人以上の世帯）'] / df22['消費支出（二人以上の世帯）'] * 100

# 地域ラベル付与
df22['地域'] = df22.index.map(region_map)
df22 = df22.reset_index()

# ── 都道府県略称 ──────────────────────────────────────────────────────────────
abbr = {
    '北海道':'北海道', '青森県':'青森', '岩手県':'岩手', '宮城県':'宮城', '秋田県':'秋田',
    '山形県':'山形', '福島県':'福島', '茨城県':'茨城', '栃木県':'栃木', '群馬県':'群馬',
    '埼玉県':'埼玉', '千葉県':'千葉', '東京都':'東京', '神奈川県':'神奈川',
    '新潟県':'新潟', '富山県':'富山', '石川県':'石川', '福井県':'福井',
    '山梨県':'山梨', '長野県':'長野', '岐阜県':'岐阜', '静岡県':'静岡', '愛知県':'愛知',
    '三重県':'三重', '滋賀県':'滋賀', '京都府':'京都', '大阪府':'大阪',
    '兵庫県':'兵庫', '奈良県':'奈良', '和歌山県':'和歌山',
    '鳥取県':'鳥取', '島根県':'島根', '岡山県':'岡山', '広島県':'広島', '山口県':'山口',
    '徳島県':'徳島', '香川県':'香川', '愛媛県':'愛媛', '高知県':'高知',
    '福岡県':'福岡', '佐賀県':'佐賀', '長崎県':'長崎', '熊本県':'熊本',
    '大分県':'大分', '宮崎県':'宮崎', '鹿児島県':'鹿児島', '沖縄県':'沖縄'
}

fig1, ax1 = plt.subplots(figsize=(10, 7))

x_col = '保育所密度'
y_col = '合計特殊出生率'
df_plot = df22[['都道府県', x_col, y_col, '地域']].dropna()

for region, grp in df_plot.groupby('地域'):
    col = region_colors.get(region, '#888')
    ax1.scatter(grp[x_col], grp[y_col], c=col, s=60, alpha=0.85,
                label=region, edgecolors='white', linewidths=0.5, zorder=3)
    for _, row in grp.iterrows():
        ax1.annotate(abbr.get(row['都道府県'], row['都道府県']),
                     (row[x_col], row[y_col]),
                     fontsize=7.5, ha='center', va='bottom', xytext=(0, 3),
                     textcoords='offset points', color='#333', zorder=4)

# 回帰線
x_arr = df_plot[x_col].values
y_arr = df_plot[y_col].values
slope, intercept, r_val, p_val, _ = stats.linregress(x_arr, y_arr)
x_line = np.linspace(x_arr.min(), x_arr.max(), 100)
ax1.plot(x_line, intercept + slope * x_line, 'k--', linewidth=1.5, alpha=0.7, zorder=2)
ax1.text(0.98, 0.04,
         f'r = {r_val:.3f}  (p = {p_val:.3f})',
         transform=ax1.transAxes, ha='right', fontsize=10,
         bbox=dict(boxstyle='round,pad=0.4', fc='white', ec='#ccc', alpha=0.9))

ax1.set_xlabel('保育所密度（施設数/千人）', fontsize=12)
ax1.set_ylabel('合計特殊出生率', fontsize=12)
ax1.set_title('図1　保育所密度と合計特殊出生率の関係（2022年・47都道府県）', fontsize=13, fontweight='bold')
ax1.legend(loc='upper left', fontsize=9, framealpha=0.9)
ax1.grid(True, alpha=0.3, zorder=1)
fig1.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2021_H5_1_fig1.png'), dpi=150, bbox_inches='tight')
plt.close(fig1)
print('fig1 saved')

years_avail = sorted(df_b['年度'].unique())

region_pop = {}
for yr in years_avail:
    dfyr = df_b[df_b['年度'] == yr].copy()
    dfyr['地域'] = dfyr['都道府県'].map(region_map)
    grp = dfyr.groupby('地域')['総人口'].sum()
    region_pop[yr] = grp

pop_df = pd.DataFrame(region_pop).T
# 2012年を100として指数化
base = pop_df.loc[2012]
pop_idx = pop_df.div(base) * 100

fig2, ax2 = plt.subplots(figsize=(10, 6))
for region in region_colors:
    if region in pop_idx.columns:
        col = region_colors[region]
        ax2.plot(pop_idx.index, pop_idx[region], color=col,
                 linewidth=2.2, marker='o', markersize=4, label=region, zorder=3)

ax2.axhline(100, color='#aaa', linestyle='--', linewidth=1, alpha=0.7)
ax2.set_xlabel('年度', fontsize=12)
ax2.set_ylabel('人口指数（2012年=100）', fontsize=12)
ax2.set_title('図2　地域別人口推移（2012〜2023年、2012年=100）', fontsize=13, fontweight='bold')
ax2.set_xticks(years_avail)
ax2.set_xticklabels([str(y) for y in years_avail], rotation=45, fontsize=9)
ax2.legend(loc='lower left', fontsize=9, framealpha=0.9)
ax2.grid(True, alpha=0.3, zorder=1)
fig2.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2021_H5_1_fig2.png'), dpi=150, bbox_inches='tight')
plt.close(fig2)
print('fig2 saved')

fig3, ax3 = plt.subplots(figsize=(9, 5))

var_labels = {
    '保育所密度':    '保育所密度\n（施設数/千人）',
    '合計特殊出生率': '合計特殊出生率',
    '転入率':        '転入率\n（転入数/人口 %）',
    '転出率':        '転出率\n（転出数/人口 %）',
    '住居費割合':    '住居費割合\n（住居費/消費支出 %）',
}
labels = [var_labels.get(v, v) for v in coef_df['変数']]
colors_bar = ['#e05c5c' if sig else '#aaaaaa' for sig in coef_df['有意']]
values = coef_df['標準化係数'].values

bars = ax3.barh(labels, values, color=colors_bar, edgecolor='white', height=0.6)
ax3.axvline(0, color='#333', linewidth=1.2)

# p値ラベル
for i, (bar, row) in enumerate(zip(bars, coef_df.itertuples())):
    sig_str = '***' if row.p値 < 0.001 else '**' if row.p値 < 0.01 else '*' if row.p値 < 0.05 else 'n.s.'
    x_pos = row.標準化係数
    ha = 'left' if x_pos >= 0 else 'right'
    offset = 0.01 if x_pos >= 0 else -0.01
    ax3.text(x_pos + offset, i, f'  {sig_str}  (p={row.p値:.3f})',
             va='center', ha=ha, fontsize=9)

# 凡例パッチ
import matplotlib.patches as mpatches
sig_patch = mpatches.Patch(color='#e05c5c', label='有意（p<0.05）')
ns_patch  = mpatches.Patch(color='#aaaaaa', label='非有意（n.s.）')
ax3.legend(handles=[sig_patch, ns_patch], loc='lower right', fontsize=9)

ax3.set_xlabel('標準化偏回帰係数 (β)', fontsize=12)
ax3.set_title(f'図3　重回帰分析：標準化偏回帰係数（目的変数：人口増加率, N={len(df_ols)}, Adj-R²={r2_adj:.3f}）',
              fontsize=11, fontweight='bold')
ax3.grid(True, axis='x', alpha=0.3)
fig3.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2021_H5_1_fig3.png'), dpi=150, bbox_inches='tight')
plt.close(fig3)
print('fig3 saved')

reg_cols = ['保育所密度', '合計特殊出生率', '転入率', '転出率', '住居費割合']
y_col_ols = '人口増加率'

df_ols = df22[[y_col_ols] + reg_cols].dropna().copy()

# 標準化
df_std = (df_ols - df_ols.mean()) / df_ols.std()
Y = df_std[y_col_ols]
X = sm.add_constant(df_std[reg_cols])
model = sm.OLS(Y, X).fit()

print('\n=== OLS 重回帰結果（標準化） ===')
print(model.summary())

coef_df = pd.DataFrame({
    '変数': reg_cols,
    '標準化係数': model.params[reg_cols].values,
    'p値': model.pvalues[reg_cols].values,
    '有意': model.pvalues[reg_cols].values < 0.05
})
coef_df = coef_df.sort_values('標準化係数')
print('\n係数表:')
print(coef_df.to_string(index=False))

# R²・調整済みR²
r2     = model.rsquared
r2_adj = model.rsquared_adj
print(f'\nR² = {r2:.4f}, Adj-R² = {r2_adj:.4f}, N = {len(df_ols)}')

rank_df = df22[['都道府県', '人口増加率', '地域']].dropna().copy()
rank_df = rank_df.sort_values('人口増加率', ascending=True).reset_index(drop=True)
rank_df['略称'] = rank_df['都道府県'].map(abbr)
rank_df['色'] = rank_df.apply(
    lambda row: ('#e05c5c' if row['人口増加率'] >= 0 else region_colors.get(row['地域'], '#888')),
    axis=1
)

fig4, ax4 = plt.subplots(figsize=(11, 10))
bars4 = ax4.barh(rank_df['略称'], rank_df['人口増加率'],
                 color=rank_df['色'], edgecolor='white', height=0.7)
ax4.axvline(0, color='#333', linewidth=1.2)

# 数値ラベル
for bar, val in zip(bars4, rank_df['人口増加率']):
    ha = 'left' if val >= 0 else 'right'
    offset = 0.1 if val >= 0 else -0.1
    ax4.text(val + offset, bar.get_y() + bar.get_height()/2,
             f'{val:+.1f}%', va='center', ha=ha, fontsize=8)

# 地域カラー凡例
legend_elements = [
    mpatches.Patch(facecolor='#e05c5c', label='人口増加（正）'),
] + [
    mpatches.Patch(facecolor=region_colors[r], label=r)
    for r in region_colors
]
ax4.legend(handles=legend_elements, loc='lower right', fontsize=8, framealpha=0.9)

ax4.set_xlabel('人口増加率（2012→2022年、%）', fontsize=12)
ax4.set_title('図4　都道府県別 人口増加率ランキング（2012→2022年）', fontsize=13, fontweight='bold')
ax4.grid(True, axis='x', alpha=0.3)
fig4.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2021_H5_1_fig4.png'), dpi=150, bbox_inches='tight')
plt.close(fig4)
print('fig4 saved')

print('\n===== HTML埋め込み用統計サマリー =====')
print(f'分析対象: 47都道府県 (2022年断面)')
print(f'人口増加率 (2012→2022) 平均: {df22["人口増加率"].mean():.2f}%')
print(f'人口増加率 最大: {df22["人口増加率"].max():.2f}% ({df22.loc[df22["人口増加率"].idxmax(), "都道府県"]})')
print(f'人口増加率 最小: {df22["人口増加率"].min():.2f}% ({df22.loc[df22["人口増加率"].idxmin(), "都道府県"]})')
print(f'保育所密度 平均: {df22["保育所密度"].mean():.4f} 施設/千人')
print(f'合計特殊出生率 全国平均 (2022): {df22["合計特殊出生率"].mean():.2f}')
print(f'純移動率 平均: {df22["純移動率"].mean():.4f}%')
pos_pref = (df22['人口増加率'] > 0).sum()
print(f'人口増加都道府県数: {pos_pref}都道府県')
print(f'\nOLS結果 R²={r2:.4f}, Adj-R²={r2_adj:.4f}')
print(f'有意な変数: {coef_df[coef_df["有意"]]["変数"].tolist()}')

# 散布図の相関
r_nursery, p_nursery = stats.pearsonr(
    df_plot['保育所密度'].values, df_plot['合計特殊出生率'].values)
print(f'\n保育所密度×合計特殊出生率: r={r_nursery:.3f}, p={p_nursery:.3f}')

print('\n===== 完了 =====')