生活保護・社会福祉支出と
地域経済格差の分析

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

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)

# 2022年断面
df22 = df_b[df_b['年度'] == 2022].copy()
assert len(df22) == 47, f"47都道府県が必要ですが {len(df22)} 行のみ取得"

# ── 派生変数の計算 ──────────────────────────────────────
# 高齢化率（%）
df22['高齢化率'] = df22['65歳以上人口'] / df22['総人口'] * 100
# 転出率（1000人あたり転出者数）
df22['転出率'] = df22['転出者数（日本人移動者）'] / df22['総人口'] * 1000
# 求人倍率代理（有効求人数 / 有効求職者数）
df22['求人倍率'] = df22['月間有効求人数（一般）'] / df22['月間有効求職者数（一般）']
# 保育所充足率（定員/在所児数）
df22['保育充足率'] = df22['保育所等定員数'] / df22['保育所等在所児数']
# 1人あたり保健医療費（目的変数として使用）
df22['保健医療費_1人'] = df22['保健医療費（二人以上の世帯）']
# 消費支出（生活水準の代理）
df22['消費支出'] = df22['消費支出（二人以上の世帯）']
# 待機児童率
df22['待機児童率'] = df22['保育所等利用待機児童数'] / df22['保育所等定員数'] * 100

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

df22['地域'] = df22['都道府県'].map(region_map)

# ── 分析変数の確認 ──────────────────────────────────────
print("=== 2022年 基本統計量 ===")
vars_show = ['高齢化率', '転出率', '求人倍率', '保健医療費_1人', '消費支出']
print(df22[vars_show].describe().round(2))

# ── 図1: 保健医療費/世帯 都道府県別ランキング（地域色分け棒グラフ） ─
print("\n=== 図1: 保健医療費ランキング ===")
df_rank = df22[['都道府県', '保健医療費_1人', '地域']].sort_values('保健医療費_1人', ascending=True)

fig, ax = plt.subplots(figsize=(10, 11))
colors = [region_colors[r] for r in df_rank['地域']]
bars = ax.barh(df_rank['都道府県'], df_rank['保健医療費_1人'], color=colors, edgecolor='white', linewidth=0.4)

# 値ラベル
for bar, val in zip(bars, df_rank['保健医療費_1人']):
    ax.text(val + 100, bar.get_y() + bar.get_height()/2,
            f'{int(val):,}', va='center', ha='left', fontsize=7.5)

ax.set_xlabel('保健医療費（円/世帯・月）', fontsize=11)
ax.set_title('都道府県別 保健医療費（2022年）\n地域別色分け', fontsize=13, fontweight='bold', pad=10)
ax.set_xlim(0, df_rank['保健医療費_1人'].max() * 1.20)
ax.tick_params(axis='y', labelsize=8.5)

# 凡例
from matplotlib.patches import Patch
legend_elements = [Patch(facecolor=c, label=r) for r, c in region_colors.items()]
ax.legend(handles=legend_elements, loc='lower right', fontsize=9, title='地域', title_fontsize=9)
ax.grid(axis='x', linestyle='--', alpha=0.4)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)

plt.tight_layout()
out1 = os.path.join(FIG_DIR, '2020_H4_fig1.png')
plt.savefig(out1, dpi=150, bbox_inches='tight')
plt.close()

print(f"  保存: {out1}")

# ── 図2: 高齢化率 vs 保健医療費 散布図（地域色分け・回帰線・ラベル） ──
print("\n=== 図2: 高齢化率 vs 保健医療費 ===")
r, p = stats.pearsonr(df22['高齢化率'], df22['保健医療費_1人'])
print(f"  r = {r:.3f}, p = {p:.4f}")

fig, ax = plt.subplots(figsize=(10, 7))
for region, grp in df22.groupby('地域'):
    ax.scatter(grp['高齢化率'], grp['保健医療費_1人'],
               color=region_colors[region], label=region,
               s=65, zorder=3, edgecolors='white', linewidth=0.6)

# ラベル
for _, row in df22.iterrows():
    ax.annotate(row['都道府県'], (row['高齢化率'], row['保健医療費_1人']),
                fontsize=6.5, xytext=(3, 2), textcoords='offset points', color='#333')

# 回帰直線
x_line = np.linspace(df22['高齢化率'].min(), df22['高齢化率'].max(), 100)
slope, intercept = np.polyfit(df22['高齢化率'], df22['保健医療費_1人'], 1)
ax.plot(x_line, slope * x_line + intercept, 'k--', linewidth=1.5,
        label=f'回帰線 (r={r:.3f}{"*" if p < 0.05 else ""})')

