都道府県別スポーツ活動・体力と健康指標の関係
重回帰による規定要因分析

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

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

def get_short_name(full_name):
    for suffix in ['都', '道', '府', '県']:
        if full_name.endswith(suffix):
            return full_name[:-1]
    return full_name

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

# ── 分析用データセット作成（2023年断面） ─────────────────────────
latest_yr = 2023
df_cross = df_b[df_b['年度'] == latest_yr].copy()

df_cross['short_name'] = df_cross['都道府県'].apply(get_short_name)
df_cross['region']     = df_cross['short_name'].map(region_map_raw)

# プロキシ変数の計算
df_cross['高齢化率']      = df_cross['65歳以上人口'] / df_cross['総人口'] * 100
df_cross['消費支出_万']   = df_cross['消費支出（二人以上の世帯）'] / 10000
df_cross['教養娯楽費_千'] = df_cross['教養娯楽費（二人以上の世帯）'] / 1000
df_cross['光熱水道費_千'] = df_cross['光熱・水道費（二人以上の世帯）'] / 1000

y_col = '保健医療費（二人以上の世帯）'
analysis_cols = [y_col, '消費支出_万', '教養娯楽費_千', '高齢化率', '光熱水道費_千', 'region']
df_cross = df_cross.dropna(subset=analysis_cols).copy()

print(f"分析サンプル数: {len(df_cross)} 都道府県（{latest_yr}年）")
print()
print("=== 基本統計量 ===")
print(df_cross[[y_col, '消費支出_万', '教養娯楽費_千', '高齢化率', '光熱水道費_千']].describe().round(2))
print()

national_avg = df_cross[y_col].mean()
print(f"保健医療費（全国平均）: {national_avg:.0f} 円/月")
print()

# Top/Bottom 5 都道府県
df_sorted_rank = df_cross.sort_values(y_col, ascending=False)
print("=== 保健医療費 上位5都道府県 ===")
print(df_sorted_rank[['short_name', y_col]].head(5).to_string(index=False))
print("\n=== 保健医療費 下位5都道府県 ===")
print(df_sorted_rank[['short_name', y_col]].tail(5).to_string(index=False))
print()

df_sorted = df_cross.sort_values(y_col, ascending=False).copy()

fig1, ax1 = plt.subplots(figsize=(14, 7))
colors_bar = [region_colors[r] for r in df_sorted['region']]
ax1.bar(range(len(df_sorted)), df_sorted[y_col] / 1000,
        color=colors_bar, edgecolor='white', linewidth=0.5, alpha=0.88)

ax1.axhline(national_avg / 1000, color='#333333', linewidth=1.8,
            linestyle='--', label=f'全国平均: {national_avg/1000:.1f} 千円/月')

ax1.set_xticks(range(len(df_sorted)))
ax1.set_xticklabels(df_sorted['short_name'], rotation=55, ha='right', fontsize=8.5)
ax1.set_ylabel('保健医療費（千円/月）', fontsize=12)
ax1.set_title('図1：都道府県別 保健医療費（2023年）\n'
              '〈健康意識・健康状態の代理指標：二人以上の世帯の月次保健医療支出〉',
              fontsize=13, fontweight='bold', pad=14)
ax1.set_xlim(-0.8, len(df_sorted) - 0.2)
ax1.grid(axis='y', alpha=0.3, linewidth=0.8)
ax1.spines['top'].set_visible(False)
ax1.spines['right'].set_visible(False)

legend_elements = [Patch(facecolor=region_colors[r], label=r, alpha=0.88)
                   for r in region_order]
legend_elements.append(
    plt.Line2D([0], [0], color='#333333', linewidth=1.8, linestyle='--',
               label=f'全国平均: {national_avg/1000:.1f} 千円/月')
)
ax1.legend(handles=legend_elements, fontsize=9, loc='upper right', ncol=2)

plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2018_H4_fig1.png'), bbox_inches='tight')
plt.close(fig1)
print("fig1 saved.")

x2 = df_cross['消費支出_万'].values
y2 = df_cross[y_col].values / 1000  # 千円
r2, p2 = stats.pearsonr(x2, y2)

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

