女性就労と少子化
女性就業率・育児環境と合計特殊出生率の関係

df_raw = pd.read_csv(DATA_B, encoding='cp932', header=0)
# 行0はコード行（英字）、行1はラベル行（日本語）、行2以降がデータ
data = df_raw.iloc[1:].copy()
data.columns = df_raw.columns

# 列名マッピング（コード → 意味）
# SSDSE-B-2026  = 年度
# Code          = 地域コード
# Prefecture    = 都道府県
# A4103         = 合計特殊出生率
# A1302         = 15〜64歳人口（総数）
# A130202       = 15〜64歳人口（女）
# A1101         = 総人口
# A1303         = 65歳以上人口
# A9101         = 婚姻件数
# J2503         = 保育所等数
# L3221         = 消費支出（二人以上の世帯）

num_cols = ['A4103', 'A1302', 'A130202', 'A1101', 'A1303', 'A9101', 'J2503', 'L3221']
for c in num_cols:
    data[c] = pd.to_numeric(data[c], errors='coerce')
data['SSDSE-B-2026'] = pd.to_numeric(data['SSDSE-B-2026'], errors='coerce')

# 年度・都道府県・国全体コードを整理
data = data[data['Code'] != 'Z00000'].copy()   # 全国計除外
data = data.dropna(subset=num_cols + ['SSDSE-B-2026']).copy()
data = data[(data['SSDSE-B-2026'] >= 2012) & (data['SSDSE-B-2026'] <= 2023)].copy()

print("=" * 60)
print(f"■ 読み込み完了: {len(data)}件（{data['SSDSE-B-2026'].nunique()}年 × {data['Prefecture'].nunique()}都道府県）")
print("=" * 60)

# ── 特徴量エンジニアリング ────────────────────────────────────────
data['TFR']             = data['A4103'].astype(float)
data['女性就業率_代理'] = data['A130202'] / data['A1302'] * 100   # 15-64歳女性/15-64歳総人口
data['保育所密度']      = data['J2503'] / data['A1101'] * 10000   # 保育所等数/総人口×10000
data['高齢化率']        = data['A1303'] / data['A1101'] * 100     # 65歳以上/総人口×100
data['婚姻率']          = data['A9101'] / data['A1101'] * 1000    # 婚姻件数/総人口×1000
data['消費支出_log']    = np.log(data['L3221'].clip(lower=1))     # log(消費支出)

data['年度'] = data['SSDSE-B-2026'].astype(int)
data['都道府県'] = data['Prefecture']

vars_stat = ['TFR', '女性就業率_代理', '保育所密度', '高齢化率', '婚姻率', '消費支出_log']

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

region_order  = ['北海道・東北', '関東甲信越', '東海・北陸', '近畿', '中国・四国', '九州・沖縄']
region_colors = ['#1565C0', '#C62828', '#2E7D32', '#E65100', '#6A1B9A', '#00695C']

# 2022年断面データ（地域列を含む）
df22 = data[data['年度'] == 2022].copy()

print("\n【変数の記述統計（2022年断面）】")
print(df22[vars_stat].describe().round(3))

print("\n図1を生成中 ...")

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

tfr_region = (
    data[data['地域'].isin(region_order)]
    .groupby(['年度', '地域'])['TFR']
    .mean()
    .reset_index()
)

for reg, col in zip(region_order, region_colors):
    subset = tfr_region[tfr_region['地域'] == reg].sort_values('年度')
    ax1.plot(subset['年度'], subset['TFR'], marker='o', markersize=4,
             linewidth=1.8, color=col, label=reg, alpha=0.9)

ax1.set_xlabel('年度', fontsize=12)
ax1.set_ylabel('合計特殊出生率（TFR）', fontsize=12)
ax1.set_title('合計特殊出生率の時系列推移（地域別平均）\nSSDSE-B 2012〜2023年 実データ',
              fontsize=13, fontweight='bold')
ax1.legend(fontsize=9, loc='upper right', ncol=2)
ax1.set_xlim(2011.5, 2023.5)
ax1.set_xticks(range(2012, 2024))
ax1.xaxis.set_tick_params(rotation=45)
ax1.grid(True, alpha=0.3)
ax1.axvline(2020, color='gray', linestyle=':', linewidth=1.2, alpha=0.6)
ax1.text(2020.1, ax1.get_ylim()[0] + 0.03, 'COVID-19', fontsize=8, color='gray')

