""" 教育用再現コード: 2025年 統計データ分析コンペティション 審査員奨励賞 [高校生の部] ================================================================= 論文タイトル：消滅可能性自治体を救うには 著者：早稲田大学系属 早稲田実業学校 【分析概要】 データ：SSDSE-A-2025.csv（市区町村データ） 人口5,000〜20,000人の小規模自治体 523自治体 目的変数：15〜64歳女性人口比率（働き手女性の地域定着指標） Step1. 無相関検定（各説明変数との相関） Step2. 重回帰分析（バックワード変数選択：全p<0.05まで除外） Step3. 箱ひげ図（外れ値IQR法） Step4. 事例分析（川島町） 【変数】（人口1万人あたり換算） 保育所数, 幼稚園数, 小学校数, 中学校数, 製造業従業者割合, 農業従業者割合, 女性就業率（15〜64歳女性対）， 一般病院数, 一般診療所数, 農家比率（販売農家） 【データ出典】 SSDSE-A-2025.csv: 社会・人口統計体系（市区町村データ） 統計センターより公表の実データ 【データサイエンス学習ポイント】 1. 無相関検定と相関係数の解釈 2. バックワード変数選択の手順（p値基準） 3. 箱ひげ図によるデータの分布把握 4. 自治体レベルデータの特性（地域差の大きさ） ================================================================= """ # ============================================================ # 【データの準備】実行前に以下のデータファイルを用意してください # # 必要ファイル: # ・SSDSE-A-2025.csv # → data/raw/SSDSE-A-2025.csv に配置 # # ダウンロード先: # https://www.nstac.go.jp/use/literacy/ssdse/ # （SSDSE-A（社会・人口統計体系 都道府県時系列） の CSV をダウンロード） # # フォルダ配置（プロジェクトルートからの相対パス）: # code/ ← このスクリプトの場所 # data/raw/ ← CSV ファイルをここに配置 # html/figures/ ← 図の出力先（自動生成） # # 実行方法（ファイルを一切編集せず実行可能）: # python3 code/2025_H5_3_shorei.py # ============================================================ 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 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) # ================================================================ # ■ 実データ読み込み（SSDSE-A-2025: 市区町村データ） # ================================================================ df_a_raw = pd.read_csv( os.path.join(DATA_DIR, 'SSDSE-A-2025.csv'), encoding='cp932', header=1 ) df_a = df_a_raw.iloc[1:].copy() df_a.columns = df_a_raw.iloc[0].values # 数値変換 num_cols = [ '総人口', '総人口（女）', '15歳未満人口（女）', '15～64歳人口（女）', '65歳以上人口（女）', '75歳以上人口（女）', '就業者数', '就業者数（女）', '第１次産業就業者数', '第２次産業就業者数', '第３次産業就業者数', '従業者数（民営）（製造業）', '保育所等数（基本票）', '幼稚園数', '小学校数', '中学校数', '高等学校数', '一般病院数', '一般診療所数', '医師数', '農家数（販売農家）', '農家数（自給的農家）', ] for c in num_cols: if c in df_a.columns: df_a[c] = pd.to_numeric(df_a[c], errors='coerce') # 人口5,000〜20,000人の自治体に絞る（消滅可能性自治体の規模） df = df_a[ (df_a['総人口'] >= 5000) & (df_a['総人口'] <= 20000) ].copy().reset_index(drop=True) print("=" * 60) print(f"■ 対象自治体数: N={len(df)}（人口5,000〜20,000人）") print("=" * 60) # ── 特徴量エンジニアリング ───────────────────────────────────── pop = df['総人口'].values.clip(1) pop_f = df['総人口（女）'].values.clip(1) emp = df['就業者数'].values.clip(1) # 目的変数: 15〜64歳女性人口比率（%）— 消滅可能性に直結する指標 y_vals = df['15～64歳人口（女）'].values / pop_f * 100 # 説明変数（人口1万人あたり等） nursery = df['保育所等数（基本票）'].values / pop * 10000 # 保育所数（万人対） kindergarten = df['幼稚園数'].values / pop * 10000 # 幼稚園数（万人対） elementary = df['小学校数'].values / pop * 10000 # 小学校数（万人対） junior_high = df['中学校数'].values / pop * 10000 # 中学校数（万人対） mfg_ratio = df['従業者数（民営）（製造業）'].values / emp * 100 # 製造業従業者割合（%） agri_ratio = df['第１次産業就業者数'].values / emp * 100 # 農業従業者割合（%） female_emp_rate = df['就業者数（女）'].values / df['15～64歳人口（女）'].values.clip(1) * 100 # 女性就業率 hospital = df['一般病院数'].values / pop * 10000 # 一般病院数（万人対） clinic = df['一般診療所数'].values / pop * 10000 # 一般診療所数（万人対） farmer_rate = df['農家数（販売農家）'].values / pop * 10000 # 販売農家数（万人対） VAR_NAMES = [ '保育所数', '幼稚園数', '小学校数', '中学校数', '製造業従業者割合', '農業従業者割合', '女性就業率', '一般病院数', '一般診療所数', '農家比率' ] X_raw = np.column_stack([ nursery, kindergarten, elementary, junior_high, mfg_ratio, agri_ratio, female_emp_rate, hospital, clinic, farmer_rate ]) # 欠損を含む行を除外 valid = np.all(np.isfinite(X_raw), axis=1) & np.isfinite(y_vals) X_raw = X_raw[valid] y_vals = y_vals[valid] df_valid = df[valid].reset_index(drop=True) N = len(y_vals) print(f"有効自治体数（欠損除外後）: N={N}") # ── 川島町（埼玉県）のインデックスを特定 ───────────────────────────── kawashima_mask = df_valid['市区町村'].str.contains('川島', na=False) kawashima_idx = df_valid[kawashima_mask].index.tolist() has_kawashima = len(kawashima_idx) > 0 if has_kawashima: kaw_idx = kawashima_idx[0] print(f"川島町: index={kaw_idx}, 15-64歳女性比率={y_vals[kaw_idx]:.2f}%") print(f" 保育所数={nursery[valid][kaw_idx]:.2f}, 製造業割合={mfg_ratio[valid][kaw_idx]:.2f}%") else: print("川島町は対象外（人口範囲外または欠損）") df_main = pd.DataFrame(X_raw, columns=VAR_NAMES) df_main['若年女性比率'] = y_vals print() print(df_main[['若年女性比率'] + VAR_NAMES].describe().round(2)) # IQR外れ値チェック Q1 = df_main['若年女性比率'].quantile(0.25) Q3 = df_main['若年女性比率'].quantile(0.75) IQR = Q3 - Q1 lower_iqr = Q1 - 1.5 * IQR upper_iqr = Q3 + 1.5 * IQR mask_iqr = (df_main['若年女性比率'] >= lower_iqr) & (df_main['若年女性比率'] <= upper_iqr) n_out = (~mask_iqr).sum() print(f"\n IQR法外れ値: {n_out}件（IQR={IQR:.2f}, 範囲=[{lower_iqr:.2f}, {upper_iqr:.2f}]）") # ================================================================ # ■ Step1. 無相関検定 # ================================================================ print("\n" + "=" * 60) print("■ Step1. 無相関検定（ピアソン相関）") print("=" * 60) print(f"\n {'変数':<18} {'r':>8} {'p値':>10} {'有意':>6}") print(" " + "-" * 46) for name, col in zip(VAR_NAMES, X_raw.T): r, p = stats.pearsonr(col, y_vals) sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else 'n.s.' print(f" {name:<18} {r:>8.4f} {p:>10.4f} {sig:>6}") # ================================================================ # ■ Step2. バックワード変数選択（p<0.05） # ================================================================ print("\n" + "=" * 60) print("■ Step2. バックワード変数選択（閾値 p < 0.05）") print("=" * 60) current_vars = list(range(len(VAR_NAMES))) step = 0 while True: X_cur = sm.add_constant(X_raw[:, current_vars]) model = sm.OLS(y_vals, X_cur).fit() pvals = np.asarray(model.pvalues[1:]) # 定数項を除く（numpy配列に変換） max_p_idx = int(np.argmax(pvals)) max_p = pvals[max_p_idx] if max_p < 0.05: break removed_var_idx = current_vars[max_p_idx] step += 1 print(f" Step{step}: {VAR_NAMES[removed_var_idx]} を除外（p={max_p:.4f}）") current_vars.pop(max_p_idx) selected_names = [VAR_NAMES[i] for i in current_vars] final_model = sm.OLS(y_vals, sm.add_constant(X_raw[:, current_vars])).fit() print(f"\n 最終選択変数: {selected_names}") print(f" R² = {final_model.rsquared:.4f}") print(final_model.summary2()) # ================================================================ # ■ Step3. 外れ値処理後の再推定 # ================================================================ print("\n" + "=" * 60) print("■ Step3. 外れ値除外（IQR×1.5）後の再推定") print("=" * 60) X_clean = X_raw[mask_iqr.values][:, current_vars] y_clean = y_vals[mask_iqr.values] clean_model = sm.OLS(y_clean, sm.add_constant(X_clean)).fit() print(f" 除外後 N={mask_iqr.sum()}, R² = {clean_model.rsquared:.4f}") # ================================================================ # ■ 図の生成（4枚） # ================================================================ # ──────────────────────────────────────────────────────────────── # 図1 (fig1_corr): 相関係数棒グラフ（無相関検定） # ──────────────────────────────────────────────────────────────── fig1, ax1 = plt.subplots(figsize=(9, 5)) corrs = [stats.pearsonr(X_raw[:, i], y_vals)[0] for i in range(len(VAR_NAMES))] pvals2 = [stats.pearsonr(X_raw[:, i], y_vals)[1] for i in range(len(VAR_NAMES))] cols2 = ['#C62828' if p < 0.05 else '#9E9E9E' for p in pvals2] y_pos2 = np.arange(len(VAR_NAMES)) ax1.barh(y_pos2, corrs, color=cols2, alpha=0.75, edgecolor='white', height=0.6) ax1.axvline(0, color='gray', linestyle='--', linewidth=1.0) ax1.set_yticks(y_pos2) ax1.set_yticklabels(VAR_NAMES, fontsize=10) ax1.set_xlabel('ピアソン相関係数', fontsize=11) ax1.set_title('各変数と若年女性比率の相関係数\n（無相関検定: 赤=p<0.05）', fontsize=12, fontweight='bold') ax1.invert_yaxis() ax1.grid(axis='x', alpha=0.3) for i, (r, p) in enumerate(zip(corrs, pvals2)): sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else '' ax1.text(r + (0.005 if r >= 0 else -0.005), i, f' {r:.3f}{sig}', va='center', ha='left' if r >= 0 else 'right', fontsize=8.5) plt.tight_layout() fig1.savefig(os.path.join(FIG_DIR, '2025_H5_3_fig1_corr.png'), bbox_inches='tight', dpi=150) plt.close(fig1) print("\n図1保存: 2025_H5_3_fig1_corr.png") # ──────────────────────────────────────────────────────────────── # 図2 (fig2_backward): バックワード選択後の回帰係数 # ──────────────────────────────────────────────────────────────── fig2, ax2 = plt.subplots(figsize=(9, max(4, len(current_vars) * 0.8))) n_sel = len(current_vars) coefs3 = np.asarray(final_model.params[1:]) ses3 = np.asarray(final_model.bse[1:]) pvals3 = np.asarray(final_model.pvalues[1:]) cols3 = ['#C62828' if p < 0.01 else '#FF8F00' if p < 0.05 else '#9E9E9E' for p in pvals3] y_pos3 = np.arange(n_sel) ax2.barh(y_pos3, coefs3, color=cols3, alpha=0.75, edgecolor='white', height=0.6) ax2.errorbar(coefs3, y_pos3, xerr=1.96*ses3, fmt='none', color='#333', capsize=4, linewidth=1.2) ax2.axvline(0, color='gray', linestyle='--', linewidth=1.0) ax2.set_yticks(y_pos3) ax2.set_yticklabels(selected_names, fontsize=10) ax2.set_xlabel('回帰係数（±1.96SE）', fontsize=11) ax2.set_title(f'バックワード変数選択後の重回帰係数\n' f'R²={final_model.rsquared:.3f}（全p<0.05）, N={N}自治体', 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.01'), Patch(color='#FF8F00', alpha=0.75, label='p<0.05')], fontsize=9) plt.tight_layout() fig2.savefig(os.path.join(FIG_DIR, '2025_H5_3_fig2_backward.png'), bbox_inches='tight', dpi=150) plt.close(fig2) print("図2保存: 2025_H5_3_fig2_backward.png") # ──────────────────────────────────────────────────────────────── # 図3 (fig3_box): 若年女性比率の分布（箱ひげ図 + ヒストグラム） # ──────────────────────────────────────────────────────────────── fig3, axes3 = plt.subplots(1, 2, figsize=(12, 5)) fig3.suptitle(f'若年女性比率（15〜64歳女性/総女性人口）の分布\n' f'人口5,000〜20,000人 {N}自治体（SSDSE-A 実データ）', fontsize=12, fontweight='bold') ax3a = axes3[0] bp = ax3a.boxplot(df_main['若年女性比率'], vert=True, patch_artist=True, boxprops=dict(facecolor='#BBDEFB', color='#1565C0'), medianprops=dict(color='#C62828', linewidth=2), flierprops=dict(marker='o', markerfacecolor='#FF8F00', markersize=4)) ax3a.axhline(upper_iqr, color='#C62828', linestyle='--', linewidth=1.2, alpha=0.7, label=f'IQR上限={upper_iqr:.1f}%') ax3a.axhline(lower_iqr, color='#C62828', linestyle='--', linewidth=1.2, alpha=0.7, label=f'IQR下限={lower_iqr:.1f}%') ax3a.set_ylabel('若年女性比率（%）', fontsize=11) ax3a.set_title('箱ひげ図（IQR×1.5外れ値）', fontsize=11, fontweight='bold') ax3a.legend(fontsize=9) ax3a.grid(axis='y', alpha=0.3) if has_kawashima: kaw_y = y_vals[kaw_idx] ax3a.scatter([1], [kaw_y], color='#2E7D32', s=80, zorder=5) ax3a.text(1.15, kaw_y, '←川島町', va='center', fontsize=9, color='#2E7D32') ax3b = axes3[1] ax3b.hist(df_main['若年女性比率'], bins=35, color='#1565C0', alpha=0.7, edgecolor='white') mean_y = df_main['若年女性比率'].mean() ax3b.axvline(mean_y, color='#C62828', linestyle='--', linewidth=1.5, label=f'平均={mean_y:.2f}%') ax3b.set_xlabel('若年女性比率（%）', fontsize=11) ax3b.set_ylabel('自治体数', fontsize=11) ax3b.set_title('ヒストグラム', fontsize=11, fontweight='bold') ax3b.legend(fontsize=9) ax3b.grid(axis='y', alpha=0.3) plt.tight_layout() fig3.savefig(os.path.join(FIG_DIR, '2025_H5_3_fig3_box.png'), bbox_inches='tight', dpi=150) plt.close(fig3) print("図3保存: 2025_H5_3_fig3_box.png") # ──────────────────────────────────────────────────────────────── # 図4 (fig4_case): 事例分析 — 保育所数・製造業割合 vs 若年女性比率 # ──────────────────────────────────────────────────────────────── fig4, axes4 = plt.subplots(1, 2, figsize=(13, 5)) fig4.suptitle('保育・雇用環境と若年女性比率の関係（実データ）', fontsize=13, fontweight='bold') # (a) 保育所数 vs 若年女性比率 ax4a = axes4[0] x_nursery = df_main['保育所数'].values ax4a.scatter(x_nursery, y_vals, color='#1565C0', s=20, alpha=0.4, label='各自治体', zorder=2) if has_kawashima: ax4a.scatter([x_nursery[kaw_idx]], [y_vals[kaw_idx]], color='#C62828', s=150, zorder=5, label='川島町', marker='*') z4 = np.polyfit(x_nursery, y_vals, 1) xs4 = np.linspace(x_nursery.min(), x_nursery.max(), 100) ax4a.plot(xs4, np.poly1d(z4)(xs4), 'r-', linewidth=1.5, alpha=0.7) r4, p4 = stats.pearsonr(x_nursery, y_vals) ax4a.set_xlabel('保育所数（人口1万人あたり）', fontsize=11) ax4a.set_ylabel('若年女性比率（%）', fontsize=11) ax4a.set_title(f'保育所数との関係 r={r4:.3f} (p={p4:.4f})', fontsize=11, fontweight='bold') ax4a.legend(fontsize=9) ax4a.grid(True, alpha=0.3) # (b) 製造業従業者割合 vs 若年女性比率 ax4b = axes4[1] x_mfg = df_main['製造業従業者割合'].values ax4b.scatter(x_mfg, y_vals, color='#1565C0', s=20, alpha=0.4, label='各自治体', zorder=2) if has_kawashima: ax4b.scatter([x_mfg[kaw_idx]], [y_vals[kaw_idx]], color='#C62828', s=150, zorder=5, label='川島町', marker='*') z4b = np.polyfit(x_mfg, y_vals, 1) xs4b = np.linspace(x_mfg.min(), x_mfg.max(), 100) ax4b.plot(xs4b, np.poly1d(z4b)(xs4b), 'r-', linewidth=1.5, alpha=0.7) r4b, p4b = stats.pearsonr(x_mfg, y_vals) ax4b.set_xlabel('製造業従業者割合（%）', fontsize=11) ax4b.set_ylabel('若年女性比率（%）', fontsize=11) ax4b.set_title(f'製造業従業者割合との関係 r={r4b:.3f} (p={p4b:.4f})', fontsize=11, fontweight='bold') ax4b.legend(fontsize=9) ax4b.grid(True, alpha=0.3) plt.tight_layout() fig4.savefig(os.path.join(FIG_DIR, '2025_H5_3_fig4_case.png'), bbox_inches='tight', dpi=150) plt.close(fig4) print("図4保存: 2025_H5_3_fig4_case.png") print("\n全図の生成完了（4枚）") print(f" fig1_corr.png : 無相関検定（相関係数棒グラフ）") print(f" fig2_backward.png: バックワード変数選択後の回帰係数") print(f" fig3_box.png : 箱ひげ図・ヒストグラム") print(f" fig4_case.png : 事例分析（川島町・散布図）")