都道府県別食料費支出の地域格差分析
消費構造と社会経済要因の重回帰

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)

# 全国行を除外（都道府県のみ）
df_b = df_b[df_b['地域コード'] != 'R00000'].copy()

print("=== SSDSE-B-2026 columns ===")
print(df_b.columns.tolist())
print(f"\nShape: {df_b.shape}")
print(f"Years: {sorted(df_b['年度'].unique())}")

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

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

# 数値型変換
num_cols = ['食料費（二人以上の世帯）', '消費支出（二人以上の世帯）',
            '光熱・水道費（二人以上の世帯）', '住居費（二人以上の世帯）',
            '保健医療費（二人以上の世帯）', '総人口', '65歳以上人口']
for c in num_cols:
    df_b[c] = pd.to_numeric(df_b[c], errors='coerce')

latest = df_b[df_b['年度'] == 2022].copy()
latest['高齢化率'] = latest['65歳以上人口'] / latest['総人口'] * 100
latest['食料費率'] = (latest['食料費（二人以上の世帯）'] /
                     latest['消費支出（二人以上の世帯）'] * 100)

# 全国平均（SSDSE-B には全国行がないため 47都道府県平均）
national_avg = latest['食料費（二人以上の世帯）'].mean()
national_rate = latest['食料費率'].mean()

print(f"\n=== 基本統計（2022年）===")
print(f"全国平均 食料費: {national_avg:.0f}円/月")
print(f"全国平均 食料費率: {national_rate:.1f}%")
print(f"\n上位5都道府県（食料費）:")
top5 = latest.nlargest(5, '食料費（二人以上の世帯）')[['都道府県', '食料費（二人以上の世帯）', '地域']]
print(top5.to_string(index=False))
print(f"\n下位5都道府県（食料費）:")
bot5 = latest.nsmallest(5, '食料費（二人以上の世帯）')[['都道府県', '食料費（二人以上の世帯）', '地域']]
print(bot5.to_string(index=False))

print("\n=== 地域別 食料費平均（2022年） ===")
region_stats = latest.groupby('地域')['食料費（二人以上の世帯）'].agg(['mean', 'std'])
region_stats.columns = ['平均', '標準偏差']
region_stats['変動係数(%)'] = region_stats['標準偏差'] / region_stats['平均'] * 100
print(region_stats.round(0).to_string())

print("\n=== 時系列変化（COVID前後比較）===")
pre_covid  = df_b[df_b['年度'] == 2019].groupby('地域')['食料費（二人以上の世帯）'].mean()
post_covid = df_b[df_b['年度'] == 2021].groupby('地域')['食料費（二人以上の世帯）'].mean()
change = (post_covid - pre_covid) / pre_covid * 100
print("COVID前後（2019→2021）食料費変化率:")
print(change.round(1).to_string())

print("\nDONE: 2019_H5_2_shorei")

fig1, ax1 = plt.subplots(figsize=(14, 8))

sorted_latest = latest.sort_values('食料費（二人以上の世帯）', ascending=False).reset_index(drop=True)
colors_bar = [region_colors[r] for r in sorted_latest['地域']]

bars = ax1.bar(range(len(sorted_latest)),
               sorted_latest['食料費（二人以上の世帯）'],
               color=colors_bar, width=0.85, alpha=0.88, edgecolor='white', linewidth=0.5)

# 全国平均線
ax1.axhline(y=national_avg, color='#333333', linewidth=1.8, linestyle='--',
            label=f'全国平均: {national_avg:,.0f}円')

ax1.set_xticks(range(len(sorted_latest)))
ax1.set_xticklabels(sorted_latest['都道府県'].str.replace('県', '').str.replace('府', '').str.replace('都', '').str.replace('道', ''),
                     rotation=90, fontsize=7.5)
ax1.set_ylabel('食料費（円/月）', fontsize=11)
ax1.set_title('都道府県別 食料費（二人以上の世帯）ランキング\n2022年 SSDSE-B-2026', fontsize=13, fontweight='bold', pad=14)
ax1.set_ylim(55000, 97000)
ax1.yaxis.set_major_formatter(plt.FuncFormatter(lambda x, _: f'{x/1000:.0f}千'))
ax1.grid(axis='y', alpha=0.3, linewidth=0.6)

# 凡例（地域別）
from matplotlib.patches import Patch
legend_handles = [Patch(facecolor=region_colors[r], label=r) for r in region_order]
legend_handles.append(plt.Line2D([0], [0], color='#333333', linewidth=1.8,
                                  linestyle='--', label=f'全国平均'))