ax.set_xlabel('高齢化率（%）', fontsize=11)
ax.set_ylabel('保健医療費（円/世帯・月）', fontsize=11)
ax.set_title('高齢化率と保健医療費の関係（2022年, n=47）', fontsize=13, fontweight='bold')
ax.legend(fontsize=9, loc='upper left')
ax.grid(linestyle='--', alpha=0.4)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)

# 統計注記
p_str = f"p={p:.4f}" if p >= 0.0001 else "p<0.0001"
ax.text(0.98, 0.04, f'r = {r:.3f}, {p_str}',
        transform=ax.transAxes, ha='right', fontsize=10,
        bbox=dict(boxstyle='round,pad=0.4', facecolor='#fff9c4', edgecolor='#f9a825'))

plt.tight_layout()
out2 = os.path.join(FIG_DIR, '2020_H4_fig2.png')
plt.savefig(out2, dpi=150, bbox_inches='tight')
plt.close()

print(f"  保存: {out2}")

# ── 図3: 相関ヒートマップ ──────────────────────────────────
print("\n=== 図3: 相関ヒートマップ ===")
corr_vars = {
    '保健医療費':  '保健医療費_1人',
    '消費支出':    '消費支出',
    '高齢化率':    '高齢化率',
    '転出率':      '転出率',
    '求人倍率':    '求人倍率',
    '保育充足率':  '保育充足率',
    '待機児童率':  '待機児童率',
}
df_corr = df22[[v for v in corr_vars.values()]].rename(columns={v: k for k, v in corr_vars.items()})
corr_matrix = df_corr.corr()

# p値行列
n = len(df_corr)
p_matrix = pd.DataFrame(np.ones((len(corr_matrix), len(corr_matrix))),
                         index=corr_matrix.index, columns=corr_matrix.columns)
for c1 in corr_matrix.columns:
    for c2 in corr_matrix.columns:
        if c1 != c2:
            _, p_val = stats.pearsonr(df_corr[c1].dropna(), df_corr[c2].dropna())
            p_matrix.loc[c1, c2] = p_val

fig, ax = plt.subplots(figsize=(8, 7))
import matplotlib.colors as mcolors
cmap = plt.cm.RdBu_r
im = ax.imshow(corr_matrix.values, cmap=cmap, vmin=-1, vmax=1, aspect='auto')
plt.colorbar(im, ax=ax, shrink=0.82, label='Pearson r')

labels = list(corr_matrix.columns)
ax.set_xticks(range(len(labels)))
ax.set_yticks(range(len(labels)))
ax.set_xticklabels(labels, fontsize=10, rotation=30, ha='right')
ax.set_yticklabels(labels, fontsize=10)

for i in range(len(labels)):
    for j in range(len(labels)):
        val = corr_matrix.iloc[i, j]
        p_val = p_matrix.iloc[i, j]
        star = '**' if p_val < 0.01 else ('*' if p_val < 0.05 else '')
        txt_color = 'white' if abs(val) > 0.6 else 'black'
        ax.text(j, i, f'{val:.2f}{star}', ha='center', va='center',
                fontsize=8.5, color=txt_color, fontweight='bold' if star else 'normal')

ax.set_title('福祉・経済指標の相関ヒートマップ（2022年）\n* p<0.05  ** p<0.01',
             fontsize=12, fontweight='bold', pad=10)
plt.tight_layout()
out3 = os.path.join(FIG_DIR, '2020_H4_fig3.png')
plt.savefig(out3, dpi=150, bbox_inches='tight')
plt.close()

print(f"  保存: {out3}")
print(corr_matrix.round(3))

# ── 重回帰分析 ────────────────────────────────────────────
print("\n=== 重回帰分析 ===")
X_cols_raw = ['高齢化率', '転出率', '求人倍率', '消費支出']
X_cols_disp = ['高齢化率', '転出率', '求人倍率', '消費支出']
y_col = '保健医療費_1人'

df_reg = df22[X_cols_raw + [y_col]].dropna()
print(f"  回帰サンプル数: {len(df_reg)}")

# 標準化
df_std = df_reg.copy()
for col in X_cols_raw:
    df_std[col] = (df_reg[col] - df_reg[col].mean()) / df_reg[col].std()
df_std[y_col] = (df_reg[y_col] - df_reg[y_col].mean()) / df_reg[y_col].std()

X = sm.add_constant(df_std[X_cols_raw])
y = df_std[y_col]
result = sm.OLS(y, X).fit()
print(result.summary())

