都道府県ごとの学力の差

df_raw = pd.read_csv(
    os.path.join(DATA_DIR, 'SSDSE-B-2026.csv'),
    encoding='cp932', header=None
)

# 行0がコード行・行1がラベル行・行2以降がデータ
df_raw.columns = df_raw.iloc[0]
df_all = df_raw.iloc[2:].copy().reset_index(drop=True)

# 都道府県コード（R\d{5}形式）に絞る
mask_pref = df_all['Code'].str.match(r'^R\d{5}$', na=False)
df_pref = df_all[mask_pref].copy()

# 2022年度のみ抽出
df_pref['年度_str'] = df_pref['SSDSE-B-2026'].astype(str).str.strip()
df = df_pref[df_pref['年度_str'] == '2022'].copy().reset_index(drop=True)

print("=" * 60)
print(f"■ 対象都道府県数: N={len(df)}（2022年度）")
print("=" * 60)

# ── 必要列を数値変換 ─────────────────────────────────────────────
NUM_COLS = [
    'E4602', 'E4601',            # 大学進学者数・高校卒業者数
    'L322108', 'L3221',          # 教育費・消費支出
    'E2401', 'E2501',            # 小学教員数・児童数
    'E3401', 'E3501',            # 中学教員数・生徒数
    'C5401',                     # 住宅地標準価格
    'A1303', 'A1101',            # 65歳以上人口・総人口
    'G7101',                     # 延べ宿泊者数
    'B4101',                     # 年平均気温
]
for c in NUM_COLS:
    df[c] = pd.to_numeric(df[c], errors='coerce')

# ── 目的変数 ─────────────────────────────────────────────────────
df['大学進学率'] = df['E4602'] / df['E4601'] * 100

# ── 説明変数 ─────────────────────────────────────────────────────
df['教育費割合']     = df['L322108'] / df['L3221'] * 100
df['消費支出']       = df['L3221']
df['小学教員児童比'] = df['E2401'] / df['E2501']
df['中学教員生徒比'] = df['E3401'] / df['E3501']
df['住宅地価格']     = df['C5401']
df['高齢化率']       = df['A1303'] / df['A1101'] * 100
df['宿泊者per capita'] = df['G7101'] / df['A1101']
df['年平均気温']     = df['B4101']

VAR_NAMES = [
    '教育費割合', '消費支出', '小学教員児童比', '中学教員生徒比',
    '住宅地価格', '高齢化率', '宿泊者per capita', '年平均気温'
]
VAR_LABELS = {
    '教育費割合':     '教育費割合（%）',
    '消費支出':       '消費支出（円/月）',
    '小学教員児童比': '小学教員/児童比',
    '中学教員生徒比': '中学教員/生徒比',
    '住宅地価格':     '住宅地価格（千円/m²）',
    '高齢化率':       '高齢化率（%）',
    '宿泊者per capita': '宿泊者per capita',
    '年平均気温':     '年平均気温（℃）',
}

# 都道府県名
pref_names = df['Prefecture'].values

# 有効データ確認
y = df['大学進学率'].values
X_df = df[VAR_NAMES].copy()

valid_mask = (
    np.isfinite(y) &
    X_df.apply(lambda col: np.isfinite(col)).all(axis=1).values
)
df_valid = df[valid_mask].reset_index(drop=True)
y = df_valid['大学進学率'].values
X_df = df_valid[VAR_NAMES].copy()
pref_names = df_valid['Prefecture'].values

N = len(y)
print(f"有効都道府県数（欠損除外後）: N={N}")

print()
print(df_valid[['Prefecture', '大学進学率'] + VAR_NAMES].describe().round(3).to_string())

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

sort_idx = np.argsort(y)
y_sorted = y[sort_idx]
names_sorted = pref_names[sort_idx]

top10_threshold = np.sort(y)[-10]
colors_bar = ['#C62828' if v >= top10_threshold else '#1565C0' for v in y_sorted]
alphas_bar = [0.95 if v >= top10_threshold else 0.65 for v in y_sorted]

bars = ax1.barh(
    np.arange(N), y_sorted,
    color=colors_bar,
    alpha=0.75,
    edgecolor='white',
    height=0.7
)
for bar, alpha in zip(bars, alphas_bar):
    bar.set_alpha(alpha)

ax1.set_yticks(np.arange(N))
ax1.set_yticklabels(names_sorted, fontsize=9)
ax1.set_xlabel('大学進学率（%）', fontsize=12)
ax1.set_title('都道府県別 大学進学率（2022年度）\n（赤：上位10都道府県）',
              fontsize=13, fontweight='bold', pad=12)
