食料費支出と健康
都道府県別食生活指標の回帰分析

df_raw = pd.read_csv(DATA_B, encoding='cp932', header=1)

# 都道府県のみ（地域コードが R + 5桁数字）
pref_mask = df_raw['地域コード'].str.match(r'^R\d{5}$', na=False)
df_all = df_raw[pref_mask].copy().reset_index(drop=True)

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

# ── 派生変数を計算 ─────────────────────────────────────────────
df_all['食料費割合']    = df_all['食料費（二人以上の世帯）'] / df_all['消費支出（二人以上の世帯）'] * 100
df_all['保健医療費_千円'] = df_all['保健医療費（二人以上の世帯）'] / 1000
df_all['高齢化率']      = df_all['65歳以上人口'] / df_all['総人口'] * 100
df_all['消費支出_万円']  = df_all['消費支出（二人以上の世帯）'] / 10000
df_all['医療機関密度']   = df_all['一般診療所数'] / df_all['総人口'] * 10000

print("=" * 65)
print(f"■ 全データ: {len(df_all)}行（47都道府県 × 12年間）")
print("=" * 65)

# 2022年断面
df2022 = df_all[df_all['年度'] == 2022].copy().reset_index(drop=True)
print(f"  2022年データ: {len(df2022)}都道府県")

# 欠損除外（2022年）
key_vars = ['食料費割合', '保健医療費_千円', '高齢化率', '消費支出_万円', '医療機関密度']
df2022_clean = df2022.dropna(subset=key_vars).reset_index(drop=True)
N = len(df2022_clean)
print(f"  欠損除外後: N={N}都道府県")
print()
print(df2022_clean[['都道府県'] + key_vars].describe().round(2))

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

df_all['地域ブロック']   = df_all['都道府県'].map(region_map)
df2022_clean['地域ブロック'] = df2022_clean['都道府県'].map(region_map)

# ── 時系列用データ（全年度）────────────────────────────────────
df_ts = df_all.dropna(subset=['食料費割合', '地域ブロック']).copy()

# 相関確認（2022年）
print("\n■ 2022年 主要変数間の相関（vs 保健医療費_千円）")
for v in ['食料費割合', '高齢化率', '消費支出_万円', '医療機関密度']:
    r, p = stats.pearsonr(df2022_clean[v], df2022_clean['保健医療費_千円'])
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else 'n.s.'
    print(f"  {v:<18} r={r:>7.4f}  p={p:.4f}  {sig}")

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

for region in region_order:
    df_r = df_ts[df_ts['地域ブロック'] == region].groupby('年度')['食料費割合'].mean().sort_index()
    if len(df_r) == 0:
        continue
    ax1.plot(df_r.index, df_r.values,
             color=region_colors[region], linewidth=2.0,
             marker='o', markersize=5, label=region)

ax1.set_xlabel('年度', fontsize=12)
ax1.set_ylabel('食料費割合（%）\n[食料費 / 消費支出 × 100]', fontsize=11)
ax1.set_title('食料費割合の地域別推移（2012〜2023年）\n'
              '6地域ブロック平均値（SSDSE-B実データ）',
              fontsize=13, fontweight='bold')
ax1.legend(fontsize=10, loc='upper right')
ax1.grid(True, alpha=0.3)
ax1.set_xticks(range(2012, 2024, 2))

# 全国平均を太線で重ね書き
df_nat = df_ts.groupby('年度')['食料費割合'].mean().sort_index()
ax1.plot(df_nat.index, df_nat.values,
         color='black', linewidth=2.5, linestyle='--',
         marker='D', markersize=5, label='全国平均', zorder=10)
ax1.legend(fontsize=10, loc='upper right')

plt.tight_layout()
fig1_path = os.path.join(FIG_DIR, '2022_H5_4_fig1_timeseries.png')
fig1.savefig(fig1_path, bbox_inches='tight', dpi=150)
plt.close(fig1)
print(f"\n図1保存: {os.path.basename(fig1_path)}")

fig2, ax2 = plt.subplots(figsize=(12, 8))

for region in region_order:
    mask = df2022_clean['地域ブロック'] == region
    ax2.scatter(df2022_clean.loc[mask, '食料費割合'],
                df2022_clean.loc[mask, '保健医療費_千円'],
                color=region_colors[region], s=60, alpha=0.8,
                label=region, zorder=3)

# 都道府県ラベル
for _, row in df2022_clean.iterrows():
    ax2.annotate(row['都道府県'],
                 (row['食料費割合'], row['保健医療費_千円']),
                 fontsize=7.5, alpha=0.85,
                 xytext=(3, 3), textcoords='offset points')

# 回帰直線（単相関）
x_vals = df2022_clean['食料費割合'].values
y_vals_s = df2022_clean['保健医療費_千円'].values
z = np.polyfit(x_vals, y_vals_s, 1)
xs = np.linspace(x_vals.min() - 0.5, x_vals.max() + 0.5, 200)
ax2.plot(xs, np.poly1d(z)(xs), color='gray', linewidth=1.8,
         linestyle='--', alpha=0.7, zorder=2)

r2, p2 = stats.pearsonr(x_vals, y_vals_s)
ax2.set_xlabel('食料費割合（%）[食料費 / 消費支出 × 100]', fontsize=12)
ax2.set_ylabel('保健医療費（千円/月）', fontsize=12)
ax2.set_title(f'食料費割合 vs 保健医療費（2022年・47都道府県）\n'
              f'Pearson r = {r2:.3f}（p = {p2:.4f}）',
              fontsize=13, fontweight='bold')
