"""
2021_H1_daijin.py
都道府県別教育環境と学力格差の要因分析
総務大臣賞（高校生の部）2021年度
"""

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

# 2022年断面
df = df_b[df_b['年度'] == 2022].copy().reset_index(drop=True)
assert len(df) == 47, f"都道府県数が47ではありません: {len(df)}"

# ─── 変数構築 ─────────────────────────────────────────────────
# 大学進学率（%）: 高校卒業者のうち大学・短大進学者 / 高校卒業者数
df['大学進学率'] = df['高等学校卒業者のうち進学者数'] / df['高等学校卒業者数'] * 100

# 教育費（円/月, 二人以上世帯）
df['教育費'] = df['教育費（二人以上の世帯）'].astype(float)

# 人口1万人あたり学校数（幼小中高合計）
df['学校数_万人あたり'] = (
    (df['幼稚園数'] + df['小学校数'] + df['中学校数'] + df['高等学校数'])
    / df['総人口'] * 10000
)

# 消費支出（円/月）: 経済力の代理
df['消費支出'] = df['消費支出（二人以上の世帯）'].astype(float)

# 合計特殊出生率
df['出生率'] = df['合計特殊出生率'].astype(float)

# 15歳未満人口比率（%）
df['子供比率'] = df['15歳未満人口'] / df['総人口'] * 100

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

# 分析用データフレーム（欠損除去）
vars_use = ['大学進学率', '教育費', '学校数_万人あたり', '消費支出', '出生率', '子供比率']
df_ana = df[['都道府県', '地域'] + vars_use].dropna().copy()
print(f"分析対象都道府県数: {len(df_ana)}")

# ─── 重回帰分析（標準化OLS）────────────────────────────────────
Y_col  = '大学進学率'
X_cols = ['教育費', '学校数_万人あたり', '消費支出', '出生率', '子供比率']

Y  = df_ana[Y_col]
X  = df_ana[X_cols]

# 標準化
Y_z = (Y - Y.mean()) / Y.std()
X_z = (X - X.mean()) / X.std()

X_z_sm = sm.add_constant(X_z)
model  = sm.OLS(Y_z, X_z_sm).fit()
print(model.summary())

# 標準化偏回帰係数（切片除く）
coef_df = pd.DataFrame({
    '変数':  X_cols,
    '係数':  model.params[X_cols].values,
    'CI_lo': model.conf_int().loc[X_cols, 0].values,
    'CI_hi': model.conf_int().loc[X_cols, 1].values,
    'p値':   model.pvalues[X_cols].values,
})
coef_df = coef_df.sort_values('係数')
print("\n標準化偏回帰係数:")
print(coef_df)

# 相関行列（ピアソン）
corr_mat = df_ana[vars_use].corr()
print("\n相関行列:")
print(corr_mat.round(3))

# R²・調整済みR²
r2     = model.rsquared
r2_adj = model.rsquared_adj
print(f"\nR² = {r2:.3f}, 調整済みR² = {r2_adj:.3f}")

# 都道府県別進学率ランキング
df_rank = df_ana.sort_values('大学進学率', ascending=False).reset_index(drop=True)
top10    = df_rank.head(10)
bot10    = df_rank.tail(10)
print("\n上位10都道府県:")
print(top10[['都道府県', '大学進学率', '地域']].to_string(index=False))
print("\n下位10都道府県:")
print(bot10[['都道府県', '大学進学率', '地域']].to_string(index=False))

# 教育費と大学進学率の相関係数
r_edu, p_edu = stats.pearsonr(df_ana['教育費'], df_ana['大学進学率'])
r_sho, p_sho = stats.pearsonr(df_ana['消費支出'], df_ana['大学進学率'])
print(f"\n教育費×大学進学率: r={r_edu:.3f}, p={p_edu:.4f}")
print(f"消費支出×大学進学率: r={r_sho:.3f}, p={p_sho:.4f}")

# ─── 図1: 散布図（教育費 vs 大学進学率）──────────────────────────
fig1, ax1 = plt.subplots(figsize=(10, 7))

for region, grp in df_ana.groupby('地域'):
    ax1.scatter(
        grp['教育費'], grp['大学進学率'],
        c=region_colors[region], label=region, s=60, zorder=3, alpha=0.85
    )

