"""
教育用再現コード: 2024年 統計データ分析コンペティション 統計数理賞（高校生）
=================================================================
論文タイトル：子供の体力・運動能力
著者：大河内花音（愛知県立一宮高等学校）

【分析概要】
  データ：SSDSE-B-2026（政府統計の総合窓口）2022年度

  Step1. 相関分析（全47都道府県）
         - 教育指標（教員児童比・教員生徒比）と高校進学率の相関
         - 多重共線性の確認（小学教員児童比 vs 中学教員生徒比, 男子率 vs 女子率）
  Step2. 多重共線性の偏相関係数による解決・VIF確認
  Step3. 重回帰分析（4グループ）
         - 小学校水準×男子多都道府県（小5男子相当）
         - 小学校水準×女子多都道府県（小5女子相当）
         - 中学校水準×男子多都道府県（中2男子相当）
         - 中学校水準×女子多都道府県（中2女子相当）
  Step4. 都道府県別 教育関与度スコアランキング分析

【変数設計】（SSDSE-B-2026より）
  目的変数：高校進学率 = 高等学校卒業者のうち進学者数 / 高等学校卒業者数 × 100
  説明変数：
    小学教員児童比  = 小学校教員数 / 小学校児童数 × 100  （小5学習環境の代理）
    中学教員生徒比  = 中学校教員数 / 中学校生徒数 × 100  （中2学習環境の代理）
    男子率          = 15歳未満人口（男）/ 総人口 × 100
    女子率          = 15歳未満人口（女）/ 総人口 × 100   ← 男子率と高相関（多重共線性）
    高齢化率        = 65歳以上人口 / 総人口 × 100
    消費支出        = 消費支出（二人以上の世帯）/ 1000 千円
    気温            = 年平均気温 (℃)
    降水日数        = 降水日数（年間）(日)

【データサイエンス学習ポイント】
  1. 偏相関係数による多重共線性の解決（小学/中学教員比 r=0.922、男子/女子率 r=0.997）
  2. 4サブグループ（小5男・小5女・中2男・中2女）の比較分析
  3. VIFと偏相関の違いと使い分け
  4. 都道府県別 教育関与度スコアによる順位分析

【データ】SSDSE-B-2026（実公的データ）—合成データは一切使用しない
=================================================================
"""

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


import pandas as pd
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import statsmodels.api as sm
from statsmodels.stats.outliers_influence import variance_inflation_factor
from scipy.stats import pearsonr
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_PATH = 'data/raw/SSDSE-B-2026.csv'
os.makedirs(FIG_DIR, exist_ok=True)

# ================================================================
# ■ Step 0. データ読み込み（SSDSE-B-2026, 2022年度, 47都道府県）
# ================================================================

print("=" * 65)
print("■ データ読み込み（SSDSE-B-2026 / 2022年度 / 47都道府県）")
print("=" * 65)

df_b = pd.read_csv(DATA_PATH, encoding='cp932', header=1)
df_b = df_b[df_b['地域コード'].str.match(r'^R\d{5}$', na=False)]
df = df_b[df_b['年度'] == 2022].copy().reset_index(drop=True)

assert len(df) == 47, f"都道府県数が{len(df)}です（47が必要）"
print(f"読み込み完了: {len(df)} 都道府県, {len(df.columns)} 変数")

# ================================================================
# ■ Step 1. 変数の構築（派生変数）
# ================================================================

print("\n" + "=" * 65)
print("■ 変数の構築（派生変数）")
print("=" * 65)

# 目的変数
df['高校進学率'] = (
    df['高等学校卒業者のうち進学者数'] / df['高等学校卒業者数'] * 100
)

# 説明変数（小・中学校の学習環境プロキシ）
df['小学教員児童比'] = df['小学校教員数'] / df['小学校児童数'] * 100   # 大きいほど少人数
df['中学教員生徒比'] = df['中学校教員数'] / df['中学校生徒数'] * 100   # 大きいほど少人数
df['中学進学率']     = (
    df['中学校卒業者のうち進学者数'] / df['中学校卒業者数'] * 100
)

# 人口構成
df['男子率']   = df['15歳未満人口（男）'] / df['総人口'] * 100
df['女子率']   = df['15歳未満人口（女）'] / df['総人口'] * 100
df['高齢化率'] = df['65歳以上人口'] / df['総人口'] * 100

