都道府県別交通事故・安全指標の規定要因
経済・都市化・医療アクセスの複合分析

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 matplotlib.patches import Patch

plt.rcParams['font.family'] = 'Hiragino Sans'
plt.rcParams['axes.unicode_minus'] = False
plt.rcParams['figure.dpi'] = 150

FIG_DIR = 'html/figures'
DATA_B  = 'data/raw/SSDSE-B-2026.csv'
DATA_E  = 'data/raw/SSDSE-E-2026.csv'
os.makedirs(FIG_DIR, exist_ok=True)

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_b['年度'] = df_b['年度'].astype(int)

print("SSDSE-B columns:", df_b.columns.tolist())
print("Years available:", sorted(df_b['年度'].unique()))

# 最新年（2023年）データ抽出
latest_year = df_b['年度'].max()
d = df_b[df_b['年度'] == latest_year].copy().reset_index(drop=True)
print(f"\n{latest_year}年 都道府県数: {len(d)}")

# SSDSE-E（横断面・都道府県別）
df_e_raw = pd.read_csv(DATA_E, encoding='cp932', header=0)
col_names_e = df_e_raw.iloc[1].tolist()
df_e = df_e_raw.iloc[2:].copy()
df_e.columns = col_names_e
df_e = df_e[df_e['地域コード'].str.match(r'R\d{5}', na=False)].copy()
df_e = df_e[df_e['地域コード'] != 'R00000'].copy().reset_index(drop=True)

def normalize_pref(name):
    name = str(name).strip()
    if name == '北海道':
        return '北海道'
    for s in ['都', '道', '府', '県']:
        if name.endswith(s):
            return name[:-1]
    return name

d['pref_key'] = d['都道府県'].apply(normalize_pref)
df_e['pref_key'] = df_e['都道府県'].apply(normalize_pref)

print("\nSSDSE-B 都道府県 sample:", d['pref_key'].tolist()[:5])
print("SSDSE-E 都道府県 sample:", df_e['pref_key'].tolist()[:5])

df_e['可住地面積_num'] = pd.to_numeric(df_e['可住地面積'], errors='coerce')
df_e['総人口_e_num']  = pd.to_numeric(df_e['総人口'], errors='coerce')

# 説明変数5: 人口密度（人/km²）　可住地面積は ha 単位 → / 100 で km²
df_e['人口密度'] = df_e['総人口_e_num'] / (df_e['可住地面積_num'] / 100)

region_map = {
    '北海道': '北海道・東北', '青森': '北海道・東北', '岩手': '北海道・東北', '宮城': '北海道・東北',
    '秋田': '北海道・東北', '山形': '北海道・東北', '福島': '北海道・東北',
    '茨城': '関東', '栃木': '関東', '群馬': '関東', '埼玉': '関東', '千葉': '関東', '東京': '関東', '神奈川': '関東',
    '新潟': '中部', '富山': '中部', '石川': '中部', '福井': '中部', '山梨': '中部',
    '長野': '中部', '岐阜': '中部', '静岡': '中部', '愛知': '中部',
    '三重': '近畿', '滋賀': '近畿', '京都': '近畿', '大阪': '近畿', '兵庫': '近畿', '奈良': '近畿', '和歌山': '近畿',
    '鳥取': '中国・四国', '島根': '中国・四国', '岡山': '中国・四国', '広島': '中国・四国',
    '山口': '中国・四国', '徳島': '中国・四国', '香川': '中国・四国', '愛媛': '中国・四国', '高知': '中国・四国',
    '福岡': '九州・沖縄', '佐賀': '九州・沖縄', '長崎': '九州・沖縄', '熊本': '九州・沖縄',
    '大分': '九州・沖縄', '宮崎': '九州・沖縄', '鹿児島': '九州・沖縄', '沖縄': '九州・沖縄'
}
region_colors = {
    '北海道・東北': '#4e9af1', '関東': '#e05c5c', '中部': '#f0a500',
    '近畿': '#5cb85c', '中国・四国': '#9b59b6', '九州・沖縄': '#f39c12'
}

df_merge['地域'] = df_merge['pref_key'].map(region_map)
df_merge['地域色'] = df_merge['地域'].map(region_colors)

vars_analysis = ['死亡率', '高齢化率', '有効求人倍率', '人口密度', '消費支出_万円', '人口10万対病院数']
var_labels = {
    '死亡率':         '死亡率\n(‰)',
    '高齢化率':       '高齢化率\n(%)',
    '有効求人倍率':    '有効求人\n倍率',
    '人口密度':       '人口密度\n(人/km²)',
    '消費支出_万円':  '消費支出\n(万円/月)',
    '人口10万対病院数': '病院数\n(10万対)',
}

