"""
2021_H5_2_shorei.py
「少子化と女性就業・育児支援の関係：都道府県パネルデータ分析」
2021年度 統計データ分析コンペ 審査員奨励賞（高校生の部）
教育用再現スクリプト — SSDSE-B-2026.csv 実データ使用
合成データ（np.random.seed 等）は一切使用しない
"""

# ============================================================
# 【データの準備】実行前に以下のデータファイルを用意してください
#
#   必要ファイル:
#     ・SSDSE-B-2026.csv
#       → data/raw/SSDSE-B-2026.csv に配置
#
#   ダウンロード先:
#     https://www.nstac.go.jp/use/literacy/ssdse/
#     （SSDSE-B（社会・人口統計体系 都道府県データ） の CSV をダウンロード）
#
#   フォルダ配置（プロジェクトルートからの相対パス）:
#     code/                ← このスクリプトの場所
#     data/raw/            ← CSV ファイルをここに配置
#     html/figures/        ← 図の出力先（自動生成）
#
#   実行方法（ファイルを一切編集せず実行可能）:
#     python3 code/2021_H5_2_shorei.py
# ============================================================


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)

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

# ── 派生変数作成（全年度）───────────────────────────────────────────────────
df_b['地域区分']    = df_b['都道府県'].map(region_map)
df_b['婚姻率']      = df_b['婚姻件数']       / df_b['総人口']        * 1000
df_b['保育所数_千人'] = df_b['保育所等数']    / df_b['15歳未満人口']  * 1000
df_b['保育定員率']  = df_b['保育所等定員数']  / df_b['15歳未満人口']  * 100
df_b['高齢化率']    = df_b['65歳以上人口']    / df_b['総人口']        * 100
df_b['就職率']      = df_b['就職件数（一般）'] / df_b['総人口']        * 10000

# ── 2022年断面データ ─────────────────────────────────────────────────────────
df22 = df_b[df_b['年度'] == 2022].copy()

# ============================================================
# 図1: 地域別 TFR 時系列推移 2012〜2023
# ============================================================
fig1, ax1 = plt.subplots(figsize=(10, 6))

for reg in region_order:
    sub = df_b[df_b['地域区分'] == reg].groupby('年度')['合計特殊出生率'].mean()
    ax1.plot(sub.index, sub.values,
             color=region_colors[reg], marker='o', markersize=4,
             linewidth=2, label=reg)

ax1.axhline(2.07, color='gray', linestyle='--', linewidth=1.2, label='人口置換水準（2.07）')
ax1.set_xlabel('年度', fontsize=12)
ax1.set_ylabel('合計特殊出生率（TFR）', fontsize=12)
ax1.set_title('地域別 合計特殊出生率の推移（2012〜2023年）', fontsize=14, fontweight='bold')
ax1.set_xticks(range(2012, 2024))
ax1.tick_params(axis='x', rotation=45)
ax1.legend(fontsize=10, loc='upper right')
ax1.set_ylim(0.9, 2.3)
ax1.grid(axis='y', alpha=0.3)
fig1.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2021_H5_2_fig1.png'))
plt.close(fig1)
print('fig1 saved')

# ============================================================
# 図2: 散布図 — 保育定員率 vs TFR（地域色分け・都道府県ラベル・回帰線）
# 保育定員率 = 保育所等定員数 / 15歳未満人口 × 100（育児支援の充実度）
# ============================================================
fig2, ax2 = plt.subplots(figsize=(10, 8))

for reg in region_order:
    sub = df22[df22['地域区分'] == reg]
    ax2.scatter(sub['保育定員率'], sub['合計特殊出生率'],
                color=region_colors[reg], s=65, zorder=3, label=reg, alpha=0.85)
    for _, row in sub.iterrows():
        pref = row['都道府県'].replace('県', '').replace('府', '').replace('都', '').replace('道', '')
        ax2.annotate(pref,
                     (row['保育定員率'], row['合計特殊出生率']),
                     fontsize=7.5, ha='left', va='bottom',
                     xytext=(2, 2), textcoords='offset points', color='#444')

# 回帰直線
sub_clean = df22[['保育定員率', '合計特殊出生率']].dropna()
x_vals = sub_clean['保育定員率'].values
y_vals = sub_clean['合計特殊出生率'].values
slope, intercept, r_val2, p_val2, _ = stats.linregress(x_vals, y_vals)
x_range = np.linspace(x_vals.min(), x_vals.max(), 100)
label_r = f'回帰線 (r={r_val2:.3f}, p<0.01)' if p_val2 < 0.01 else f'回帰線 (r={r_val2:.3f}, p={p_val2:.3f})'
ax2.plot(x_range, intercept + slope * x_range,
         color='#333', linewidth=1.5, linestyle='-', zorder=2, label=label_r)

