味覚から眺める地域別「ふるさとの味」

df_b_raw = pd.read_csv(DATA_B, encoding='cp932', header=1)
df_b = df_b_raw[df_b_raw['地域コード'].str.match(r'^R\d{5}$', na=False)].copy()

print("=" * 65)
print("■ SSDSE-B-2026 都道府県データ読み込み")
print("=" * 65)
print(f"  行数: {len(df_b)}, 年度範囲: {df_b['年度'].min()}〜{df_b['年度'].max()}")
print(f"  都道府県数（年平均）: {df_b['年度'].value_counts().median():.0f}")

df_b['pref_num'] = df_b['地域コード'].str.extract(r'R(\d{2})').astype(int)

# 地域区分の定義
def assign_region(n):
    if n == 1:
        return '北海道'
    elif 2 <= n <= 7:
        return '東北'
    elif 8 <= n <= 14:
        return '関東'
    elif 15 <= n <= 23:
        return '中部'
    elif 24 <= n <= 30:
        return '関西'
    elif 31 <= n <= 39:
        return '中国・四国'
    elif 40 <= n <= 47:
        return '九州・沖縄'
    else:
        return 'その他'

df_b['地域'] = df_b['pref_num'].apply(assign_region)

df_b['食料費割合'] = df_b['食料費（二人以上の世帯）'] / df_b['消費支出（二人以上の世帯）'] * 100

# 光熱・水道費割合（調理・生活エネルギー使用の代理変数）
df_b['光熱・水道費割合'] = df_b['光熱・水道費（二人以上の世帯）'] / df_b['消費支出（二人以上の世帯）'] * 100

# 教養娯楽費割合（生活スタイルの指標）
df_b['教養娯楽費割合'] = df_b['教養娯楽費（二人以上の世帯）'] / df_b['消費支出（二人以上の世帯）'] * 100

print(f"\n  食料費割合の基本統計:")
print(df_b['食料費割合'].describe().round(2))

# 関東・関西フラグ
KANTO_PREFS  = list(range(8, 15))   # R08〜R14（茨城〜神奈川）
KANSAI_PREFS = list(range(24, 31))  # R24〜R30（三重〜和歌山）

df_b['地域2'] = df_b['pref_num'].apply(
    lambda n: '関東' if n in KANTO_PREFS else ('関西' if n in KANSAI_PREFS else 'その他')
)

latest_year = df_b['年度'].max()
df_latest = df_b[df_b['年度'] == latest_year].copy()
print(f"\n  最新年度: {latest_year}年 ({len(df_latest)}都道府県)")

print("\n" + "=" * 65)
print("■ Kruskal-Wallis検定：6地域別 食料費割合の差")
print("=" * 65)

REGION_ORDER = ['北海道', '東北', '関東', '中部', '関西', '中国・四国', '九州・沖縄']

# 全年度データを使って地域別に集計
region_groups = [
    df_b[df_b['地域'] == r]['食料費割合'].dropna().values
    for r in REGION_ORDER
]

kw_stat, kw_p = stats.kruskal(*region_groups)
print(f"  H統計量={kw_stat:.4f}, p値={kw_p:.6f}")
sig_kw = '*** (p<0.001)' if kw_p < 0.001 else '** (p<0.01)' if kw_p < 0.01 else '* (p<0.05)' if kw_p < 0.05 else 'n.s.'
print(f"  判定: {sig_kw}")

# 地域別の記述統計
print(f"\n  {'地域':<12} {'N':>5} {'中央値':>8} {'平均':>8} {'SD':>7}")
print("  " + "-" * 45)
for r, grp in zip(REGION_ORDER, region_groups):
    print(f"  {r:<12} {len(grp):>5} {np.median(grp):>8.2f} {np.mean(grp):>8.2f} {np.std(grp):>7.2f}")

print("\n" + "=" * 65)
print("■ Mann-Whitney U検定：関東 vs 関西（全年度）")
print("=" * 65)

kanto_vals  = df_b[df_b['地域2'] == '関東']['食料費割合'].dropna().values
kansai_vals = df_b[df_b['地域2'] == '関西']['食料費割合'].dropna().values

