合計特殊出生率の変動要因分析

df = pd.DataFrame()
df['都道府県'] = df_b['都道府県']

# TFR（合計特殊出生率, 2022年）
df['TFR'] = pd.to_numeric(df_b['合計特殊出生率'], errors='coerce')

# 婚姻率 = 婚姻件数 / 15〜64歳人口 × 1000
df['婚姻率'] = (
    pd.to_numeric(df_b['婚姻件数'], errors='coerce') /
    pd.to_numeric(df_b['15～64歳人口'], errors='coerce') * 1000
)

# 高齢化率 = 65歳以上人口 / 総人口 × 100
df['高齢化率'] = (
    pd.to_numeric(df_b['65歳以上人口'], errors='coerce') /
    pd.to_numeric(df_b['総人口'], errors='coerce') * 100
)

# 大学進学率 = 高等学校卒業者のうち進学者数 / 高等学校卒業者数 × 100
df['大学進学率'] = (
    pd.to_numeric(df_b['高等学校卒業者のうち進学者数'], errors='coerce') /
    pd.to_numeric(df_b['高等学校卒業者数'], errors='coerce') * 100
)

# 保育所数（千対）= 保育所等数 / 総人口 × 1000
df['保育所数千対'] = (
    pd.to_numeric(df_b['保育所等数'], errors='coerce') /
    pd.to_numeric(df_b['総人口'], errors='coerce') * 1000
)

# 年平均気温（℃）
df['年平均気温'] = pd.to_numeric(df_b['年平均気温'], errors='coerce')

# 1人当たり県民所得（万円, SSDSE-E 2021年基準）
income_map = dict(zip(
    ssdse_e['都道府県'],
    pd.to_numeric(ssdse_e['1人当たり県民所得（平成27年基準）'], errors='coerce')
))
df['1人当たり県民所得'] = df['都道府県'].map(income_map)

# 人口密度（人/km²）= 総人口 / 総面積(ha÷100)
area_map = dict(zip(
    ssdse_e['都道府県'],
    pd.to_numeric(ssdse_e['総面積（北方地域及び竹島を除く）'], errors='coerce') / 100
))
df['面積_km2'] = df['都道府県'].map(area_map)
df['人口密度'] = pd.to_numeric(df_b['総人口'], errors='coerce') / df['面積_km2']
df['ln人口密度'] = np.log(df['人口密度'])

# 3区分の設定
df['区分'] = 'その他'
df.loc[df['都道府県'] == '東京都', '区分'] = '東京都'
df.loc[(df['人口密度'] < 150) & (df['都道府県'] != '東京都'), '区分'] = '過疎化進行県'

df = df.dropna().reset_index(drop=True)
N = len(df)

VAR_NAMES = ['婚姻率', '1人当たり県民所得', '高齢化率', 'ln人口密度',
             '保育所数千対', '大学進学率', '年平均気温']

print("=" * 60)
print("■ 記述統計（都道府県別, 2022年）N =", N)
print("=" * 60)
print(df[['TFR'] + VAR_NAMES].describe().round(3))

print("\n3区分の内訳:")
for g in ['東京都', '過疎化進行県', 'その他']:
    prefs_g = df[df['区分'] == g]['都道府県'].tolist()
    print(f"  {g}（{len(prefs_g)}）: {prefs_g}")

print("\n" + "=" * 60)
print("■ Step3. 3区分別TFR比較")
print("=" * 60)
for g in ['東京都', 'その他', '過疎化進行県']:
    sub = df[df['区分'] == g]['TFR']
    print(f"  {g}: N={len(sub)}, 平均TFR={sub.mean():.3f}, 中央値={sub.median():.3f}")

fig1, ax1 = plt.subplots(figsize=(10, 10))
sorted_idx = np.argsort(tfr)
vals_s = tfr[sorted_idx]
prefs_s = df['都道府県'].values[sorted_idx]
grp_s   = df['区分'].values[sorted_idx]
bar_cols = [color_map[g] for g in grp_s]

ax1.barh(range(N), vals_s, color=bar_cols, alpha=0.8, edgecolor='white', height=0.75)
ax1.set_yticks(range(N))
ax1.set_yticklabels(prefs_s, fontsize=8)
ax1.axvline(tfr.mean(), color='gray', linestyle='--', linewidth=1.2)
ax1.axvline(2.07, color='black', linestyle=':', linewidth=1.2, alpha=0.6)
ax1.set_xlabel('合計特殊出生率（TFR）', fontsize=12)
ax1.set_title('都道府県別TFR（2022年）\n3区分（過疎化進行県・東京都・その他）', fontsize=12, fontweight='bold')
ax1.grid(axis='x', alpha=0.3)
legend_els = [Patch(color=c, alpha=0.8, label=g) for g, c in color_map.items()]
legend_els += [
    plt.Line2D([0], [0], color='gray', linestyle='--', label=f'全国平均={tfr.mean():.3f}'),
    plt.Line2D([0], [0], color='black', linestyle=':', label='人口置換水準=2.07')
]
ax1.legend(handles=legend_els, fontsize=9)
plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2025_H5_6_fig1_rank.png'), bbox_inches='tight', dpi=150)
plt.close(fig1)
print("\n図1保存: 2025_H5_6_fig1_rank.png")