ax1.legend(handles=legend_handles, loc='upper right', fontsize=8.5,
           framealpha=0.92, ncol=2)

# 変動係数
cv = latest['食料費（二人以上の世帯）'].std() / latest['食料費（二人以上の世帯）'].mean() * 100
ax1.text(0.02, 0.04, f'変動係数 CV = {cv:.1f}%', transform=ax1.transAxes,
         fontsize=9, color='#555', bbox=dict(boxstyle='round,pad=0.3', facecolor='#f5f5f5', alpha=0.8))

ax1.spines[['top', 'right']].set_visible(False)
plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2019_H5_2_fig1.png'), bbox_inches='tight')
plt.close(fig1)
print("\n[OK] fig1 saved")

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

for region in region_order:
    sub = latest[latest['地域'] == region]
    ax2.scatter(sub['消費支出（二人以上の世帯）'],
                sub['食料費（二人以上の世帯）'],
                color=region_colors[region], s=60, alpha=0.85,
                edgecolors='white', linewidth=0.7, label=region, zorder=3)

# ラベル（重ならないように選択）
for _, row in latest.iterrows():
    pref = row['都道府県'].replace('県', '').replace('府', '').replace('都', '').replace('道', '')
    ax2.annotate(pref,
                 (row['消費支出（二人以上の世帯）'], row['食料費（二人以上の世帯）']),
                 fontsize=6.5, xytext=(3, 3), textcoords='offset points',
                 alpha=0.8, color='#444')

# 回帰直線
x_vals = latest['消費支出（二人以上の世帯）']
y_vals = latest['食料費（二人以上の世帯）']
mask = x_vals.notna() & y_vals.notna()
slope, intercept, r_val, p_val, _ = stats.linregress(x_vals[mask], y_vals[mask])
xline = np.linspace(x_vals.min(), x_vals.max(), 100)
ax2.plot(xline, slope * xline + intercept, color='#333', linewidth=1.8,
         linestyle='--', alpha=0.7, label=f'回帰直線')

print(f"\nPearson r (食料費 vs 消費支出): r={r_val:.3f}, p={p_val:.4f}")

ax2.set_xlabel('消費支出（円/月）', fontsize=11)
ax2.set_ylabel('食料費（円/月）', fontsize=11)
ax2.set_title(f'食料費 vs 消費支出 （2022年・47都道府県）\nr = {r_val:.3f}, p {"< 0.001" if p_val < 0.001 else f"= {p_val:.4f}"}',
              fontsize=12, fontweight='bold', pad=12)

ax2.xaxis.set_major_formatter(plt.FuncFormatter(lambda x, _: f'{x/10000:.0f}万'))
ax2.yaxis.set_major_formatter(plt.FuncFormatter(lambda x, _: f'{x/1000:.0f}千'))
ax2.grid(alpha=0.25, linewidth=0.6)
ax2.legend(loc='upper left', fontsize=8.5, framealpha=0.92)
ax2.spines[['top', 'right']].set_visible(False)

# r, p テキスト
ax2.text(0.97, 0.05,
         f'r = {r_val:.3f}\np < 0.001',
         transform=ax2.transAxes, fontsize=10, ha='right', va='bottom',
         bbox=dict(boxstyle='round,pad=0.5', facecolor='#fff3e0', alpha=0.9))

plt.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2019_H5_2_fig2.png'), bbox_inches='tight')
plt.close(fig2)
print("[OK] fig2 saved")

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

# COVID期間
ax3.axvspan(2020, 2021, color='#cccccc', alpha=0.35, label='COVID-19 期（2020-2021）')

# 地域別年平均
ts_data = df_b.groupby(['年度', '地域'])['食料費（二人以上の世帯）'].mean().reset_index()

for region in region_order:
    sub = ts_data[ts_data['地域'] == region].sort_values('年度')
    ax3.plot(sub['年度'], sub['食料費（二人以上の世帯）'],
             color=region_colors[region], linewidth=2.2,
             marker='o', markersize=5, label=region)

ax3.set_xlabel('年度', fontsize=11)
ax3.set_ylabel('食料費（円/月）', fontsize=11)
ax3.set_title('地域別 食料費の時系列推移（2012〜2022年）\n6地域別平均・SSDSE-B-2026',
              fontsize=12, fontweight='bold', pad=12)
