"""
教育用再現コード: 2022年 統計データ分析コンペティション 審査員奨励賞 [高校生の部]
=================================================================
論文タイトル：食料費支出と健康：都道府県別食生活指標の回帰分析
受賞：審査員奨励賞（高校生の部）

【分析概要】
  データ：SSDSE-B-2026.csv（都道府県別パネルデータ, 2012〜2023年度）
  対象：47都道府県 × 12年間（時系列・横断面）

  食料費支出（エンゲル係数類似指標）と健康関連アウトカム（保健医療費）の
  関連を相関分析・OLS回帰・時系列・地域比較で多角的に分析する。

  分析の流れ
  Fig1. 食料費割合の地域別時系列（2012〜2023年）
  Fig2. 食料費割合 vs 保健医療費 散布図（2022年・47都道府県ラベル付き）
  Fig3. OLS回帰係数プロット（保健医療費の決定要因）
  Fig4. 地域ブロック別 食料費割合の箱ひげ図比較

【変数設計】
  食料費割合（%）       = 食料費（二人以上の世帯）/ 消費支出（二人以上の世帯）× 100
  保健医療費_千円        = 保健医療費（二人以上の世帯）/ 1000
  高齢化率（%）         = 65歳以上人口 / 総人口 × 100
  消費支出_万円          = 消費支出（二人以上の世帯）/ 10000
  医療機関密度（/万人）  = 一般診療所数 / 総人口 × 10000

【データ出典】
  SSDSE-B-2026.csv: 社会・人口統計体系（都道府県データ）
  統計数理研究所・統計センター公表の実データ

【データサイエンス学習ポイント】
  1. エンゲル係数と食料費割合（計算式と消費水準との関係）
  2. 健康の代理変数（医療費を健康指標として使う根拠と限界）
  3. 外れ値の検出と対処（箱ひげ図・残差プロット）
  4. 食生活政策への統計的示唆
=================================================================
"""

# ============================================================
# 【データの準備】実行前に以下のデータファイルを用意してください
#
#   必要ファイル:
#     ・SSDSE-B-2026.csv
#       → data/raw/SSDSE-B-2026.csv に配置
#
#   ダウンロード先:
#     https://www.nstac.go.jp/use/literacy/ssdse/
#     （SSDSE-B（社会・人口統計体系 都道府県データ） の CSV をダウンロード）
#
#   フォルダ配置（プロジェクトルートからの相対パス）:
#     code/                ← このスクリプトの場所
#     data/raw/            ← CSV ファイルをここに配置
#     html/figures/        ← 図の出力先（自動生成）
#
#   実行方法（ファイルを一切編集せず実行可能）:
#     python3 code/2022_H5_4_shorei.py
# ============================================================


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)

# ================================================================
# ■ 実データ読み込み（SSDSE-B-2026: 都道府県データ）
# ================================================================
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}")

# ================================================================
# ■ OLS回帰（保健医療費_千円 の決定要因）
# ================================================================
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:]

# ================================================================
# ■ 図の生成（4枚）
# ================================================================

# ────────────────────────────────────────────────────────────────
# 図1: 食料費割合の地域別時系列（2012〜2023年）
# ────────────────────────────────────────────────────────────────
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)}")

# ────────────────────────────────────────────────────────────────
# 図2: 食料費割合 vs 保健医療費 散布図（2022年・都道府県ラベル付き）
# ────────────────────────────────────────────────────────────────
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)}")

# ────────────────────────────────────────────────────────────────
# 図3: OLS回帰係数プロット（保健医療費の決定要因）
# ────────────────────────────────────────────────────────────────
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)}")

# ────────────────────────────────────────────────────────────────
# 図4: 地域別 食料費割合の箱ひげ図（6地域ブロック比較）
# ────────────────────────────────────────────────────────────────
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)}")