print("\n=== 記述統計 ===")
print(df_merge[vars_analysis].describe().round(3).to_string())

# Pearson相関行列
print("\n=== Pearson相関行列 ===")
corr_mat = df_merge[vars_analysis].corr()
print(corr_mat.round(3).to_string())

df_sorted = df_merge.sort_values('死亡率', ascending=False).reset_index(drop=True)
colors_bar = [region_colors[r] for r in df_sorted['地域']]
mean_val = df_merge['死亡率'].mean()

fig1, ax1 = plt.subplots(figsize=(15, 7))
ax1.bar(range(len(df_sorted)), df_sorted['死亡率'],
        color=colors_bar, width=0.7, edgecolor='white', linewidth=0.4)
ax1.axhline(mean_val, color='#333', linewidth=2.0, linestyle='--', zorder=5)

ax1.set_xticks(range(len(df_sorted)))
ax1.set_xticklabels(df_sorted['pref_key'], rotation=90, fontsize=8.5)
ax1.set_ylabel('死亡率（人口千人あたり）', fontsize=11)
ax1.set_title(f'都道府県別 死亡率ランキング（{latest_year}年・47都道府県）\n'
              '高齢化・医療アクセス・経済水準の複合指標', fontsize=13, fontweight='bold', pad=12)

ax1.set_xlim(-0.7, len(df_sorted) - 0.3)
ax1.set_ylim(0, df_sorted['死亡率'].max() * 1.15)
ax1.spines['top'].set_visible(False)
ax1.spines['right'].set_visible(False)
ax1.yaxis.grid(True, alpha=0.3)
ax1.set_axisbelow(True)

# 全国平均ラベル
ax1.text(len(df_sorted) - 1, mean_val + 0.2, f'全国平均\n{mean_val:.2f}‰',
         ha='right', va='bottom', fontsize=9, color='#333')

# 地域凡例
legend_els = [Patch(facecolor=c, label=r) for r, c in region_colors.items()]
legend_els.append(plt.Line2D([0], [0], color='#333', linewidth=2.0,
                              linestyle='--', label=f'全国平均 {mean_val:.2f}‰'))
ax1.legend(handles=legend_els, fontsize=9, loc='upper right', framealpha=0.9)

plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2019_H5_1_fig1.png'), bbox_inches='tight')
plt.close(fig1)
print("\nFig1 saved.")

fig2, ax2 = plt.subplots(figsize=(11, 7))

for region, color in region_colors.items():
    mask = df_merge['地域'] == region
    ax2.scatter(df_merge.loc[mask, '有効求人倍率'],
                df_merge.loc[mask, '死亡率'],
                color=color, s=65, alpha=0.85,
                edgecolors='white', linewidths=0.5, label=region, zorder=3)

# 都道府県ラベル
for _, row in df_merge.iterrows():
    ax2.annotate(row['pref_key'],
                 (row['有効求人倍率'], row['死亡率']),
                 fontsize=6.5, ha='left', va='bottom',
                 xytext=(2, 2), textcoords='offset points', color='#333')

# 回帰直線
x_sc = df_merge['有効求人倍率'].values
y_sc = df_merge['死亡率'].values
r_val, p_val = stats.pearsonr(x_sc, y_sc)
slope, intercept, _, _, _ = stats.linregress(x_sc, y_sc)
x_line = np.linspace(x_sc.min() - 0.05, x_sc.max() + 0.05, 200)
ax2.plot(x_line, slope * x_line + intercept,
         color='#333', linewidth=1.8, linestyle='--', zorder=2)

p_str = f'p = {p_val:.4f}' if p_val >= 0.001 else 'p < 0.001'
ax2.text(0.05, 0.95, f'r = {r_val:.3f}\n{p_str}',
         transform=ax2.transAxes, fontsize=11, va='top', ha='left',
         bbox=dict(boxstyle='round,pad=0.4', facecolor='white',
                   edgecolor='#aaa', alpha=0.9))

ax2.set_xlabel('有効求人倍率（月間有効求人数 / 月間有効求職者数）', fontsize=11)
ax2.set_ylabel('死亡率（人口千人あたり）', fontsize=11)
ax2.set_title(f'有効求人倍率 vs 死亡率（{latest_year}年・47都道府県）\n'
              '経済的活力と健康・安全指標の地域差', fontsize=13, fontweight='bold', pad=12)

legend_els2 = [Patch(facecolor=c, label=r) for r, c in region_colors.items()]
ax2.legend(handles=legend_els2, fontsize=9, loc='lower right', framealpha=0.9)
ax2.spines['top'].set_visible(False)
ax2.spines['right'].set_visible(False)