ax1.axvline(y.mean(), color='#FF8F00', linestyle='--', linewidth=1.5,
            label=f'全国平均 {y.mean():.1f}%')
ax1.grid(axis='x', alpha=0.3)
ax1.legend(fontsize=10)

for i, (v, name) in enumerate(zip(y_sorted, names_sorted)):
    ax1.text(v + 0.3, i, f'{v:.1f}%', va='center', fontsize=7.5,
             color='#C62828' if v >= top10_threshold else '#555')

plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2023_H5_2_fig1_bar.png'), bbox_inches='tight', dpi=150)
plt.close(fig1)
print("\n図1保存: 2023_H5_2_fig1_bar.png")

fig2, ax2 = plt.subplots(figsize=(9, 5))

corr_list = [corr_vals[n] for n in VAR_NAMES]
pval_list = [pval_vals[n] for n in VAR_NAMES]
colors2 = ['#C62828' if p < 0.05 else '#9E9E9E' for p in pval_list]
y_pos2 = np.arange(len(VAR_NAMES))

ax2.barh(y_pos2, corr_list, color=colors2, alpha=0.75, edgecolor='white', height=0.6)
ax2.axvline(0, color='gray', linestyle='--', linewidth=1.0)
ax2.set_yticks(y_pos2)
ax2.set_yticklabels(VAR_NAMES, fontsize=10)
ax2.set_xlabel('ピアソン相関係数', fontsize=11)
ax2.set_title('各説明変数と大学進学率の相関係数\n（無相関検定: 赤=p<0.05）',
              fontsize=12, fontweight='bold')
ax2.invert_yaxis()
ax2.grid(axis='x', alpha=0.3)

from matplotlib.patches import Patch
ax2.legend(
    handles=[
        Patch(color='#C62828', alpha=0.75, label='有意（p<0.05）'),
        Patch(color='#9E9E9E', alpha=0.75, label='非有意'),
    ],
    fontsize=9
)

for i, (r, p) in enumerate(zip(corr_list, pval_list)):
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else ''
    offset = 0.01 if r >= 0 else -0.01
    ha = 'left' if r >= 0 else 'right'
    ax2.text(r + offset, i, f' {r:.3f}{sig}', va='center', ha=ha, fontsize=8.5)

plt.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2023_H5_2_fig2_corr.png'), bbox_inches='tight', dpi=150)
plt.close(fig2)
print("図2保存: 2023_H5_2_fig2_corr.png")

fig3, ax3 = plt.subplots(figsize=(9, max(4, len(selected_cols) * 0.9)))

colors3 = ['#C62828' if p < 0.01 else '#FF8F00' if p < 0.05 else '#9E9E9E'
           for p in beta_pvals]
y_pos3 = np.arange(len(selected_cols))

ax3.barh(y_pos3, beta_coefs, color=colors3, alpha=0.75, edgecolor='white', height=0.6)
ax3.errorbar(beta_coefs, y_pos3, xerr=1.96 * beta_ses,
             fmt='none', color='#333', capsize=4, linewidth=1.2)
ax3.axvline(0, color='gray', linestyle='--', linewidth=1.0)
ax3.set_yticks(y_pos3)
ax3.set_yticklabels(selected_cols, fontsize=10)
ax3.set_xlabel('標準化係数 β（±1.96SE）', fontsize=11)
ax3.set_title(
    f'バックワードAIC選択後の重回帰 — 標準化係数\n'
    f'R²={final_model.rsquared:.3f}  Adj.R²={final_model.rsquared_adj:.3f}  N={N}都道府県',
    fontsize=12, fontweight='bold'
)
ax3.invert_yaxis()
ax3.grid(axis='x', alpha=0.3)

ax3.legend(
    handles=[
        Patch(color='#C62828', alpha=0.75, label='p<0.01'),
        Patch(color='#FF8F00', alpha=0.75, label='p<0.05'),
        Patch(color='#9E9E9E', alpha=0.75, label='n.s.'),
    ],
    fontsize=9
)

for i, (b, p) in enumerate(zip(beta_coefs, beta_pvals)):
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else ''
    offset = 0.01 if b >= 0 else -0.01
    ha = 'left' if b >= 0 else 'right'
    ax3.text(b + offset, i, f' {b:.3f}{sig}', va='center', ha=ha, fontsize=8.5)

plt.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2023_H5_2_fig3_beta.png'), bbox_inches='tight', dpi=150)
plt.close(fig3)
print("図3保存: 2023_H5_2_fig3_beta.png")

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

