都道府県別社会保障給付費と
生活保護受給率の規定要因分析

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(3) ===")
print(df_b.head(3))
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['月間有効求人数（一般）'] / df_cross['月間有効求職者数（一般）']
df_cross['病院数_per10万'] = df_cross['一般病院数'] / df_cross['総人口'] * 100000
df_cross['消費支出_万'] = df_cross['消費支出（二人以上の世帯）'] / 10000

# 目的変数：保健医療費（代理変数）
y_col = '保健医療費（二人以上の世帯）'
df_cross = df_cross.dropna(subset=[y_col, '高齢化率', '有効求人倍率', '病院数_per10万', '消費支出_万', 'region'])

print(f"分析サンプル数: {len(df_cross)} 都道府県（{latest_yr}年）")
print()
print("=== 基本統計量 ===")
stats_cols = [y_col, '高齢化率', '有効求人倍率', '病院数_per10万', '消費支出_万']
print(df_cross[stats_cols].describe().round(2))
print()

# ── 全国平均 ──────────────────────────────────────────────────────
national_avg = df_cross[y_col].mean()
print(f"保健医療費（全国平均）: {national_avg:.0f} 円/月")
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']]
bars = 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.legend(fontsize=11, loc='upper right')
ax1.grid(axis='y', alpha=0.3, linewidth=0.8)
ax1.spines['top'].set_visible(False)
ax1.spines['right'].set_visible(False)

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

plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2019_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=(9, 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 = np.polyfit(x2, y2, 1)
xr = np.linspace(x2.min() - 0.5, x2.max() + 0.5, 200)
ax2.plot(xr, m * xr + b, 'k-', linewidth=1.8, alpha=0.7, label=f'回帰直線')

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'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, '2019_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)
print(f"\nKruskal-Wallis: H={kw_stat:.3f}, p={kw_p:.4f}")

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)
kw_p_str = f'{kw_p:.4f}' if kw_p >= 0.0001 else f'{kw_p:.2e}'
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数表示
for i, (region, grp) in enumerate(zip(region_order, groups)):
    ax3.text(i + 1, ax3.get_ylim()[0] + 0.2,
             f'n={len(grp)}', ha='center', va='bottom', fontsize=9, color='#555')

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

X_vars = ['高齢化率', '有効求人倍率', '病院数_per10万', '消費支出_万']
var_labels = {
    '高齢化率': '高齢化率（%）',
    '有効求人倍率': '有効求人倍率',
    '病院数_per10万': '病院数（人口10万対）',
    '消費支出_万': '消費支出（万円/月）'
}

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

# 標準化係数・信頼区間・p値
coefs    = ols_result.params[1:]      # 定数項除く
ci       = ols_result.conf_int().iloc[1:]  # 95%CI
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=(9, 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.025 if c >= 0 else -0.025
    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=11)
ax4.set_xlabel('標準化偏回帰係数（β）', fontsize=12)
ax4.set_title(f'図4：重回帰分析 — 標準化偏回帰係数（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, '2019_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"高齢化率との相関: 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):
    print(f"  {var_labels[v]:20s}: β={c:+.4f}  {sig_star(p)}  (p={p:.4f})")

print("\nDONE: 2019_H4_katsuyo")