# 都道府県ラベル
for _, row in df_ana.iterrows():
    ax1.annotate(
        row['都道府県'],
        (row['教育費'], row['大学進学率']),
        fontsize=6.5, ha='left', va='bottom',
        xytext=(2, 2), textcoords='offset points', color='#333333'
    )

# 回帰線
x_reg = np.linspace(df_ana['教育費'].min(), df_ana['教育費'].max(), 200)
slope, intercept, rval, pval, _ = stats.linregress(df_ana['教育費'], df_ana['大学進学率'])
ax1.plot(x_reg, slope * x_reg + intercept, color='#333', linewidth=1.5,
         linestyle='--', zorder=2, label=f'回帰線 r={rval:.2f}')

ax1.set_xlabel('教育費（円/月, 二人以上世帯）', fontsize=12)
ax1.set_ylabel('大学進学率（%）', fontsize=12)
ax1.set_title('図1: 教育費と大学進学率の関係（2022年・47都道府県）', fontsize=13, pad=12)
ax1.legend(fontsize=9, loc='upper left', framealpha=0.9)
ax1.grid(True, alpha=0.3)
fig1.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2021_H1_fig1.png'), bbox_inches='tight')
plt.close(fig1)
print("fig1 saved")

# ─── 図2: 相関行列ヒートマップ ───────────────────────────────────
var_labels = {
    '大学進学率': '大学\n進学率',
    '教育費': '教育費',
    '学校数_万人あたり': '学校数\n/万人',
    '消費支出': '消費\n支出',
    '出生率': '出生率',
    '子供比率': '子供\n比率'
}
corr_disp = corr_mat.copy()
corr_disp.index   = [var_labels[v] for v in corr_disp.index]
corr_disp.columns = [var_labels[v] for v in corr_disp.columns]

fig2, ax2 = plt.subplots(figsize=(7, 6))
n = len(corr_disp)
im = ax2.imshow(corr_disp.values, cmap='RdBu_r', vmin=-1, vmax=1, aspect='auto')
plt.colorbar(im, ax=ax2, shrink=0.8)

for i in range(n):
    for j in range(n):
        val = corr_disp.values[i, j]
        color = 'white' if abs(val) > 0.6 else 'black'
        ax2.text(j, i, f'{val:.2f}', ha='center', va='center',
                 fontsize=10, color=color, fontweight='bold')

ax2.set_xticks(range(n))
ax2.set_yticks(range(n))
ax2.set_xticklabels(corr_disp.columns, fontsize=10)
ax2.set_yticklabels(corr_disp.index, fontsize=10)
ax2.set_title('図2: 教育・経済関連変数の相関行列（2022年）', fontsize=12, pad=12)
fig2.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2021_H1_fig2.png'), bbox_inches='tight')
plt.close(fig2)
print("fig2 saved")

# ─── 図3: 標準化偏回帰係数（横棒グラフ）────────────────────────────
var_name_map = {
    '教育費': '教育費（円/月）',
    '学校数_万人あたり': '学校数/万人',
    '消費支出': '消費支出（円/月）',
    '出生率': '合計特殊出生率',
    '子供比率': '15歳未満人口比率（%）'
}
coef_df['変数名'] = coef_df['変数'].map(var_name_map)
colors_bar = ['#e05c5c' if c > 0 else '#4e9af1' for c in coef_df['係数']]
sig_mark   = ['*' if p < 0.05 else '' for p in coef_df['p値']]

fig3, ax3 = plt.subplots(figsize=(8, 5))
bars = ax3.barh(
    coef_df['変数名'], coef_df['係数'],
    color=colors_bar, alpha=0.8, edgecolor='white', linewidth=0.5
)
# 95%CI
for i, (_, row) in enumerate(coef_df.iterrows()):
    ax3.plot(
        [row['CI_lo'], row['CI_hi']], [i, i],
        color='#333', linewidth=1.5, zorder=5
    )
    ax3.plot([row['CI_lo']], [i], '|', color='#333', markersize=7)
    ax3.plot([row['CI_hi']], [i], '|', color='#333', markersize=7)

