保育所整備と女性就業率が
合計特殊出生率に与える影響の
パネルデータ分析

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)

# 2012〜2022年度に絞る
df_b = df_b[df_b['年度'].between(2012, 2022)].copy()
df_b = df_b.reset_index(drop=True)

df_b['TFR'] = df_b['合計特殊出生率'].astype(float)

# 保育所数（人口千人当たり）
df_b['保育所千人当たり'] = df_b['保育所等数'] / df_b['総人口'] * 1000

# 女性就業率の代理変数：就職件数（一般）/ 15〜64歳女性人口 × 1000
# （労働市場での女性活動強度の代理指標）
df_b['女性就業活動率'] = df_b['就職件数（一般）'] / df_b['15～64歳人口（女）'] * 1000

# 婚姻率（人口千人当たり）
df_b['婚姻率'] = df_b['婚姻件数'] / df_b['総人口'] * 1000

# 消費支出（所得水準の代理変数）：対数変換
df_b['log消費支出'] = np.log(df_b['消費支出（二人以上の世帯）'].astype(float))

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

# 保育所等定員充足率（待機児童割合の逆）
df_b['保育所定員率'] = df_b['保育所等定員数'] / (df_b['保育所等在所児数'] + 1) * 100

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

region_colors = {
    '北海道東北': '#4e9af1',
    '関東': '#e05c5c',
    '中部': '#f0a500',
    '近畿': '#5cb85c',
    '中国四国': '#9b59b6',
    '九州沖縄': '#f39c12',
}

# 分析対象カラムの欠損除去
analysis_cols = ['TFR', '保育所千人当たり', '女性就業活動率', '婚姻率', 'log消費支出']
df_b = df_b.dropna(subset=analysis_cols).copy()

print(f"分析データ：{df_b['都道府県'].nunique()}都道府県 × {df_b['年度'].nunique()}年度 = {len(df_b)}観測値")
print(f"年度範囲：{df_b['年度'].min()}〜{df_b['年度'].max()}")

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

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

# 地域別平均
for region, color in region_colors.items():
    region_data = df_b[df_b['地域'] == region].groupby('年度')['TFR'].mean()
    ax.plot(region_data.index, region_data.values,
            color=color, linewidth=2, marker='o', markersize=4,
            label=region, alpha=0.85)

# 全国平均
national_avg = df_b.groupby('年度')['TFR'].mean()
ax.plot(national_avg.index, national_avg.values,
        color='#1a1a2e', linewidth=3, marker='D', markersize=6,
        label='全国平均', linestyle='--', zorder=5)

ax.set_xlabel('年度', fontsize=12)
ax.set_ylabel('合計特殊出生率', fontsize=12)
ax.set_title('都道府県別・合計特殊出生率の時系列変化（2012〜2022年）', fontsize=13, fontweight='bold')
ax.legend(loc='upper right', fontsize=10, framealpha=0.9)
ax.set_xticks(years)
ax.tick_params(axis='x', rotation=45)
ax.grid(axis='y', alpha=0.3, linestyle='--')
ax.set_ylim(1.0, 2.2)

# 注釈：2020年コロナ
ax.axvline(x=2020, color='gray', linestyle=':', alpha=0.5, linewidth=1.5)
ax.text(2020.1, 1.05, 'COVID-19', fontsize=9, color='gray')

plt.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2019_U4_fig1.png'), dpi=150, bbox_inches='tight')
plt.close(fig)
print("fig1 saved")

df_2022 = df_b[df_b['年度'] == 2022].copy()

fig, ax = plt.subplots(figsize=(9, 7))

for region, color in region_colors.items():
    sub = df_2022[df_2022['地域'] == region]
    ax.scatter(sub['女性就業活動率'], sub['TFR'],
               color=color, s=70, alpha=0.85, label=region, zorder=3, edgecolors='white', linewidths=0.5)

# 回帰直線
x_arr = df_2022['女性就業活動率'].values
y_arr = df_2022['TFR'].values
mask = np.isfinite(x_arr) & np.isfinite(y_arr)
x_arr, y_arr = x_arr[mask], y_arr[mask]