ax2.legend(fontsize=10, loc='upper right')
ax2.grid(True, alpha=0.3)

plt.tight_layout()
fig2_path = os.path.join(FIG_DIR, '2022_H5_4_fig2_scatter.png')
fig2.savefig(fig2_path, bbox_inches='tight', dpi=150)
plt.close(fig2)
print(f"図2保存: {os.path.basename(fig2_path)}")

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

n_vars = len(ols_varnames)
y_pos  = np.arange(n_vars)
bar_colors = ['#C62828' if p < 0.01 else '#FF8F00' if p < 0.05 else '#9E9E9E'
              for p in pvals]

ax3.barh(y_pos, coefs, color=bar_colors, alpha=0.78, edgecolor='white', height=0.55)
ax3.errorbar(coefs, y_pos, xerr=1.96 * ses,
             fmt='none', color='#333333', capsize=5, linewidth=1.4)
ax3.axvline(0, color='gray', linestyle='--', linewidth=1.0)
ax3.set_yticks(y_pos)
ax3.set_yticklabels(ols_varnames, fontsize=11)
ax3.set_xlabel('回帰係数（±1.96×SE）', fontsize=12)
ax3.set_title(f'保健医療費（千円）の決定要因 — OLS回帰係数\n'
              f'R² = {ols_model.rsquared:.3f}、N = {N}都道府県（2022年）',
              fontsize=12, fontweight='bold')
ax3.invert_yaxis()
ax3.grid(axis='x', alpha=0.3)

legend_handles = [
    Patch(color='#C62828', alpha=0.78, label='p < 0.01（有意）'),
    Patch(color='#FF8F00', alpha=0.78, label='p < 0.05（有意）'),
    Patch(color='#9E9E9E', alpha=0.78, label='n.s.（非有意）'),
]
ax3.legend(handles=legend_handles, fontsize=9, loc='lower right')

# p値ラベル
for i, (c, p) in enumerate(zip(coefs, pvals)):
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else ''
    offset = 0.05 if c >= 0 else -0.05
    ax3.text(c + offset, i, sig, va='center',
             ha='left' if c >= 0 else 'right', fontsize=10, fontweight='bold')

plt.tight_layout()
fig3_path = os.path.join(FIG_DIR, '2022_H5_4_fig3_coef.png')
fig3.savefig(fig3_path, bbox_inches='tight', dpi=150)
plt.close(fig3)
print(f"図3保存: {os.path.basename(fig3_path)}")

import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from matplotlib.patches import Patch
import statsmodels.api as sm
from scipy import stats
import warnings
warnings.filterwarnings('ignore')

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

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

y_ols = df2022_clean['保健医療費_千円'].values
X_ols = df2022_clean[['食料費割合', '高齢化率', '消費支出_万円', '医療機関密度']].values
X_ols_c = sm.add_constant(X_ols)
ols_model = sm.OLS(y_ols, X_ols_c).fit()

print("\n■ OLS回帰結果（目的変数: 保健医療費_千円）")
print(ols_model.summary2())

ols_varnames = ['食料費割合', '高齢化率', '消費支出_万円', '医療機関密度']
coefs  = ols_model.params[1:]
ses    = ols_model.bse[1:]
pvals  = ols_model.pvalues[1:]

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

# 2022年の地域別データ
box_data   = []
box_labels = []
for region in region_order:
    vals = df2022_clean.loc[df2022_clean['地域ブロック'] == region, '食料費割合'].dropna().values
    box_data.append(vals)
    box_labels.append(f'{region}\n(n={len(vals)})')

bp = ax4.boxplot(
    box_data, patch_artist=True, notch=False,
    medianprops=dict(color='black', linewidth=2.0),
    flierprops=dict(marker='o', markersize=5, alpha=0.6),
    whiskerprops=dict(linewidth=1.4),
    capprops=dict(linewidth=1.4),
)

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

# 全国平均線
national_mean = df2022_clean['食料費割合'].mean()
ax4.axhline(national_mean, color='black', linestyle='--',
            linewidth=1.5, label=f'全国平均 {national_mean:.1f}%')

ax4.set_xticklabels(box_labels, fontsize=10)
ax4.set_ylabel('食料費割合（%）', fontsize=12)
ax4.set_title('地域ブロック別 食料費割合の分布（2022年）\n'
              '箱ひげ図：中央値・四分位範囲・外れ値（SSDSE-B実データ）',
              fontsize=13, fontweight='bold')
ax4.grid(axis='y', alpha=0.3)
ax4.legend(fontsize=10)

plt.tight_layout()
fig4_path = os.path.join(FIG_DIR, '2022_H5_4_fig4_boxplot.png')
fig4.savefig(fig4_path, bbox_inches='tight', dpi=150)
plt.close(fig4)
print(f"図4保存: {os.path.basename(fig4_path)}")

print("\n" + "=" * 65)
print("■ 全図の生成完了（4枚）")
print("=" * 65)
print(f"  fig1_timeseries.png : 食料費割合の地域別時系列（2012〜2023年）")
print(f"  fig2_scatter.png    : 食料費割合 vs 保健医療費 散布図（2022年）")
print(f"  fig3_coef.png       : OLS回帰係数プロット（保健医療費の決定要因）")
print(f"  fig4_boxplot.png    : 地域別 食料費割合 箱ひげ図（6ブロック比較）")
print(f"  保存先: {os.path.abspath(FIG_DIR)}")