ax3.set_xticks(sorted(df_b['年度'].unique()))
ax3.set_xticklabels(sorted(df_b['年度'].unique()), rotation=45, fontsize=9)
ax3.yaxis.set_major_formatter(plt.FuncFormatter(lambda x, _: f'{x/1000:.0f}千'))
ax3.grid(alpha=0.25, linewidth=0.6)
ax3.legend(loc='upper left', fontsize=9, framealpha=0.92, ncol=2)
ax3.spines[['top', 'right']].set_visible(False)

# COVID注釈
ax3.text(2020.5, ax3.get_ylim()[0] * 1.005 + ax3.get_ylim()[1] * 0.005,
         'COVID-19', fontsize=8.5, ha='center', color='#666', style='italic')

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

y_reg = latest['食料費（二人以上の世帯）']
X_reg = latest[['消費支出（二人以上の世帯）', '高齢化率', '光熱・水道費（二人以上の世帯）']].copy()
mask_reg = X_reg.notna().all(axis=1) & y_reg.notna()
X_reg_c = X_reg[mask_reg]
y_reg_c = y_reg[mask_reg]

# 標準化
X_std = (X_reg_c - X_reg_c.mean()) / X_reg_c.std()
y_std = (y_reg_c - y_reg_c.mean()) / y_reg_c.std()
X_sm = sm.add_constant(X_std)
model = sm.OLS(y_std, X_sm).fit()

coefs = model.params.drop('const')
cis   = model.conf_int().drop('const')
pvals = model.pvalues.drop('const')

var_labels = ['消費支出\n（二人以上世帯）', '高齢化率\n（65歳以上比率）', '光熱・水道費\n（二人以上世帯）']

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

print(f"\n=== 重回帰結果（標準化係数）===")
print(f"R² = {model.rsquared:.3f}, Adj.R² = {model.rsquared_adj:.3f}")
for vn, coef, pv in zip(var_labels, coefs.values, pvals.values):
    print(f"  {vn.replace(chr(10),' ')}: β={coef:.3f}, p={pv:.4f} {sig_mark(pv)}")

fig4, ax4 = plt.subplots(figsize=(9, 5))

colors_coef = ['#e05c5c' if c > 0 else '#4e9af1' for c in coefs.values]
y_pos = range(len(coefs))

bars4 = ax4.barh(y_pos, coefs.values,
                  color=colors_coef, alpha=0.85,
                  edgecolor='white', linewidth=0.8, height=0.55)

# 95% CI
for i, (lo, hi) in enumerate(zip(cis[0].values, cis[1].values)):
    ax4.plot([lo, hi], [i, i], color='#333', linewidth=2.5, zorder=4)
    ax4.plot([lo, lo], [i - 0.12, i + 0.12], color='#333', linewidth=2)
    ax4.plot([hi, hi], [i - 0.12, i + 0.12], color='#333', linewidth=2)

# 有意水準マーク
for i, (coef, pv) in enumerate(zip(coefs.values, pvals.values)):
    mark = sig_mark(pv)
    offset = 0.02
    x_pos = coef + offset if coef >= 0 else coef - offset
    ha = 'left' if coef >= 0 else 'right'
    color = '#c0392b' if pv < 0.05 else '#888'
    ax4.text(x_pos, i, mark, va='center', ha=ha, fontsize=13, color=color, fontweight='bold')

ax4.set_yticks(list(y_pos))
ax4.set_yticklabels(var_labels, fontsize=10)
ax4.axvline(x=0, color='#333', linewidth=1.2, linestyle='-', alpha=0.6)
ax4.set_xlabel('標準化偏回帰係数 (β)', fontsize=11)
ax4.set_title(f'重回帰分析：食料費の決定要因（標準化偏回帰係数 + 95%CI）\nN=47都道府県, R²={model.rsquared:.3f}, 2022年',
              fontsize=11, fontweight='bold', pad=14)
ax4.set_xlim(-0.85, 0.95)
ax4.grid(axis='x', alpha=0.3, linewidth=0.6)
ax4.spines[['top', 'right']].set_visible(False)

# 凡例説明
from matplotlib.patches import Patch
legend_e = [Patch(facecolor='#e05c5c', label='正の効果（食料費を増加）'),
            Patch(facecolor='#4e9af1', label='負の効果（食料費を減少）')]
ax4.legend(handles=legend_e, loc='lower right', fontsize=9, framealpha=0.9)

# 有意水準凡例
ax4.text(0.01, -0.12,
         '有意水準: *** p<0.001  ** p<0.01  * p<0.05  n.s. 非有意',
         transform=ax4.transAxes, fontsize=8.5, color='#555')

plt.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2019_H5_2_fig4.png'), bbox_inches='tight')
plt.close(fig4)
print("[OK] fig4 saved")