健康寿命の地域差：
生活習慣・医療環境の重回帰分析

df_raw = pd.read_csv(DATA_B, encoding='cp932', header=1)
df_b = df_raw[df_raw['地域コード'].str.match(r'^R\d{5}$', na=False)].copy()

# 使用列の数値変換
num_cols = [
    '総人口', '65歳以上人口',
    '消費支出（二人以上の世帯）',
    '食料費（二人以上の世帯）',
    '保健医療費（二人以上の世帯）',
    '一般病院数',
]
for c in num_cols:
    df_b[c] = pd.to_numeric(df_b[c], errors='coerce')

df_b['年度'] = pd.to_numeric(df_b['年度'], errors='coerce')

df_b['保健医療費比率']   = df_b['保健医療費（二人以上の世帯）'] / df_b['消費支出（二人以上の世帯）'] * 100
df_b['食料費割合']       = df_b['食料費（二人以上の世帯）']   / df_b['消費支出（二人以上の世帯）'] * 100
df_b['高齢化率']         = df_b['65歳以上人口'] / df_b['総人口'] * 100
df_b['医療機関密度']     = df_b['一般病院数']   / df_b['総人口'] * 10000
df_b['消費支出_log']     = np.log(df_b['消費支出（二人以上の世帯）'])

# 地域ブロック付与
df_b['地域ブロック'] = df_b['都道府県'].map(REGION_MAP)

years = sorted(df_b['年度'].dropna().unique())

print("=" * 60)
print("■ データ概要")
print("=" * 60)
print(f"  総観測数  : {len(df_b)} 行 ({df_b['都道府県'].nunique()} 都道府県 × {len(years)} 年度)")
print(f"  期間      : {int(years[0])}〜{int(years[-1])} 年度")
print()

df22 = df_b[df_b['年度'] == 2022].dropna(
    subset=['保健医療費比率', '高齢化率', '食料費割合', '医療機関密度', '消費支出_log']
).copy()

print("■ 2022年の基本統計")
for col in ['保健医療費比率', '高齢化率', '食料費割合', '医療機関密度']:
    m = df22[col].mean()
    s = df22[col].std()
    print(f"  {col:<14}: mean={m:.2f}, std={s:.2f}")
print()

fig1, ax1 = plt.subplots(figsize=(12, 6))

for region in REGION_ORDER:
    df_reg = df_b[df_b['地域ブロック'] == region].copy()
    ts = df_reg.groupby('年度')['保健医療費比率'].mean().sort_index()
    ax1.plot(ts.index, ts.values,
             color=REGION_COLORS[region], linewidth=2.2,
             marker='o', markersize=4.5,
             label=region)

ax1.set_xlabel('年度', fontsize=12)
ax1.set_ylabel('保健医療費比率（保健医療費 / 消費支出 × 100）', fontsize=11)
ax1.set_title('保健医療費比率の時系列推移（地域ブロック別平均）\n'
              'SSDSE-B 2012〜2023年度, 二人以上世帯',
              fontsize=13, fontweight='bold')
ax1.set_xticks(years)
ax1.legend(fontsize=9.5, loc='upper left', ncol=2)
ax1.grid(True, alpha=0.3, axis='y')
ax1.yaxis.set_major_formatter(plt.FuncFormatter(lambda x, _: f'{x:.1f}%'))
plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2022_H5_2_fig1_ts_health_ratio.png'),
             bbox_inches='tight', dpi=150)
plt.close(fig1)
print("図1保存: 2022_H5_2_fig1_ts_health_ratio.png")

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

for region in REGION_ORDER:
    sub = df22[df22['地域ブロック'] == region]
    ax2.scatter(sub['高齢化率'], sub['保健医療費比率'],
                color=REGION_COLORS[region], s=65, alpha=0.82,
                label=region, zorder=3)

