都道府県別学力格差の要因分析
社会経済環境と教育支出の影響

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)

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()

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 保存完了")

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 保存完了")

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 保存完了")

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()

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()

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")