"""
教育用再現コード: 2022年 統計データ分析コンペティション 総務大臣賞 [高校生の部]
=================================================================
論文タイトル：体力と基礎学力の関係：都道府県別教育指標の相関・回帰分析
受賞：総務大臣賞（高校生の部）— 最優秀賞

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

  ※ 体力・学力の直接データは都道府県別統計センターデータには含まれないため、
    利用可能な教育関連指標（教育費・学校設備・進学率等）を代理指標として分析

  分析の流れ
  1. 相関ヒートマップ：教育関連指標間の相関構造
  2. 散布図：教育費 vs 大学進学率（2022年）
  3. OLS回帰：大学進学率の決定要因
  4. Ward法クラスタリング：教育環境による都道府県分類

【主要変数】
  大学進学率 = 高等学校卒業者のうち進学者数 / 高等学校卒業者数 × 100
  教育費（二人以上の世帯）
  学校師弟比（小）= 小学校児童数 / 小学校教員数
  学校師弟比（中）= 中学校生徒数 / 中学校教員数
  保育所密度 = 保育所等数 / 総人口 × 10000
  高齢化率 = 65歳以上人口 / 総人口 × 100

【データ出典】
  SSDSE-B-2026.csv: 社会・人口統計体系（都道府県データ）

【データサイエンス学習ポイント】
  1. 教育指標の代理変数設計（間接指標による分析アプローチ）
  2. Pearson相関係数と統計的有意性の検定
  3. 多変量OLS回帰と標準化係数の比較
  4. Ward法による教育環境グループの分類
=================================================================
"""

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


import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from sklearn.preprocessing import StandardScaler
from scipy.cluster.hierarchy import dendrogram, linkage
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_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)
df_2022 = df_b[df_b['年度'] == 2022].copy()

# 変数生成
df_2022['大学進学率'] = df_2022['高等学校卒業者のうち進学者数'] / df_2022['高等学校卒業者数'].replace(0, np.nan) * 100
df_2022['教育費_万円'] = df_2022['教育費（二人以上の世帯）'] / 10000
df_2022['師弟比_小学'] = df_2022['小学校児童数'] / df_2022['小学校教員数'].replace(0, np.nan)
df_2022['師弟比_中学'] = df_2022['中学校生徒数'] / df_2022['中学校教員数'].replace(0, np.nan)
df_2022['保育所密度'] = df_2022['保育所等数'] / df_2022['総人口'] * 10000
df_2022['高齢化率'] = df_2022['65歳以上人口'] / df_2022['総人口'] * 100
df_2022['消費支出_log'] = np.log(df_2022['消費支出（二人以上の世帯）'].clip(lower=1))

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

# Fig1: 相関ヒートマップ
analysis_vars = ['大学進学率', '教育費_万円', '師弟比_小学', '師弟比_中学', '保育所密度', '高齢化率']
df_corr = df_2022[analysis_vars].dropna()
corr = df_corr.corr()

fig, ax = plt.subplots(figsize=(8, 6))
im = ax.imshow(corr, cmap='RdBu_r', vmin=-1, vmax=1)
ax.set_xticks(range(len(analysis_vars)))
ax.set_yticks(range(len(analysis_vars)))
ax.set_xticklabels(analysis_vars, rotation=45, ha='right', fontsize=9)
ax.set_yticklabels(analysis_vars, fontsize=9)
for i in range(len(analysis_vars)):
    for j in range(len(analysis_vars)):
        val = corr.iloc[i, j]
        ax.text(j, i, f'{val:.2f}', ha='center', va='center',
                fontsize=8, color='white' if abs(val) > 0.5 else 'black')
plt.colorbar(im, ax=ax, label='Pearson相関係数')
ax.set_title('教育指標間の相関ヒートマップ（2022年, N=47）', fontsize=13, fontweight='bold')
plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2022_H1_fig1_corr.png'), dpi=150, bbox_inches='tight')
plt.close()
print("Fig1 saved")

# Fig2: 散布図 教育費 vs 大学進学率
df_sc = df_2022[['都道府県', '教育費_万円', '大学進学率', '地域']].dropna()
fig, ax = plt.subplots(figsize=(9, 6))
for reg, grp in df_sc.groupby('地域'):
    ax.scatter(grp['教育費_万円'], grp['大学進学率'],
               color=region_colors.get(reg, 'gray'), label=reg, s=70, alpha=0.8)