fig2, axes2 = plt.subplots(2, 2, figsize=(12, 9))
fig2.suptitle('TFRと主要説明変数の散布図（2022年 SSDSE-B/E）', fontsize=13, fontweight='bold')

key_pairs = [
    (VAR_NAMES.index('婚姻率'), '婚姻率（‰）'),
    (VAR_NAMES.index('ln人口密度'), 'ln(人口密度)'),
    (VAR_NAMES.index('保育所数千対'), '保育所数（千対）'),
    (VAR_NAMES.index('大学進学率'), '大学進学率（%）'),
]
groups = df['区分'].values

for ax, (idx, xlabel) in zip(axes2.flat, key_pairs):
    xv = X_raw[:, idx]
    for g, col in color_map.items():
        mask = groups == g
        ax.scatter(xv[mask], tfr[mask], color=col, s=60, alpha=0.75,
                   edgecolors='white', linewidth=0.5, label=g, zorder=3)
    z = np.polyfit(xv, tfr, 1)
    xs = np.linspace(xv.min() - 0.1, xv.max() + 0.1, 100)
    ax.plot(xs, np.poly1d(z)(xs), 'k--', linewidth=1.5, alpha=0.6)
    r, p = stats.pearsonr(xv, tfr)
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else 'n.s.'
    ax.set_xlabel(xlabel, fontsize=10)
    ax.set_ylabel('TFR', fontsize=10)
    ax.set_title(f'r={r:.3f} {sig}', fontsize=11, fontweight='bold')
    ax.grid(True, alpha=0.3)
    if idx == 0:
        ax.legend(fontsize=7)

plt.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2025_H5_6_fig2_scatter.png'), bbox_inches='tight', dpi=150)
plt.close(fig2)
print("図2保存: 2025_H5_6_fig2_scatter.png")

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
from sklearn.preprocessing import StandardScaler
from matplotlib.patches import Patch
import warnings
warnings.filterwarnings('ignore')
import os

plt.rcParams['font.family'] = 'Hiragino Sans'
plt.rcParams['axes.unicode_minus'] = False
plt.rcParams['figure.dpi'] = 150

BASE_DIR = os.path.join(_script_dir, '..')
FIG_DIR = os.path.join(BASE_DIR, 'html', 'figures')
DATA_DIR = os.path.join(BASE_DIR, 'data', 'raw')

ssdse_b = pd.read_csv(os.path.join(DATA_DIR, 'SSDSE-B-2026.csv'), encoding='cp932', header=1)
df_b = ssdse_b[
    (ssdse_b['年度'] == 2022) &
    ssdse_b['地域コード'].str.match(r'^R\d{5}$', na=False)
].copy().reset_index(drop=True)

# SSDSE-E: 都道府県別クロスセクション（1人当たり県民所得, 総面積）
ssdse_e_raw = pd.read_csv(os.path.join(DATA_DIR, 'SSDSE-E-2026.csv'), encoding='cp932', header=1)
ssdse_e = ssdse_e_raw.iloc[1:].copy()
ssdse_e.columns = ssdse_e_raw.iloc[0].values
ssdse_e = ssdse_e[ssdse_e['都道府県'] != '全国'].reset_index(drop=True)

print("\n" + "=" * 60)
print("■ Step1. 相関分析（TFR vs 各説明変数）")
print("=" * 60)
print(f"\n  {'変数':<16} {'r':>8} {'p値':>10} {'有意':>6}")
print("  " + "-" * 44)
X_raw = df[VAR_NAMES].values
tfr = df['TFR'].values
for i, name in enumerate(VAR_NAMES):
    r, p = stats.pearsonr(X_raw[:, i], tfr)
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else 'n.s.'
    print(f"  {name:<16} {r:>8.4f} {p:>10.4f} {sig:>6}")

print("\n" + "=" * 60)
print("■ Step2. 重回帰分析（全説明変数）")
print("=" * 60)
X_ols = sm.add_constant(X_raw)
model = sm.OLS(tfr, X_ols).fit()
print(model.summary2())

color_map = {'過疎化進行県': '#E65100', '東京都': '#1565C0', 'その他': '#2E7D32'}

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

