"""
教育用再現コード: 2022年 統計データ分析コンペティション 審査員奨励賞 [高校生の部]
=======================================================================
論文タイトル：少子化と保育環境：都道府県別出生率と保育所整備の関係分析
受賞：審査員奨励賞（高校生の部）

【分析概要】
  データ：SSDSE-B-2026.csv（都道府県別データ, 2012〜2023年度）

  変数設計：
    TFR（合計特殊出生率）= 目的変数
    保育所密度 = 保育所等数 / 総人口 × 10000（人口1万人あたり保育所数）
    保育所待機児童率 = 保育所等利用待機児童数 / 保育所等定員数 × 100
    高齢化率 = 65歳以上人口 / 総人口 × 100
    女性就業代理 = 15〜64歳人口（女）/ 15〜64歳人口 × 100
    消費支出_log = log(消費支出（二人以上の世帯）)

  分析の流れ：
    Fig1: TFRの都道府県別時系列推移（地域別色分け）
    Fig2: 保育所密度 vs TFRの散布図（2022年、都道府県ラベル付き）
    Fig3: OLS回帰係数プロット（TFRの決定要因）
    Fig4: 保育所待機児童率の都道府県別ランキング（2022年）

【データ出典】
  SSDSE-B-2026.csv: 社会・人口統計体系データセット（都道府県データ）
  提供元：統計数理研究所 https://www.ism.ac.jp/

【データサイエンス学習ポイント】
  1. 相関分析：ピアソン相関係数とp値による変数間関係の定量化
  2. OLS重回帰：statsmodelsによる最小二乗法と係数の解釈
  3. 代理変数設計：直接観測できない概念を測定可能な変数で代替する手法
  4. 時系列・パネルデータ：複数年・複数地域の縦断データの可視化
=======================================================================
"""

# ============================================================
# 【データの準備】実行前に以下のデータファイルを用意してください
#
#   必要ファイル:
#     ・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/2022_H5_1_shorei.py
# ============================================================


import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from scipy import stats
import statsmodels.api as sm
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)

# ── データ読み込み ──────────────────────────────────────────────────────────
df_raw = pd.read_csv(DATA_B, encoding='cp932', header=1)

# 都道府県行のみ抽出 (地域コード R#####)
df = df_raw[df_raw['地域コード'].str.match(r'^R\d{5}$', na=False)].copy()
df['年度'] = df['年度'].astype(int)

# ── 変数構築 ──────────────────────────────────────────────────────────────
df['TFR'] = pd.to_numeric(df['合計特殊出生率'], errors='coerce')
df['保育所密度'] = (
    pd.to_numeric(df['保育所等数'], errors='coerce') /
    pd.to_numeric(df['総人口'], errors='coerce') * 10000
)
df['保育所待機児童率'] = (
    pd.to_numeric(df['保育所等利用待機児童数'], errors='coerce') /
    pd.to_numeric(df['保育所等定員数'], errors='coerce') * 100
)
df['高齢化率'] = (
    pd.to_numeric(df['65歳以上人口'], errors='coerce') /
    pd.to_numeric(df['総人口'], errors='coerce') * 100
)
df['女性就業代理'] = (
    pd.to_numeric(df['15～64歳人口（女）'], errors='coerce') /
    pd.to_numeric(df['15～64歳人口'], errors='coerce') * 100
)
df['消費支出_log'] = np.log(
    pd.to_numeric(df['消費支出（二人以上の世帯）'], errors='coerce')
)

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

region_colors = {
    '北海道':   '#1565C0',
    '東北':     '#00838F',
    '関東':     '#E65100',
    '中部':     '#558B2F',
    '近畿':     '#6A1B9A',
    '中国':     '#AD1457',
    '四国':     '#37474F',
    '九州・沖縄': '#F9A825',
}

