都道府県別死亡率の規定要因分析
高齢化・医療アクセス・経済格差の複合効果

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

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

# 図の保存先フォルダを用意（プロジェクトルートから実行する前提）
os.makedirs('html/figures', exist_ok=True)

# CSV読み込み → 47都道府県 → 最新年度 に絞り込み
df_b = pd.read_csv('data/raw/SSDSE-B-2026.csv', 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)
latest_year = df_b['年度'].max()
df_cross = df_b[df_b['年度'] == latest_year].copy().reset_index(drop=True)
print(f"最新年度: {latest_year}, n={len(df_cross)}")

# 目的変数：死亡率（人口千人当たり）
df_cross['死亡率'] = df_cross['死亡数'] / df_cross['総人口'] * 1000

# 高齢化率（65歳以上人口割合 %）
df_cross['高齢化率'] = df_cross['65歳以上人口'] / df_cross['総人口'] * 100

# 医療アクセス指標：人口10万人当たり一般病院数
df_cross['病院数_10万対'] = df_cross['一般病院数'] / df_cross['総人口'] * 100000

# 保健医療費（世帯）
df_cross['保健医療費'] = df_cross['保健医療費（二人以上の世帯）']

# 有効求人倍率（月間有効求人数 / 月間有効求職者数）
df_cross['有効求人倍率'] = df_cross['月間有効求人数（一般）'] / df_cross['月間有効求職者数（一般）']

# 消費支出
df_cross['消費支出'] = df_cross['消費支出（二人以上の世帯）']

# 都道府県名の「県・都・道・府」を省略（表示用）
def shorten_pref(name):
    for suffix in ['道', '都', '府', '県']:
        if name.endswith(suffix) and len(name) > len(suffix):
            return name[:-len(suffix)]
    return name

df_cross['pref_short'] = df_cross['都道府県'].apply(shorten_pref)

print("\n変数サマリー:")
print(df_cross[['都道府県', '死亡率', '高齢化率', '病院数_10万対',
                '有効求人倍率', '保健医療費']].describe())

region_colors = {
    '北海道・東北': '#4e9af1',
    '関東':       '#e05c5c',
    '中部':       '#f0a500',
    '近畿':       '#5cb85c',
    '中国・四国':  '#9b59b6',
    '九州・沖縄':  '#f39c12',
}

region_map_full = {
    '北海道': '北海道・東北', '青森県': '北海道・東北', '岩手県': '北海道・東北',
    '宮城県': '北海道・東北', '秋田県': '北海道・東北', '山形県': '北海道・東北', '福島県': '北海道・東北',
    '茨城県': '関東', '栃木県': '関東', '群馬県': '関東', '埼玉県': '関東',
    '千葉県': '関東', '東京都': '関東', '神奈川県': '関東',
    '新潟県': '中部', '富山県': '中部', '石川県': '中部', '福井県': '中部', '山梨県': '中部',
    '長野県': '中部', '岐阜県': '中部', '静岡県': '中部', '愛知県': '中部',
    '三重県': '近畿', '滋賀県': '近畿', '京都府': '近畿', '大阪府': '近畿',
    '兵庫県': '近畿', '奈良県': '近畿', '和歌山県': '近畿',
    '鳥取県': '中国・四国', '島根県': '中国・四国', '岡山県': '中国・四国',
    '広島県': '中国・四国', '山口県': '中国・四国', '徳島県': '中国・四国',
    '香川県': '中国・四国', '愛媛県': '中国・四国', '高知県': '中国・四国',
    '福岡県': '九州・沖縄', '佐賀県': '九州・沖縄', '長崎県': '九州・沖縄',
    '熊本県': '九州・沖縄', '大分県': '九州・沖縄', '宮崎県': '九州・沖縄',
    '鹿児島県': '九州・沖縄', '沖縄県': '九州・沖縄',
}
df_cross['地域'] = df_cross['都道府県'].map(region_map_full)
df_cross['color'] = df_cross['地域'].map(region_colors)

df_sorted = df_cross.sort_values('死亡率', ascending=False).reset_index(drop=True)
national_avg = df_cross['死亡率'].mean()

fig, ax = plt.subplots(figsize=(14, 7))
ax.bar(range(len(df_sorted)), df_sorted['死亡率'],
       color=df_sorted['color'], edgecolor='white', linewidth=0.5)