slope, intercept, r_value, p_value, std_err = stats.linregress(x_arr, y_arr)
x_line = np.linspace(x_arr.min(), x_arr.max(), 100)
y_line = slope * x_line + intercept
ax.plot(x_line, y_line, color='#333333', linewidth=2, linestyle='--', zorder=4)

# 相関係数・p値
ax.text(0.05, 0.95,
        f'r = {r_value:.3f}\np = {p_value:.4f}',
        transform=ax.transAxes, fontsize=11, verticalalignment='top',
        bbox=dict(boxstyle='round,pad=0.5', facecolor='white', edgecolor='#cccccc', alpha=0.9))

# 都道府県ラベル（上位・下位）
label_prefs = ['沖縄', '東京', '秋田', '北海道', '宮崎', '島根', '愛知', '大阪']
for _, row in df_2022[df_2022['都道府県'].isin(label_prefs)].iterrows():
    if np.isfinite(row['女性就業活動率']) and np.isfinite(row['TFR']):
        ax.annotate(row['都道府県'],
                    (row['女性就業活動率'], row['TFR']),
                    xytext=(4, 4), textcoords='offset points', fontsize=8.5, color='#333333')

ax.set_xlabel('女性就業活動率（就職件数/15〜64歳女性人口×1000）', fontsize=11)
ax.set_ylabel('合計特殊出生率（2022年）', fontsize=11)
ax.set_title('女性就業活動率と合計特殊出生率の関係（2022年・都道府県別）', fontsize=12, fontweight='bold')
ax.legend(loc='lower right', fontsize=9, framealpha=0.9)
ax.grid(alpha=0.3, linestyle='--')

plt.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2019_U4_fig2.png'), dpi=150, bbox_inches='tight')
plt.close(fig)
print("fig2 saved")

fig, ax = plt.subplots(figsize=(9, 7))

for region, color in region_colors.items():
    sub = df_2022[df_2022['地域'] == region]
    ax.scatter(sub['保育所千人当たり'], sub['TFR'],
               color=color, s=70, alpha=0.85, label=region, zorder=3, edgecolors='white', linewidths=0.5)

# 回帰直線
x_arr2 = df_2022['保育所千人当たり'].values
y_arr2 = df_2022['TFR'].values
mask2 = np.isfinite(x_arr2) & np.isfinite(y_arr2)
x2, y2 = x_arr2[mask2], y_arr2[mask2]

slope2, intercept2, r2, p2, se2 = stats.linregress(x2, y2)
x_line2 = np.linspace(x2.min(), x2.max(), 100)
y_line2 = slope2 * x_line2 + intercept2
ax.plot(x_line2, y_line2, color='#333333', linewidth=2, linestyle='--', zorder=4)

ax.text(0.05, 0.95,
        f'r = {r2:.3f}\np = {p2:.4f}',
        transform=ax.transAxes, fontsize=11, verticalalignment='top',
        bbox=dict(boxstyle='round,pad=0.5', facecolor='white', edgecolor='#cccccc', alpha=0.9))

# 都道府県ラベル
label_prefs2 = ['沖縄', '東京', '秋田', '福井', '長崎', '鳥取', '大阪', '神奈川']
for _, row in df_2022[df_2022['都道府県'].isin(label_prefs2)].iterrows():
    if np.isfinite(row['保育所千人当たり']) and np.isfinite(row['TFR']):
        ax.annotate(row['都道府県'],
                    (row['保育所千人当たり'], row['TFR']),
                    xytext=(4, 4), textcoords='offset points', fontsize=8.5, color='#333333')

ax.set_xlabel('保育所等数（人口千人当たり）', fontsize=11)
ax.set_ylabel('合計特殊出生率（2022年）', fontsize=11)
ax.set_title('保育所数（人口千人当たり）と合計特殊出生率の関係（2022年）', fontsize=12, fontweight='bold')
ax.legend(loc='lower right', fontsize=9, framealpha=0.9)
ax.grid(alpha=0.3, linestyle='--')

plt.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2019_U4_fig3.png'), dpi=150, bbox_inches='tight')
plt.close(fig)
print("fig3 saved")