# 経済・気候
df['消費支出'] = df['消費支出（二人以上の世帯）'] / 1000   # 千円単位
df['気温']     = df['年平均気温']
df['降水日数'] = df['降水日数（年間）']

# 全説明変数（相関分析用: 多重共線性ペアを含む）
EXPLAIN_FULL = [
    '小学教員児童比', '中学教員生徒比', '男子率', '女子率',
    '高齢化率', '消費支出', '気温', '降水日数'
]

# 回帰分析用（多重共線性ペアから片方を除く: 女子率を除く）
EXPLAIN_REG = [
    '小学教員児童比', '中学教員生徒比', '男子率',
    '高齢化率', '消費支出', '気温', '降水日数'
]

TARGET = '高校進学率'

# 変数ラベル（表示用）
LABEL_MAP = {
    '小学教員児童比': '小学教員/児童比(%)',
    '中学教員生徒比': '中学教員/生徒比(%)',
    '男子率':         '15歳未満男子率(%)',
    '女子率':         '15歳未満女子率(%)',
    '高齢化率':       '高齢化率(%)',
    '消費支出':       '消費支出(千円)',
    '気温':           '年平均気温(℃)',
    '降水日数':       '降水日数(日)',
}

print("目的変数 高校進学率:")
print(f"  平均 {df[TARGET].mean():.2f}%, 標準偏差 {df[TARGET].std():.2f}%")
print(f"  最小 {df[TARGET].min():.2f}%（{df.loc[df[TARGET].idxmin(), '都道府県']}）")
print(f"  最大 {df[TARGET].max():.2f}%（{df.loc[df[TARGET].idxmax(), '都道府県']}）")

# ================================================================
# ■ Step 2. 相関分析と多重共線性の確認
# ================================================================

print("\n" + "=" * 65)
print("■ 相関分析（高校進学率との相関）")
print("=" * 65)

corr_results = {}
for v in EXPLAIN_FULL:
    r, p = pearsonr(df[v], df[TARGET])
    corr_results[v] = (r, p)
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else ''
    print(f"  {LABEL_MAP[v]:<22} r={r:+.3f}  p={p:.4f} {sig}")

# 多重共線性ペアの確認
print("\n【多重共線性の確認】")
r_pair1, p_pair1 = pearsonr(df['小学教員児童比'], df['中学教員生徒比'])
r_pair2, p_pair2 = pearsonr(df['男子率'], df['女子率'])
print(f"  小学教員児童比 × 中学教員生徒比: r={r_pair1:.3f} (p={p_pair1:.4f}) ← 高相関")
print(f"  男子率 × 女子率:                 r={r_pair2:.3f} (p={p_pair2:.4f}) ← 高相関")
print("  → 偏相関係数による解決と VIF 確認が必要")

# 偏相関（制御変数あり）
def partial_corr(x, y, z):
    """z を統制した x-y の偏相関係数を返す"""
    if z.ndim == 1:
        z = z.reshape(-1, 1)
    res_x = sm.OLS(x, sm.add_constant(z)).fit().resid
    res_y = sm.OLS(y, sm.add_constant(z)).fit().resid
    r, p = pearsonr(res_x, res_y)
    return r, p

z_ctrl = df[['高齢化率', '消費支出']].values
r_p1, p_p1 = partial_corr(
    df['小学教員児童比'].values,
    df['中学教員生徒比'].values,
    z_ctrl
)
r_p2, p_p2 = partial_corr(
    df['男子率'].values,
    df['女子率'].values,
    z_ctrl
)
print(f"\n  偏相関（高齢化率・消費支出を統制）:")
print(f"    小学教員児童比 × 中学教員生徒比: 偏相関 r={r_p1:.3f} (p={p_p1:.4f})")
print(f"    男子率 × 女子率:                 偏相関 r={r_p2:.3f} (p={p_p2:.4f})")

# ================================================================
# ■ Step 3. VIF 確認（多重共線性ペア除外後）
# ================================================================

print("\n" + "=" * 65)
print("■ VIF 確認（回帰分析用変数セット）")
print("=" * 65)

