""" 教育用再現コード: 2022年 統計データ分析コンペティション 審査員奨励賞 [大学生・一般の部] ================================================================= 論文タイトル：建物火災発生率予測モデルの構築および発生率の決定要因分析 受賞区分：審査員奨励賞（大学生・一般の部） 【分析概要】 データ：SSDSE-B-2026.csv（都道府県別・2022年度） N=47 目的変数：着工建築物密度（着工建築物数 / 総人口） ── 建物活動率（建物火災リスクの代理変数） 説明変数： 消費支出（二人以上の世帯） → 都市部生活水準 高齢化率（65歳以上人口/総人口）→ 高齢者比率 年平均気温 → 気候要因 人口密度（総人口 / 都道府県面積） → 都市化・密集度 一般病院密度（一般病院数/総人口） → 緊急対応密度の代理変数 Step1. 相関ヒートマップ（主要変数） Step2. OLS重回帰（標準化係数プロット） Step3. RandomForest 変数重要度 Step4. 着工建築物密度 vs 病院密度 散布図（都道府県ラベル） 【データ出典】 SSDSE-B-2026.csv: 社会・人口統計体系（都道府県データ） 統計センターより公表の実データ（2022年度 N=47都道府県） 【データサイエンス学習ポイント】 1. 代理変数（Proxy Variable）の設計 2. OLS回帰の標準化係数による効果量比較 3. RandomForest の変数重要度（MDI: Mean Decrease Impurity） 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/2022_U5_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 sklearn.ensemble import RandomForestRegressor from sklearn.preprocessing import StandardScaler 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_B = 'data/raw/SSDSE-B-2026.csv' os.makedirs(FIG_DIR, exist_ok=True) # ================================================================ # ■ 実データ読み込み（SSDSE-B-2026: 都道府県データ） # ================================================================ 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_2022 = df_b[df_b['年度'] == 2022].copy() # N=47 print("=" * 60) print(f"■ 都道府県データ 2022年度: N={len(df_2022)}") print("=" * 60) # ── 都道府県面積（km²）── 国土地理院 全国都道府県市区町村別面積調（令和5年） area_km2 = { '北海道': 83424, '青森県': 9646, '岩手県': 15275, '宮城県': 7282, '秋田県': 11638, '山形県': 9323, '福島県': 13784, '茨城県': 6097, '栃木県': 6408, '群馬県': 6362, '埼玉県': 3798, '千葉県': 5158, '東京都': 2194, '神奈川県': 2416, '新潟県': 12584, '富山県': 4248, '石川県': 4186, '福井県': 4191, '山梨県': 4465, '長野県': 13562, '岐阜県': 10621, '静岡県': 7777, '愛知県': 5173, '三重県': 5774, '滋賀県': 4017, '京都府': 4612, '大阪府': 1905, '兵庫県': 8401, '奈良県': 3691, '和歌山県': 4725, '鳥取県': 3507, '島根県': 6708, '岡山県': 7114, '広島県': 8480, '山口県': 6114, '徳島県': 4147, '香川県': 1877, '愛媛県': 5676, '高知県': 7104, '福岡県': 4987, '佐賀県': 2441, '長崎県': 4130, '熊本県': 7409, '大分県': 6341, '宮崎県': 7735, '鹿児島県': 9187,'沖縄県': 2281, } df_2022 = df_2022.copy() df_2022['面積_km2'] = df_2022['都道府県'].map(area_km2) # ── 特徴量エンジニアリング ───────────────────────────────────── pop = df_2022['総人口'].values.astype(float).clip(1) pop65 = df_2022['65歳以上人口'].values.astype(float) hosp = df_2022['一般病院数'].values.astype(float) kouji = df_2022['着工建築物数'].values.astype(float) temp = df_2022['年平均気温'].values.astype(float) spend = df_2022['消費支出（二人以上の世帯）'].values.astype(float) area = df_2022['面積_km2'].values.astype(float).clip(1) # 目的変数: 着工建築物密度（建物火災リスク代理変数） y_vals = kouji / pop * 10000 # 人口1万人あたり着工建築物数 # 説明変数 aging_rate = pop65 / pop * 100 # 高齢化率（%） hospital_den = hosp / pop * 10000 # 一般病院密度（人口1万人対） pop_density = pop / area # 人口密度（人/km²） temp_vals = temp # 年平均気温（℃） spend_vals = spend / 1000 # 消費支出（千円） VAR_NAMES = ['消費支出', '高齢化率', '年平均気温', '人口密度', '病院密度'] pref_names = df_2022['都道府県'].values X_raw = np.column_stack([spend_vals, aging_rate, temp_vals, pop_density, hospital_den]) # 欠損値チェック・除外 valid = np.all(np.isfinite(X_raw), axis=1) & np.isfinite(y_vals) X_raw = X_raw[valid] y_vals = y_vals[valid] pref_names = pref_names[valid] N = len(y_vals) print(f"有効都道府県数: N={N}") print(f"\n 目的変数（着工建築物密度）: mean={y_vals.mean():.2f}, std={y_vals.std():.2f}") print(f" 範囲: [{y_vals.min():.2f}, {y_vals.max():.2f}]") # ── 基本統計 ─────────────────────────────────────────────────── df_main = pd.DataFrame(X_raw, columns=VAR_NAMES) df_main['着工建築物密度'] = y_vals df_main['都道府県'] = pref_names print("\n 基本統計:") print(df_main[['着工建築物密度'] + VAR_NAMES].describe().round(2)) # ================================================================ # ■ Step1. OLS重回帰（標準化係数） # ================================================================ print("\n" + "=" * 60) print("■ Step1. OLS重回帰分析（標準化係数）") print("=" * 60) scaler = StandardScaler() X_scaled = scaler.fit_transform(X_raw) y_scaled = (y_vals - y_vals.mean()) / y_vals.std() model_ols = sm.OLS(y_scaled, sm.add_constant(X_scaled)).fit() print(model_ols.summary2()) std_coefs = model_ols.params[1:] # 定数項を除く std_ses = model_ols.bse[1:] std_pvals = model_ols.pvalues[1:] print(f"\n R² = {model_ols.rsquared:.4f}, adj.R² = {model_ols.rsquared_adj:.4f}") print(f"\n {'変数':<14} {'標準化係数':>10} {'SE':>8} {'p値':>10} {'有意':>6}") print(" " + "-" * 52) for name, b, se, p in zip(VAR_NAMES, std_coefs, std_ses, std_pvals): sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else 'n.s.' print(f" {name:<14} {b:>10.4f} {se:>8.4f} {p:>10.4f} {sig:>6}") # ================================================================ # ■ Step2. RandomForest変数重要度 # ================================================================ print("\n" + "=" * 60) print("■ Step2. RandomForest 変数重要度（MDI）") print("=" * 60) rf = RandomForestRegressor(n_estimators=500, max_features='sqrt', random_state=42, n_jobs=-1) rf.fit(X_raw, y_vals) importances = rf.feature_importances_ order_idx = np.argsort(importances)[::-1] print(f"\n OOB R²（参考）: RandomForestは交差検証が推奨") print(f"\n {'順位':>4} {'変数':<14} {'重要度':>10}") print(" " + "-" * 32) for rank, idx in enumerate(order_idx, 1): print(f" {rank:>4} {VAR_NAMES[idx]:<14} {importances[idx]:>10.4f}") # ================================================================ # ■ 図の生成（4枚） # ================================================================ # ──────────────────────────────────────────────────────────────── # 図1 (fig1_corr): 相関ヒートマップ（主要変数） # ──────────────────────────────────────────────────────────────── all_col_names = ['着工建築物密度'] + VAR_NAMES df_corr = df_main[all_col_names].copy() corr_mat = df_corr.corr() fig1, ax1 = plt.subplots(figsize=(8, 6.5)) im = ax1.imshow(corr_mat.values, cmap='RdBu_r', vmin=-1, vmax=1, aspect='auto') plt.colorbar(im, ax=ax1, shrink=0.8, label='ピアソン相関係数') ax1.set_xticks(range(len(all_col_names))) ax1.set_yticks(range(len(all_col_names))) ax1.set_xticklabels(all_col_names, rotation=35, ha='right', fontsize=10) ax1.set_yticklabels(all_col_names, fontsize=10) for i in range(len(all_col_names)): for j in range(len(all_col_names)): val = corr_mat.values[i, j] color = 'white' if abs(val) > 0.5 else 'black' ax1.text(j, i, f'{val:.2f}', ha='center', va='center', fontsize=9.5, color=color, fontweight='bold') ax1.set_title('相関ヒートマップ（主要変数・2022年度 都道府県別 N=47）', fontsize=12, fontweight='bold', pad=14) plt.tight_layout() fig1.savefig(os.path.join(FIG_DIR, '2022_U5_3_fig1_corr.png'), bbox_inches='tight', dpi=150) plt.close(fig1) print("\n図1保存: 2022_U5_3_fig1_corr.png") # ──────────────────────────────────────────────────────────────── # 図2 (fig2_importance): RF変数重要度棒グラフ # ──────────────────────────────────────────────────────────────── sorted_names = [VAR_NAMES[i] for i in order_idx] sorted_imp = importances[order_idx] colors_rf = ['#C62828' if imp >= sorted_imp[0] * 0.7 else '#FF8F00' if imp >= sorted_imp[0] * 0.4 else '#1565C0' for imp in sorted_imp] fig2, ax2 = plt.subplots(figsize=(8, 5)) bars = ax2.barh(range(len(sorted_names)), sorted_imp[::-1], color=colors_rf[::-1], alpha=0.80, edgecolor='white', height=0.6) ax2.set_yticks(range(len(sorted_names))) ax2.set_yticklabels(sorted_names[::-1], fontsize=11) ax2.set_xlabel('変数重要度（MDI: 不純度の平均減少量）', fontsize=11) ax2.set_title('RandomForest 変数重要度\n（建物火災リスク代理変数 = 着工建築物密度）', fontsize=12, fontweight='bold') ax2.axvline(0, color='gray', linewidth=0.8) ax2.grid(axis='x', alpha=0.3) for i, (imp, idx) in enumerate(zip(sorted_imp[::-1], order_idx[::-1])): ax2.text(imp + 0.003, i, f'{imp:.4f}', va='center', fontsize=9) plt.tight_layout() fig2.savefig(os.path.join(FIG_DIR, '2022_U5_3_fig2_importance.png'), bbox_inches='tight', dpi=150) plt.close(fig2) print("図2保存: 2022_U5_3_fig2_importance.png") # ──────────────────────────────────────────────────────────────── # 図3 (fig3_coef): OLS標準化係数プロット # ──────────────────────────────────────────────────────────────── sig_colors = [] for p in std_pvals: if p < 0.01: sig_colors.append('#C62828') elif p < 0.05: sig_colors.append('#FF8F00') else: sig_colors.append('#9E9E9E') fig3, ax3 = plt.subplots(figsize=(8, 5)) y_pos3 = np.arange(len(VAR_NAMES)) ax3.barh(y_pos3, std_coefs, color=sig_colors, alpha=0.78, edgecolor='white', height=0.6) ax3.errorbar(std_coefs, y_pos3, xerr=1.96 * std_ses, fmt='none', color='#333', capsize=4, linewidth=1.2) ax3.axvline(0, color='gray', linestyle='--', linewidth=1.0) ax3.set_yticks(y_pos3) ax3.set_yticklabels(VAR_NAMES, fontsize=11) ax3.set_xlabel('標準化回帰係数（±1.96 SE）', fontsize=11) ax3.set_title( f'OLS重回帰 標準化係数（N=47都道府県）\n' f'R²={model_ols.rsquared:.3f}, adj.R²={model_ols.rsquared_adj:.3f}', fontsize=12, fontweight='bold') ax3.invert_yaxis() ax3.grid(axis='x', alpha=0.3) from matplotlib.patches import Patch legend_items = [ Patch(color='#C62828', alpha=0.78, label='p<0.01'), Patch(color='#FF8F00', alpha=0.78, label='p<0.05'), Patch(color='#9E9E9E', alpha=0.78, label='n.s.'), ] ax3.legend(handles=legend_items, fontsize=9, loc='lower right') for i, (b, p) in enumerate(zip(std_coefs, std_pvals)): sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else '' offset = 0.02 if b >= 0 else -0.02 ha = 'left' if b >= 0 else 'right' ax3.text(b + offset, i, f'{b:.3f}{sig}', va='center', ha=ha, fontsize=8.5) plt.tight_layout() fig3.savefig(os.path.join(FIG_DIR, '2022_U5_3_fig3_coef.png'), bbox_inches='tight', dpi=150) plt.close(fig3) print("図3保存: 2022_U5_3_fig3_coef.png") # ──────────────────────────────────────────────────────────────── # 図4 (fig4_scatter): 着工建築物密度 vs 病院密度 散布図（都道府県ラベル） # ──────────────────────────────────────────────────────────────── x_hosp = df_main['病院密度'].values x_kouji = df_main['着工建築物密度'].values # 相関 r_sc, p_sc = stats.pearsonr(x_hosp, x_kouji) fig4, ax4 = plt.subplots(figsize=(11, 8)) sc = ax4.scatter(x_hosp, x_kouji, c=df_main['高齢化率'].values, cmap='RdYlBu_r', s=60, alpha=0.75, zorder=3, edgecolors='gray', linewidths=0.5) cbar4 = plt.colorbar(sc, ax=ax4, shrink=0.75, pad=0.02) cbar4.set_label('高齢化率（%）', fontsize=10) # 回帰直線 z4 = np.polyfit(x_hosp, x_kouji, 1) xs4 = np.linspace(x_hosp.min(), x_hosp.max(), 200) ax4.plot(xs4, np.poly1d(z4)(xs4), 'r--', linewidth=1.5, alpha=0.7, label=f'回帰直線 (r={r_sc:.3f}, p={p_sc:.3f})') # 都道府県ラベル for i, pref in enumerate(pref_names): label_short = pref.replace('都', '').replace('道', '').replace('府', '').replace('県', '') ax4.annotate(label_short, (x_hosp[i], x_kouji[i]), textcoords='offset points', xytext=(3, 3), fontsize=7.5, alpha=0.85) ax4.set_xlabel('一般病院密度（人口1万人あたり）', fontsize=12) ax4.set_ylabel('着工建築物密度（人口1万人あたり）', fontsize=12) ax4.set_title('着工建築物密度 vs 一般病院密度\n（都道府県別 2022年度, 色=高齢化率）', fontsize=13, fontweight='bold') ax4.legend(fontsize=10) ax4.grid(True, alpha=0.25) plt.tight_layout() fig4.savefig(os.path.join(FIG_DIR, '2022_U5_3_fig4_scatter.png'), bbox_inches='tight', dpi=150) plt.close(fig4) print("図4保存: 2022_U5_3_fig4_scatter.png") print("\n" + "=" * 60) print("全図の生成完了（4枚）") print(f" 2022_U5_3_fig1_corr.png : 相関ヒートマップ") print(f" 2022_U5_3_fig2_importance.png: RF変数重要度棒グラフ") print(f" 2022_U5_3_fig3_coef.png : OLS標準化係数プロット") print(f" 2022_U5_3_fig4_scatter.png : 着工建築物密度 vs 病院密度 散布図") print("=" * 60)