"""
2019_H2_yushu.py
都道府県別学力格差の要因分析：社会経済環境と教育支出の影響
優秀賞（高校生部門）2019年度（令和元年度）

代理変数について：
  SSDSE-B には全国学力・学習状況調査の得点データは収録されていない。
  本分析では「高等学校卒業者のうち進学者数 / 高等学校卒業者数」を
  大学進学率として計算し、学力・学習達成度の代理変数として使用する。
  大学進学率は家庭の社会経済的背景や教育投資と強く連動することが
  先行研究で示されており（文部科学省 教育指標の国際比較 等）、
  学力格差の地域差を反映する合理的な代理変数とみなせる。
"""

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

print("=== 列名一覧 ===")
print(df_b.columns.tolist())
print()

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

df_b['地域'] = df_b['都道府県'].map(region_map)

# ----------------------------------------------------------------
# 変数の計算（2023年断面）
# ----------------------------------------------------------------
df_2023 = df_b[df_b['年度'] == 2023].copy()

# 目的変数: 大学進学率（学力の代理変数）
df_2023['大学進学率'] = (
    pd.to_numeric(df_2023['高等学校卒業者のうち進学者数'], errors='coerce') /
    pd.to_numeric(df_2023['高等学校卒業者数'], errors='coerce') * 100
)

# 説明変数候補
df_2023['教育費'] = pd.to_numeric(df_2023['教育費（二人以上の世帯）'], errors='coerce')
df_2023['消費支出'] = pd.to_numeric(df_2023['消費支出（二人以上の世帯）'], errors='coerce')
df_2023['高齢化率'] = (
    pd.to_numeric(df_2023['65歳以上人口'], errors='coerce') /
    pd.to_numeric(df_2023['総人口'], errors='coerce') * 100
)
df_2023['少子化率'] = (
    pd.to_numeric(df_2023['15歳未満人口'], errors='coerce') /
    pd.to_numeric(df_2023['総人口'], errors='coerce') * 100
)
df_2023['log_消費支出'] = np.log(df_2023['消費支出'])
df_2023['log_教育費'] = np.log(df_2023['教育費'])

# 都道府県略称
def abbrev(name):
    return name.replace('県','').replace('都','').replace('府','').replace('道','')

df_2023['都道府県略'] = df_2023['都道府県'].apply(abbrev)

# 欠損除去
df_ana = df_2023.dropna(subset=['大学進学率','教育費','消費支出','高齢化率','少子化率']).copy()
df_ana = df_ana.reset_index(drop=True)

print("=== 2023年 分析データ ===")
print(f"観測数: {len(df_ana)} 都道府県")
print(f"大学進学率 平均: {df_ana['大学進学率'].mean():.2f}%")
print(f"大学進学率 最大: {df_ana['大学進学率'].max():.2f}% ({df_ana.loc[df_ana['大学進学率'].idxmax(),'都道府県']})")
print(f"大学進学率 最小: {df_ana['大学進学率'].min():.2f}% ({df_ana.loc[df_ana['大学進学率'].idxmin(),'都道府県']})")
print()

# ----------------------------------------------------------------
# 相関分析
# ----------------------------------------------------------------
r_edu, p_edu = stats.pearsonr(df_ana['教育費'], df_ana['大学進学率'])
r_con, p_con = stats.pearsonr(df_ana['消費支出'], df_ana['大学進学率'])
r_age, p_age = stats.pearsonr(df_ana['高齢化率'], df_ana['大学進学率'])
r_kid, p_kid = stats.pearsonr(df_ana['少子化率'], df_ana['大学進学率'])
print("=== 相関分析 ===")
print(f"教育費   vs 大学進学率: r={r_edu:.4f}, p={p_edu:.4f}")
print(f"消費支出 vs 大学進学率: r={r_con:.4f}, p={p_con:.4f}")
print(f"高齢化率 vs 大学進学率: r={r_age:.4f}, p={p_age:.4f}")
print(f"少子化率 vs 大学進学率: r={r_kid:.4f}, p={p_kid:.4f}")
print()