# 標準化
df_z = df.copy()
for v in EXPLAIN_FULL + [TARGET]:
    mu, sd = df[v].mean(), df[v].std()
    df_z[v + '_z'] = (df[v] - mu) / sd

Z_FULL = [v + '_z' for v in EXPLAIN_FULL]
Z_REG  = [v + '_z' for v in EXPLAIN_REG]

# VIF: フルセット（多重共線性あり）
X_full = df_z[Z_FULL].values
vif_full = {}
print("VIF（フルセット — 多重共線性ペアを含む）:")
for i, v in enumerate(EXPLAIN_FULL):
    vf = variance_inflation_factor(X_full, i)
    vif_full[v] = vf
    flag = ' ★高VIF' if vf > 5 else ''
    print(f"  {LABEL_MAP[v]:<22} VIF={vf:.2f}{flag}")

# VIF: 回帰分析用（女子率除外）
X_reg = df_z[Z_REG].values
vif_reg = {}
print("\nVIF（回帰分析用 — 女子率除外後）:")
for i, v in enumerate(EXPLAIN_REG):
    vf = variance_inflation_factor(X_reg, i)
    vif_reg[v] = vf
    flag = ' ★高VIF' if vf > 5 else ''
    print(f"  {LABEL_MAP[v]:<22} VIF={vf:.2f}{flag}")

# ================================================================
# ■ Step 4. 重回帰分析（4グループ）
# ================================================================

print("\n" + "=" * 65)
print("■ 重回帰分析（4グループ）")
print("=" * 65)

# サブグループ定義
# 小学校水準 vs 中学校水準: 教員比の中央値で分割
# 男子率 vs 女子率優勢: 男子率の中央値で分割
elem_med = df['小学教員児童比'].median()
jhs_med  = df['中学教員生徒比'].median()
boy_med  = df['男子率'].median()

mask_elem_hi = df['小学教員児童比'] >= elem_med  # 小学校水準が高い都道府県
mask_jhs_hi  = df['中学教員生徒比'] >= jhs_med   # 中学校水準が高い都道府県
mask_boy_hi  = df['男子率'] >= boy_med            # 男子率が高い都道府県

GROUP_MASKS = {
    '小5男子型': mask_elem_hi &  mask_boy_hi,
    '小5女子型': mask_elem_hi & ~mask_boy_hi,
    '中2男子型': mask_jhs_hi  &  mask_boy_hi,
    '中2女子型': mask_jhs_hi  & ~mask_boy_hi,
}
GROUP_LABELS = list(GROUP_MASKS.keys())

reg_results = {}
for label, mask in GROUP_MASKS.items():
    sub = df_z[mask].reset_index(drop=True)
    n_sub = len(sub)
    if n_sub < len(Z_REG) + 2:
        # サブグループが小さすぎる場合は変数を絞る
        z_use = Z_REG[:4]
    else:
        z_use = Z_REG
    X_r = sm.add_constant(sub[z_use])
    y_r = sub[TARGET + '_z']
    reg = sm.OLS(y_r, X_r).fit(cov_type='HC1')
    reg_results[label] = (reg, z_use, n_sub)
    print(f"  {label} (n={n_sub}): R²={reg.rsquared:.3f}, "
          f"adj.R²={reg.rsquared_adj:.3f}")

# ================================================================
# ■ Step 5. 教育関与度スコアと都道府県ランキング
# ================================================================

print("\n" + "=" * 65)
print("■ 都道府県別 教育関与度スコアとランキング")
print("=" * 65)

# 教育関与度スコア: 標準化した教員児童比と消費支出の合成
df['教育関与度スコア'] = (
    df_z['小学教員児童比_z'] * 0.4
    + df_z['中学教員生徒比_z'] * 0.4
    + df_z['消費支出_z'] * 0.2
)

df_rank = df[['都道府県', '教育関与度スコア', '高校進学率',
              '小学教員児童比', '中学教員生徒比', '消費支出']].copy()
df_rank = df_rank.sort_values('教育関与度スコア', ascending=False).reset_index(drop=True)
df_rank['順位'] = range(1, 48)