u_stat, u_p = stats.mannwhitneyu(kanto_vals, kansai_vals, alternative='two-sided')
print(f"  U統計量={u_stat:.1f}, p値={u_p:.6f}")
sig_u = '*** (p<0.001)' if u_p < 0.001 else '** (p<0.01)' if u_p < 0.01 else '* (p<0.05)' if u_p < 0.05 else 'n.s.'
print(f"  判定: {sig_u}")
print(f"  関東 中央値={np.median(kanto_vals):.2f}%, 平均={np.mean(kanto_vals):.2f}%")
print(f"  関西 中央値={np.median(kansai_vals):.2f}%, 平均={np.mean(kansai_vals):.2f}%")

# 最新年度の関東・関西比較
kanto_latest  = df_latest[df_latest['地域2'] == '関東']
kansai_latest = df_latest[df_latest['地域2'] == '関西']
metrics = {
    '食料費割合（%）':    ('食料費割合', 2),
    '光熱・水道費割合（%）': ('光熱・水道費割合', 2),
    '教養娯楽費割合（%）': ('教養娯楽費割合', 2),
}
print(f"\n  {latest_year}年 関東 vs 関西 比較:")
print(f"  {'指標':<22} {'関東':>10} {'関西':>10}")
print("  " + "-" * 45)
for label, (col, nd) in metrics.items():
    v_k = kanto_latest[col].mean()
    v_ks = kansai_latest[col].mean()
    print(f"  {label:<22} {v_k:>10.{nd}f} {v_ks:>10.{nd}f}")

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

bp_data = [grp for grp in region_groups]
bp = ax1.boxplot(
    bp_data,
    labels=REGION_ORDER,
    patch_artist=True,
    medianprops=dict(color='black', linewidth=2),
    flierprops=dict(marker='o', markersize=4, alpha=0.6),
    boxprops=dict(linewidth=1.2),
    whiskerprops=dict(linewidth=1.2),
    capprops=dict(linewidth=1.5),
)

for patch, region in zip(bp['boxes'], REGION_ORDER):
    c = REGION_COLORS[region]
    patch.set_facecolor(c)
    patch.set_alpha(0.65)

for flier, region in zip(bp['fliers'], REGION_ORDER):
    flier.set(markerfacecolor=REGION_COLORS[region], alpha=0.6)

ax1.set_ylabel('食料費割合（食料費 / 消費支出，%）', fontsize=12)
ax1.set_xlabel('地域', fontsize=12)
ax1.set_title(
    f'6地域別 食料費割合の分布（2012〜{latest_year}年）\n'
    f'Kruskal-Wallis: H={kw_stat:.2f}, {sig_kw}',
    fontsize=13, fontweight='bold'
)
ax1.grid(axis='y', alpha=0.3)

# 各地域の中央値をテキスト表示
for i, grp in enumerate(region_groups, start=1):
    med = np.median(grp)
    ax1.text(i, med + 0.1, f'{med:.1f}%', ha='center', va='bottom', fontsize=9,
             fontweight='bold', color=REGION_COLORS[REGION_ORDER[i-1]])

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

fig2, axes2 = plt.subplots(1, 3, figsize=(13, 5))
fig2.suptitle(
    f'関東 vs 関西 消費支出内訳の比較（{latest_year}年，各都道府県平均）',
    fontsize=13, fontweight='bold'
)

metric_list = [
    ('食料費割合', '食料費割合\n（食料費 / 消費支出，%）',  '#C62828', '#6A1B9A'),
    ('光熱・水道費割合', '光熱・水道費割合\n（光熱・水道費 / 消費支出，%）', '#E65100', '#9C27B0'),
    ('教養娯楽費割合', '教養娯楽費割合\n（教養娯楽費 / 消費支出，%）', '#1565C0', '#00695C'),
]