plt.tight_layout()
fig1_path = os.path.join(FIG_DIR, '2022_H5_10_fig1_timeseries.png')
fig1.savefig(fig1_path, bbox_inches='tight', dpi=150)
plt.close(fig1)
print(f"図1保存: {fig1_path}")

print("図2を生成中 ...")

fig2, ax2 = plt.subplots(figsize=(11, 7))

scatter_colors = [region_colors[region_order.index(reg)] if reg in region_order else '#888888'
                  for reg in df22_clean['地域']]

ax2.scatter(df22_clean['女性就業率_代理'], df22_clean['TFR'],
            c=scatter_colors, s=60, alpha=0.8, edgecolors='white', linewidths=0.5, zorder=3)

# 都道府県ラベル（略称）
label_map = {
    '北海道': '北海', '青森県': '青森', '岩手県': '岩手', '宮城県': '宮城', '秋田県': '秋田',
    '山形県': '山形', '福島県': '福島', '茨城県': '茨城', '栃木県': '栃木', '群馬県': '群馬',
    '埼玉県': '埼玉', '千葉県': '千葉', '東京都': '東京', '神奈川県': '神奈', '新潟県': '新潟',
    '富山県': '富山', '石川県': '石川', '福井県': '福井', '山梨県': '山梨', '長野県': '長野',
    '岐阜県': '岐阜', '静岡県': '静岡', '愛知県': '愛知', '三重県': '三重', '滋賀県': '滋賀',
    '京都府': '京都', '大阪府': '大阪', '兵庫県': '兵庫', '奈良県': '奈良', '和歌山県': '和歌',
    '鳥取県': '鳥取', '島根県': '島根', '岡山県': '岡山', '広島県': '広島', '山口県': '山口',
    '徳島県': '徳島', '香川県': '香川', '愛媛県': '愛媛', '高知県': '高知', '福岡県': '福岡',
    '佐賀県': '佐賀', '長崎県': '長崎', '熊本県': '熊本', '大分県': '大分', '宮崎県': '宮崎',
    '鹿児島県': '鹿児', '沖縄県': '沖縄',
}
for _, row in df22_clean.iterrows():
    short = label_map.get(row['都道府県'], row['都道府県'][:2])
    ax2.annotate(short,
                 xy=(row['女性就業率_代理'], row['TFR']),
                 xytext=(3, 3), textcoords='offset points',
                 fontsize=6.5, alpha=0.85)

# 回帰直線
x_emp = df22_clean['女性就業率_代理'].astype(float).values
y_tfr = df22_clean['TFR'].astype(float).values
r2, p2 = stats.pearsonr(x_emp, y_tfr)
z2 = np.polyfit(x_emp, y_tfr, 1)
xs2 = np.linspace(x_emp.min(), x_emp.max(), 100)
ax2.plot(xs2, np.poly1d(z2)(xs2), color='#C62828', linewidth=1.5,
         linestyle='--', alpha=0.8, label=f'回帰直線 r={r2:.3f} (p={p2:.3f})')

# 凡例（地域別）
from matplotlib.lines import Line2D
handles = [Line2D([0], [0], marker='o', color='w', markerfacecolor=col,
                  markersize=8, label=reg)
           for reg, col in zip(region_order, region_colors)]
handles.append(Line2D([0], [0], color='#C62828', linestyle='--', linewidth=1.5, label='回帰直線'))
ax2.legend(handles=handles, fontsize=8, loc='upper left')

ax2.set_xlabel('女性就業率代理（15〜64歳女性/15〜64歳総人口 ×100, %）', fontsize=11)
ax2.set_ylabel('合計特殊出生率（TFR）', fontsize=11)
ax2.set_title('女性就業率代理 vs 合計特殊出生率（2022年, 47都道府県）\nSSDSE-B 実データ',
              fontsize=12, fontweight='bold')
ax2.grid(True, alpha=0.3)

plt.tight_layout()
fig2_path = os.path.join(FIG_DIR, '2022_H5_10_fig2_scatter.png')
fig2.savefig(fig2_path, bbox_inches='tight', dpi=150)
plt.close(fig2)
print(f"図2保存: {fig2_path}")

print("図3を生成中 ...")