ax.axhline(national_avg, color='black', linestyle='--', linewidth=1.5,
           label=f'全国平均 {national_avg:.2f}‰')

ax.set_xticks(range(len(df_sorted)))
ax.set_xticklabels(df_sorted['pref_short'], rotation=90, fontsize=8)
ax.set_ylabel('死亡率（人口千人当たり）', fontsize=12)
ax.set_title(f'図1：都道府県別死亡率ランキング（{latest_year}年）',
             fontsize=13, fontweight='bold')
ax.set_ylim(0, df_sorted['死亡率'].max() * 1.15)

legend_handles = [plt.Rectangle((0, 0), 1, 1, color=c, label=r)
                  for r, c in region_colors.items()]
legend_handles.append(plt.Line2D([0], [0], color='black', linestyle='--',
                                 linewidth=1.5, label=f'全国平均 {national_avg:.2f}‰'))
ax.legend(handles=legend_handles, loc='upper right', fontsize=9, ncol=2)
ax.spines['top'].set_visible(False); ax.spines['right'].set_visible(False)
plt.tight_layout()
fig.savefig('html/figures/2018_H1_fig1.png', dpi=150, bbox_inches='tight')
plt.close()

def gini(arr):
    v = np.sort(arr)
    n = len(v)
    cumv = np.cumsum(v)
    return (n + 1 - 2 * np.sum(cumv) / cumv[-1]) / n

health_vals = df_cross['保健医療費'].dropna().values
g_health = gini(health_vals)

hosp_vals = df_cross['病院数_10万対'].dropna().values
g_hosp = gini(hosp_vals)

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# 左右2つのサブプロットの設定をリストにまとめる
plot_configs = [
    (axes[0], health_vals, '保健医療費（世帯月額・円）', g_health, '#1565C0'),
    (axes[1], hosp_vals,   '一般病院数（人口10万人当たり）', g_hosp,   '#00695C'),
]

for ax, vals, label, g, color in plot_configs:
    sorted_vals = np.sort(vals)
    cum_pop = np.linspace(0, 1, len(sorted_vals) + 1)
    cum_val = np.concatenate([[0], np.cumsum(sorted_vals) / sorted_vals.sum()])
    ax.fill_between(cum_pop, cum_pop, cum_val, alpha=0.3, color=color, label='格差面積')
    ax.plot(cum_pop, cum_val, color=color, linewidth=2.5, label=f'Lorenz曲線 (Gini={g:.3f})')
    ax.plot([0, 1], [0, 1], 'k--', linewidth=1.5, label='完全平等線')
    ax.set_xlabel('都道府県の累積人口割合'); ax.set_ylabel('累積シェア')
    ax.set_title(f'{label}\nGini係数 = {g:.3f}', fontweight='bold')
    ax.legend(fontsize=9); ax.set_xlim(0, 1); ax.set_ylim(0, 1); ax.grid(alpha=0.3)
    ax.spines['top'].set_visible(False); ax.spines['right'].set_visible(False)

plt.tight_layout()
fig.savefig('html/figures/2018_H1_fig2.png', dpi=150, bbox_inches='tight')
plt.close()
print(f"保健医療費Gini={g_health:.3f}, 病院数Gini={g_hosp:.3f}")

x_var = df_cross['有効求人倍率'].values
y_var = df_cross['死亡率'].values

r_val, p_val = stats.pearsonr(x_var, y_var)
slope, intercept, _, _, _ = stats.linregress(x_var, y_var)
x_line = np.linspace(x_var.min(), x_var.max(), 100)
y_line = slope * x_line + intercept

fig, ax = plt.subplots(figsize=(11, 8))
for _, row in df_cross.iterrows():
    ax.scatter(row['有効求人倍率'], row['死亡率'],
               color=row['color'], s=60, zorder=3, alpha=0.85,
               edgecolors='white', linewidth=0.5)
    ax.annotate(row['pref_short'],
                xy=(row['有効求人倍率'], row['死亡率']),
                xytext=(3, 3), textcoords='offset points',
                fontsize=7, color='#333333')

ax.plot(x_line, y_line, 'k-', linewidth=1.8, alpha=0.7,
        label=f'回帰直線 (r={r_val:.3f}, p={p_val:.4f})')
ax.set_xlabel('有効求人倍率', fontsize=12)
ax.set_ylabel('死亡率（人口千人当たり）', fontsize=12)
ax.set_title(f'図3：有効求人倍率と死亡率（{latest_year}年）  r={r_val:.3f}, p={p_val:.4f}',
             fontweight='bold')