# 回帰直線（全体）
x_all = df22['高齢化率'].values
y_all = df22['保健医療費比率'].values
m_fit, b_fit, r_val, p_val, _ = stats.linregress(x_all, y_all)
x_line = np.linspace(x_all.min(), x_all.max(), 100)
ax2.plot(x_line, m_fit * x_line + b_fit,
         color='#333333', linewidth=1.8, linestyle='--',
         label=f'回帰直線 (r={r_val:+.3f}, p={p_val:.3f})', zorder=2)

# 代表的な都道府県にラベル
for _, row in df22.iterrows():
    pref = row['都道府県']
    if pref in ['秋田県', '東京都', '沖縄県', '高知県', '神奈川県', '愛知県']:
        ax2.annotate(pref,
                     xy=(row['高齢化率'], row['保健医療費比率']),
                     xytext=(5, 3), textcoords='offset points',
                     fontsize=8.5, color='#333')

ax2.set_xlabel('高齢化率（65歳以上人口 / 総人口 × 100）', fontsize=12)
ax2.set_ylabel('保健医療費比率（保健医療費 / 消費支出 × 100）', fontsize=11)
ax2.set_title('高齢化率 vs 保健医療費比率（2022年度, 47都道府県）',
              fontsize=13, fontweight='bold')
ax2.legend(fontsize=9, loc='upper left')
ax2.grid(True, alpha=0.3)
plt.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2022_H5_2_fig2_scatter_aging.png'),
             bbox_inches='tight', dpi=150)
plt.close(fig2)
print("図2保存: 2022_H5_2_fig2_scatter_aging.png")

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

# 標準化係数を計算（比較のため）
df22_std = df22[x_cols + [y_col]].copy()
for c in df22_std.columns:
    df22_std[c] = (df22_std[c] - df22_std[c].mean()) / df22_std[c].std()

X_std = sm.add_constant(df22_std[x_cols])
model_std = sm.OLS(df22_std[y_col], X_std).fit()

vars_plot = x_cols
coefs     = model_std.params[vars_plot].values
ci_low    = model_std.conf_int().loc[vars_plot, 0].values
ci_high   = model_std.conf_int().loc[vars_plot, 1].values
pvals     = model_std.pvalues[vars_plot].values

colors_bar = []
for p in pvals:
    if p < 0.001:
        colors_bar.append('#C62828')
    elif p < 0.05:
        colors_bar.append('#E65100')
    else:
        colors_bar.append('#9E9E9E')

y_pos = np.arange(len(vars_plot))
ax3.barh(y_pos, coefs, color=colors_bar, alpha=0.82, height=0.55, zorder=3)
ax3.errorbar(coefs, y_pos,
             xerr=[coefs - ci_low, ci_high - coefs],
             fmt='none', color='#333', capsize=5, linewidth=1.5, zorder=4)
ax3.axvline(0, color='black', linewidth=1.2, linestyle='-')
ax3.set_yticks(y_pos)
ax3.set_yticklabels(vars_plot, fontsize=11)
ax3.set_xlabel('標準化回帰係数（95%信頼区間）', fontsize=12)
ax3.set_title('OLS重回帰：標準化回帰係数プロット（2022年度）\n'
              '目的変数：保健医療費比率',
              fontsize=13, fontweight='bold')

# 凡例パッチ
import matplotlib.patches as mpatches
patches = [
    mpatches.Patch(color='#C62828', label='p < 0.001 ***'),
    mpatches.Patch(color='#E65100', label='p < 0.05  *'),
    mpatches.Patch(color='#9E9E9E', label='n.s.'),
]
ax3.legend(handles=patches, fontsize=9, loc='lower right')

# p値アノテーション
for i, (coef, p) in enumerate(zip(coefs, pvals)):
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else 'n.s.'
    ax3.text(ci_high[i] + 0.005, i, f' {sig}', va='center', fontsize=10)

ax3.grid(True, alpha=0.3, axis='x')
plt.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2022_H5_2_fig3_coef.png'),
             bbox_inches='tight', dpi=150)
plt.close(fig3)
print("図3保存: 2022_H5_2_fig3_coef.png")

import os
import warnings
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