# 係数・信頼区間
coef_df = pd.DataFrame({
    '変数': X_cols_disp,
    '標準化係数': result.params[X_cols_raw].values,
    'p値': result.pvalues[X_cols_raw].values,
    'CI下限': result.conf_int().loc[X_cols_raw, 0].values,
    'CI上限': result.conf_int().loc[X_cols_raw, 1].values,
})
print("\n標準化偏回帰係数:")
print(coef_df.round(4))

# ── 図4: 標準化偏回帰係数プロット ──────────────────────────
print("\n=== 図4: 標準化偏回帰係数プロット ===")
coef_sorted = coef_df.sort_values('標準化係数')

fig, ax = plt.subplots(figsize=(9, 5))
y_pos = range(len(coef_sorted))
bar_colors = ['#e05c5c' if v > 0 else '#4e9af1' for v in coef_sorted['標準化係数']]
bars = ax.barh(y_pos, coef_sorted['標準化係数'], color=bar_colors,
               edgecolor='white', linewidth=0.5, height=0.55)

# 信頼区間（エラーバー）
for i, (_, row) in enumerate(coef_sorted.iterrows()):
    ax.plot([row['CI下限'], row['CI上限']], [i, i], 'k-', linewidth=2.2, zorder=4)
    ax.plot([row['CI下限'], row['CI上限']], [i, i], 'k|', markersize=10, zorder=4)

# 値ラベル
for bar, (_, row) in zip(bars, coef_sorted.iterrows()):
    x_pos = row['標準化係数']
    star = '**' if row['p値'] < 0.01 else ('*' if row['p値'] < 0.05 else ' (n.s.)')
    offset = 0.005 if x_pos >= 0 else -0.005
    ha = 'left' if x_pos >= 0 else 'right'
    ax.text(x_pos + offset, bar.get_y() + bar.get_height()/2,
            f'{x_pos:.3f}{star}', va='center', ha=ha, fontsize=9.5)

ax.axvline(0, color='black', linewidth=0.8)
ax.set_yticks(y_pos)
ax.set_yticklabels(coef_sorted['変数'], fontsize=11)
ax.set_xlabel('標準化偏回帰係数（β）', fontsize=11)
ax.set_title(f'保健医療費の決定要因（標準化偏回帰係数）\n'
             f'目的変数: 保健医療費/世帯   Adj. R² = {result.rsquared_adj:.3f}',
             fontsize=12, fontweight='bold', pad=10)
ax.set_xlim(coef_sorted['CI下限'].min() * 1.3,
            coef_sorted['CI上限'].max() * 1.3)
ax.grid(axis='x', linestyle='--', alpha=0.4)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)

note = '* p<0.05  ** p<0.01  エラーバー: 95%信頼区間'
ax.text(0.99, 0.02, note, transform=ax.transAxes, fontsize=9, ha='right', color='#555')

plt.tight_layout()
out4 = os.path.join(FIG_DIR, '2020_H4_fig4.png')
plt.savefig(out4, dpi=150, bbox_inches='tight')
plt.close()

print(f"  保存: {out4}")

# ── 地域別集計（HTML用） ─────────────────────────────────
print("\n=== 地域別集計 ===")
region_stats = df22.groupby('地域').agg(
    保健医療費平均=('保健医療費_1人', 'mean'),
    高齢化率平均=('高齢化率', 'mean'),
    消費支出平均=('消費支出', 'mean'),
    転出率平均=('転出率', 'mean'),
).round(1)
print(region_stats)

# ── 要約統計 ──────────────────────────────────────────────
print("\n=== 要約統計 ===")
print(f"保健医療費 全国平均: {df22['保健医療費_1人'].mean():.0f} 円/世帯")
print(f"保健医療費 最大: {df22['保健医療費_1人'].max():.0f} ({df22.loc[df22['保健医療費_1人'].idxmax(), '都道府県']})")
print(f"保健医療費 最小: {df22['保健医療費_1人'].min():.0f} ({df22.loc[df22['保健医療費_1人'].idxmin(), '都道府県']})")
print(f"高齢化率 全国平均: {df22['高齢化率'].mean():.1f}%")
print(f"高齢化率 最大: {df22['高齢化率'].max():.1f}% ({df22.loc[df22['高齢化率'].idxmax(), '都道府県']})")
print(f"R²: {result.rsquared:.3f}, Adj.R²: {result.rsquared_adj:.3f}")
print(f"F統計量: {result.fvalue:.3f}, p値: {result.f_pvalue:.4f}")
print(f"高齢化率 vs 保健医療費 r={r:.3f}, p={p:.4f}")

print("\n=== 全図生成完了 ===")
for f in [out1, out2, out3, out4]:
    exists = os.path.exists(f)
    print(f"  {'OK' if exists else 'NG'}: {os.path.basename(f)}")