plt.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2019_H5_1_fig2.png'), bbox_inches='tight')
plt.close(fig2)
print("Fig2 saved.")

from matplotlib.colors import TwoSlopeNorm

n_v = len(vars_analysis)
pval_mat = np.ones((n_v, n_v))
for i in range(n_v):
    for j in range(n_v):
        if i != j:
            _, p_ij = stats.pearsonr(df_merge[vars_analysis[i]],
                                     df_merge[vars_analysis[j]])
            pval_mat[i, j] = p_ij

corr_arr = corr_mat.values
norm = TwoSlopeNorm(vmin=-1, vcenter=0, vmax=1)

fig3, ax3 = plt.subplots(figsize=(9, 7))
im = ax3.imshow(corr_arr, cmap='RdBu_r', norm=norm, aspect='auto')

tick_labels = [var_labels[v] for v in vars_analysis]
ax3.set_xticks(range(n_v))
ax3.set_yticks(range(n_v))
ax3.set_xticklabels(tick_labels, fontsize=9)
ax3.set_yticklabels(tick_labels, fontsize=9)

for i in range(n_v):
    for j in range(n_v):
        r_ij  = corr_arr[i, j]
        p_ij  = pval_mat[i, j]
        sig_mark = ('***' if p_ij < 0.001 else
                    '**'  if p_ij < 0.01  else
                    '*'   if p_ij < 0.05  else '')
        text_color = 'white' if abs(r_ij) > 0.6 else 'black'
        ax3.text(j, i, f'{r_ij:.2f}{sig_mark}',
                 ha='center', va='center', fontsize=9,
                 color=text_color,
                 fontweight='bold' if sig_mark else 'normal')

plt.colorbar(im, ax=ax3, shrink=0.8, label='Pearson r')
ax3.set_title(f'変数間Pearson相関行列（{latest_year}年・47都道府県）\n'
              '*p<0.05  **p<0.01  ***p<0.001',
              fontsize=12, fontweight='bold', pad=14)
plt.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2019_H5_1_fig3.png'), bbox_inches='tight')
plt.close(fig3)
print("Fig3 saved.")

d['総人口_num']      = pd.to_numeric(d['総人口'], errors='coerce')
d['65歳以上_num']    = pd.to_numeric(d['65歳以上人口'], errors='coerce')
d['死亡数_num']      = pd.to_numeric(d['死亡数'], errors='coerce')
d['有効求人数_num']  = pd.to_numeric(d['月間有効求人数（一般）'], errors='coerce')
d['有効求職者数_num'] = pd.to_numeric(d['月間有効求職者数（一般）'], errors='coerce')
d['病院数_num']      = pd.to_numeric(d['一般病院数'], errors='coerce')
d['消費支出_num']    = pd.to_numeric(d['消費支出（二人以上の世帯）'], errors='coerce')

# 目的変数: 死亡率（人口千人あたり）
d['死亡率'] = d['死亡数_num'] / d['総人口_num'] * 1000

# 説明変数1: 高齢化率（%）
d['高齢化率'] = d['65歳以上_num'] / d['総人口_num'] * 100

# 説明変数2: 有効求人倍率
d['有効求人倍率'] = d['有効求人数_num'] / d['有効求職者数_num']

# 説明変数3: 人口10万対病院数（医療アクセス）
d['人口10万対病院数'] = d['病院数_num'] / d['総人口_num'] * 100000

# 説明変数4: 消費支出（万円換算）
d['消費支出_万円'] = d['消費支出_num'] / 10000

df_merge = d[['pref_key', '都道府県', '死亡率', '高齢化率', '有効求人倍率',
              '人口10万対病院数', '消費支出_万円']].merge(
    df_e[['pref_key', '人口密度']],
    on='pref_key', how='inner'
)

print(f"\nマージ後データ件数: {len(df_merge)}")
print(df_merge[['pref_key', '死亡率', '高齢化率', '有効求人倍率',
                '人口密度', '消費支出_万円', '人口10万対病院数']].head(10).to_string())

# 欠損値確認・除去
print("\n欠損値:")
print(df_merge.isnull().sum())
df_merge = df_merge.dropna().reset_index(drop=True)
print(f"欠損除去後: {len(df_merge)} 都道府県")

y = df_merge['死亡率'].values
X_vars = ['高齢化率', '有効求人倍率', '人口密度', '消費支出_万円', '人口10万対病院数']
X = df_merge[X_vars].values

# 標準化（標準化偏回帰係数のため）
y_std = (y - y.mean()) / y.std()
X_std = (X - X.mean(axis=0)) / X.std(axis=0)