# ─────────────────────────────────────────────────────────────────────────────
# Fig 1: TFRの都道府県別時系列推移（地域別色分け）
# ─────────────────────────────────────────────────────────────────────────────
fig1, ax1 = plt.subplots(figsize=(12, 7))

for region, color in region_colors.items():
    prefs_in_region = [p for p, r in region_map.items() if r == region]
    for pref in prefs_in_region:
        subset = df[df['都道府県'] == pref].sort_values('年度')
        if subset.empty or subset['TFR'].isna().all():
            continue
        ax1.plot(subset['年度'], subset['TFR'], color=color, alpha=0.55,
                 linewidth=1.2, zorder=2)

# 全国平均（各年の中央値）
yearly_med = df.groupby('年度')['TFR'].median()
ax1.plot(yearly_med.index, yearly_med.values, color='black', linewidth=2.5,
         linestyle='--', zorder=5, label='全国中央値')

# 凡例（地方区分）
from matplotlib.lines import Line2D
legend_handles = [
    Line2D([0], [0], color=c, linewidth=2.5, label=r)
    for r, c in region_colors.items()
]
legend_handles.append(
    Line2D([0], [0], color='black', linewidth=2.5, linestyle='--', label='全国中央値')
)
ax1.legend(handles=legend_handles, loc='lower left', fontsize=9,
           ncol=2, framealpha=0.85)

# 注目年をハイライト
ax1.axvline(2022, color='red', linewidth=1.2, linestyle=':', alpha=0.7, zorder=3)
ax1.text(2022.05, ax1.get_ylim()[1] * 0.98, '2022年', color='red',
         fontsize=9, va='top')

ax1.set_xlabel('年度', fontsize=12)
ax1.set_ylabel('合計特殊出生率（TFR）', fontsize=12)
ax1.set_title('Fig 1: 都道府県別 合計特殊出生率の時系列推移（地域別色分け）',
              fontsize=13, fontweight='bold', pad=12)
ax1.set_xlim(2012, 2023)
ax1.set_xticks(range(2012, 2024))
ax1.grid(axis='y', linestyle='--', alpha=0.4)
ax1.spines['top'].set_visible(False)
ax1.spines['right'].set_visible(False)

plt.tight_layout()
fig1_path = os.path.join(FIG_DIR, '2022_H5_1_fig1_tfr_timeseries.png')
fig1.savefig(fig1_path, dpi=150, bbox_inches='tight')
plt.close(fig1)
print(f'Saved: {fig1_path}')

# ─────────────────────────────────────────────────────────────────────────────
# Fig 2: 保育所密度 vs TFR の散布図（2022年、都道府県ラベル付き）
# ─────────────────────────────────────────────────────────────────────────────
df22 = df[df['年度'] == 2022].copy().dropna(subset=['TFR', '保育所密度'])

fig2, ax2 = plt.subplots(figsize=(12, 8))

for region, color in region_colors.items():
    mask = df22['地方'] == region
    ax2.scatter(df22.loc[mask, '保育所密度'], df22.loc[mask, 'TFR'],
                color=color, s=60, alpha=0.85, zorder=3, label=region,
                edgecolors='white', linewidths=0.5)

# 回帰直線
x_arr = df22['保育所密度'].values
y_arr = df22['TFR'].values
slope, intercept, r_val, p_val, _ = stats.linregress(x_arr, y_arr)
x_line = np.linspace(x_arr.min(), x_arr.max(), 200)
ax2.plot(x_line, intercept + slope * x_line, color='#C62828', linewidth=1.8,
         linestyle='--', zorder=2,
         label=f'回帰直線 (r={r_val:.3f}, p={p_val:.3f})')

# 都道府県ラベル（全47都道府県）
for _, row in df22.iterrows():
    label = row['都道府県'].replace('県', '').replace('府', '').replace('都', '').replace('道', '')
    ax2.annotate(label,
                 xy=(row['保育所密度'], row['TFR']),
                 xytext=(3, 3), textcoords='offset points',
                 fontsize=6.5, color='#333333', zorder=4)