print("\n上位10都道府県（教育関与度スコア）:")
for _, row in df_rank.head(10).iterrows():
    print(f"  {int(row['順位']):2d}位 {row['都道府県']:<5} スコア={row['教育関与度スコア']:+.3f}  "
          f"高校進学率={row['高校進学率']:.1f}%")

print("\n下位10都道府県:")
for _, row in df_rank.tail(10).iterrows():
    print(f"  {int(row['順位']):2d}位 {row['都道府県']:<5} スコア={row['教育関与度スコア']:+.3f}  "
          f"高校進学率={row['高校進学率']:.1f}%")

# ================================================================
# ■ 図の生成（4枚）
# ================================================================

print("\n" + "=" * 65)
print("■ 図の生成（4枚）")
print("=" * 65)

# ─── 図1: 相関行列ヒートマップ ─────────────────────────────────
print("図1: 相関行列ヒートマップを作成中...")

HEAT_VARS = EXPLAIN_FULL + [TARGET]
HEAT_LABELS = [LABEL_MAP.get(v, v) for v in EXPLAIN_FULL] + ['高校進学率(%)']

corr_matrix = df[HEAT_VARS].corr()

fig1, (ax1a, ax1b) = plt.subplots(1, 2, figsize=(16, 6))
fig1.suptitle('教育指標と高校進学率の相関分析（SSDSE-B-2026 / 2022年度 / 47都道府県）',
              fontsize=12, fontweight='bold')

# 左: 相関行列ヒートマップ
im = ax1a.imshow(corr_matrix.values, cmap='RdBu_r', vmin=-1, vmax=1, aspect='auto')
plt.colorbar(im, ax=ax1a, shrink=0.8, label='Pearson r')
ax1a.set_xticks(range(len(HEAT_VARS)))
ax1a.set_xticklabels(HEAT_LABELS, rotation=45, ha='right', fontsize=8)
ax1a.set_yticks(range(len(HEAT_VARS)))
ax1a.set_yticklabels(HEAT_LABELS, fontsize=8)
ax1a.set_title('相関行列\n（赤:正の相関, 青:負の相関）', fontsize=10, fontweight='bold')
for i in range(len(HEAT_VARS)):
    for j in range(len(HEAT_VARS)):
        val = corr_matrix.values[i, j]
        color = 'white' if abs(val) > 0.6 else 'black'
        ax1a.text(j, i, f'{val:.2f}', ha='center', va='center', fontsize=7, color=color)

# 右: 高校進学率との相関係数棒グラフ
corr_with_target = [(v, corr_results[v][0], corr_results[v][1]) for v in EXPLAIN_FULL]
corr_with_target.sort(key=lambda x: x[1])
labels_sorted = [LABEL_MAP[v] for v, _, _ in corr_with_target]
vals_sorted = [r for _, r, _ in corr_with_target]
pvals_sorted = [p for _, _, p in corr_with_target]

bar_colors = []
for v, r, p in corr_with_target:
    if v in ('男子率', '女子率'):
        bar_colors.append('#FF8F00')  # 多重共線性ペアはオレンジ
    elif r < 0:
        bar_colors.append('#1565C0')
    else:
        bar_colors.append('#E53935')

bars = ax1b.barh(range(len(corr_with_target)), vals_sorted,
                 color=bar_colors, alpha=0.85, edgecolor='white')
ax1b.set_yticks(range(len(corr_with_target)))
ax1b.set_yticklabels(labels_sorted, fontsize=9)
ax1b.axvline(0, color='black', linewidth=1.0)
ax1b.set_xlabel('Pearson 相関係数（高校進学率との相関）', fontsize=10)
ax1b.set_title('各変数と高校進学率との相関係数\n（オレンジ: 多重共線性ペア）',
               fontsize=10, fontweight='bold')
ax1b.grid(axis='x', alpha=0.3)

for i, (v, r, p) in enumerate(corr_with_target):
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else ''
    if sig:
        offset = 0.02 if r > 0 else -0.02
        ha = 'left' if r > 0 else 'right'
        ax1b.text(r + offset, i, sig, va='center', ha=ha, fontsize=10, fontweight='bold')
    ax1b.text(r + (0.01 if r >= 0 else -0.01), i,
              f'{r:+.3f}', va='center', ha='left' if r >= 0 else 'right',
              fontsize=7.5, color='#333')