ax2.set_xlabel('保育定員率（保育所等定員数 / 15歳未満人口 × 100）', fontsize=11)
ax2.set_ylabel('合計特殊出生率（TFR）', fontsize=12)
ax2.set_title('保育定員率と合計特殊出生率の関係（2022年、47都道府県）', fontsize=13, fontweight='bold')
ax2.legend(fontsize=9, loc='upper left')
ax2.grid(alpha=0.25)
fig2.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2021_H5_2_fig2.png'))
plt.close(fig2)
print('fig2 saved')

# ============================================================
# 図3: 相関行列ヒートマップ
# ============================================================
vars_corr = {
    '合計特殊出生率': df22['合計特殊出生率'],
    '保育定員率':     df22['保育定員率'],
    '保育所数/千人':   df22['保育所数_千人'],
    '就職率':         df22['就職率'],
    '高齢化率':       df22['高齢化率'],
    '消費支出(万円)':  df22['消費支出（二人以上の世帯）'] / 10000,
}
df_corr = pd.DataFrame(vars_corr).dropna()
corr_mat = df_corr.corr()
n_var = len(corr_mat)

# p値行列
pval_mat = pd.DataFrame(np.ones((n_var, n_var)),
                         index=corr_mat.index, columns=corr_mat.columns)
for i, c1 in enumerate(df_corr.columns):
    for j, c2 in enumerate(df_corr.columns):
        if i != j:
            _, p = stats.pearsonr(df_corr[c1], df_corr[c2])
            pval_mat.iloc[i, j] = p

fig3, ax3 = plt.subplots(figsize=(8, 7))
cmap = plt.cm.RdBu_r
im = ax3.imshow(corr_mat.values, cmap=cmap, vmin=-1, vmax=1, aspect='auto')
plt.colorbar(im, ax=ax3, shrink=0.8, label='Pearson r')

ax3.set_xticks(range(n_var))
ax3.set_yticks(range(n_var))
ax3.set_xticklabels(corr_mat.columns, rotation=30, ha='right', fontsize=10)
ax3.set_yticklabels(corr_mat.index, fontsize=10)

for i in range(n_var):
    for j in range(n_var):
        r = corr_mat.iloc[i, j]
        p = pval_mat.iloc[i, j]
        star = '***' if p < 0.001 else ('**' if p < 0.01 else ('*' if p < 0.05 else ''))
        txt = f'{r:.2f}{star}'
        color = 'white' if abs(r) > 0.5 else 'black'
        ax3.text(j, i, txt, ha='center', va='center', fontsize=9, color=color)

ax3.set_title('相関行列（2022年断面、N=47）\n*** p<0.001  ** p<0.01  * p<0.05', fontsize=12, fontweight='bold')
fig3.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2021_H5_2_fig3.png'))
plt.close(fig3)
print('fig3 saved')

# ============================================================
# 図4: 標準化偏回帰係数プロット（OLS）
# ============================================================
X_vars = ['保育定員率', '保育所数/千人', '就職率', '高齢化率', '消費支出(万円)']
df_ols = df_corr[['合計特殊出生率'] + X_vars].dropna().copy()
df_std = (df_ols - df_ols.mean()) / df_ols.std()

y_std = df_std['合計特殊出生率'].values
X_std = sm.add_constant(df_std[X_vars].values)
model = sm.OLS(y_std, X_std).fit()

coefs = np.array(model.params[1:])   # 定数除く
pvals = np.array(model.pvalues[1:])
conf_arr = model.conf_int()
if hasattr(conf_arr, 'values'):
    conf_arr = conf_arr.values
conf  = conf_arr[1:]  # shape (n_var, 2)

# ソート（係数の絶対値で昇順）
order = np.argsort(np.abs(coefs))
coefs_s = coefs[order]
pvals_s = pvals[order]
conf_s  = conf[order]
labels_s = [X_vars[i] for i in order]
colors_s = ['#e05c5c' if p < 0.05 else '#aaaaaa' for p in pvals_s]

fig4, ax4 = plt.subplots(figsize=(8, 5))
y_pos = np.arange(len(labels_s))
ax4.barh(y_pos, coefs_s, color=colors_s, edgecolor='white', height=0.55)