scaler = StandardScaler()
X_std = scaler.fit_transform(X_raw)
model_std = sm.OLS(tfr, sm.add_constant(X_std)).fit()
coefs = model_std.params[1:]
ses   = model_std.bse[1:]
pvals = model_std.pvalues[1:]

bar_cols3 = ['#C62828' if p < 0.05 else '#1565C0' if p < 0.1 else '#9E9E9E' for p in pvals]
yp = np.arange(len(VAR_NAMES))
ax3.barh(yp, coefs, color=bar_cols3, alpha=0.75, edgecolor='white', height=0.6)
ax3.errorbar(coefs, yp, xerr=1.96 * ses, fmt='none', color='#333', capsize=4, linewidth=1.2)
ax3.axvline(0, color='gray', linestyle='--', linewidth=1.0)
ax3.set_yticks(yp)
ax3.set_yticklabels(VAR_NAMES, fontsize=11)
ax3.set_xlabel('標準化回帰係数（±1.96SE）', fontsize=11)
ax3.set_title(f'TFR重回帰の標準化係数（データ：SSDSE-B/E 2022年）\nR²={model.rsquared:.3f} N={N}都道府県', fontsize=12, fontweight='bold')
ax3.invert_yaxis()
ax3.grid(axis='x', alpha=0.3)
legend_els3 = [Patch(color='#C62828', alpha=0.75, label='p<0.05'),
               Patch(color='#1565C0', alpha=0.75, label='p<0.10'),
               Patch(color='#9E9E9E', alpha=0.75, label='非有意')]
ax3.legend(handles=legend_els3, fontsize=9)
plt.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2025_H5_6_fig3_coef.png'), bbox_inches='tight', dpi=150)
plt.close(fig3)
print("図3保存: 2025_H5_6_fig3_coef.png")

fig4, axes4 = plt.subplots(1, 2, figsize=(13, 5))
fig4.suptitle('少子化の2タイプ：東京都 vs 過疎化進行県', fontsize=13, fontweight='bold')

group_order = ['東京都', 'その他', '過疎化進行県']
group_data  = [df[df['区分'] == g]['TFR'].values for g in group_order]
group_cols  = [color_map[g] for g in group_order]

ax4a = axes4[0]
bp4 = ax4a.boxplot(group_data,
                   labels=['東京都\n（高コスト）', 'その他', '過疎化進行県\n（若年層流出）'],
                   patch_artist=True, medianprops=dict(color='white', linewidth=2))
for patch, col in zip(bp4['boxes'], group_cols):
    patch.set_facecolor(col)
    patch.set_alpha(0.65)
ax4a.set_ylabel('TFR（合計特殊出生率）', fontsize=11)
ax4a.set_title('3区分別TFR分布（2022年）', fontsize=11, fontweight='bold')
ax4a.axhline(2.07, color='black', linestyle=':', linewidth=1.2, alpha=0.6, label='人口置換水準')
ax4a.axhline(tfr.mean(), color='gray', linestyle='--', linewidth=1.0, label=f'全国平均={tfr.mean():.3f}')
ax4a.legend(fontsize=8)
ax4a.grid(axis='y', alpha=0.3)

# 3区分の指標比較（全国平均比）
ax4b = axes4[1]
indicators = ['婚姻率', '保育所数千対', '大学進学率', 'TFR']
x_idx = np.arange(len(indicators))
width = 0.25
grp_colors_bar = [color_map[g] for g in group_order]
for i, (g, col) in enumerate(zip(group_order, grp_colors_bar)):
    sub = df[df['区分'] == g]
    vals_g = []
    for var in indicators:
        if var == 'TFR':
            vals_g.append(sub['TFR'].mean())
        else:
            vals_g.append(sub[var].mean() / df[var].mean())
    ax4b.bar(x_idx + i * width, vals_g, width, label=g, color=col, alpha=0.75, edgecolor='white')
ax4b.set_xticks(x_idx + width)
ax4b.set_xticklabels(['婚姻率', '保育所数', '大学進学率', 'TFR'], fontsize=10)
ax4b.set_ylabel('全国平均比（TFRは実値）', fontsize=10)
ax4b.set_title('3区分の特徴比較\n（各指標の全国平均比）', fontsize=11, fontweight='bold')
ax4b.axhline(1.0, color='gray', linestyle='--', linewidth=1.0, alpha=0.6, label='全国平均=1')
ax4b.legend(fontsize=8)
ax4b.grid(axis='y', alpha=0.3)

plt.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2025_H5_6_fig4_compare.png'), bbox_inches='tight', dpi=150)
plt.close(fig4)
print("図4保存: 2025_H5_6_fig4_compare.png")

print("\n全図の生成完了（4枚）")
print(f"\nデータ出典: SSDSE-B-2026.csv（2022年度）, SSDSE-E-2026.csv（統計センター）")