coef_names = ['女性就業率_代理', '保育所密度', '高齢化率', '婚姻率', '消費支出_log']
coef_labels = [
    '女性就業率代理\n(15-64歳女性比率)',
    '保育所密度\n(保育所数/万人)',
    '高齢化率\n(65歳以上/総人口)',
    '婚姻率\n(婚姻件数/千人)',
    '消費支出(log)\n(二人以上世帯)',
]

coefs3 = np.asarray(ols_model.params[1:])   # 定数項を除く
ses3   = np.asarray(ols_model.bse[1:])
pvals3 = np.asarray(ols_model.pvalues[1:])

bar_colors3 = []
for p in pvals3:
    if p < 0.01:
        bar_colors3.append('#C62828')    # 赤: p<0.01
    elif p < 0.05:
        bar_colors3.append('#FF8F00')    # オレンジ: p<0.05
    elif p < 0.1:
        bar_colors3.append('#1565C0')    # 青: p<0.10
    else:
        bar_colors3.append('#9E9E9E')    # 灰: n.s.

fig3, ax3 = plt.subplots(figsize=(9, 5))
y_pos3 = np.arange(len(coef_names))
bars3 = ax3.barh(y_pos3, coefs3, color=bar_colors3, alpha=0.80,
                 edgecolor='white', height=0.6)
ax3.errorbar(coefs3, y_pos3, xerr=1.96 * ses3,
             fmt='none', color='#333333', capsize=4, linewidth=1.2, zorder=5)
ax3.axvline(0, color='gray', linestyle='--', linewidth=1.0)
ax3.set_yticks(y_pos3)
ax3.set_yticklabels(coef_labels, fontsize=10)
ax3.set_xlabel('標準化なし回帰係数（±1.96 SE）', fontsize=11)
ax3.set_title(f'OLS回帰係数プロット（目的変数: TFR）\nR²={ols_model.rsquared:.3f}, N=47都道府県（2022年）',
              fontsize=12, fontweight='bold')
ax3.invert_yaxis()
ax3.grid(axis='x', alpha=0.3)

from matplotlib.patches import Patch
legend_elements = [
    Patch(facecolor='#C62828', alpha=0.80, label='p < 0.01'),
    Patch(facecolor='#FF8F00', alpha=0.80, label='p < 0.05'),
    Patch(facecolor='#1565C0', alpha=0.80, label='p < 0.10'),
    Patch(facecolor='#9E9E9E', alpha=0.80, label='n.s.'),
]
ax3.legend(handles=legend_elements, fontsize=9, loc='lower right')

# p値注記
for i, (coef, se, pval) in enumerate(zip(coefs3, ses3, pvals3)):
    sig = '***' if pval < 0.001 else '**' if pval < 0.01 else '*' if pval < 0.05 else '†' if pval < 0.10 else ''
    offset = 0.003 if coef >= 0 else -0.003
    ax3.text(coef + 1.96 * se + offset, i, f' {sig}',
             va='center', ha='left', fontsize=10, color='black')

plt.tight_layout()
fig3_path = os.path.join(FIG_DIR, '2022_H5_10_fig3_coef.png')
fig3.savefig(fig3_path, bbox_inches='tight', dpi=150)
plt.close(fig3)
print(f"図3保存: {fig3_path}")

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
import warnings
warnings.filterwarnings('ignore')

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

import os
FIG_DIR = 'html/figures'
DATA_B  = 'data/raw/SSDSE-B-2026.csv'
os.makedirs(FIG_DIR, exist_ok=True)

print("\n" + "=" * 60)
print("■ OLS回帰（2022年断面, N=47都道府県）")
print("=" * 60)

df22_clean = df22.dropna(subset=vars_stat).copy()
X_vars = ['女性就業率_代理', '保育所密度', '高齢化率', '婚姻率', '消費支出_log']
X_ols = sm.add_constant(df22_clean[X_vars].astype(float))
y_ols = df22_clean['TFR'].astype(float)

ols_model = sm.OLS(y_ols, X_ols).fit()
print(ols_model.summary2())

# 相関行列（2022年）
print("\n【相関行列（2022年, N=47）】")
corr_vars = ['TFR', '女性就業率_代理', '保育所密度', '婚姻率', '高齢化率']
corr_df = df22_clean[corr_vars].astype(float).corr()
print(corr_df.round(3))

print("図4を生成中 ...")

heat_vars  = ['TFR', '女性就業率_代理', '保育所密度', '婚姻率', '高齢化率']
heat_labels = [
    'TFR\n(合計特殊出生率)',
    '女性就業率代理\n(15-64歳)',
    '保育所密度\n(保育所/万人)',
    '婚姻率\n(件/千人)',
    '高齢化率\n(65歳以上%)',
]