for _, row in df_cross.iterrows():
    c = region_colors[row['region']]
    ax2.scatter(row['消費支出_万'], row[y_col] / 1000,
                color=c, s=72, alpha=0.85, zorder=3,
                edgecolors='white', linewidths=0.8)
    ax2.annotate(row['short_name'],
                 (row['消費支出_万'], row[y_col] / 1000),
                 xytext=(3, 3), textcoords='offset points',
                 fontsize=6.5, color='#333', zorder=4)

m, b_int = np.polyfit(x2, y2, 1)
xr = np.linspace(x2.min() - 0.5, x2.max() + 0.5, 200)
ax2.plot(xr, m * xr + b_int, 'k-', linewidth=1.8, alpha=0.7)

p2_str = f'{p2:.4f}' if p2 >= 0.0001 else f'{p2:.2e}'
ax2.set_xlabel('消費支出（万円/月）', fontsize=12)
ax2.set_ylabel('保健医療費（千円/月）', fontsize=12)
ax2.set_title(f'図2：消費支出 vs 保健医療費（2023年、47都道府県）\n'
              f'Pearson r = {r2:.3f},  p = {p2_str}',
              fontsize=13, fontweight='bold', pad=12)

legend_elem2 = [Patch(facecolor=region_colors[r], label=r, alpha=0.85)
                for r in region_order]
ax2.legend(handles=legend_elem2, fontsize=9, loc='upper left',
           framealpha=0.85, ncol=2)
ax2.grid(alpha=0.25, linewidth=0.8)
ax2.spines['top'].set_visible(False)
ax2.spines['right'].set_visible(False)

plt.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2018_H4_fig2.png'), bbox_inches='tight')
plt.close(fig2)
print(f"fig2 saved.  r={r2:.3f}, p={p2_str}")

groups = [df_cross[df_cross['region'] == r][y_col].dropna().values
          for r in region_order]
groups_nonempty = [g for g in groups if len(g) > 0]
kw_stat, kw_p = kruskal(*groups_nonempty)
kw_p_str = f'{kw_p:.4f}' if kw_p >= 0.0001 else f'{kw_p:.2e}'
print(f"\nKruskal-Wallis: H={kw_stat:.3f}, p={kw_p_str}")

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

bp_data = [g / 1000 for g in groups]
bp = ax3.boxplot(bp_data, patch_artist=True, notch=False,
                 medianprops=dict(color='black', linewidth=2),
                 whiskerprops=dict(linewidth=1.4),
                 capprops=dict(linewidth=1.4),
                 flierprops=dict(marker='o', markersize=5, alpha=0.6))

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

ax3.set_xticks(range(1, len(region_order) + 1))
ax3.set_xticklabels(region_order, rotation=20, ha='right', fontsize=10)
ax3.set_ylabel('保健医療費（千円/月）', fontsize=12)
ax3.set_title(f'図3：地域ブロック別 保健医療費の分布（2023年）\n'
              f'Kruskal-Wallis検定: H = {kw_stat:.2f},  p = {kw_p_str}',
              fontsize=13, fontweight='bold', pad=12)
ax3.grid(axis='y', alpha=0.3, linewidth=0.8)
ax3.spines['top'].set_visible(False)
ax3.spines['right'].set_visible(False)

# n数表示（ylimを先に取得）
y3_lim = ax3.get_ylim()
for i, (region, grp) in enumerate(zip(region_order, groups)):
    ax3.text(i + 1, y3_lim[0] + (y3_lim[1] - y3_lim[0]) * 0.01,
             f'n={len(grp)}', ha='center', va='bottom', fontsize=9, color='#555')

plt.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2018_H4_fig3.png'), bbox_inches='tight')
plt.close(fig3)
print("fig3 saved.")

X_vars = ['消費支出_万', '教養娯楽費_千', '高齢化率', '光熱水道費_千']
var_labels = {
    '消費支出_万':   '消費支出（万円/月）\n〈所得水準の代理〉',
    '教養娯楽費_千': '教養娯楽費（千円/月）\n〈余暇・スポーツ支出の代理〉',
    '高齢化率':      '高齢化率（%）\n〈年齢構成・体力水準の代理〉',
    '光熱水道費_千': '光熱・水道費（千円/月）\n〈生活様式の指標〉',
}

