建物火災発生率予測モデルの構築および
発生率の決定要因分析

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))

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")

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")

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")

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)

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}")

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}")

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)