自己肯定感を高める要因となる
共通特徴の探求

df_b_raw = pd.read_csv(os.path.join(DATA_DIR, 'SSDSE-B-2026.csv'),
                       encoding='cp932', header=1)
df_b = df_b_raw[
    (df_b_raw['年度'] == 2022) &
    df_b_raw['地域コード'].str.match(r'^R\d{5}$', na=False)
].copy().reset_index(drop=True)

df_e_raw = pd.read_csv(os.path.join(DATA_DIR, 'SSDSE-E-2026.csv'),
                       encoding='cp932', header=None)
col_e = df_e_raw.iloc[2].tolist()
df_e = df_e_raw.iloc[3:].copy()
df_e.columns = col_e
df_e = df_e[df_e['都道府県'] != '全国'].reset_index(drop=True)

df_d_raw = pd.read_csv(os.path.join(DATA_DIR, 'SSDSE-D-2023.csv'),
                       encoding='cp932', header=1)
df_d = df_d_raw[
    (df_d_raw['男女の別'] == '0_総数') &
    df_d_raw['地域コード'].str.match(r'^R\d{5}$', na=False) &
    (df_d_raw['地域コード'] != 'R00000')
].copy().reset_index(drop=True)

def to_num(series):
    return pd.to_numeric(series, errors='coerce')

PREFS = df_b['都道府県'].tolist()
N = len(PREFS)

# 目的変数: 大学進学率
y_raw = to_num(df_b['高等学校卒業者のうち進学者数']) / to_num(df_b['高等学校卒業者数']) * 100
y = y_raw.values

# 説明変数
pop_total   = to_num(df_b['総人口'])
pop_65up    = to_num(df_b['65歳以上人口'])
pop_1564    = to_num(df_b['15～64歳人口'])
婚姻件数     = to_num(df_b['婚姻件数'])
保育所数     = to_num(df_b['保育所等数'])
気温         = to_num(df_b['年平均気温'])

所得         = to_num(df_e['1人当たり県民所得（平成27年基準）'])
面積ha       = to_num(df_e['総面積（北方地域及び竹島を除く）'])

睡眠時間     = to_num(df_d['睡眠'])           # 分/日
趣味娯楽時間 = to_num(df_d['趣味・娯楽'])     # 分/日（SNS・動画の代替指標）
通勤通学時間 = to_num(df_d['通勤・通学'])     # 分/日

高齢化率     = pop_65up / pop_total * 100
保育所数千対 = 保育所数 / pop_total * 1000
婚姻率       = 婚姻件数 / pop_1564 * 1000
ln_所得      = np.log(所得)
面積km2      = 面積ha / 100
人口密度     = pop_total / 面積km2
ln_人口密度  = np.log(人口密度)

VAR_NAMES = [
    '睡眠時間(分/日)', '趣味・娯楽時間(分/日)', '通勤通学時間(分/日)',
    'ln(県民所得)', '高齢化率', '保育所数千対', '婚姻率',
    '年平均気温', 'ln(人口密度)',
]

X_raw = np.column_stack([
    睡眠時間.values, 趣味娯楽時間.values, 通勤通学時間.values,
    ln_所得.values, 高齢化率.values, 保育所数千対.values,
    婚姻率.values, 気温.values, ln_人口密度.values,
])

df = pd.DataFrame(X_raw, columns=VAR_NAMES)
df['大学進学率'] = y
df['都道府県'] = PREFS

print("=" * 60)
print("■ 記述統計")
print("=" * 60)
print(df[['大学進学率'] + VAR_NAMES].describe().round(2))

fig1, ax1 = plt.subplots(figsize=(11, 9))
all_vars = VAR_NAMES + ['大学進学率']
corr_mat = df[all_vars].corr().values
im = ax1.imshow(corr_mat, cmap='RdBu_r', vmin=-1, vmax=1, aspect='auto')
plt.colorbar(im, ax=ax1, fraction=0.046, pad=0.04)
ax1.set_xticks(range(len(all_vars)))
ax1.set_yticks(range(len(all_vars)))
ax1.set_xticklabels(all_vars, rotation=45, ha='right', fontsize=9)
ax1.set_yticklabels(all_vars, fontsize=9)
for i in range(len(all_vars)):
    for j in range(len(all_vars)):
        ax1.text(j, i, f'{corr_mat[i,j]:.2f}', ha='center', va='center',
                 fontsize=7, color='white' if abs(corr_mat[i,j]) > 0.5 else 'black')
ax1.set_title(
    '変数間の相関行列\n（目的変数：大学進学率, SSDSE-B 2022年度）',
    fontsize=12, fontweight='bold'
)
plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2025_H5_2_fig1_corr.png'), bbox_inches='tight', dpi=150)
plt.close(fig1)
print("\n図1保存: 2025_H5_2_fig1_corr.png")

