"""
教育用再現コード: 2025年 統計データ分析コンペティション 審査員奨励賞 [高校生の部]
=================================================================
論文タイトル：合計特殊出生率の変動要因分析
著者：京華高等学校

【データ】
  SSDSE-B-2026.csv（2022年度 47都道府県パネルデータ）
  SSDSE-E-2026.csv（1人当たり県民所得, 総面積）

【分析概要】
  目的変数：TFR（合計特殊出生率, 2022年）
  Step1. 相関分析（TFR vs 各説明変数）
  Step2. 重回帰分析（全変数 → 有意変数の確認）
  Step3. 3区分（東京都・過疎化進行県・その他）での比較

【変数】
  婚姻率, 1人あたり県民所得, 高齢化率, ln(人口密度),
  保育所数（千対）, 大学進学率（高校卒業者のうち進学者比率）, 年平均気温
=================================================================
"""

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


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
from sklearn.preprocessing import StandardScaler
from matplotlib.patches import Patch
import warnings
warnings.filterwarnings('ignore')
import os

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

BASE_DIR = os.path.join(_script_dir, '..')
FIG_DIR = os.path.join(BASE_DIR, 'html', 'figures')
DATA_DIR = os.path.join(BASE_DIR, 'data', 'raw')

# ================================================================
# ■ データ読み込み
# ================================================================
# SSDSE-B: 都道府県別パネルデータ（2022年）
ssdse_b = pd.read_csv(os.path.join(DATA_DIR, 'SSDSE-B-2026.csv'), encoding='cp932', header=1)
df_b = ssdse_b[
    (ssdse_b['年度'] == 2022) &
    ssdse_b['地域コード'].str.match(r'^R\d{5}$', na=False)
].copy().reset_index(drop=True)

# SSDSE-E: 都道府県別クロスセクション（1人当たり県民所得, 総面積）
ssdse_e_raw = pd.read_csv(os.path.join(DATA_DIR, 'SSDSE-E-2026.csv'), encoding='cp932', header=1)
ssdse_e = ssdse_e_raw.iloc[1:].copy()
ssdse_e.columns = ssdse_e_raw.iloc[0].values
ssdse_e = ssdse_e[ssdse_e['都道府県'] != '全国'].reset_index(drop=True)

# ================================================================
# ■ 変数の構築
# ================================================================
df = pd.DataFrame()
df['都道府県'] = df_b['都道府県']

# TFR（合計特殊出生率, 2022年）
df['TFR'] = pd.to_numeric(df_b['合計特殊出生率'], errors='coerce')

# 婚姻率 = 婚姻件数 / 15〜64歳人口 × 1000
df['婚姻率'] = (
    pd.to_numeric(df_b['婚姻件数'], errors='coerce') /
    pd.to_numeric(df_b['15～64歳人口'], errors='coerce') * 1000
)

# 高齢化率 = 65歳以上人口 / 総人口 × 100
df['高齢化率'] = (
    pd.to_numeric(df_b['65歳以上人口'], errors='coerce') /
    pd.to_numeric(df_b['総人口'], errors='coerce') * 100
)

# 大学進学率 = 高等学校卒業者のうち進学者数 / 高等学校卒業者数 × 100
df['大学進学率'] = (
    pd.to_numeric(df_b['高等学校卒業者のうち進学者数'], errors='coerce') /
    pd.to_numeric(df_b['高等学校卒業者数'], errors='coerce') * 100
)

# 保育所数（千対）= 保育所等数 / 総人口 × 1000
df['保育所数千対'] = (
    pd.to_numeric(df_b['保育所等数'], errors='coerce') /
    pd.to_numeric(df_b['総人口'], errors='coerce') * 1000
)

# 年平均気温（℃）
df['年平均気温'] = pd.to_numeric(df_b['年平均気温'], errors='coerce')

# 1人当たり県民所得（万円, SSDSE-E 2021年基準）
income_map = dict(zip(
    ssdse_e['都道府県'],
    pd.to_numeric(ssdse_e['1人当たり県民所得（平成27年基準）'], errors='coerce')
))
df['1人当たり県民所得'] = df['都道府県'].map(income_map)