# ----------------------------------------------------------------
# バックワード AIC 変数選択
# ----------------------------------------------------------------
def backward_aic(y, X_df):
    cols = list(X_df.columns)
    while len(cols) > 0:
        X = sm.add_constant(X_df[cols].astype(float))
        current_aic = sm.OLS(y, X).fit().aic
        drop_col = None
        for c in cols:
            remaining = [x for x in cols if x != c]
            if len(remaining) == 0:
                break
            X_try = sm.add_constant(X_df[remaining].astype(float))
            aic_try = sm.OLS(y, X_try).fit().aic
            if aic_try < current_aic:
                current_aic = aic_try
                drop_col = c
        if drop_col is None:
            break
        cols.remove(drop_col)
        print(f"  除外: {drop_col} → AIC={current_aic:.2f}")
    return cols

y_reg = df_ana['大学進学率']
X_cands = df_ana[['教育費', '消費支出', '高齢化率', '少子化率']].copy()

print("=== バックワード AIC 変数選択 ===")
# フル AIC
X_full = sm.add_constant(X_cands.astype(float))
full_aic = sm.OLS(y_reg, X_full).fit().aic
print(f"フルモデル AIC: {full_aic:.2f}")

final_cols = backward_aic(y_reg, X_cands)
print(f"最終変数: {final_cols}")

# 最終モデル
X_final = sm.add_constant(df_ana[final_cols].astype(float))
model_final = sm.OLS(y_reg, X_final).fit()

print()
print("=== 最終モデル結果 ===")
print(model_final.summary())
print(f"R²     = {model_final.rsquared:.4f}")
print(f"Adj R² = {model_final.rsquared_adj:.4f}")
print(f"AIC    = {model_final.aic:.4f}")
print()
print("変数ごとの係数・p値:")
for var in final_cols:
    print(f"  {var}: β={model_final.params[var]:.4f}, p={model_final.pvalues[var]:.4f}")
print()

# 標準化偏回帰係数
std_betas = {}
for var in final_cols:
    beta_std = model_final.params[var] * df_ana[var].std() / df_ana['大学進学率'].std()
    std_betas[var] = beta_std
    print(f"  {var}: 標準化β = {beta_std:.4f}")
print()

# ----------------------------------------------------------------
# 図1: 都道府県別大学進学率ランキング（横棒グラフ）
# ----------------------------------------------------------------
df_sorted = df_ana.sort_values('大学進学率', ascending=True).copy()
national_avg = df_ana['大学進学率'].mean()

fig1, ax1 = plt.subplots(figsize=(10, 13))

bar_colors = [region_colors[r] for r in df_sorted['地域']]
bars = ax1.barh(df_sorted['都道府県略'], df_sorted['大学進学率'],
                color=bar_colors, alpha=0.85, height=0.75)

ax1.axvline(national_avg, color='#333333', linewidth=2, linestyle='--',
            label=f'全国平均 {national_avg:.1f}%', zorder=5)

ax1.set_xlabel('大学進学率（%）', fontsize=13)
ax1.set_title('都道府県別 大学進学率ランキング（2023年）\n（高校卒業者のうち大学等進学者の割合）',
              fontsize=14, fontweight='bold', pad=12)
ax1.set_xlim(35, 85)

# 凡例（地域色）
from matplotlib.patches import Patch
legend_elems = [Patch(facecolor=region_colors[r], label=r, alpha=0.85) for r in regions_order]
legend_elems.append(plt.Line2D([0],[0], color='#333333', linewidth=2, linestyle='--',
                                label=f'全国平均 {national_avg:.1f}%'))
ax1.legend(handles=legend_elems, loc='lower right', fontsize=9, framealpha=0.9)

ax1.grid(True, axis='x', alpha=0.3)
ax1.tick_params(axis='y', labelsize=9)
plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2019_H2_fig1.png'), dpi=150, bbox_inches='tight')
plt.close(fig1)
print("図1 保存完了")

# ----------------------------------------------------------------
# 図2: 教育費 vs 大学進学率 散布図
# ----------------------------------------------------------------
fig2, ax2 = plt.subplots(figsize=(11, 8))

for region in regions_order:
    d = df_ana[df_ana['地域'] == region]
    ax2.scatter(d['教育費'], d['大学進学率'], color=region_colors[region],
                s=65, alpha=0.85, label=region, zorder=3)

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

