都道府県別健康寿命の決定要因
社会経済・医療アクセス要因の重回帰分析

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 scipy.stats import kruskal

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'
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("=== df_b.columns ===")
print(df_b.columns.tolist())
print()
print("=== df_b.head() ===")
print(df_b.head())

YEAR = 2019
df_cross = df_b[df_b['年度'] == YEAR].copy().reset_index(drop=True)
print(f"\n=== {YEAR}年度 都道府県数: {len(df_cross)} ===")

df_cross['死亡率'] = df_cross['死亡数'] / df_cross['総人口'] * 1000

# 説明変数
# (1) 高齢化率（%）― 高齢者が多いほど死亡率は上がる可能性
df_cross['高齢化率'] = df_cross['65歳以上人口'] / df_cross['総人口'] * 100

# (2) 消費支出（万円/世帯月額）― 社会経済的地位の代理変数
df_cross['消費支出_万円'] = df_cross['消費支出（二人以上の世帯）'] / 10000

# (3) 保健医療費（万円/世帯月額）― 医療へのアクセス・支出
df_cross['保健医療費_万円'] = df_cross['保健医療費（二人以上の世帯）'] / 10000

# (4) 病院数（千人当たり）― 医療インフラアクセス
df_cross['病院数_千人'] = df_cross['一般病院数'] / df_cross['総人口'] * 1000

# (5) 食料費（万円/世帯月額）― 食生活・栄養の代理変数
df_cross['食料費_万円'] = df_cross['食料費（二人以上の世帯）'] / 10000

print("\n=== 主要変数の記述統計 ===")
key_vars = ['死亡率', '高齢化率', '消費支出_万円', '保健医療費_万円', '病院数_千人', '食料費_万円']
print(df_cross[key_vars].describe().round(3))

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

region_order = ['北海道・東北', '関東', '中部', '近畿', '中国・四国', '九州・沖縄']
region_colors = {
    '北海道・東北': '#2196F3',
    '関東':       '#4CAF50',
    '中部':       '#FF9800',
    '近畿':       '#9C27B0',
    '中国・四国':  '#F44336',
    '九州・沖縄':  '#009688'
}

heat_vars = ['死亡率', '高齢化率', '消費支出_万円', '保健医療費_万円', '病院数_千人', '食料費_万円']
heat_labels = ['死亡率\n(千人当)', '高齢化率\n(%)', '消費支出\n(万円)', '保健医療費\n(万円)', '病院数\n(千人当)', '食料費\n(万円)']
n_vars_heat = len(heat_vars)

corr_matrix = df_cross[heat_vars].corr()
pval_matrix = pd.DataFrame(np.ones((n_vars_heat, n_vars_heat)),
                            index=heat_vars, columns=heat_vars)
for i, v1 in enumerate(heat_vars):
    for j, v2 in enumerate(heat_vars):
        if i != j:
            r, p = stats.pearsonr(df_cross[v1].dropna(), df_cross[v2].dropna())
            pval_matrix.loc[v1, v2] = p

fig1, ax1 = plt.subplots(figsize=(8, 6.5))
im = ax1.imshow(corr_matrix.values, cmap='RdBu_r', vmin=-1, vmax=1, aspect='auto')
plt.colorbar(im, ax=ax1, shrink=0.8, label='Pearson相関係数 r')

ax1.set_xticks(range(n_vars_heat))
ax1.set_yticks(range(n_vars_heat))
ax1.set_xticklabels(heat_labels, fontsize=10)
ax1.set_yticklabels(heat_labels, fontsize=10)

for i in range(n_vars_heat):
    for j in range(n_vars_heat):
        r_val = corr_matrix.values[i, j]
        p_val = pval_matrix.values[i, j]
        stars = '***' if p_val < 0.001 else '**' if p_val < 0.01 else '*' if p_val < 0.05 else ''
        txt = f"{r_val:.2f}{stars}"
        color = 'white' if abs(r_val) > 0.5 else 'black'
        ax1.text(j, i, txt, ha='center', va='center', fontsize=10,
                 color=color, fontweight='bold' if stars else 'normal')

ax1.set_title('図1：主要変数間のPearson相関係数ヒートマップ\n（2019年度 都道府県別，N=47）\n*p<0.05, **p<0.01, ***p<0.001',
              fontsize=11, pad=12)
plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2019_U5_2_fig1.png'), bbox_inches='tight')
plt.close(fig1)
print("\n図1 保存完了")

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

for region in region_order:
    mask = df_cross['地域'] == region
    ax2.scatter(df_cross.loc[mask, '消費支出_万円'],
                df_cross.loc[mask, '死亡率'],
                c=region_colors[region], label=region, s=60, alpha=0.85, zorder=3)

# 都道府県ラベル
for _, row in df_cross.iterrows():
    ax2.annotate(row['都道府県'].replace('県', '').replace('都', '').replace('府', '').replace('道', ''),
                 (row['消費支出_万円'], row['死亡率']),
                 fontsize=6.5, ha='left', va='bottom',
                 xytext=(2, 2), textcoords='offset points', color='#333333')

# 回帰直線
x_plot = df_cross['消費支出_万円'].dropna()
y_plot = df_cross['死亡率'].dropna()
idx_common = df_cross[['消費支出_万円', '死亡率']].dropna().index
x_c = df_cross.loc[idx_common, '消費支出_万円']
y_c = df_cross.loc[idx_common, '死亡率']
r_val, p_val = stats.pearsonr(x_c, y_c)
slope, intercept, _, _, _ = stats.linregress(x_c, y_c)
x_line = np.linspace(x_c.min(), x_c.max(), 100)
ax2.plot(x_line, slope * x_line + intercept, 'k--', linewidth=1.5, alpha=0.7)

stars = '***' if p_val < 0.001 else '**' if p_val < 0.01 else '*' if p_val < 0.05 else ''
ax2.text(0.05, 0.95, f'r = {r_val:.3f}{stars}\np = {p_val:.4f}\nN = 47',
         transform=ax2.transAxes, fontsize=11,
         verticalalignment='top',
         bbox=dict(boxstyle='round', facecolor='white', alpha=0.8))

ax2.set_xlabel('消費支出（万円/月，二人以上の世帯）', fontsize=12)
ax2.set_ylabel('死亡率（千人当たり）', fontsize=12)
ax2.set_title('図2：消費支出と死亡率の散布図\n（2019年度 都道府県別，地域別色分け）', fontsize=12)
ax2.legend(loc='upper right', fontsize=9, title='地域', title_fontsize=9)
ax2.grid(True, alpha=0.3)
fig2.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2019_U5_2_fig2.png'), bbox_inches='tight')
plt.close(fig2)
print("図2 保存完了")

fig3, ax3 = plt.subplots(figsize=(9, 6))

groups = [df_cross.loc[df_cross['地域'] == reg, '死亡率'].dropna().values
          for reg in region_order]
stat_kw, p_kw = kruskal(*groups)

bp = ax3.boxplot(groups, patch_artist=True, notch=False,
                 medianprops=dict(color='black', linewidth=2),
                 whiskerprops=dict(linewidth=1.5),
                 capprops=dict(linewidth=1.5))

for patch, region in zip(bp['boxes'], region_order):
    patch.set_facecolor(region_colors[region])
    patch.set_alpha(0.7)

# 個別データ点をオーバーレイ
for i, (region, grp) in enumerate(zip(region_order, groups), start=1):
    jitter = np.random.default_rng(seed=42 + i).uniform(-0.15, 0.15, len(grp))
    ax3.scatter([i + j for j in jitter], grp,
                color=region_colors[region], s=30, alpha=0.8, zorder=3)

ax3.set_xticks(range(1, len(region_order) + 1))
ax3.set_xticklabels(region_order, fontsize=9, rotation=15)
ax3.set_ylabel('死亡率（千人当たり）', fontsize=12)
ax3.set_title('図3：地域別死亡率の比較（箱ひげ図）\n（2019年度，N=47）', fontsize=12)
ax3.text(0.98, 0.97,
         f'Kruskal-Wallis: H={stat_kw:.2f}, p={p_kw:.4f}',
         transform=ax3.transAxes, ha='right', va='top', fontsize=10,
         bbox=dict(boxstyle='round', facecolor='lightyellow', alpha=0.9))
ax3.grid(True, axis='y', alpha=0.3)
fig3.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2019_U5_2_fig3.png'), bbox_inches='tight')
plt.close(fig3)
print("図3 保存完了")