for ax, (col, ylabel, c_kanto, c_kansai) in zip(axes2, metric_list):
    v_kanto  = kanto_latest[col].mean()
    v_kansai = kansai_latest[col].mean()
    se_kanto  = kanto_latest[col].sem()
    se_kansai = kansai_latest[col].sem()

    bars = ax.bar(['関東', '関西'],
                  [v_kanto, v_kansai],
                  color=[c_kanto, c_kansai], alpha=0.75,
                  edgecolor='white', linewidth=1.2, width=0.5)
    ax.errorbar(['関東', '関西'], [v_kanto, v_kansai],
                yerr=[1.96*se_kanto, 1.96*se_kansai],
                fmt='none', color='#333', capsize=5, linewidth=1.5)

    # 値のラベル
    for bar, val in zip(bars, [v_kanto, v_kansai]):
        ax.text(bar.get_x() + bar.get_width() / 2,
                bar.get_height() + 0.05,
                f'{val:.2f}%', ha='center', va='bottom',
                fontsize=11, fontweight='bold')

    # Mann-Whitney検定
    col_kanto_series  = kanto_latest[col].dropna().values
    col_kansai_series = kansai_latest[col].dropna().values
    if len(col_kanto_series) >= 2 and len(col_kansai_series) >= 2:
        _, p_tmp = stats.mannwhitneyu(col_kanto_series, col_kansai_series, alternative='two-sided')
        sig_tmp = '***' if p_tmp < 0.001 else '**' if p_tmp < 0.01 else '*' if p_tmp < 0.05 else 'n.s.'
        ax.set_title(f'Mann-Whitney: {sig_tmp}', fontsize=10)

    ax.set_ylabel(ylabel, fontsize=9)
    ax.set_ylim(0, max(v_kanto, v_kansai) * 1.25)
    ax.grid(axis='y', alpha=0.3)
    ax.tick_params(axis='x', labelsize=12)

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

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

years = sorted(df_b['年度'].unique())

# 全国平均
national_ts = df_b.groupby('年度')['食料費割合'].mean()

# 関東平均
kanto_ts  = df_b[df_b['地域2'] == '関東'].groupby('年度')['食料費割合'].mean()
kansai_ts = df_b[df_b['地域2'] == '関西'].groupby('年度')['食料費割合'].mean()

# 関東・関西の信頼区間（SE×1.96）
kanto_sem  = df_b[df_b['地域2'] == '関東'].groupby('年度')['食料費割合'].sem()
kansai_sem = df_b[df_b['地域2'] == '関西'].groupby('年度')['食料費割合'].sem()

ax3.fill_between(kanto_ts.index,
                 kanto_ts - 1.96 * kanto_sem,
                 kanto_ts + 1.96 * kanto_sem,
                 color='#C62828', alpha=0.12)
ax3.fill_between(kansai_ts.index,
                 kansai_ts - 1.96 * kansai_sem,
                 kansai_ts + 1.96 * kansai_sem,
                 color='#6A1B9A', alpha=0.12)

ax3.plot(national_ts.index, national_ts.values,
         color='#455A64', linewidth=2.0, linestyle='--',
         marker='o', markersize=5, label='全国平均')
ax3.plot(kanto_ts.index, kanto_ts.values,
         color='#C62828', linewidth=2.5,
         marker='s', markersize=6, label='関東（7都県平均）')
ax3.plot(kansai_ts.index, kansai_ts.values,
         color='#6A1B9A', linewidth=2.5,
         marker='^', markersize=6, label='関西（7府県平均）')

ax3.set_xlabel('年度', fontsize=12)
ax3.set_ylabel('食料費割合（食料費 / 消費支出，%）', fontsize=12)
ax3.set_title(
    '食料費割合の時系列推移：関東・関西・全国平均（SSDSE-B 実データ）',
    fontsize=13, fontweight='bold'
)
ax3.set_xticks(years)
ax3.set_xticklabels([str(y) for y in years], rotation=45, fontsize=9)
ax3.legend(fontsize=11)
ax3.grid(alpha=0.3)

# 直近年の差を注釈
diff_latest = kanto_ts[latest_year] - kansai_ts[latest_year]
arrow_x = latest_year
ax3.annotate(
    f'差={diff_latest:+.2f}pp',
    xy=(arrow_x, kanto_ts[latest_year]),
    xytext=(arrow_x - 1.5, kanto_ts[latest_year] + 0.3),
    fontsize=9, color='#333',
    arrowprops=dict(arrowstyle='->', color='#333', lw=1.2)
)

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

import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
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)

print("\n" + "=" * 65)
print("■ 相関分析：住宅地価格 vs 食料費割合（最新年度）")
print("=" * 65)

scatter_df = df_latest[['都道府県', '地域', '食料費割合',
                          '標準価格（平均価格）（住宅地）']].dropna()