# 係数値とp値マーク
for i, (_, row) in enumerate(coef_df.iterrows()):
    offset = 0.005
    sign = 1 if row['係数'] >= 0 else -1
    ax3.text(
        row['CI_hi'] + offset if row['係数'] >= 0 else row['CI_lo'] - offset,
        i, f'{row["係数"]:.3f}{sig_mark[i]}',
        va='center', ha='left' if row['係数'] >= 0 else 'right',
        fontsize=9, color='#333'
    )

ax3.axvline(0, color='#333', linewidth=0.8, linestyle='-')
ax3.set_xlabel('標準化偏回帰係数（β）', fontsize=11)
ax3.set_title('図3: 標準化偏回帰係数と95%信頼区間\n（目的変数：大学進学率）', fontsize=12, pad=10)
ax3.set_xlim(coef_df['CI_lo'].min() - 0.15, coef_df['CI_hi'].max() + 0.2)
note = f'R²={r2:.3f}  調整済みR²={r2_adj:.3f}   *p<0.05'
ax3.text(0.98, 0.02, note, transform=ax3.transAxes, fontsize=9,
         ha='right', va='bottom', color='#555')
ax3.grid(axis='x', alpha=0.3)
fig3.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2021_H1_fig3.png'), bbox_inches='tight')
plt.close(fig3)
print("fig3 saved")

# ─── 図4: 都道府県別大学進学率ランキング（上位・下位10県）──────────
top10_disp = df_rank.head(10).copy()
bot10_disp = df_rank.tail(10).copy()
rank_df = pd.concat([top10_disp, bot10_disp]).reset_index(drop=True)
rank_df['色'] = rank_df['地域'].map(region_colors)

fig4, ax4 = plt.subplots(figsize=(10, 6))
x_pos = list(range(10)) + [v + 11.5 for v in range(10)]
colors_rank = [region_colors[r] for r in rank_df['地域']]
bars4 = ax4.bar(x_pos, rank_df['大学進学率'], color=colors_rank, alpha=0.85, edgecolor='white')

# 値ラベル
for bar, val in zip(bars4, rank_df['大学進学率']):
    ax4.text(bar.get_x() + bar.get_width() / 2, bar.get_height() + 0.3,
             f'{val:.1f}%', ha='center', va='bottom', fontsize=7.5, color='#333')

# x軸ラベル
ax4.set_xticks(x_pos)
ax4.set_xticklabels(rank_df['都道府県'], rotation=45, ha='right', fontsize=9)

# 区切り線
ax4.axvline(10.25, color='#aaa', linewidth=1.5, linestyle='--')
ax4.text(4.5, ax4.get_ylim()[1] * 0.97, '上位10都道府県', ha='center',
         fontsize=10, color='#c0392b', fontweight='bold')
ax4.text(16, ax4.get_ylim()[1] * 0.97, '下位10都道府県', ha='center',
         fontsize=10, color='#2980b9', fontweight='bold')

# 凡例
from matplotlib.patches import Patch
legend_elements = [Patch(facecolor=c, label=r) for r, c in region_colors.items()]
ax4.legend(handles=legend_elements, fontsize=8, loc='upper right',
           framealpha=0.9, ncol=2)

ax4.set_ylabel('大学進学率（%）', fontsize=11)
ax4.set_title('図4: 都道府県別大学進学率ランキング（2022年）', fontsize=13, pad=10)
ax4.set_ylim(0, rank_df['大学進学率'].max() * 1.12)
ax4.grid(axis='y', alpha=0.3)
fig4.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2021_H1_fig4.png'), bbox_inches='tight')
plt.close(fig4)
print("fig4 saved")

print("\n=== 統計サマリー ===")
print(f"大学進学率 全国平均: {df_ana['大学進学率'].mean():.1f}%")
print(f"大学進学率 最高: {df_ana.loc[df_ana['大学進学率'].idxmax(), '都道府県']} "
      f"{df_ana['大学進学率'].max():.1f}%")
print(f"大学進学率 最低: {df_ana.loc[df_ana['大学進学率'].idxmin(), '都道府県']} "
      f"{df_ana['大学進学率'].min():.1f}%")
print(f"教育費 全国平均: {df_ana['教育費'].mean():.0f}円/月")
print(f"出生率 全国平均: {df_ana['出生率'].mean():.2f}")

# 地域別平均
print("\n地域別大学進学率平均:")
print(df_ana.groupby('地域')['大学進学率'].mean().sort_values(ascending=False).round(1))