# 凡例
from matplotlib.patches import Patch
handles = [Patch(color='#E53935', alpha=0.85, label='正の相関'),
           Patch(color='#1565C0', alpha=0.85, label='負の相関'),
           Patch(color='#FF8F00', alpha=0.85, label='多重共線性ペア')]
ax1b.legend(handles=handles, fontsize=8, loc='lower right')

plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2024_H3_fig1_corr.png'), bbox_inches='tight', dpi=150)
plt.close(fig1)
print("  → 2024_H3_fig1_corr.png 保存完了")

# ─── 図2: 4グループの重回帰係数比較 ────────────────────────────
print("図2: 4グループの重回帰係数比較を作成中...")

fig2, axes2 = plt.subplots(2, 2, figsize=(14, 10))
fig2.suptitle('4グループ別 重回帰分析結果（目的変数: 高校進学率の標準化スコア）',
              fontsize=12, fontweight='bold')
axes2_flat = axes2.flatten()
colors4 = ['#1565C0', '#E53935', '#43A047', '#FF8F00']

for idx, (label, (reg, z_use, n_sub)) in enumerate(reg_results.items()):
    ax = axes2_flat[idx]
    color = colors4[idx]
    var_names = [v.replace('_z', '') for v in z_use]
    coefs = [reg.params.get(zv, 0) for zv in z_use]
    ses   = [reg.bse.get(zv, 0) for zv in z_use]
    pvals = [reg.pvalues.get(zv, 1) for zv in z_use]
    display_labels = [LABEL_MAP.get(v, v) for v in var_names]

    sorted_idx_r = sorted(range(len(coefs)), key=lambda i: coefs[i])
    bar_colors_r = [color if pvals[i] < 0.05 else '#BDBDBD' for i in sorted_idx_r]

    ax.barh(range(len(z_use)),
            [coefs[i] for i in sorted_idx_r],
            xerr=[1.96 * ses[i] for i in sorted_idx_r],
            color=bar_colors_r,
            alpha=0.85, edgecolor='white', capsize=3,
            error_kw={'elinewidth': 1.2, 'ecolor': '#555'})
    ax.set_yticks(range(len(z_use)))
    ax.set_yticklabels([display_labels[i] for i in sorted_idx_r], fontsize=8.5)
    ax.axvline(0, color='black', linewidth=1.0)
    ax.set_xlabel('標準化回帰係数（±95%CI）', fontsize=9)
    ax.set_title(
        f'{label}  (n={n_sub})\n'
        f'R²={reg.rsquared:.3f}  adj.R²={reg.rsquared_adj:.3f}',
        fontsize=10, fontweight='bold', color=color
    )
    ax.grid(axis='x', alpha=0.3)
    # 有意な係数にラベル
    for rank_i, orig_i in enumerate(sorted_idx_r):
        if pvals[orig_i] < 0.05:
            c = coefs[orig_i]
            offset = 0.02 if c >= 0 else -0.02
            ha = 'left' if c >= 0 else 'right'
            sig_str = ('***' if pvals[orig_i] < 0.001
                       else '**' if pvals[orig_i] < 0.01 else '*')
            ax.text(c + offset, rank_i, sig_str, va='center', ha=ha,
                    fontsize=9, fontweight='bold', color=color)

plt.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2024_H3_fig2_coef.png'), bbox_inches='tight', dpi=150)
plt.close(fig2)
print("  → 2024_H3_fig2_coef.png 保存完了")

# ─── 図3: VIF棒グラフ（フルセット vs 縮約セット）───────────────
print("図3: VIF棒グラフを作成中...")

fig3, (ax3a, ax3b) = plt.subplots(1, 2, figsize=(14, 6))
fig3.suptitle('多重共線性の確認：VIF と 偏相関係数',
              fontsize=12, fontweight='bold')

# 左: VIF フルセット vs 縮約セット
all_vars = EXPLAIN_FULL
vif_full_vals = [vif_full[v] for v in all_vars]
colors_vif = ['#E53935' if vif_full[v] > 5 else '#43A047' for v in all_vars]
labels_vif = [LABEL_MAP[v] for v in all_vars]

y_pos = range(len(all_vars))
bars3a = ax3a.barh(y_pos, [min(v, 30) for v in vif_full_vals],
                   color=colors_vif, alpha=0.85, edgecolor='white')