print("\n" + "=" * 60)
print("■ Step1. ピアソン相関係数（目的変数：大学進学率）")
print("=" * 60)
print(f"\n  {'変数':<20} {'r':>8} {'p値':>10} {'有意':>6}")
print("  " + "-" * 48)
corr_vals = {}
pval_vals = {}
for name in VAR_NAMES:
    r, p = stats.pearsonr(X_df[name].values, y)
    corr_vals[name] = r
    pval_vals[name] = p
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else 'n.s.'
    print(f"  {name:<20} {r:>8.4f} {p:>10.4f} {sig:>6}")

def backward_aic(X_df, y):
    """バックワードAIC変数選択。AICが改善しなくなるまで変数を除外する。"""
    cols = list(X_df.columns)
    step = 0
    print("\n" + "=" * 60)
    print("■ Step2. バックワードAIC変数選択")
    print("=" * 60)
    while len(cols) > 1:
        full_model = sm.OLS(y, sm.add_constant(X_df[cols])).fit()
        full_aic = full_model.aic
        aic_drop = {}
        for c in cols:
            sub = [x for x in cols if x != c]
            m = sm.OLS(y, sm.add_constant(X_df[sub])).fit()
            aic_drop[c] = m.aic
        best_drop = min(aic_drop, key=aic_drop.get)
        if aic_drop[best_drop] < full_aic:
            step += 1
            print(f"  Step{step}: '{best_drop}' を除外"
                  f"（AIC: {full_aic:.3f} → {aic_drop[best_drop]:.3f}）")
            cols.remove(best_drop)
        else:
            break
    return cols

selected_cols = backward_aic(X_df, y)
print(f"\n  最終選択変数: {selected_cols}")

# ── 最終モデルの推定 ──────────────────────────────────────────────
X_final = sm.add_constant(X_df[selected_cols])
final_model = sm.OLS(y, X_final).fit()
print(f"\n  R² = {final_model.rsquared:.4f}")
print(f"  Adj. R² = {final_model.rsquared_adj:.4f}")
print(f"  AIC = {final_model.aic:.3f}")
print()
print(final_model.summary2())

# ── 標準化係数（βeta） ──────────────────────────────────────────
y_std = (y - y.mean()) / y.std()
X_sel_std = X_df[selected_cols].copy()
for c in selected_cols:
    X_sel_std[c] = (X_sel_std[c] - X_sel_std[c].mean()) / X_sel_std[c].std()

std_model = sm.OLS(y_std, sm.add_constant(X_sel_std)).fit()
beta_coefs = std_model.params[1:].values  # 定数項除く
beta_ses   = std_model.bse[1:].values
beta_pvals = std_model.pvalues[1:].values

print("\n標準化係数（β）:")
for name, b, p in zip(selected_cols, beta_coefs, beta_pvals):
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else 'n.s.'
    print(f"  {name:<20}: β={b:.4f}  p={p:.4f}  {sig}")

fig4, ax4 = plt.subplots(figsize=(11, 7))

x4 = X_df['教育費割合'].values
r4, p4 = stats.pearsonr(x4, y)

ax4.scatter(x4, y, color='#1565C0', s=50, alpha=0.7, zorder=3, edgecolors='white', linewidth=0.5)

# 回帰直線
z4 = np.polyfit(x4, y, 1)
xs4 = np.linspace(x4.min() - 0.05, x4.max() + 0.05, 200)
ax4.plot(xs4, np.poly1d(z4)(xs4), color='#C62828', linewidth=1.8, alpha=0.85,
         label=f'回帰直線  r={r4:.3f}（p={p4:.4f}）')

# 都道府県ラベル
for xi, yi, name in zip(x4, y, pref_names):
    # 短縮ラベル
    label = name.replace('都', '').replace('道', '').replace('府', '').replace('県', '')
    ax4.text(xi + 0.01, yi + 0.15, label, fontsize=7, color='#333',
             ha='left', va='bottom', alpha=0.85)

ax4.set_xlabel('教育費割合（消費支出に占める割合 %）', fontsize=12)
ax4.set_ylabel('大学進学率（%）', fontsize=12)
ax4.set_title('教育費割合と大学進学率の関係（2022年度・47都道府県）',
              fontsize=13, fontweight='bold')
ax4.legend(fontsize=10)
ax4.grid(True, alpha=0.3)

plt.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2023_H5_2_fig4_scatter.png'), bbox_inches='tight', dpi=150)
plt.close(fig4)
print("図4保存: 2023_H5_2_fig4_scatter.png")

print("\n全図の生成完了（4枚）")
print("  fig1_bar.png    : 都道府県別大学進学率（横棒グラフ）")
print("  fig2_corr.png   : 各変数との相関係数棒グラフ")
print("  fig3_beta.png   : 最終OLSモデル標準化係数")
print("  fig4_scatter.png: 教育費割合 vs 大学進学率 散布図")