fig2, axes2 = plt.subplots(1, 2, figsize=(13, 5))
fig2.suptitle('AICステップワイズ変数選択と最終モデル', fontsize=13, fontweight='bold')

# (a) AIC推移
ax2a = axes2[0]
init_aic, _ = calc_aic(y, X_raw)
aic_steps = [init_aic] + [h[1] for h in history]
labels = ['全変数'] + [f'−{h[0][:6]}' for h in history]
ax2a.plot(range(len(aic_steps)), aic_steps, 'o-', color='#1565C0', linewidth=2, markersize=8)
ax2a.set_xticks(range(len(aic_steps)))
ax2a.set_xticklabels(labels, rotation=35, ha='right', fontsize=8)
ax2a.set_ylabel('AIC', fontsize=11)
ax2a.set_title('後退法AICの推移\n（低いほど良いモデル）', fontsize=11, fontweight='bold')
ax2a.grid(True, alpha=0.3)
if aic_steps:
    min_idx = aic_steps.index(min(aic_steps))
    ax2a.axvline(min_idx, color='red', linestyle='--', linewidth=1.2,
                 alpha=0.7, label='最小AIC')
    ax2a.legend(fontsize=9)

# (b) 最終モデル係数（標準化）
ax2b = axes2[1]
ss = StandardScaler()
X_std = ss.fit_transform(X_final)
std_model = sm.OLS(y, sm.add_constant(X_std)).fit()
std_coefs = std_model.params[1:]
std_pvals = std_model.pvalues[1:]
bar_cols = ['#C62828' if p < 0.05 else '#1565C0' if p < 0.1 else '#9E9E9E' for p in std_pvals]
y_pos = np.arange(len(current_vars))
ax2b.barh(y_pos, std_coefs, color=bar_cols, alpha=0.75, edgecolor='white', height=0.6)
ax2b.axvline(0, color='gray', linestyle='--', linewidth=1.0)
ax2b.set_yticks(y_pos)
ax2b.set_yticklabels(selected_vars, fontsize=9)
ax2b.set_xlabel('標準化回帰係数', fontsize=11)
ax2b.set_title(f'最終モデルの係数（標準化）\nR²={final_model.rsquared:.3f}',
               fontsize=11, fontweight='bold')
ax2b.invert_yaxis()
ax2b.grid(axis='x', alpha=0.3)
from matplotlib.patches import Patch
legend_els = [Patch(color='#C62828', alpha=0.75, label='p<0.05'),
              Patch(color='#1565C0', alpha=0.75, label='p<0.10'),
              Patch(color='#9E9E9E', alpha=0.75, label='非有意')]
ax2b.legend(handles=legend_els, fontsize=8, loc='lower right')

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

fig3, axes3 = plt.subplots(2, 2, figsize=(11, 9))
fig3.suptitle('主要変数と大学進学率の散布図（SSDSE実データ, 47都道府県）',
              fontsize=13, fontweight='bold')

# 相関の絶対値上位4変数を選ぶ
sorted_corr = sorted(corr_results, key=lambda x: abs(x[1]), reverse=True)
top4_names = [r[0] for r in sorted_corr[:4]]
top4_idx   = [VAR_NAMES.index(n) for n in top4_names]

for ax, idx in zip(axes3.flat, top4_idx):
    xv = X_raw[:, idx]
    r, p = stats.pearsonr(xv, y)
    ax.scatter(xv, y, color='#1565C0', s=50, alpha=0.65,
               edgecolors='white', linewidth=0.5)
    z = np.polyfit(xv, y, 1)
    xs = np.linspace(xv.min(), xv.max(), 100)
    ax.plot(xs, np.poly1d(z)(xs), 'r-', linewidth=1.8, alpha=0.8)
    ax.set_xlabel(VAR_NAMES[idx], fontsize=10)
    ax.set_ylabel('大学進学率（%）', fontsize=10)
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else 'n.s.'
    ax.set_title(f'r = {r:.3f} {sig}', fontsize=11, fontweight='bold')
    ax.grid(True, alpha=0.3)

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

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
from sklearn.preprocessing import StandardScaler
import warnings
warnings.filterwarnings('ignore')

plt.rcParams['font.family'] = 'Hiragino Sans'
plt.rcParams['axes.unicode_minus'] = False
plt.rcParams['figure.dpi'] = 150

BASE_DIR = os.path.join(_script_dir, '..')
FIG_DIR  = os.path.join(BASE_DIR, 'html', 'figures')
DATA_DIR = os.path.join(BASE_DIR, 'data', 'raw')

