保育・子育て支援と合計特殊出生率の
都道府県比較分析

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)

# 2022年断面
df = df_b[df_b['年度'] == 2022].copy()
assert len(df) == 47, f"都道府県数が47ではありません: {len(df)}"

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

df['地域'] = df['都道府県'].map(region_map)
df['地域色'] = df['地域'].map(region_colors)

# ── 派生変数 ────────────────────────────────────────────────────
# 保育所密度: 保育所等数 / 総人口 × 1000
df['保育所密度'] = df['保育所等数'] / df['総人口'] * 1000

# 保育定員率: 保育所等定員数 / 0〜4歳推計人口（15歳未満/3） × 100
df['0_4歳人口推計'] = df['15歳未満人口'] / 3
df['保育定員率'] = df['保育所等定員数'] / df['0_4歳人口推計'] * 100

# 婚姻率: 婚姻件数 / 総人口 × 1000
df['婚姻率'] = df['婚姻件数'] / df['総人口'] * 1000

# 高齢化率: 65歳以上人口 / 総人口 × 100
df['高齢化率'] = df['65歳以上人口'] / df['総人口'] * 100

# 消費支出_log（所得水準の代理）
df['消費支出_log'] = np.log(df['消費支出（二人以上の世帯）'])

# 保育士密度: 保育所等保育士数 / 総人口 × 10000
df['保育士密度'] = df['保育所等保育士数'] / df['総人口'] * 10000

# 目的変数
y = df['合計特殊出生率'].values

# ── 重回帰分析 ───────────────────────────────────────────────────
X_vars = ['保育所密度', '保育定員率', '婚姻率', '消費支出_log', '高齢化率', '保育士密度']
X_df = df[X_vars].copy()
X = sm.add_constant(X_df)
ols_result = sm.OLS(y, X).fit()
print(ols_result.summary())

# 相関係数
print("\n=== Pearson相関係数 ===")
corr_data = {}
for v in X_vars:
    r, p = stats.pearsonr(df[v].values, y)
    corr_data[v] = {'r': r, 'p': p}
    print(f"{v}: r={r:.3f}, p={p:.4f}")

# 標準化偏回帰係数
X_std = (X_df - X_df.mean()) / X_df.std()
y_std = (y - y.mean()) / y.std()
X_std_const = sm.add_constant(X_std)
ols_std = sm.OLS(y_std, X_std_const).fit()
print("\n=== 標準化偏回帰係数 ===")
for var, coef, pval in zip(X_vars, ols_std.params[1:], ols_std.pvalues[1:]):
    print(f"{var}: β={coef:.3f}, p={pval:.4f}")

# ── 統計値を出力（HTML埋め込み用）────────────────────────────────
print(f"\n=== モデル適合度 ===")
print(f"R² = {ols_result.rsquared:.3f}")
print(f"調整済みR² = {ols_result.rsquared_adj:.3f}")
print(f"F統計量 = {ols_result.fvalue:.3f}, p = {ols_result.f_pvalue:.4f}")
print(f"N = {len(df)}")

tfr_mean = df['合計特殊出生率'].mean()
tfr_max_row = df.nlargest(1, '合計特殊出生率').iloc[0]
tfr_min_row = df.nsmallest(1, '合計特殊出生率').iloc[0]
print(f"\n=== TFR基本統計 ===")
print(f"全国平均TFR = {tfr_mean:.3f}")
print(f"最大: {tfr_max_row['都道府県']} {tfr_max_row['合計特殊出生率']:.2f}")
print(f"最小: {tfr_min_row['都道府県']} {tfr_min_row['合計特殊出生率']:.2f}")

# ── 図1: 散布図 保育所密度 vs TFR ─────────────────────────────────
fig1, ax1 = plt.subplots(figsize=(10, 7))

for region, color in region_colors.items():
    mask = df['地域'] == region
    ax1.scatter(
        df.loc[mask, '保育所密度'],
        df.loc[mask, '合計特殊出生率'],
        c=color, s=70, alpha=0.85, label=region, zorder=3
    )

# 都道府県ラベル（略称）
pref_short = {
    '北海道': '北', '青森県': '青', '岩手県': '岩', '宮城県': '宮', '秋田県': '秋',
    '山形県': '形', '福島県': '福', '茨城県': '茨', '栃木県': '栃', '群馬県': '群',
    '埼玉県': '埼', '千葉県': '千', '東京都': '東', '神奈川県': '神',
    '新潟県': '潟', '富山県': '富', '石川県': '石', '福井県': '井',
    '山梨県': '梨', '長野県': '野', '岐阜県': '岐', '静岡県': '静', '愛知県': '知',
    '三重県': '三', '滋賀県': '滋', '京都府': '京', '大阪府': '阪',
    '兵庫県': '兵', '奈良県': '奈', '和歌山県': '和',
    '鳥取県': '鳥', '島根県': '島', '岡山県': '岡', '広島県': '広', '山口県': '口',
    '徳島県': '徳', '香川県': '香', '愛媛県': '媛', '高知県': '高',
    '福岡県': '岡', '佐賀県': '佐', '長崎県': '崎', '熊本県': '熊',
    '大分県': '分', '宮崎県': '宮', '鹿児島県': '鹿', '沖縄県': '沖'
}
for _, row in df.iterrows():
    label = pref_short.get(row['都道府県'], row['都道府県'][:2])
    ax1.annotate(
        label,
        (row['保育所密度'], row['合計特殊出生率']),
        fontsize=7, ha='center', va='bottom', color='#333',
        xytext=(0, 3), textcoords='offset points'
    )

