""" 教育用再現コード: 2025年 統計データ分析コンペティション 審査員奨励賞 [高校生の部] ================================================================= 論文タイトル：自己肯定感を高める要因となる共通特徴の探求 著者：慶應義塾湘南藤沢高等部 【データ出典】 目的変数：大学進学率 ＝ 高等学校卒業者のうち進学者数 / 高等学校卒業者数 × 100 出典: SSDSE-B-2026.csv（政府統計SSDSE, 2022年度） ※全国学力・学習状況調査（自己肯定感割合）は都道府県別の集計値が非公開のため、 教育的達成・向上意欲の代理指標として大学進学率を採用。 説明変数: SSDSE-D-2023.csv（社会生活基本調査 2021年度, 単位: 分/日） - 睡眠時間（総数） - 趣味・娯楽時間（SNS・動画視聴等の代替指標） - 通勤・通学時間 SSDSE-B-2026.csv（2022年度） - 高齢化率 = 65歳以上人口 / 総人口 × 100 - 保育所数千対 = 保育所等数 / 総人口 × 1,000 - 婚姻率 = 婚姻件数 / 15～64歳人口 × 1,000 - 年平均気温 SSDSE-E-2026.csv - 1人当たり県民所得（平成27年基準, 万円） - ln(人口密度) = log(総人口 / 総面積[km²]) 【分析概要】 Step1. 相関分析（ピアソン相関 + 無相関検定） Step2. AICステップワイズ変数選択（後退法） Step3. 重回帰分析（最終モデル） Step4. 外れ値除外（IQR×2法）と再推定 【データサイエンス学習ポイント】 1. ピアソン相関と無相関検定（scipy.stats.pearsonr） 2. AIC（赤池情報量基準）による変数選択 3. IQR法による外れ値検出 4. 重回帰分析の解釈（標準化係数） ================================================================= """ # ============================================================ # 【データの準備】実行前に以下のデータファイルを用意してください # # 必要ファイル: # ・SSDSE-B-2026.csv # → data/raw/SSDSE-B-2026.csv に配置 # ・SSDSE-D-2023.csv # → data/raw/SSDSE-D-2023.csv に配置 # ・SSDSE-E-2026.csv # → data/raw/SSDSE-E-2026.csv に配置 # # ダウンロード先: # https://www.nstac.go.jp/use/literacy/ssdse/ # （SSDSE-B（社会・人口統計体系 都道府県データ） の CSV をダウンロード） # （SSDSE-D（社会・人口統計体系 都道府県の指標） の CSV をダウンロード） # （SSDSE-E（社会・人口統計体系 都道府県の指標2） の CSV をダウンロード） # # フォルダ配置（プロジェクトルートからの相対パス）: # code/ ← このスクリプトの場所 # data/raw/ ← CSV ファイルをここに配置 # html/figures/ ← 図の出力先（自動生成） # # 実行方法（ファイルを一切編集せず実行可能）: # python3 code/2025_H5_2_shorei.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 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') # ================================================================ # ■ SSDSE-B-2026 読み込み（2022年度, 47都道府県） # ================================================================ 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) # ================================================================ # ■ SSDSE-E-2026 読み込み（47都道府県） # ================================================================ 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) # ================================================================ # ■ SSDSE-D-2023 読み込み（社会生活基本調査, 総数, 47都道府県） # ================================================================ 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)) # ================================================================ # ■ Step1. 相関分析（無相関検定） # ================================================================ 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}") # ================================================================ # ■ Step2. AICステップワイズ変数選択（後退法） # ================================================================ 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}") # ================================================================ # ■ Step3. 最終重回帰モデル # ================================================================ 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}") # ================================================================ # ■ Step4. 外れ値除外（IQR×2法） # ================================================================ 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}") # ================================================================ # ■ 図の生成（4枚） # ================================================================ os.makedirs(FIG_DIR, exist_ok=True) # ──────────────────────────────────────────────────────────────── # 図1: 相関行列ヒートマップ # ──────────────────────────────────────────────────────────────── 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") # ──────────────────────────────────────────────────────────────── # 図2: AICステップワイズ選択の推移 + 最終モデル係数 # ──────────────────────────────────────────────────────────────── 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") # ──────────────────────────────────────────────────────────────── # 図3: 主要変数との散布図（2×2） # ──────────────────────────────────────────────────────────────── 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") # ──────────────────────────────────────────────────────────────── # 図4: IQR外れ値診断 + 外れ値除外前後のR²比較 # ──────────────────────────────────────────────────────────────── 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枚）")