# 標準化モデル
X_std_const = sm.add_constant(X_std)
model_std = sm.OLS(y_std, X_std_const).fit()

print("\n=== 重回帰分析（標準化） ===")
print(model_std.summary())

# 非標準化モデル（Cook's distance 等）
X_const = sm.add_constant(X)
model = sm.OLS(y, X_const).fit()

print("\n=== 重回帰分析（非標準化） ===")
print(model.summary())

# 標準化偏回帰係数の抽出
betas      = model_std.params[1:]
beta_se    = model_std.bse[1:]
beta_pvals = model_std.pvalues[1:]

print("\n=== 標準化偏回帰係数 ===")
for i, v in enumerate(X_vars):
    sig = ('***' if beta_pvals[i] < 0.001 else
           '**'  if beta_pvals[i] < 0.01  else
           '*'   if beta_pvals[i] < 0.05  else 'n.s.')
    print(f"  {v}: β={betas[i]:.4f}, SE={beta_se[i]:.4f}, p={beta_pvals[i]:.4f} {sig}")

print(f"\nR² = {model_std.rsquared:.4f}")
print(f"Adj.R² = {model_std.rsquared_adj:.4f}")

# Cook's distance
influence = model.get_influence()
cooks_d = influence.cooks_distance[0]
print("\n=== Cook's Distance（上位5） ===")
top5_cooks = pd.DataFrame({'都道府県': df_merge['pref_key'].values, 'Cooks_D': cooks_d})
print(top5_cooks.nlargest(5, 'Cooks_D').to_string(index=False))

# 上位・下位 5 都道府県（死亡率）
print("\n=== 死亡率 上位5都道府県 ===")
print(df_merge[['pref_key', '死亡率']].nlargest(5, '死亡率').to_string(index=False))
print("\n=== 死亡率 下位5都道府県 ===")
print(df_merge[['pref_key', '死亡率']].nsmallest(5, '死亡率').to_string(index=False))

x_labels = ['高齢化率', '有効求人倍率', '人口密度\n(都市化)', '消費支出\n(万円/月)', '病院数\n(10万対)']
ci_95 = 1.96 * beta_se

fig4, ax4 = plt.subplots(figsize=(9, 6))
y_pos = np.arange(len(X_vars))[::-1]

bar_colors = []
for i, pv in enumerate(beta_pvals):
    if pv < 0.05:
        bar_colors.append('#e05c5c' if betas[i] > 0 else '#4e9af1')
    else:
        bar_colors.append('#cccccc')

ax4.barh(y_pos, betas, xerr=ci_95,
         color=bar_colors, height=0.55,
         edgecolor='white', linewidth=0.5,
         error_kw=dict(ecolor='#555', capsize=4, lw=1.5))
ax4.axvline(0, color='#333', linewidth=1.2)

ax4.set_yticks(y_pos)
ax4.set_yticklabels(x_labels, fontsize=10)
ax4.set_xlabel('標準化偏回帰係数（β）', fontsize=11)
ax4.set_title('重回帰分析：標準化偏回帰係数と95%信頼区間\n'
              '目的変数：死亡率（安全・健康指標）',
              fontsize=13, fontweight='bold', pad=12)

# 有意水準マーク
for i, (b, pv, ci) in enumerate(zip(betas, beta_pvals, ci_95)):
    sig_mark = ('***' if pv < 0.001 else
                '**'  if pv < 0.01  else
                '*'   if pv < 0.05  else 'n.s.')
    offset = 0.025
    x_text = b + ci + offset if b >= 0 else b - ci - offset
    ha = 'left' if b >= 0 else 'right'
    ax4.text(x_text, y_pos[i], sig_mark,
             ha=ha, va='center', fontsize=9,
             color='#c00' if pv < 0.05 else '#888')

legend_els4 = [
    Patch(facecolor='#e05c5c', label='正の有意効果 (p<0.05)'),
    Patch(facecolor='#4e9af1', label='負の有意効果 (p<0.05)'),
    Patch(facecolor='#cccccc', label='非有意 (n.s.)'),
]
ax4.legend(handles=legend_els4, fontsize=9, loc='lower right')
ax4.text(0.02, 0.02,
         f'R²={model_std.rsquared:.3f}  Adj.R²={model_std.rsquared_adj:.3f}  N={len(df_merge)}',
         transform=ax4.transAxes, fontsize=9, color='#555')
ax4.spines['top'].set_visible(False)
ax4.spines['right'].set_visible(False)

plt.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2019_H5_1_fig4.png'), bbox_inches='tight')
plt.close(fig4)
print("Fig4 saved.")

print("\nDONE: 2019_H5_1_shorei")