ax2.set_xlabel('保育所密度（保育所等数 / 総人口 × 10,000）', fontsize=12)
ax2.set_ylabel('合計特殊出生率（TFR）', fontsize=12)
ax2.set_title('Fig 2: 保育所密度 vs 合計特殊出生率（2022年、都道府県別）',
              fontsize=13, fontweight='bold', pad=12)
ax2.legend(fontsize=9, loc='upper right', framealpha=0.85)
ax2.grid(linestyle='--', alpha=0.35)
ax2.spines['top'].set_visible(False)
ax2.spines['right'].set_visible(False)

# 統計情報を右下に表示
info_text = f'N=47都道府県\nr = {r_val:.3f}\np = {p_val:.3f}'
ax2.text(0.02, 0.03, info_text, transform=ax2.transAxes,
         fontsize=10, verticalalignment='bottom',
         bbox=dict(boxstyle='round', facecolor='#FFF9C4', alpha=0.85))

plt.tight_layout()
fig2_path = os.path.join(FIG_DIR, '2022_H5_1_fig2_scatter_nursery_tfr.png')
fig2.savefig(fig2_path, dpi=150, bbox_inches='tight')
plt.close(fig2)
print(f'Saved: {fig2_path}')

# ─────────────────────────────────────────────────────────────────────────────
# Fig 3: OLS回帰係数プロット（TFRの決定要因）
# ─────────────────────────────────────────────────────────────────────────────
df22_reg = df22.dropna(subset=['TFR', '保育所密度', '保育所待機児童率',
                                '高齢化率', '女性就業代理', '消費支出_log'])

VARS = {
    '保育所密度\n(人口1万対)':    '保育所密度',
    '保育所\n待機児童率(%)':      '保育所待機児童率',
    '高齢化率(%)':               '高齢化率',
    '女性就業\n代理指標(%)':      '女性就業代理',
    '消費支出\n(log)':            '消費支出_log',
}

X_cols = list(VARS.values())
X_data = df22_reg[X_cols].copy()

# 標準化（係数の比較のため）
X_std = (X_data - X_data.mean()) / X_data.std()
X_sm  = sm.add_constant(X_std)
y_sm  = df22_reg['TFR'].values

ols_model = sm.OLS(y_sm, X_sm).fit()

coef_df = pd.DataFrame({
    '変数':  list(VARS.keys()),
    '係数':  ols_model.params[1:],
    'se':    ols_model.bse[1:],
    'p値':   ols_model.pvalues[1:],
})
coef_df = coef_df.sort_values('係数')

colors_bar = ['#C62828' if c < 0 else '#1565C0' for c in coef_df['係数']]
sig_alpha  = [1.0 if p < 0.05 else 0.45 for p in coef_df['p値']]

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

bars = ax3.barh(coef_df['変数'], coef_df['係数'],
                color=colors_bar, alpha=0.85, edgecolor='white',
                linewidth=0.8, height=0.55, zorder=3)

# 透明度によるp値表示（アルファを個別設定は困難なのでエラーバーで表現）
ax3.errorbar(coef_df['係数'], coef_df['変数'],
             xerr=1.96 * coef_df['se'],
             fmt='none', color='#333333', capsize=5, linewidth=1.5, zorder=4)

# 有意ラベル
for i, (_, row) in enumerate(coef_df.iterrows()):
    sig = '***' if row['p値'] < 0.001 else ('**' if row['p値'] < 0.01
          else ('*' if row['p値'] < 0.05 else 'n.s.'))
    x_pos = row['係数'] + 1.96 * row['se'] + 0.003
    ax3.text(x_pos, row['変数'], sig, va='center', fontsize=10,
             color='#C62828' if 'n.s.' not in sig else '#888888')