warnings.filterwarnings('ignore')

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

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

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

REGION_ORDER  = ['北海道・東北', '関東', '中部', '近畿', '中国・四国', '九州・沖縄']
REGION_COLORS = {
    '北海道・東北': '#1565C0',
    '関東':         '#E65100',
    '中部':         '#2E7D32',
    '近畿':         '#6A1B9A',
    '中国・四国':   '#00838F',
    '九州・沖縄':   '#C62828',
}

y_col  = '保健医療費比率'
x_cols = ['高齢化率', '食料費割合', '医療機関密度', '消費支出_log']

print("■ 相関分析（2022年, n=47）")
corr_results = {}
for xc in x_cols:
    r, p = stats.pearsonr(df22[xc].dropna(), df22[y_col][df22[xc].notna()])
    corr_results[xc] = {'r': r, 'p': p}
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else 'n.s.'
    print(f"  {xc:<16}: r = {r:+.3f}  p = {p:.4f}  {sig}")
print()

X_ols = sm.add_constant(df22[x_cols])
model = sm.OLS(df22[y_col], X_ols).fit()

print("■ OLS 重回帰結果（2022年）")
print(f"  R² = {model.rsquared:.4f},  調整済みR² = {model.rsquared_adj:.4f},  n = {len(df22)}")
print(f"  F統計量 = {model.fvalue:.2f},  p(F) = {model.f_pvalue:.4f}")
print()
print(f"  {'変数':<16}  {'係数':>8}  {'標準誤差':>8}  {'t値':>7}  {'p値':>8}  有意性")
print("  " + "-" * 62)
for nm, coef, se, t, p in zip(
        model.params.index,
        model.params.values,
        model.bse.values,
        model.tvalues.values,
        model.pvalues.values):
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else 'n.s.'
    print(f"  {nm:<16}  {coef:>8.4f}  {se:>8.4f}  {t:>7.3f}  {p:>8.4f}  {sig}")
print()

fig4, ax4 = plt.subplots(figsize=(11, 6))

data_by_region = []
labels_region  = []
colors_box     = []

for region in REGION_ORDER:
    vals = df_b[df_b['地域ブロック'] == region]['医療機関密度'].dropna().values
    data_by_region.append(vals)
    labels_region.append(region)
    colors_box.append(REGION_COLORS[region])

bp = ax4.boxplot(data_by_region,
                 labels=labels_region,
                 patch_artist=True,
                 notch=False,
                 widths=0.5,
                 medianprops=dict(color='white', linewidth=2.5),
                 whiskerprops=dict(linewidth=1.5),
                 capprops=dict(linewidth=1.5),
                 flierprops=dict(marker='o', markersize=4, alpha=0.6))

for patch, color in zip(bp['boxes'], colors_box):
    patch.set_facecolor(color)
    patch.set_alpha(0.78)

for i, (vals, region) in enumerate(zip(data_by_region, labels_region)):
    jitter = np.random.default_rng(42).uniform(-0.18, 0.18, size=len(vals))
    ax4.scatter(np.full(len(vals), i + 1) + jitter, vals,
                color=REGION_COLORS[region], alpha=0.45, s=20, zorder=3)

ax4.set_xlabel('地域ブロック', fontsize=12)
ax4.set_ylabel('医療機関密度（一般病院数 / 総人口 × 10,000）', fontsize=11)
ax4.set_title('地域ブロック別 医療機関密度の分布\n（SSDSE-B 2012〜2023年度, 47都道府県×12年分）',
              fontsize=13, fontweight='bold')
ax4.grid(True, alpha=0.3, axis='y')
plt.xticks(rotation=15, ha='right')
plt.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2022_H5_2_fig4_boxplot_density.png'),
             bbox_inches='tight', dpi=150)
plt.close(fig4)
print("図4保存: 2022_H5_2_fig4_boxplot_density.png")

print()
print("=" * 60)
print("■ 全図の生成完了")
print(f"  保存先: {FIG_DIR}")
print("=" * 60)