社会増減が合計特殊出生率に与える影響

import os, numpy as np, pandas as pd
import matplotlib; matplotlib.use('Agg')
import matplotlib.pyplot as plt
from matplotlib.patches import Patch
import statsmodels.api as sm

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'
YEAR    = 2017          # 原論文の断面年（SSDSE2020B の 2017 年に対応）
os.makedirs(FIG_DIR, exist_ok=True)

# ── SSDSE-B の読み込み（1行目=英字コード, 2行目=日本語列名） ──
df = pd.read_csv(DATA_B, encoding='cp932', header=0, skiprows=[1])
df = df.rename(columns={'SSDSE-B-2026': '年度', 'Code': '地域コード', 'Prefecture': '都道府県'})
df = df[df['地域コード'].str.match(r'^R\d{5}', na=False)].copy()   # 47都道府県のみ
df['年度'] = df['年度'].astype(int)

# ── 2017年断面（原論文の分析年） ──
d = df[df['年度'] == YEAR].copy().reset_index(drop=True)
d['合計特殊出生率'] = d['A4103']                              # A4103: TFR
d['学生数per1000']  = (d['E6301'] + d['E6302']) / d['A1101'] * 1000  # (短大+大学)/総人口×1000

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

r = d.sort_values('合計特殊出生率', ascending=True).reset_index(drop=True)
mean_tfr = d['合計特殊出生率'].mean()

fig1, ax1 = plt.subplots(figsize=(10, 12))
colors_bar = [region_colors[g] for g in r['地域']]
ax1.barh(r['都道府県'], r['合計特殊出生率'], color=colors_bar, edgecolor='white', linewidth=0.5)
# 京都府を強調
kyoto_pos = r.index[r['都道府県'] == '京都府'][0]
ax1.barh('京都府', r.loc[kyoto_pos, '合計特殊出生率'], color='#1a3e72', edgecolor='white')
ax1.axvline(mean_tfr, color='black', linestyle='--', linewidth=1.2, alpha=0.7)
ax1.text(mean_tfr + 0.01, 1, f'47都道府県\n単純平均 {mean_tfr:.2f}', fontsize=9, va='bottom', ha='left')
ax1.set_xlabel('合計特殊出生率（TFR）', fontsize=12)
ax1.set_title('都道府県別 合計特殊出生率ランキング（2017年）\n【実再現：SSDSE-B-2026 の 2017 年断面】',
              fontsize=13, fontweight='bold', pad=12)
ax1.set_xlim(1.0, 2.0)
ax1.tick_params(axis='y', labelsize=8.5)
legend_elements = [Patch(facecolor=v, label=k) for k, v in region_colors.items()]
ax1.legend(handles=legend_elements, loc='lower right', fontsize=9, framealpha=0.8)
plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2020_U5_4_fig1.png'), bbox_inches='tight')
plt.close(fig1)
print('[OK] fig1 saved')

X = sm.add_constant(d['学生数per1000'])
model = sm.OLS(d['合計特殊出生率'], X).fit()
b0, b1 = model.params['const'], model.params['学生数per1000']

fig2, ax2 = plt.subplots(figsize=(9, 7))
for region, grp in d.groupby('地域'):
    ax2.scatter(grp['学生数per1000'], grp['合計特殊出生率'],
                color=region_colors[region], label=region, s=60, zorder=3, alpha=0.85)
for name in ['京都府', '東京都', '沖縄県']:
    row = d[d['都道府県'] == name].iloc[0]
    ax2.annotate(name, xy=(row['学生数per1000'], row['合計特殊出生率']),
                 fontsize=10, fontweight='bold', ha='center', va='bottom',
                 xytext=(0, 6), textcoords='offset points')
xline = np.linspace(d['学生数per1000'].min(), d['学生数per1000'].max(), 100)
ax2.plot(xline, b0 + b1 * xline, 'k-', linewidth=1.6, alpha=0.75,
         label=f'回帰線 (係数={b1:.4f}, R²={model.rsquared:.3f})')
ax2.set_xlabel('人口1000人あたり学生数（短大＋大学）／千人', fontsize=12)
ax2.set_ylabel('合計特殊出生率（TFR）', fontsize=12)
ax2.set_title('人口あたり学生数と TFR の関係（2017年, 47都道府県）\n【実再現：SSDSE-B-2026】',
              fontsize=13, fontweight='bold', pad=10)
ax2.legend(fontsize=9, loc='upper right', framealpha=0.85, ncol=2)
ax2.grid(True, alpha=0.3)
plt.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2020_U5_4_fig2.png'), bbox_inches='tight')

ages   = ['15-19歳', '20-24歳', '25-29歳', '30-34歳', '35-39歳', '40-44歳', '45-49歳']
coef   = [-0.039, -0.012, -0.003, 0.116, 0.150, -0.065, 0.088]   # 偏回帰係数（報告値）
ci_lo  = [-0.064, -0.027, -0.044, 0.020, -0.005, -0.297, -0.172] # 95%C.I.下限（報告値）
ci_hi  = [-0.015,  0.002,  0.039, 0.213, 0.304,  0.167,  0.349]  # 95%C.I.上限（報告値）
sig    = ['**', '†', 'n.s.', '*', '†', 'n.s.', 'n.s.']