ax3a.axvline(5, color='#E53935', linestyle='--', linewidth=2, label='VIF=5（警戒線）')
ax3a.axvline(3, color='#FF8F00', linestyle='--', linewidth=1.5, label='VIF=3（注意線）')
ax3a.set_yticks(y_pos)
ax3a.set_yticklabels(labels_vif, fontsize=9)
ax3a.set_xlabel('VIF（分散拡大係数）', fontsize=11)
ax3a.set_title('VIF確認（フルセット）\n赤: 多重共線性ペア（VIF>5）', fontsize=10, fontweight='bold')
ax3a.legend(fontsize=8)
ax3a.grid(axis='x', alpha=0.3)
ax3a.invert_yaxis()
ax3a.set_xlim(0, 35)
for i, v in enumerate(vif_full_vals):
    label_v = f'{v:.1f}' if v <= 30 else f'{v:.0f}★'
    ax3a.text(min(v, 30) + 0.3, i, label_v, va='center', fontsize=8)

# 右: VIF 縮約セット + 偏相関まとめ
vif_reg_vals = [vif_reg[v] for v in EXPLAIN_REG]
labels_vif_r = [LABEL_MAP[v] for v in EXPLAIN_REG]
colors_vif_r = ['#E53935' if vif_reg[v] > 5 else '#FF8F00' if vif_reg[v] > 3 else '#43A047'
                for v in EXPLAIN_REG]

y_pos_r = range(len(EXPLAIN_REG))
ax3b.barh(y_pos_r, vif_reg_vals, color=colors_vif_r, alpha=0.85, edgecolor='white')
ax3b.axvline(5, color='#E53935', linestyle='--', linewidth=2, label='VIF=5（警戒線）')
ax3b.axvline(3, color='#FF8F00', linestyle='--', linewidth=1.5, label='VIF=3（注意線）')
ax3b.set_yticks(y_pos_r)
ax3b.set_yticklabels(labels_vif_r, fontsize=9)
ax3b.set_xlabel('VIF（分散拡大係数）', fontsize=11)
ax3b.set_title('VIF確認（女子率除外後）\nVIFが低下し多重共線性が解消', fontsize=10, fontweight='bold')
ax3b.legend(fontsize=8)
ax3b.grid(axis='x', alpha=0.3)
ax3b.invert_yaxis()
ax3b.set_xlim(0, 8)
for i, v in enumerate(vif_reg_vals):
    ax3b.text(v + 0.05, i, f'{v:.2f}', va='center', fontsize=9)

# 右下に偏相関情報を注釈
note_text = (
    f"【偏相関による多重共線性の定量化】\n"
    f"高齢化率・消費支出を統制した偏相関:\n"
    f"  小学教員比 × 中学教員比: r={r_p1:.3f}"
    + (" ***" if p_p1 < 0.001 else " **" if p_p1 < 0.01 else " *" if p_p1 < 0.05 else "") + "\n"
    f"  男子率 × 女子率: r={r_p2:.3f}"
    + (" ***" if p_p2 < 0.001 else " **" if p_p2 < 0.01 else " *" if p_p2 < 0.05 else "")
)
ax3b.text(0.02, 0.05, note_text, transform=ax3b.transAxes,
          fontsize=8, verticalalignment='bottom',
          bbox=dict(boxstyle='round', facecolor='#FFF9C4', alpha=0.8))

plt.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2024_H3_fig3_vif.png'), bbox_inches='tight', dpi=150)
plt.close(fig3)
print("  → 2024_H3_fig3_vif.png 保存完了")

# ─── 図4: 都道府県別 教育関与度スコアランキング ─────────────────
print("図4: 都道府県別 教育関与度スコアランキングを作成中...")

fig4, (ax4a, ax4b) = plt.subplots(1, 2, figsize=(14, 11))
fig4.suptitle('都道府県別 教育関与度スコアと高校進学率ランキング\n（SSDSE-B-2026 / 2022年度）',
              fontsize=12, fontweight='bold')

N = len(df_rank)

# 左: 教育関与度スコアランキング（横棒）
short_names = (df_rank['都道府県']
               .str.replace('県', '').str.replace('府', '')
               .str.replace('都', '').str.replace('道', ''))