# 人口密度（人/km²）= 総人口 / 総面積(ha÷100)
area_map = dict(zip(
    ssdse_e['都道府県'],
    pd.to_numeric(ssdse_e['総面積（北方地域及び竹島を除く）'], errors='coerce') / 100
))
df['面積_km2'] = df['都道府県'].map(area_map)
df['人口密度'] = pd.to_numeric(df_b['総人口'], errors='coerce') / df['面積_km2']
df['ln人口密度'] = np.log(df['人口密度'])

# 3区分の設定
df['区分'] = 'その他'
df.loc[df['都道府県'] == '東京都', '区分'] = '東京都'
df.loc[(df['人口密度'] < 150) & (df['都道府県'] != '東京都'), '区分'] = '過疎化進行県'

df = df.dropna().reset_index(drop=True)
N = len(df)

VAR_NAMES = ['婚姻率', '1人当たり県民所得', '高齢化率', 'ln人口密度',
             '保育所数千対', '大学進学率', '年平均気温']

print("=" * 60)
print("■ 記述統計（都道府県別, 2022年）N =", N)
print("=" * 60)
print(df[['TFR'] + VAR_NAMES].describe().round(3))

print("\n3区分の内訳:")
for g in ['東京都', '過疎化進行県', 'その他']:
    prefs_g = df[df['区分'] == g]['都道府県'].tolist()
    print(f"  {g}（{len(prefs_g)}）: {prefs_g}")

# ================================================================
# ■ Step1. 相関分析（TFRと各変数）
# ================================================================
print("\n" + "=" * 60)
print("■ Step1. 相関分析（TFR vs 各説明変数）")
print("=" * 60)
print(f"\n  {'変数':<16} {'r':>8} {'p値':>10} {'有意':>6}")
print("  " + "-" * 44)
X_raw = df[VAR_NAMES].values
tfr = df['TFR'].values
for i, name in enumerate(VAR_NAMES):
    r, p = stats.pearsonr(X_raw[:, i], tfr)
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else 'n.s.'
    print(f"  {name:<16} {r:>8.4f} {p:>10.4f} {sig:>6}")

# ================================================================
# ■ Step2. 重回帰分析
# ================================================================
print("\n" + "=" * 60)
print("■ Step2. 重回帰分析（全説明変数）")
print("=" * 60)
X_ols = sm.add_constant(X_raw)
model = sm.OLS(tfr, X_ols).fit()
print(model.summary2())

# ================================================================
# ■ Step3. 3区分別の比較
# ================================================================
print("\n" + "=" * 60)
print("■ Step3. 3区分別TFR比較")
print("=" * 60)
for g in ['東京都', 'その他', '過疎化進行県']:
    sub = df[df['区分'] == g]['TFR']
    print(f"  {g}: N={len(sub)}, 平均TFR={sub.mean():.3f}, 中央値={sub.median():.3f}")

# ================================================================
# ■ 図の生成（4枚）
# ================================================================
color_map = {'過疎化進行県': '#E65100', '東京都': '#1565C0', 'その他': '#2E7D32'}

# ──────────────────────────────────────────────────────────
# 図1: TFRの都道府県ランキング（3区分色分け）
# ──────────────────────────────────────────────────────────
fig1, ax1 = plt.subplots(figsize=(10, 10))
sorted_idx = np.argsort(tfr)
vals_s = tfr[sorted_idx]
prefs_s = df['都道府県'].values[sorted_idx]
grp_s   = df['区分'].values[sorted_idx]
bar_cols = [color_map[g] for g in grp_s]