print("\n" + "=" * 60)
print("■ Step1. 相関分析（ピアソン相関 + 無相関検定）")
print("=" * 60)
print(f"\n  {'変数':<22} {'相関係数':>8} {'p値':>10} {'有意':>6}")
print("  " + "-" * 50)
corr_results = []
for name, col in zip(VAR_NAMES, X_raw.T):
    r, p = stats.pearsonr(col, y)
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else 'n.s.'
    corr_results.append((name, r, p))
    print(f"  {name:<22} {r:>8.4f} {p:>10.4f} {sig:>6}")

print("\n" + "=" * 60)
print("■ Step2. AICステップワイズ後退法")
print("=" * 60)

def calc_aic(y_vec, X_mat):
    model = sm.OLS(y_vec, sm.add_constant(X_mat)).fit()
    return model.aic, model

current_vars = list(range(len(VAR_NAMES)))
best_aic, _ = calc_aic(y, X_raw[:, current_vars])
print(f"  初期AIC（全変数）: {best_aic:.2f}")

history = []
while len(current_vars) > 1:
    removed = None
    for i in current_vars:
        trial_vars = [v for v in current_vars if v != i]
        trial_aic, _ = calc_aic(y, X_raw[:, trial_vars])
        if trial_aic < best_aic:
            best_aic = trial_aic
            removed = i
    if removed is not None:
        current_vars.remove(removed)
        history.append((VAR_NAMES[removed], best_aic))
        print(f"  削除: {VAR_NAMES[removed]:<24} → AIC = {best_aic:.2f}")
    else:
        break

selected_vars = [VAR_NAMES[i] for i in current_vars]
print(f"\n  最終選択変数: {selected_vars}")

print("\n" + "=" * 60)
print("■ Step3. 最終重回帰モデル")
print("=" * 60)

X_final = X_raw[:, current_vars]
final_model = sm.OLS(y, sm.add_constant(X_final)).fit()
print(final_model.summary2())
print(f"\n  R² = {final_model.rsquared:.4f}")

print("\n" + "=" * 60)
print("■ Step4. 外れ値除外（IQR×2法）")
print("=" * 60)

Q1 = np.percentile(y, 25)
Q3 = np.percentile(y, 75)
IQR = Q3 - Q1
lower = Q1 - 2 * IQR
upper = Q3 + 2 * IQR
mask = (y >= lower) & (y <= upper)
n_out = (~mask).sum()
print(f"  IQR={IQR:.2f}, 範囲=[{lower:.2f}, {upper:.2f}]")
print(f"  除外観測数: {n_out}件")
if n_out > 0:
    print(f"  除外都道府県: {[PREFS[i] for i in np.where(~mask)[0]]}")
clean_model = sm.OLS(y[mask], sm.add_constant(X_final[mask])).fit()
print(f"  除外後 R² = {clean_model.rsquared:.4f}")

os.makedirs(FIG_DIR, exist_ok=True)

fig4, axes4 = plt.subplots(1, 2, figsize=(13, 5))
fig4.suptitle('外れ値診断（IQR×2法）と除外前後の比較', fontsize=13, fontweight='bold')

ax4a = axes4[0]
ax4a.boxplot(y, vert=True, patch_artist=True,
             boxprops=dict(facecolor='#BBDEFB', color='#1565C0'),
             medianprops=dict(color='#C62828', linewidth=2))
ax4a.axhline(upper, color='#C62828', linestyle='--', linewidth=1.5,
             label=f'上限 Q3+2×IQR={upper:.1f}%')
ax4a.axhline(lower, color='#C62828', linestyle='--', linewidth=1.5,
             label=f'下限 Q1−2×IQR={lower:.1f}%')
outlier_vals = y[~mask]
if len(outlier_vals) > 0:
    ax4a.scatter([1]*len(outlier_vals), outlier_vals,
                 color='red', s=60, zorder=5, label='外れ値')
ax4a.set_ylabel('大学進学率（%）', fontsize=11)
ax4a.set_title('IQR×2法による外れ値検出', fontsize=11, fontweight='bold')
ax4a.legend(fontsize=9)
ax4a.grid(axis='y', alpha=0.3)

ax4b = axes4[1]
r2_vals = [final_model.rsquared, clean_model.rsquared]
bar_labels = [f'全データ\n(N={N})', f'外れ値除外後\n(N={mask.sum()})']
bars = ax4b.bar(bar_labels, r2_vals,
                color=['#1565C0', '#2E7D32'], alpha=0.75, edgecolor='white', width=0.4)
ax4b.set_ylabel('決定係数 R²', fontsize=11)
ax4b.set_title('外れ値除外前後のR²比較', fontsize=11, fontweight='bold')
ax4b.set_ylim(0, 1.0)
ax4b.grid(axis='y', alpha=0.3)
for bar, val in zip(bars, r2_vals):
    ax4b.text(bar.get_x() + bar.get_width()/2, val + 0.02, f'{val:.3f}',
              ha='center', fontsize=12, fontweight='bold')

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

print("\n全図の生成完了（4枚）")