score_vals = df_rank['教育関与度スコア'].values
bar_colors4a = ['#E53935' if s >= score_vals[N // 3]
                else '#FF8F00' if s >= 0
                else '#1565C0'
                for s in score_vals]
# 反転（1位が上）
score_plot  = score_vals[::-1]
names_plot  = short_names.values[::-1]
colors_plot = bar_colors4a[::-1]

ax4a.barh(range(N), score_plot, color=colors_plot, alpha=0.85, edgecolor='white')
ax4a.set_yticks(range(N))
ax4a.set_yticklabels(names_plot, fontsize=7)
ax4a.axvline(0, color='black', linewidth=1.2)
ax4a.set_xlabel('教育関与度スコア（小学教員比×0.4 + 中学教員比×0.4 + 消費支出×0.2）',
                fontsize=8.5)
ax4a.set_title('都道府県別 教育関与度スコア\n（赤: 上位群, 青: 下位群）',
               fontsize=10, fontweight='bold')
ax4a.grid(axis='x', alpha=0.3)

from matplotlib.patches import Patch
handles4a = [Patch(color='#E53935', alpha=0.85, label='上位群（スコア高）'),
             Patch(color='#FF8F00', alpha=0.85, label='中位群'),
             Patch(color='#1565C0', alpha=0.85, label='下位群（スコア低）')]
ax4a.legend(handles=handles4a, fontsize=8, loc='lower right')

# 右: 散布図（教育関与度スコア vs 高校進学率）
scatter_colors = ['#E53935' if s >= score_vals[N // 3]
                  else '#FF8F00' if s >= 0
                  else '#1565C0'
                  for s in df['教育関与度スコア'].values]

ax4b.scatter(df['教育関与度スコア'], df['高校進学率'],
             c=scatter_colors, alpha=0.75, s=60, edgecolors='white', linewidth=0.5)

# 回帰直線
x_sc = df['教育関与度スコア'].values
y_sc = df['高校進学率'].values
slope, intercept, r_sc, p_sc, _ = __import__('scipy.stats', fromlist=['linregress']).linregress(x_sc, y_sc)
x_line = [x_sc.min(), x_sc.max()]
y_line = [slope * x + intercept for x in x_line]
ax4b.plot(x_line, y_line, 'k--', linewidth=1.5, alpha=0.7,
          label=f'回帰直線 (r={r_sc:.3f}{"***" if p_sc<0.001 else "**" if p_sc<0.01 else "*" if p_sc<0.05 else ""})')

# 上位5・下位5の都道府県名ラベル
for _, row in pd.concat([df_rank.head(5), df_rank.tail(5)]).iterrows():
    pref = row['都道府県']
    row_data = df[df['都道府県'] == pref].iloc[0]
    ax4b.annotate(
        pref.replace('県','').replace('府','').replace('都','').replace('道',''),
        (row_data['教育関与度スコア'], row_data['高校進学率']),
        fontsize=7.5, fontweight='bold',
        xytext=(4, 4), textcoords='offset points',
        color='#333'
    )

ax4b.set_xlabel('教育関与度スコア', fontsize=11)
ax4b.set_ylabel('高校進学率 (%)', fontsize=11)
ax4b.set_title('教育関与度スコア vs 高校進学率\n（都道府県別散布図）',
               fontsize=10, fontweight='bold')
ax4b.legend(fontsize=9)
ax4b.grid(True, alpha=0.3)

plt.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2024_H3_fig4_rank.png'), bbox_inches='tight', dpi=150)
plt.close(fig4)
print("  → 2024_H3_fig4_rank.png 保存完了")

print("\n" + "=" * 65)
print("全図の生成完了（4枚）")
print("=" * 65)
print(f"\n保存先: {FIG_DIR}")
print("  2024_H3_fig1_corr.png  - 相関行列ヒートマップ + 相関係数棒グラフ")
print("  2024_H3_fig2_coef.png  - 4グループの重回帰係数比較")
print("  2024_H3_fig3_vif.png   - VIF棒グラフ（フルセット vs 縮約セット）")
print("  2024_H3_fig4_rank.png  - 教育関与度スコアランキングと散布図")