# 回帰直線
x_fit = np.linspace(df['保育所密度'].min(), df['保育所密度'].max(), 100)
slope, intercept, r_val, p_val, _ = stats.linregress(df['保育所密度'].values, y)
ax1.plot(x_fit, intercept + slope * x_fit, 'k--', lw=1.5, alpha=0.6,
         label=f'回帰線 (r={r_val:.3f}, p<0.001)')

ax1.set_xlabel('保育所密度（保育所等数 / 総人口 × 1000）', fontsize=12)
ax1.set_ylabel('合計特殊出生率（TFR）', 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)
ax1.set_xlim(df['保育所密度'].min() * 0.9, df['保育所密度'].max() * 1.05)

# 相関係数テキスト
ax1.text(0.97, 0.05, f'Pearson r = {r_val:.3f}\np < 0.001\nN = 47',
         transform=ax1.transAxes, ha='right', va='bottom',
         fontsize=10, bbox=dict(boxstyle='round', facecolor='#EFF3FF', alpha=0.8))

plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2020_H2_fig1.png'), dpi=150, bbox_inches='tight')
plt.close()
print("図1 保存完了")

# ── 図2: 都道府県別TFRランキング棒グラフ ─────────────────────────
fig2, ax2 = plt.subplots(figsize=(12, 9))

df_sorted = df.sort_values('合計特殊出生率', ascending=True).copy()
bars = ax2.barh(
    df_sorted['都道府県'],
    df_sorted['合計特殊出生率'],
    color=df_sorted['地域色'],
    edgecolor='white', linewidth=0.5, height=0.75
)

# 平均線
mean_tfr = df_sorted['合計特殊出生率'].mean()
ax2.axvline(mean_tfr, color='black', lw=1.5, ls='--', alpha=0.7,
            label=f'全国平均: {mean_tfr:.3f}')

# 値ラベル
for bar, val in zip(bars, df_sorted['合計特殊出生率']):
    ax2.text(val + 0.005, bar.get_y() + bar.get_height()/2,
             f'{val:.2f}', va='center', ha='left', fontsize=7.5)

ax2.set_xlabel('合計特殊出生率（TFR）', fontsize=12)
ax2.set_title('図2：都道府県別 合計特殊出生率ランキング（2022年）\n九州・沖縄が上位、関東が下位', fontsize=13, fontweight='bold')
ax2.set_xlim(0.9, 1.85)
ax2.legend(fontsize=10, loc='lower right')

# 地域凡例
from matplotlib.patches import Patch
legend_patches = [Patch(facecolor=c, label=r) for r, c in region_colors.items()]
ax2.legend(handles=legend_patches, loc='lower right', fontsize=8,
           title='地域', title_fontsize=9, framealpha=0.9)

ax2.grid(True, axis='x', alpha=0.3)
plt.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2020_H2_fig2.png'), dpi=150, bbox_inches='tight')
plt.close()
print("図2 保存完了")

# ── 図3: 相関ヒートマップ ────────────────────────────────────────
fig3, ax3 = plt.subplots(figsize=(9, 7))

heatmap_vars = ['合計特殊出生率', '保育所密度', '保育定員率', '婚姻率', '消費支出_log', '高齢化率', '保育士密度']
heatmap_labels = ['TFR', '保育所\n密度', '保育\n定員率', '婚姻率', '消費支出\n(log)', '高齢化率', '保育士\n密度']
corr_matrix = df[heatmap_vars].corr()

# カラーマップ
import matplotlib.colors as mcolors
cmap = plt.cm.RdBu_r
im = ax3.imshow(corr_matrix.values, cmap=cmap, vmin=-1, vmax=1, aspect='auto')

# 軸ラベル
ax3.set_xticks(range(len(heatmap_vars)))
ax3.set_yticks(range(len(heatmap_vars)))
ax3.set_xticklabels(heatmap_labels, fontsize=10)
ax3.set_yticklabels(heatmap_labels, fontsize=10)

