"""
教育用再現コード: 2022年 統計データ分析コンペティション 統計活用奨励賞 [高校生の部]
=================================================================
論文タイトル：都市部と地方の教育格差：進学率・教育費の地域差分析
受賞：統計活用奨励賞（高校生の部）

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

  分析の流れ
  1. 時系列：大学進学率の地域別推移（都市部 vs 地方）
  2. 散布図：消費支出（所得代理）vs 大学進学率
  3. OLS回帰：大学進学率の決定要因
  4. 都道府県別教育格差の可視化（ランキング + 地域別格差）

【目的変数】
  大学進学率 = 高等学校卒業者のうち進学者数 / 高等学校卒業者数 × 100

【説明変数】
  教育費（二人以上世帯）(万円)
  高齢化率 = 65歳以上人口 / 総人口 × 100
  大学数（都道府県内の大学数）
  消費支出_log = log(消費支出)
  師弟比_中学 = 中学校生徒数 / 中学校教員数

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

【データサイエンス学習ポイント】
  1. 教育格差の測定方法（変動係数・基尼係数による定量化）
  2. 都市・地方分類の方法（人口密度・大学数による分類）
  3. OLS回帰の解釈（偏回帰係数と因果推論の限界）
  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_H4_katsuyo.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_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_b = df_b.sort_values(['都道府県', '年度']).reset_index(drop=True)

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

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

# Fig1: 大学進学率の時系列（地域別平均）
fig, ax = plt.subplots(figsize=(10, 5))
yearly = df_b.groupby(['年度', '地域'])['大学進学率'].mean().reset_index()
for reg, grp in yearly.groupby('地域'):
    ax.plot(grp['年度'], grp['大学進学率'], marker='o', markersize=4,
            label=reg, color=region_colors.get(reg, 'gray'))
ax.set_xlabel('年度', fontsize=12)
ax.set_ylabel('大学進学率（%）', fontsize=12)
ax.set_title('地域別 大学進学率の推移（2012〜2023年）', fontsize=14, fontweight='bold')
ax.legend(fontsize=9)
ax.grid(alpha=0.3)
ax.axhline(df_b[df_b['年度'] == 2022]['大学進学率'].mean(), color='red',
           linestyle=':', alpha=0.5, linewidth=1)
plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2022_H4_fig1_ts.png'), dpi=150, bbox_inches='tight')
plt.close()
print("Fig1 saved")

# Fig2: 散布図 消費支出 vs 大学進学率（2022年）
df_2022 = df_b[df_b['年度'] == 2022].copy()
df_sc = df_2022[['都道府県', '消費支出_log', '大学進学率', '地域']].dropna()
fig, ax = plt.subplots(figsize=(9, 6))
for reg, grp in df_sc.groupby('地域'):
    ax.scatter(grp['消費支出_log'], grp['大学進学率'],
               color=region_colors.get(reg, 'gray'), label=reg, s=65, alpha=0.8)
for _, row in df_sc.iterrows():
    ax.annotate(row['都道府県'][:2], (row['消費支出_log'], row['大学進学率']),
                fontsize=7, alpha=0.7, xytext=(3, 3), textcoords='offset points')
slope, intercept, r, p, _ = stats.linregress(df_sc['消費支出_log'], df_sc['大学進学率'])
xr = np.linspace(df_sc['消費支出_log'].min(), df_sc['消費支出_log'].max(), 100)
ax.plot(xr, slope * xr + intercept, 'k--', linewidth=1.5,
        label=f'回帰直線 (r={r:.2f}, p={p:.3f})')
ax.set_xlabel('消費支出（log）', 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_H4_fig2_scatter.png'), dpi=150, bbox_inches='tight')
plt.close()
print("Fig2 saved")

# Fig3: OLS回帰 係数プロット
xvars = ['教育費_万円', '高齢化率', '大学数', '消費支出_log', '師弟比_中学']
df_reg = df_2022[['大学進学率'] + xvars].dropna()
X = sm.add_constant(df_reg[xvars])
res = sm.OLS(df_reg['大学進学率'], X).fit()
print(res.summary())
coefs = res.params.drop('const')
ses = res.bse.drop('const')
pvals = res.pvalues.drop('const')
fig, ax = plt.subplots(figsize=(8, 5))
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('OLS回帰係数', 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_H4_fig3_ols.png'), dpi=150, bbox_inches='tight')
plt.close()
print("Fig3 saved")

# Fig4: 地域別格差の可視化（変動係数 CV）
cv_by_region = df_2022.groupby('地域')['大学進学率'].agg(['mean', 'std']).dropna()
cv_by_region['CV'] = cv_by_region['std'] / cv_by_region['mean'] * 100
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
# 地域別平均
colors_list = [region_colors.get(r, 'gray') for r in cv_by_region.index]
ax1.barh(cv_by_region.index, cv_by_region['mean'], color=colors_list, alpha=0.8)
ax1.set_xlabel('平均大学進学率（%）', fontsize=11)
ax1.set_title('地域別 平均大学進学率（2022年）', fontsize=12, fontweight='bold')
ax1.grid(axis='x', alpha=0.3)
ax1.axvline(df_2022['大学進学率'].mean(), color='red', linestyle='--', linewidth=1.5,
            label=f'全国平均: {df_2022["大学進学率"].mean():.1f}%')
ax1.legend(fontsize=9)
# 変動係数
ax2.barh(cv_by_region.index, cv_by_region['CV'], color=colors_list, alpha=0.8)
ax2.set_xlabel('変動係数（CV, %）', fontsize=11)
ax2.set_title('地域内格差（変動係数 CV）\n（大きいほど地域内格差が大きい）', fontsize=12, fontweight='bold')
ax2.grid(axis='x', alpha=0.3)
plt.suptitle('地域別 大学進学率の格差分析（2022年）', fontsize=13, fontweight='bold', y=1.02)
plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2022_H4_fig4_gap.png'), dpi=150, bbox_inches='tight')
plt.close()
print("Fig4 saved")
print("All done!")