# 回帰直線
x_vals = df_ana['教育費'].values
y_vals = df_ana['大学進学率'].values
slope, intercept, r_val, p_val, se = stats.linregress(x_vals, y_vals)
x_line = np.linspace(x_vals.min(), x_vals.max(), 200)
ax2.plot(x_line, slope * x_line + intercept, 'k--', linewidth=2, zorder=4,
         label=f'回帰直線 (r={r_val:.3f}, p={p_val:.4f})')

pval_text = f'p = {p_val:.4f}' if p_val >= 0.0001 else 'p < 0.0001'
ax2.text(0.05, 0.95,
         f'Pearson r = {r_val:.3f}\n{pval_text}',
         transform=ax2.transAxes, fontsize=11, va='top',
         bbox=dict(boxstyle='round,pad=0.4', facecolor='white', alpha=0.85))

ax2.set_xlabel('教育費（円/月、二人以上世帯）', fontsize=12)
ax2.set_ylabel('大学進学率（%）', fontsize=12)
ax2.set_title('教育費と大学進学率の関係（2023年、47都道府県）', fontsize=14, fontweight='bold')
ax2.legend(loc='upper right', fontsize=9, framealpha=0.9)
ax2.grid(True, alpha=0.3)
plt.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2019_H2_fig2.png'), dpi=150, bbox_inches='tight')
plt.close(fig2)
print("図2 保存完了")

# ----------------------------------------------------------------
# 図3: 消費支出 vs 大学進学率 散布図（外れ値注釈付き）
# ----------------------------------------------------------------
fig3, ax3 = plt.subplots(figsize=(11, 8))

for region in regions_order:
    d = df_ana[df_ana['地域'] == region]
    ax3.scatter(d['消費支出'] / 1000, d['大学進学率'], color=region_colors[region],
                s=65, alpha=0.85, label=region, zorder=3)

# 都道府県ラベル
for _, row in df_ana.iterrows():
    ax3.annotate(row['都道府県略'],
                 (row['消費支出'] / 1000, row['大学進学率']),
                 fontsize=7, ha='center', va='bottom',
                 xytext=(0, 4), textcoords='offset points', color='#333333')

# 回帰直線
x_con = df_ana['消費支出'].values / 1000
y_con = df_ana['大学進学率'].values
slope_c, intercept_c, r_c, p_c, se_c = stats.linregress(x_con, y_con)
x_line_c = np.linspace(x_con.min(), x_con.max(), 200)
ax3.plot(x_line_c, slope_c * x_line_c + intercept_c, 'k--', linewidth=2, zorder=4,
         label=f'回帰直線 (r={r_c:.3f}, p={p_c:.4f})')

pval_text_c = f'p = {p_c:.4f}' if p_c >= 0.0001 else 'p < 0.0001'
ax3.text(0.05, 0.95,
         f'Pearson r = {r_c:.3f}\n{pval_text_c}',
         transform=ax3.transAxes, fontsize=11, va='top',
         bbox=dict(boxstyle='round,pad=0.4', facecolor='white', alpha=0.85))

# 東京・神奈川の注釈
outlier_prefs = ['東京', '神奈川', '沖縄']
for _, row in df_ana[df_ana['都道府県略'].isin(outlier_prefs)].iterrows():
    ax3.annotate(f'← {row["都道府県略"]}',
                 (row['消費支出'] / 1000, row['大学進学率']),
                 fontsize=9, color='#C62828', fontweight='bold',
                 xytext=(8, 0), textcoords='offset points')

ax3.set_xlabel('消費支出（千円/月、二人以上世帯）', fontsize=12)
ax3.set_ylabel('大学進学率（%）', fontsize=12)
ax3.set_title('消費支出（所得水準代理）と大学進学率の関係（2023年、47都道府県）',
              fontsize=14, fontweight='bold')
ax3.legend(loc='upper left', fontsize=9, framealpha=0.9)
ax3.grid(True, alpha=0.3)
plt.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2019_H2_fig3.png'), dpi=150, bbox_inches='tight')
plt.close(fig3)
print("図3 保存完了")

# ----------------------------------------------------------------
# 図4: バックワード選択後の標準化偏回帰係数（横棒グラフ）
# ----------------------------------------------------------------
# 標準化偏回帰係数 + 95%CI
conf_int = model_final.conf_int()  # columns: [0, 1]