# 信頼区間（エラーバー）
for i, (ci_lo, ci_hi) in enumerate(zip(conf_s[:, 0], conf_s[:, 1])):
    ax4.plot([ci_lo, ci_hi], [i, i], color='#555', linewidth=2.5)
    ax4.plot([ci_lo], [i], marker='|', color='#555', markersize=8)
    ax4.plot([ci_hi], [i], marker='|', color='#555', markersize=8)

ax4.axvline(0, color='black', linewidth=0.8)
ax4.set_yticks(y_pos)
ax4.set_yticklabels(labels_s, fontsize=11)
ax4.set_xlabel('標準化偏回帰係数 (beta)', fontsize=12)
ax4.set_title('重回帰分析：標準化偏回帰係数（2022年断面、N=47）\n'
              '赤=p<0.05（有意）、灰色=p>=0.05（非有意）', fontsize=11, fontweight='bold')

# p値テキスト
for i, (c, p) in enumerate(zip(coefs_s, pvals_s)):
    star = '***' if p < 0.001 else ('**' if p < 0.01 else ('*' if p < 0.05 else 'n.s.'))
    offset = 0.015 if c >= 0 else -0.015
    ha = 'left' if c >= 0 else 'right'
    ax4.text(c + offset, i, star, ha=ha, va='center', fontsize=10)

ax4.grid(axis='x', alpha=0.25)
fig4.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2021_H5_2_fig4.png'))
plt.close(fig4)
print('fig4 saved')

# ============================================================
# 分析結果の出力（HTMLへ反映）
# ============================================================
print('\n===== 分析結果サマリー =====')
print(f'2022年 TFR 平均: {df22["合計特殊出生率"].mean():.3f}')
print(f'2022年 TFR 最小: {df22["合計特殊出生率"].min():.2f}'
      f' ({df22.loc[df22["合計特殊出生率"].idxmin(), "都道府県"]})')
print(f'2022年 TFR 最大: {df22["合計特殊出生率"].max():.2f}'
      f' ({df22.loc[df22["合計特殊出生率"].idxmax(), "都道府県"]})')

print(f'\n保育定員率 vs TFR: r={r_val2:.3f}, p={p_val2:.4f}')
sub_c = df22[['保育所数_千人','合計特殊出生率']].dropna()
r_n, p_n = stats.pearsonr(sub_c['保育所数_千人'], sub_c['合計特殊出生率'])
print(f'保育所数/千人 vs TFR: r={r_n:.3f}, p={p_n:.4f}')
sub_e = df22[['就職率','合計特殊出生率']].dropna()
r_e, p_e = stats.pearsonr(sub_e['就職率'], sub_e['合計特殊出生率'])
print(f'就職率 vs TFR: r={r_e:.3f}, p={p_e:.4f}')

print()
print('--- 重回帰分析 (標準化) ---')
print(f'R-squared: {model.rsquared:.4f}')
print(f'Adjusted R-squared: {model.rsquared_adj:.4f}')
for var, coef, pv in zip(X_vars, model.params[1:], model.pvalues[1:]):
    star = '***' if pv < 0.001 else ('**' if pv < 0.01 else ('*' if pv < 0.05 else 'n.s.'))
    print(f'  {var}: beta={coef:.3f}, p={pv:.4f} {star}')

print()
print('--- 地域別TFR 2022 ---')
for reg in region_order:
    sub = df22[df22['地域区分'] == reg]
    print(f'  {reg}: {sub["合計特殊出生率"].mean():.3f}')

tfr2012 = df_b[df_b['年度'] == 2012]['合計特殊出生率'].mean()
tfr2022 = df_b[df_b['年度'] == 2022]['合計特殊出生率'].mean()
print(f'\n全国平均 TFR 2012: {tfr2012:.3f}')
print(f'全国平均 TFR 2022: {tfr2022:.3f}')
print(f'変化: {tfr2022 - tfr2012:+.3f}')

print()
print('--- 上位・下位5都道府県 保育定員率 ---')
top5 = df22.nlargest(5, '保育定員率')[['都道府県','保育定員率','合計特殊出生率']]
bot5 = df22.nsmallest(5, '保育定員率')[['都道府県','保育定員率','合計特殊出生率']]
print('保育定員率 高い:')
print(top5.to_string(index=False))
print('保育定員率 低い:')
print(bot5.to_string(index=False))