df_reg = df_cross[X_vars + [y_col]].dropna().copy()

# 標準化
for col in X_vars + [y_col]:
    df_reg[col + '_z'] = (df_reg[col] - df_reg[col].mean()) / df_reg[col].std()

X_std = sm.add_constant(df_reg[[v + '_z' for v in X_vars]])
y_std = df_reg[y_col + '_z']
ols_result = sm.OLS(y_std, X_std).fit()

print("\n=== OLS重回帰結果（標準化変数） ===")
print(ols_result.summary())

coefs    = ols_result.params[1:]
ci       = ols_result.conf_int().iloc[1:]
pvals    = ols_result.pvalues[1:]
r2_val   = ols_result.rsquared
adj_r2   = ols_result.rsquared_adj

def sig_star(p):
    if p < 0.001:  return '***'
    elif p < 0.01: return '**'
    elif p < 0.05: return '*'
    else:          return 'n.s.'

labels_plot = [var_labels[v] for v in X_vars]
coef_vals   = coefs.values
ci_lo       = ci.iloc[:, 0].values
ci_hi       = ci.iloc[:, 1].values
xerr_lo     = coef_vals - ci_lo
xerr_hi     = ci_hi - coef_vals

colors4 = ['#e05c5c' if c > 0 else '#4e9af1' for c in coef_vals]

fig4, ax4 = plt.subplots(figsize=(10, 5.5))
y4_pos = np.arange(len(X_vars))

ax4.barh(y4_pos, coef_vals, xerr=[xerr_lo, xerr_hi],
         color=colors4, alpha=0.82, edgecolor='white',
         error_kw=dict(ecolor='#333', linewidth=1.5, capsize=5), height=0.55)
ax4.axvline(0, color='black', linewidth=1.2, linestyle='-')

for i, (c, p) in enumerate(zip(coef_vals, pvals.values)):
    star = sig_star(p)
    offset = 0.03 if c >= 0 else -0.03
    ha_a = 'left' if c >= 0 else 'right'
    color_s = '#c0392b' if star != 'n.s.' else '#888'
    ax4.text(c + offset, i, star, va='center', ha=ha_a,
             fontsize=12, fontweight='bold', color=color_s)

ax4.set_yticks(y4_pos)
ax4.set_yticklabels(labels_plot, fontsize=10)
ax4.set_xlabel('標準化偏回帰係数（β）', fontsize=12)
ax4.set_title(f'図4：重回帰分析 — 標準化偏回帰係数と95%信頼区間（2023年）\n'
              f'R² = {r2_val:.3f},  Adj.R² = {adj_r2:.3f}　'
              f'* p<0.05  ** p<0.01  *** p<0.001',
              fontsize=12, fontweight='bold', pad=12)
ax4.grid(axis='x', alpha=0.3, linewidth=0.8)
ax4.spines['top'].set_visible(False)
ax4.spines['right'].set_visible(False)

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

# ── 最終サマリー出力 ──────────────────────────────────────────────
print("\n=== 分析結果サマリー ===")
print(f"分析年度: {latest_yr}年  N={len(df_cross)}")
print(f"目的変数: 保健医療費（二人以上の世帯）[代理変数]")
print(f"消費支出との相関（Pearson r）: r={r2:.3f}, p={p2_str}")
print(f"Kruskal-Wallis（地域別）: H={kw_stat:.3f}, p={kw_p_str}")
print(f"重回帰 R²={r2_val:.3f}, Adj.R²={adj_r2:.3f}")
print("\n標準化偏回帰係数:")
for v, c, p in zip(X_vars, coef_vals, pvals.values):
    short_label = var_labels[v].split('\n')[0]
    print(f"  {short_label:25s}: β={c:+.4f}  {sig_star(p)}  (p={p:.4f})")

print("\nDONE: 2018_H4_katsuyo")