ax.grid(alpha=0.25)
ax.spines['top'].set_visible(False); ax.spines['right'].set_visible(False)
plt.tight_layout()
fig.savefig('html/figures/2018_H1_fig3.png', dpi=150, bbox_inches='tight')
plt.close()
print(f"r={r_val:.3f}, p={p_val:.4f}")

from sklearn.preprocessing import StandardScaler

df_reg = df_cross[['死亡率', '高齢化率', '病院数_10万対',
                   '有効求人倍率', '保健医療費', '消費支出']].dropna()

X_cols = ['高齢化率', '病院数_10万対', '有効求人倍率', '保健医療費', '消費支出']
y_col  = '死亡率'

# 標準化（各列を平均0・分散1に揃える）
scaler   = StandardScaler()
X_scaled = scaler.fit_transform(df_reg[X_cols])
y_scaled = (df_reg[y_col].values - df_reg[y_col].mean()) / df_reg[y_col].std()

# 切片項を追加してOLS推定
X_sm  = sm.add_constant(X_scaled)
model = sm.OLS(y_scaled, X_sm).fit()
print(model.summary())

betas    = model.params[1:]      # 切片を除く
pvalues  = model.pvalues[1:]
conf_int = model.conf_int()[1:]
var_labels = ['高齢化率 (%)', '病院数\n(10万対)', '有効求人倍率',
              '保健医療費\n(円)', '消費支出\n(円)']

fig, ax = plt.subplots(figsize=(9, 5))
bar_colors = ['#1565C0' if b > 0 else '#C62828' for b in betas]
y_pos = np.arange(len(betas))
ax.barh(y_pos, betas, color=bar_colors, alpha=0.8, height=0.55)

# 95% 信頼区間のエラーバー
for i, (lo, hi) in enumerate(conf_int):
    ax.plot([lo, hi], [i, i], 'k-', linewidth=2, zorder=4)
    ax.plot([lo, lo], [i - 0.12, i + 0.12], 'k-', linewidth=2)
    ax.plot([hi, hi], [i - 0.12, i + 0.12], 'k-', linewidth=2)

# 有意性マーク
for i, p in enumerate(pvalues):
    if   p < 0.001: mark = '***'
    elif p < 0.01:  mark = '**'
    elif p < 0.05:  mark = '*'
    else:           mark = 'n.s.'
    x_off = betas[i] + (0.03 if betas[i] >= 0 else -0.03)
    ha    = 'left' if betas[i] >= 0 else 'right'
    ax.text(x_off, i, mark, va='center', ha=ha, fontsize=11,
            color='#C62828' if mark != 'n.s.' else '#888888', fontweight='bold')

ax.axvline(0, color='black', linewidth=1.0)
ax.set_yticks(y_pos); ax.set_yticklabels(var_labels, fontsize=11)
ax.set_xlabel('標準化偏回帰係数（β）', fontsize=12)
ax.set_title(f'図4：重回帰の標準化偏回帰係数（n=47, {latest_year}年）',
             fontweight='bold')
ax.grid(axis='x', alpha=0.3)
ax.spines['top'].set_visible(False); ax.spines['right'].set_visible(False)
ax.text(0.98, 0.04,
        f'R² = {model.rsquared:.3f}  adj.R² = {model.rsquared_adj:.3f}',
        transform=ax.transAxes, fontsize=10, ha='right',
        bbox=dict(boxstyle='round,pad=0.3', facecolor='#EFF3FF', alpha=0.8))

plt.tight_layout()
fig.savefig('html/figures/2018_H1_fig4.png', dpi=150, bbox_inches='tight')
plt.close()

print("\n=== Pearson相関係数（vs 死亡率）===")
for col in X_cols:
    r, p = stats.pearsonr(df_cross[col].dropna(),
                          df_cross['死亡率'][df_cross[col].notna()])
    if   p < 0.001: sig = '***'
    elif p < 0.01:  sig = '**'
    elif p < 0.05:  sig = '*'
    else:           sig = 'n.s.'
    print(f"  {col:20s}: r={r:+.4f}, p={p:.4f} {sig}")

print(f"\nGini係数:")
print(f"  保健医療費      : {g_health:.3f}")
print(f"  病院数(10万対)  : {g_hosp:.3f}")