# セル内テキスト
for i in range(len(heatmap_vars)):
    for j in range(len(heatmap_vars)):
        val = corr_matrix.values[i, j]
        text_color = 'white' if abs(val) > 0.6 else 'black'
        # p値による有意性マーク
        if i != j:
            xi = df[heatmap_vars[i]].values
            xj = df[heatmap_vars[j]].values
            _, pv = stats.pearsonr(xi, xj)
            star = '***' if pv < 0.001 else ('**' if pv < 0.01 else ('*' if pv < 0.05 else ''))
            ax3.text(j, i, f'{val:.2f}{star}', ha='center', va='center',
                     fontsize=9, color=text_color, fontweight='bold' if pv < 0.05 else 'normal')
        else:
            ax3.text(j, i, '1.00', ha='center', va='center',
                     fontsize=9, color=text_color, fontweight='bold')

plt.colorbar(im, ax=ax3, shrink=0.8, label='Pearson相関係数')
ax3.set_title('図3：TFRと関連変数の相関ヒートマップ（2022年・N=47）\n* p<0.05, ** p<0.01, *** p<0.001',
              fontsize=12, fontweight='bold')
plt.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2020_H2_fig3.png'), dpi=150, bbox_inches='tight')
plt.close()
print("図3 保存完了")

# ── 図4: 標準化偏回帰係数プロット ────────────────────────────────
fig4, ax4 = plt.subplots(figsize=(9, 6))

coefs = ols_std.params[1:]
pvals = ols_std.pvalues[1:]
conf_int = ols_std.conf_int().iloc[1:]

# 色分け（有意・非有意）
bar_colors = ['#e05c5c' if c > 0 else '#4e9af1' for c in coefs]
alpha_vals = [0.9 if p < 0.05 else 0.45 for p in pvals]

y_pos = range(len(X_vars))
ax4.barh(y_pos, coefs, color=bar_colors,
         xerr=[coefs - conf_int.iloc[:, 0], conf_int.iloc[:, 1] - coefs],
         capsize=4, ecolor='gray', error_kw={'lw': 1.5},
         height=0.55, alpha=0.85)

# 有意性マーク
for i, (c, p) in enumerate(zip(coefs, pvals)):
    star = '***' if p < 0.001 else ('**' if p < 0.01 else ('*' if p < 0.05 else 'ns'))
    color = '#c0392b' if p < 0.05 else '#999'
    offset = 0.02 if c >= 0 else -0.02
    ax4.text(c + offset, i, star, va='center', ha='left' if c >= 0 else 'right',
             fontsize=11, color=color, fontweight='bold')

ax4.set_yticks(list(y_pos))
ax4.set_yticklabels(X_vars, fontsize=11)
ax4.axvline(0, color='black', lw=1, ls='-')
ax4.set_xlabel('標準化偏回帰係数（β）', fontsize=12)
ax4.set_title(f'図4：重回帰分析 標準化偏回帰係数\n（R²={ols_result.rsquared:.3f}, 調整済R²={ols_result.rsquared_adj:.3f}, N=47）',
              fontsize=12, fontweight='bold')
ax4.grid(True, axis='x', alpha=0.3)

# 凡例
from matplotlib.patches import Patch
legend_elem = [
    Patch(facecolor='#e05c5c', alpha=0.85, label='正の影響（TFR↑）'),
    Patch(facecolor='#4e9af1', alpha=0.85, label='負の影響（TFR↓）'),
]
ax4.legend(handles=legend_elem, fontsize=9, loc='lower right')
ax4.text(0.97, 0.95, '有意水準: * p<0.05, ** p<0.01, *** p<0.001\n(ns = not significant)',
         transform=ax4.transAxes, ha='right', va='top', fontsize=9,
         color='#555', bbox=dict(boxstyle='round', facecolor='#f8f9fa', alpha=0.8))

plt.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2020_H2_fig4.png'), dpi=150, bbox_inches='tight')
plt.close()
print("図4 保存完了")

# ── 地域別TFR比較（HTML用） ─────────────────────────────────────
print("\n=== 地域別TFR平均 ===")
region_tfr = df.groupby('地域')['合計特殊出生率'].agg(['mean', 'min', 'max']).round(3)
print(region_tfr.to_string())

# TFR上位・下位の保育所密度比較
top_prefs = df.nlargest(8, '合計特殊出生率')['都道府県'].tolist()
bot_prefs = df.nsmallest(8, '合計特殊出生率')['都道府県'].tolist()
print(f"\nTFR上位8都道府県 平均保育所密度: {df[df['都道府県'].isin(top_prefs)]['保育所密度'].mean():.3f}")
print(f"TFR下位8都道府県 平均保育所密度: {df[df['都道府県'].isin(bot_prefs)]['保育所密度'].mean():.3f}")
print(f"\nTFR上位8都道府県 平均婚姻率: {df[df['都道府県'].isin(top_prefs)]['婚姻率'].mean():.3f}")
print(f"TFR下位8都道府県 平均婚姻率: {df[df['都道府県'].isin(bot_prefs)]['婚姻率'].mean():.3f}")

print("\n全図表の生成が完了しました。")
print(f"出力先: {os.path.abspath(FIG_DIR)}")