y = np.arange(len(ages))[::-1]
bar_col = ['#e05c5c' if s in ('*', '**') else '#aaaaaa' for s in sig]
fig3, ax3 = plt.subplots(figsize=(8.5, 5.5))
ax3.barh(y, coef, color=bar_col, height=0.55, zorder=3)
ax3.errorbar(coef, y, xerr=[np.array(coef) - np.array(ci_lo), np.array(ci_hi) - np.array(coef)],
             fmt='none', color='black', capsize=4, linewidth=1.4, zorder=4)
ax3.axvline(0, color='black', linewidth=1.0)
ax3.set_yticks(y); ax3.set_yticklabels(ages, fontsize=11)
for yi, c, s in zip(y, coef, sig):
    ax3.text(c + (0.01 if c >= 0 else -0.01), yi, s, va='center',
             ha='left' if c >= 0 else 'right', fontsize=11, fontweight='bold',
             color='#cc0000' if s in ('*', '**') else '#666')
ax3.set_xlabel('偏回帰係数 β（TFR への影響）', fontsize=12)
ax3.set_title('都道府県別 TFR に係る重回帰の偏回帰係数\n'
              '【原論文 表4 の報告値の可視化・再計算ではない／R²=0.408, n=188】',
              fontsize=12, fontweight='bold', pad=10)
ax3.grid(True, axis='x', alpha=0.3)
legend_items = [Patch(facecolor='#e05c5c', label='p<0.05（有意）'),
                Patch(facecolor='#aaaaaa', label='†p<0.1 / n.s.')]
ax3.legend(handles=legend_items, loc='lower right', fontsize=9)
plt.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2020_U5_4_fig3.png'), bbox_inches='tight')
plt.close(fig3)
print('[OK] fig3 saved')

imp_pref  = [0.697, 0.489, 0.293, 0.490, 0.275, 0.184, 0.254]   # 図6（報告値）
imp_kyoto = [0.939, 0.469, 1.504, 0.677, 0.551, 0.496, 0.428]   # 図7（報告値）
x = np.arange(len(ages)); w = 0.4
fig4, ax4 = plt.subplots(figsize=(9.5, 5.5))
ax4.bar(x - w/2, imp_pref, w, label='都道府県別（図6）', color='#1565C0', zorder=3)
ax4.bar(x + w/2, imp_kyoto, w, label='京都府 市区町村別（図7）', color='#E65100', zorder=3)
ax4.set_xticks(x); ax4.set_xticklabels(ages, fontsize=10)
ax4.set_ylabel('変数重要度（Importance）', fontsize=12)
ax4.set_title('ランダムフォレストによる年齢階級別 転入超過率の変数重要度\n'
              '【原論文 図6・図7 の報告値の可視化・再計算ではない】',
              fontsize=12, fontweight='bold', pad=10)
ax4.legend(fontsize=10); ax4.grid(True, axis='y', alpha=0.3)
plt.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2020_U5_4_fig4.png'), bbox_inches='tight')
plt.close(fig4)
print('[OK] fig4 saved')

print('\n=== 図1・図2の基礎統計（2017年, SSDSE-B, 実再現） ===')
print(f'47都道府県 単純平均TFR: {mean_tfr:.3f}')
low5 = r[['都道府県', '合計特殊出生率']].head(5).to_string(index=False)
print('TFRが低い5府県:\n' + low5)
kr = r.reset_index(drop=True)
kyoto_rank = kr.index[kr['都道府県'] == '京都府'][0] + 1
print(f'京都府: TFR={d.loc[d["都道府県"]=="京都府","合計特殊出生率"].iloc[0]:.2f}, '
      f'低い方から {kyoto_rank}/47 位')
print(f'最高: 沖縄県 {d["合計特殊出生率"].max():.2f} ／ 最低: 東京都 {d["合計特殊出生率"].min():.2f}')

print('\n=== 図2：人口1000人あたり学生数 → TFR 単回帰（実再現） ===')
print(f'切片            : {b0:.3f}')
print(f'回帰係数（学生数）: {b1:.4f}  (t={model.tvalues["学生数per1000"]:.3f}, '
      f'p={model.pvalues["学生数per1000"]:.4f})')
print(f'決定係数 R²      : {model.rsquared:.3f}  （補正 R²={model.rsquared_adj:.3f}, n={int(model.nobs)}）')
print('→ 学生数が多い府県ほど TFR が低い、という原論文の指摘（表5）が符号・有意性・R²の水準で再現された。')
print('   （原論文 表5 の報告値: 回帰係数 -0.147**, R²=0.200, n=47。'
      '係数の絶対値は変数スケール・データ版の違いで異なる。）')