ax1.barh(range(N), vals_s, color=bar_cols, alpha=0.8, edgecolor='white', height=0.75)
ax1.set_yticks(range(N))
ax1.set_yticklabels(prefs_s, fontsize=8)
ax1.axvline(tfr.mean(), color='gray', linestyle='--', linewidth=1.2)
ax1.axvline(2.07, color='black', linestyle=':', linewidth=1.2, alpha=0.6)
ax1.set_xlabel('合計特殊出生率（TFR）', fontsize=12)
ax1.set_title('都道府県別TFR（2022年）\n3区分（過疎化進行県・東京都・その他）', fontsize=12, fontweight='bold')
ax1.grid(axis='x', alpha=0.3)
legend_els = [Patch(color=c, alpha=0.8, label=g) for g, c in color_map.items()]
legend_els += [
    plt.Line2D([0], [0], color='gray', linestyle='--', label=f'全国平均={tfr.mean():.3f}'),
    plt.Line2D([0], [0], color='black', linestyle=':', label='人口置換水準=2.07')
]
ax1.legend(handles=legend_els, fontsize=9)
plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2025_H5_6_fig1_rank.png'), bbox_inches='tight', dpi=150)
plt.close(fig1)
print("\n図1保存: 2025_H5_6_fig1_rank.png")

# ──────────────────────────────────────────────────────────
# 図2: 主要変数との散布図（2×2）
# ──────────────────────────────────────────────────────────
fig2, axes2 = plt.subplots(2, 2, figsize=(12, 9))
fig2.suptitle('TFRと主要説明変数の散布図（2022年 SSDSE-B/E）', fontsize=13, fontweight='bold')

key_pairs = [
    (VAR_NAMES.index('婚姻率'), '婚姻率（‰）'),
    (VAR_NAMES.index('ln人口密度'), 'ln(人口密度)'),
    (VAR_NAMES.index('保育所数千対'), '保育所数（千対）'),
    (VAR_NAMES.index('大学進学率'), '大学進学率（%）'),
]
groups = df['区分'].values

for ax, (idx, xlabel) in zip(axes2.flat, key_pairs):
    xv = X_raw[:, idx]
    for g, col in color_map.items():
        mask = groups == g
        ax.scatter(xv[mask], tfr[mask], color=col, s=60, alpha=0.75,
                   edgecolors='white', linewidth=0.5, label=g, zorder=3)
    z = np.polyfit(xv, tfr, 1)
    xs = np.linspace(xv.min() - 0.1, xv.max() + 0.1, 100)
    ax.plot(xs, np.poly1d(z)(xs), 'k--', linewidth=1.5, alpha=0.6)
    r, p = stats.pearsonr(xv, tfr)
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else 'n.s.'
    ax.set_xlabel(xlabel, fontsize=10)
    ax.set_ylabel('TFR', fontsize=10)
    ax.set_title(f'r={r:.3f} {sig}', fontsize=11, fontweight='bold')
    ax.grid(True, alpha=0.3)
    if idx == 0:
        ax.legend(fontsize=7)

plt.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2025_H5_6_fig2_scatter.png'), bbox_inches='tight', dpi=150)
plt.close(fig2)
print("図2保存: 2025_H5_6_fig2_scatter.png")

# ──────────────────────────────────────────────────────────
# 図3: 重回帰の標準化係数プロット
# ──────────────────────────────────────────────────────────
fig3, ax3 = plt.subplots(figsize=(9, 5))

scaler = StandardScaler()
X_std = scaler.fit_transform(X_raw)
model_std = sm.OLS(tfr, sm.add_constant(X_std)).fit()
coefs = model_std.params[1:]
ses   = model_std.bse[1:]
pvals = model_std.pvalues[1:]

bar_cols3 = ['#C62828' if p < 0.05 else '#1565C0' if p < 0.1 else '#9E9E9E' for p in pvals]
yp = np.arange(len(VAR_NAMES))
ax3.barh(yp, coefs, color=bar_cols3, alpha=0.75, edgecolor='white', height=0.6)
ax3.errorbar(coefs, yp, xerr=1.96 * ses, fmt='none', color='#333', capsize=4, linewidth=1.2)
ax3.axvline(0, color='gray', linestyle='--', linewidth=1.0)
ax3.set_yticks(yp)
ax3.set_yticklabels(VAR_NAMES, fontsize=11)
ax3.set_xlabel('標準化回帰係数（±1.96SE）', fontsize=11)
ax3.set_title(f'TFR重回帰の標準化係数（データ：SSDSE-B/E 2022年）\nR²={model.rsquared:.3f} N={N}都道府県', fontsize=12, fontweight='bold')
ax3.invert_yaxis()
ax3.grid(axis='x', alpha=0.3)
legend_els3 = [Patch(color='#C62828', alpha=0.75, label='p<0.05'),
               Patch(color='#1565C0', alpha=0.75, label='p<0.10'),
               Patch(color='#9E9E9E', alpha=0.75, label='非有意')]