r_scatter, p_scatter = stats.pearsonr(
    scatter_df['標準価格（平均価格）（住宅地）'],
    scatter_df['食料費割合']
)
print(f"  r={r_scatter:.4f}, p={p_scatter:.4f}")
sig_s = '*** (p<0.001)' if p_scatter < 0.001 else '** (p<0.01)' if p_scatter < 0.01 else '* (p<0.05)' if p_scatter < 0.05 else 'n.s.'
print(f"  判定: {sig_s}")

REGION_COLORS = {
    '北海道':    '#1565C0',
    '東北':      '#2E7D32',
    '関東':      '#C62828',
    '中部':      '#E65100',
    '関西':      '#6A1B9A',
    '中国・四国': '#00695C',
    '九州・沖縄': '#F9A825',
}

fig4, ax4 = plt.subplots(figsize=(10, 7))

for region in REGION_ORDER:
    sub = scatter_df[scatter_df['地域'] == region]
    ax4.scatter(
        sub['標準価格（平均価格）（住宅地）'] / 1000,   # 千円/m² 単位
        sub['食料費割合'],
        color=REGION_COLORS[region],
        s=80, alpha=0.75, zorder=3,
        label=region, edgecolors='white', linewidths=0.5
    )

# 全体回帰直線
x_sc = scatter_df['標準価格（平均価格）（住宅地）'].values / 1000
y_sc = scatter_df['食料費割合'].values
z_sc = np.polyfit(x_sc, y_sc, 1)
xs_sc = np.linspace(x_sc.min(), x_sc.max(), 200)
ax4.plot(xs_sc, np.poly1d(z_sc)(xs_sc),
         color='#333', linewidth=1.8, linestyle='--', zorder=2,
         label=f'回帰直線 r={r_scatter:.3f} ({sig_s})')

# 注目都道府県をラベル表示
highlight_prefs = ['東京都', '大阪府', '京都府', '神奈川県', '北海道', '沖縄県', '秋田県']
for _, row in scatter_df.iterrows():
    if row['都道府県'] in highlight_prefs:
        ax4.annotate(
            row['都道府県'],
            xy=(row['標準価格（平均価格）（住宅地）'] / 1000, row['食料費割合']),
            xytext=(6, 4), textcoords='offset points',
            fontsize=8.5, color='#333',
        )

ax4.set_xlabel('標準価格・住宅地（千円/m²）', fontsize=12)
ax4.set_ylabel('食料費割合（食料費 / 消費支出，%）', fontsize=12)
ax4.set_title(
    f'住宅地価格 vs 食料費割合の関係（{latest_year}年，47都道府県）\n'
    f'相関係数 r={r_scatter:.3f}，{sig_s}',
    fontsize=13, fontweight='bold'
)
ax4.legend(fontsize=9, loc='upper right', framealpha=0.85)
ax4.grid(alpha=0.3)

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

print("\n全図の生成完了（4枚）")
print(f"  2022_U5_4_fig1_region.png      : 6地域別 食料費割合 箱ひげ図")
print(f"  2022_U5_4_fig2_kanto_kansai.png: 関東 vs 関西 複数指標比較棒グラフ")
print(f"  2022_U5_4_fig3_ts.png          : 食料費割合の時系列推移")
print(f"  2022_U5_4_fig4_scatter.png     : 住宅地価格 vs 食料費割合 散布図")

print("\n" + "=" * 65)
print("■ 分析サマリー（論文の問いへの答え）")
print("=" * 65)
print(f"  ・6地域間の食料費割合の差 → Kruskal-Wallis {sig_kw}")
print(f"  ・関東 vs 関西（食料費割合）→ Mann-Whitney {sig_u}")
print(f"  ・関東 食料費割合平均: {np.mean(kanto_vals):.2f}%")
print(f"  ・関西 食料費割合平均: {np.mean(kansai_vals):.2f}%")
diff_mean = np.mean(kanto_vals) - np.mean(kansai_vals)
print(f"  ・差（関東－関西）: {diff_mean:+.2f}pp")
print(f"  ・住宅地価格 vs 食料費割合 相関: r={r_scatter:.3f} {sig_s}")
print("=" * 65)