# std βの信頼区間（deltaメソッドの簡易近似：CI/sd_x * sd_y）
std_beta_vals = []
std_beta_lo   = []
std_beta_hi   = []
var_labels     = []

label_map = {
    '教育費': '教育費（月額）',
    '消費支出': '消費支出（所得代理）',
    '高齢化率': '高齢化率',
    '少子化率': '少子化率（15歳未満割合）'
}

for var in final_cols:
    sd_x = df_ana[var].std()
    sd_y = df_ana['大学進学率'].std()
    beta_s = model_final.params[var] * sd_x / sd_y
    lo_s   = conf_int.loc[var, 0] * sd_x / sd_y
    hi_s   = conf_int.loc[var, 1] * sd_x / sd_y
    std_beta_vals.append(beta_s)
    std_beta_lo.append(lo_s)
    std_beta_hi.append(hi_s)
    var_labels.append(label_map.get(var, var))

# 有意性による色分け
bar_colors4 = []
for var in final_cols:
    p = model_final.pvalues[var]
    b = model_final.params[var]
    if p < 0.05 and b > 0:
        bar_colors4.append('#e05c5c')  # 正有意 → 赤
    elif p < 0.05 and b < 0:
        bar_colors4.append('#4e9af1')  # 負有意 → 青
    else:
        bar_colors4.append('#aaaaaa')  # 非有意 → グレー

xerr_lo4 = [abs(v - l) for v, l in zip(std_beta_vals, std_beta_lo)]
xerr_hi4 = [abs(h - v) for v, h in zip(std_beta_vals, std_beta_hi)]

fig4, ax4 = plt.subplots(figsize=(9, max(4, len(final_cols) * 1.4)))

bars4 = ax4.barh(var_labels, std_beta_vals, color=bar_colors4, alpha=0.85, height=0.55)
ax4.errorbar(std_beta_vals, var_labels,
             xerr=[xerr_lo4, xerr_hi4],
             fmt='none', color='#333333', capsize=5, linewidth=1.5, zorder=5)

ax4.axvline(0, color='black', linewidth=1.0)

# 凡例
from matplotlib.patches import Patch
leg_elem4 = [
    Patch(facecolor='#e05c5c', alpha=0.85, label='正方向・有意（p<0.05）'),
    Patch(facecolor='#4e9af1', alpha=0.85, label='負方向・有意（p<0.05）'),
    Patch(facecolor='#aaaaaa', alpha=0.85, label='非有意（p≥0.05）'),
]
ax4.legend(handles=leg_elem4, loc='lower right', fontsize=9, framealpha=0.9)

ax4.set_xlabel('標準化偏回帰係数（β）', fontsize=12)
ax4.set_title(f'バックワード AIC 選択後の標準化偏回帰係数と95%CI\n（目的変数：大学進学率、最終 Adj R²={model_final.rsquared_adj:.3f}）',
              fontsize=13, fontweight='bold')
ax4.grid(True, axis='x', alpha=0.3)
plt.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2019_H2_fig4.png'), dpi=150, bbox_inches='tight')
plt.close(fig4)
print("図4 保存完了")

# ----------------------------------------------------------------
# サマリー出力
# ----------------------------------------------------------------
print()
print("=== 最終モデル サマリー（HTML用）===")
print(f"分析年度: 2023年 （47都道府県断面）")
print(f"目的変数: 大学進学率（%）")
print(f"最終説明変数: {final_cols}")
print(f"R²     = {model_final.rsquared:.4f}")
print(f"Adj R² = {model_final.rsquared_adj:.4f}")
print(f"AIC    = {model_final.aic:.2f}")
print(f"全国平均大学進学率: {df_ana['大学進学率'].mean():.2f}%")
print(f"最大（{df_ana.loc[df_ana['大学進学率'].idxmax(),'都道府県']}）: {df_ana['大学進学率'].max():.2f}%")
print(f"最小（{df_ana.loc[df_ana['大学進学率'].idxmin(),'都道府県']}）: {df_ana['大学進学率'].min():.2f}%")
print()

print("DONE: 2019_H2_yushu")