ax3.legend(handles=legend_els3, fontsize=9)
plt.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2025_H5_6_fig3_coef.png'), bbox_inches='tight', dpi=150)
plt.close(fig3)
print("図3保存: 2025_H5_6_fig3_coef.png")

# ──────────────────────────────────────────────────────────
# 図4: 3区分別TFR箱ひげ図 + 指標比較
# ──────────────────────────────────────────────────────────
fig4, axes4 = plt.subplots(1, 2, figsize=(13, 5))
fig4.suptitle('少子化の2タイプ：東京都 vs 過疎化進行県', fontsize=13, fontweight='bold')

group_order = ['東京都', 'その他', '過疎化進行県']
group_data  = [df[df['区分'] == g]['TFR'].values for g in group_order]
group_cols  = [color_map[g] for g in group_order]

ax4a = axes4[0]
bp4 = ax4a.boxplot(group_data,
                   labels=['東京都\n（高コスト）', 'その他', '過疎化進行県\n（若年層流出）'],
                   patch_artist=True, medianprops=dict(color='white', linewidth=2))
for patch, col in zip(bp4['boxes'], group_cols):
    patch.set_facecolor(col)
    patch.set_alpha(0.65)
ax4a.set_ylabel('TFR（合計特殊出生率）', fontsize=11)
ax4a.set_title('3区分別TFR分布（2022年）', fontsize=11, fontweight='bold')
ax4a.axhline(2.07, color='black', linestyle=':', linewidth=1.2, alpha=0.6, label='人口置換水準')
ax4a.axhline(tfr.mean(), color='gray', linestyle='--', linewidth=1.0, label=f'全国平均={tfr.mean():.3f}')
ax4a.legend(fontsize=8)
ax4a.grid(axis='y', alpha=0.3)

# 3区分の指標比較（全国平均比）
ax4b = axes4[1]
indicators = ['婚姻率', '保育所数千対', '大学進学率', 'TFR']
x_idx = np.arange(len(indicators))
width = 0.25
grp_colors_bar = [color_map[g] for g in group_order]
for i, (g, col) in enumerate(zip(group_order, grp_colors_bar)):
    sub = df[df['区分'] == g]
    vals_g = []
    for var in indicators:
        if var == 'TFR':
            vals_g.append(sub['TFR'].mean())
        else:
            vals_g.append(sub[var].mean() / df[var].mean())
    ax4b.bar(x_idx + i * width, vals_g, width, label=g, color=col, alpha=0.75, edgecolor='white')
ax4b.set_xticks(x_idx + width)
ax4b.set_xticklabels(['婚姻率', '保育所数', '大学進学率', 'TFR'], fontsize=10)
ax4b.set_ylabel('全国平均比（TFRは実値）', fontsize=10)
ax4b.set_title('3区分の特徴比較\n（各指標の全国平均比）', fontsize=11, fontweight='bold')
ax4b.axhline(1.0, color='gray', linestyle='--', linewidth=1.0, alpha=0.6, label='全国平均=1')
ax4b.legend(fontsize=8)
ax4b.grid(axis='y', alpha=0.3)

plt.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2025_H5_6_fig4_compare.png'), bbox_inches='tight', dpi=150)
plt.close(fig4)
print("図4保存: 2025_H5_6_fig4_compare.png")

print("\n全図の生成完了（4枚）")
print(f"\nデータ出典: SSDSE-B-2026.csv（2022年度）, SSDSE-E-2026.csv（統計センター）")