ax3.axvline(0, color='black', linewidth=1.0, zorder=2)
ax3.set_xlabel('標準化回帰係数（95%CI）', fontsize=12)
ax3.set_title(f'Fig 3: TFRの決定要因 — OLS回帰係数プロット（2022年、N={len(df22_reg)}）',
              fontsize=13, fontweight='bold', pad=12)

r2_text = f'R² = {ols_model.rsquared:.3f}\nAdj.R² = {ols_model.rsquared_adj:.3f}\nN = {len(df22_reg)}'
ax3.text(0.97, 0.04, r2_text, transform=ax3.transAxes,
         fontsize=9.5, ha='right', va='bottom',
         bbox=dict(boxstyle='round', facecolor='#E3F2FD', alpha=0.9))

ax3.grid(axis='x', linestyle='--', alpha=0.4, zorder=1)
ax3.spines['top'].set_visible(False)
ax3.spines['right'].set_visible(False)

from matplotlib.patches import Patch
legend_elements = [
    Patch(facecolor='#1565C0', label='正の効果（TFR増加）'),
    Patch(facecolor='#C62828', label='負の効果（TFR減少）'),
]
ax3.legend(handles=legend_elements, loc='lower right', fontsize=9)

plt.tight_layout()
fig3_path = os.path.join(FIG_DIR, '2022_H5_1_fig3_ols_coef.png')
fig3.savefig(fig3_path, dpi=150, bbox_inches='tight')
plt.close(fig3)
print(f'Saved: {fig3_path}')
print(ols_model.summary())

# ─────────────────────────────────────────────────────────────────────────────
# Fig 4: 保育所待機児童率の都道府県別ランキング（2022年）
# ─────────────────────────────────────────────────────────────────────────────
df22_wait = df22.dropna(subset=['保育所待機児童率']).copy()
df22_wait = df22_wait.sort_values('保育所待機児童率', ascending=True)

# TFRに基づいて棒の色を決定
tfr_mid = df22_wait['TFR'].median()
bar_colors = ['#E65100' if t >= tfr_mid else '#546E7A'
              for t in df22_wait['TFR'].values]

fig4, ax4 = plt.subplots(figsize=(10, 13))

bars4 = ax4.barh(df22_wait['都道府県'], df22_wait['保育所待機児童率'],
                 color=bar_colors, alpha=0.85, edgecolor='white',
                 linewidth=0.5, height=0.72, zorder=3)

# 値ラベル
for bar, val in zip(bars4, df22_wait['保育所待機児童率']):
    if val > 0:
        ax4.text(val + 0.005, bar.get_y() + bar.get_height() / 2,
                 f'{val:.3f}%', va='center', fontsize=7.5, color='#212121')

ax4.set_xlabel('保育所待機児童率（待機児童数 / 定員数 × 100, %）', fontsize=11)
ax4.set_title('Fig 4: 保育所待機児童率 都道府県別ランキング（2022年）',
              fontsize=13, fontweight='bold', pad=12)
ax4.grid(axis='x', linestyle='--', alpha=0.4, zorder=1)
ax4.spines['top'].set_visible(False)
ax4.spines['right'].set_visible(False)

from matplotlib.patches import Patch
legend_elems4 = [
    Patch(facecolor='#E65100', label=f'TFR ≥ 中央値({tfr_mid:.2f})'),
    Patch(facecolor='#546E7A', label=f'TFR < 中央値({tfr_mid:.2f})'),
]
ax4.legend(handles=legend_elems4, loc='lower right', fontsize=9, framealpha=0.85)

plt.tight_layout()
fig4_path = os.path.join(FIG_DIR, '2022_H5_1_fig4_waitlist_ranking.png')
fig4.savefig(fig4_path, dpi=150, bbox_inches='tight')
plt.close(fig4)
print(f'Saved: {fig4_path}')

print('\n=== 完了 ===')
print(f'Fig1: {fig1_path}')
print(f'Fig2: {fig2_path}')
print(f'Fig3: {fig3_path}')
print(f'Fig4: {fig4_path}')