corr_mat = df22_clean[heat_vars].astype(float).corr()

fig4, ax4 = plt.subplots(figsize=(8, 6.5))
im = ax4.imshow(corr_mat.values, cmap='RdBu_r', vmin=-1, vmax=1, aspect='auto')

cbar = plt.colorbar(im, ax=ax4, fraction=0.046, pad=0.04)
cbar.set_label('相関係数', fontsize=10)

n_vars4 = len(heat_vars)
ax4.set_xticks(range(n_vars4))
ax4.set_yticks(range(n_vars4))
ax4.set_xticklabels(heat_labels, fontsize=8.5, rotation=0, ha='center')
ax4.set_yticklabels(heat_labels, fontsize=8.5)

# 相関係数と有意星を表示
for i in range(n_vars4):
    for j in range(n_vars4):
        val = corr_mat.iloc[i, j]
        if i == j:
            ax4.text(j, i, '1.00', ha='center', va='center',
                     fontsize=10, fontweight='bold', color='black')
        else:
            # p値計算
            x_c = df22_clean[heat_vars[j]].astype(float).values
            y_c = df22_clean[heat_vars[i]].astype(float).values
            mask_c = np.isfinite(x_c) & np.isfinite(y_c)
            _, pval_c = stats.pearsonr(x_c[mask_c], y_c[mask_c])
            sig_c = '***' if pval_c < 0.001 else '**' if pval_c < 0.01 else '*' if pval_c < 0.05 else ''
            text_color = 'white' if abs(val) > 0.6 else 'black'
            ax4.text(j, i, f'{val:.2f}{sig_c}',
                     ha='center', va='center', fontsize=9,
                     color=text_color, fontweight='bold' if abs(val) > 0.3 else 'normal')

ax4.set_title('相関ヒートマップ（2022年, N=47都道府県）\n* p<0.05, ** p<0.01, *** p<0.001',
              fontsize=12, fontweight='bold', pad=14)

plt.tight_layout()
fig4_path = os.path.join(FIG_DIR, '2022_H5_10_fig4_heatmap.png')
fig4.savefig(fig4_path, bbox_inches='tight', dpi=150)
plt.close(fig4)
print(f"図4保存: {fig4_path}")

print("\n" + "=" * 60)
print("■ 分析結果サマリー")
print("=" * 60)
print(f"  使用データ: SSDSE-B-2026（実データ）")
print(f"  対象期間: 2012〜2023年, 47都道府県パネル")
print()
print("  OLS回帰結果（2022年断面, N=47）:")
print(f"  {'変数':<20} {'係数':>8} {'SE':>8} {'p値':>10} {'有意'}")
print("  " + "-" * 55)
for name, coef, se, pval in zip(X_vars, coefs3, ses3, pvals3):
    sig = '***' if pval < 0.001 else '**' if pval < 0.01 else '*' if pval < 0.05 else '†' if pval < 0.10 else 'n.s.'
    print(f"  {name:<20} {coef:>8.4f} {se:>8.4f} {pval:>10.4f}  {sig}")
print(f"  R² = {ols_model.rsquared:.4f}, N=47")
print()
print("  相関分析（2022年, TFRとの相関）:")
for var in heat_vars[1:]:
    x_v = df22_clean[var].astype(float).values
    y_v = df22_clean['TFR'].astype(float).values
    mask_v = np.isfinite(x_v) & np.isfinite(y_v)
    r_v, p_v = stats.pearsonr(x_v[mask_v], y_v[mask_v])
    sig_v = '***' if p_v < 0.001 else '**' if p_v < 0.01 else '*' if p_v < 0.05 else 'n.s.'
    print(f"  TFR ↔ {var:<18}: r={r_v:+.3f} ({sig_v})")

print("\n全図の生成完了（4枚）")
print(f"  fig1: 2022_H5_10_fig1_timeseries.png  (TFR時系列推移)")
print(f"  fig2: 2022_H5_10_fig2_scatter.png     (女性就業率 vs TFR散布図)")
print(f"  fig3: 2022_H5_10_fig3_coef.png        (OLS回帰係数プロット)")
print(f"  fig4: 2022_H5_10_fig4_heatmap.png     (相関ヒートマップ)")