for _, row in df_sc.iterrows():
    ax.annotate(row['都道府県'][:2], (row['教育費_万円'], row['大学進学率']),
                fontsize=7, alpha=0.7, xytext=(3, 3), textcoords='offset points')
slope, intercept, r, p, _ = stats.linregress(df_sc['教育費_万円'], df_sc['大学進学率'])
xr = np.linspace(df_sc['教育費_万円'].min(), df_sc['教育費_万円'].max(), 100)
ax.plot(xr, slope * xr + intercept, 'k--', linewidth=1.5,
        label=f'回帰直線 (r={r:.2f}, p={p:.3f})')
ax.set_xlabel('教育費（万円/月）', fontsize=12)
ax.set_ylabel('大学進学率（%）', fontsize=12)
ax.set_title('教育費 vs 大学進学率（2022年, N=47）', fontsize=13, fontweight='bold')
ax.legend(fontsize=8, ncol=2)
ax.grid(alpha=0.3)
plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2022_H1_fig2_scatter.png'), dpi=150, bbox_inches='tight')
plt.close()
print("Fig2 saved")

# Fig3: OLS回帰 係数プロット
xvars = ['教育費_万円', '師弟比_小学', '師弟比_中学', '保育所密度', '高齢化率']
df_reg = df_2022[['大学進学率'] + xvars].dropna()
X = sm.add_constant(df_reg[xvars])
res = sm.OLS(df_reg['大学進学率'], X).fit()
print(res.summary())

# 標準化係数
from sklearn.preprocessing import StandardScaler as SS
sc = SS()
X_std = pd.DataFrame(sc.fit_transform(df_reg[xvars]), columns=xvars)
y_arr = df_reg['大学進学率'].values
y_std = pd.Series((y_arr - y_arr.mean()) / y_arr.std())
res_std = sm.OLS(y_std, sm.add_constant(X_std)).fit()

fig, ax = plt.subplots(figsize=(8, 5))
coefs = res_std.params.drop('const')
ses = res_std.bse.drop('const')
pvals = res.pvalues.drop('const')
colors_c = ['#e05c5c' if p < 0.05 else '#888888' for p in pvals]
ax.barh(range(len(coefs)), coefs, xerr=1.96 * ses, color=colors_c, alpha=0.8,
        error_kw={'elinewidth': 1.5, 'capsize': 4})
ax.set_yticks(range(len(coefs)))
ax.set_yticklabels(coefs.index, fontsize=10)
ax.axvline(0, color='black', linewidth=0.8)
ax.set_xlabel('標準化回帰係数（β）', fontsize=12)
ax.set_title(f'大学進学率の決定要因（標準化OLS, N=47）\nR²={res.rsquared:.3f}（赤=p<0.05）', fontsize=12, fontweight='bold')
ax.grid(axis='x', alpha=0.3)
plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2022_H1_fig3_ols.png'), dpi=150, bbox_inches='tight')
plt.close()
print("Fig3 saved")

# Fig4: Ward法クラスタリング
cluster_vars = ['大学進学率', '教育費_万円', '師弟比_小学', '保育所密度', '高齢化率']
df_clust = df_2022[['都道府県'] + cluster_vars].dropna().set_index('都道府県')
scaler = StandardScaler()
X_sc = scaler.fit_transform(df_clust)
Z = linkage(X_sc, method='ward')
fig, ax = plt.subplots(figsize=(12, 6))
dendrogram(Z, labels=df_clust.index.tolist(), ax=ax,
           color_threshold=Z[-3, 2], leaf_rotation=90, leaf_font_size=8)
ax.set_title('Ward法クラスタリング：教育環境による都道府県分類（2022年）', fontsize=13, fontweight='bold')
ax.set_xlabel('都道府県', fontsize=11)
ax.set_ylabel('距離（Ward法）', fontsize=11)
ax.grid(axis='y', alpha=0.3)
plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2022_H1_fig4_cluster.png'), dpi=150, bbox_inches='tight')
plt.close()
print("Fig4 saved")
print("All done!")