"""
2018_H4_katsuyo.py
都道府県別スポーツ活動・体力と健康指標の関係：重回帰による規定要因分析
統計活用奨励賞（高校生部門）2018年度

【プロキシ変数の説明】
SSDSE-B-2026には直接のスポーツ・体力データが含まれないため、
以下の代理変数を使用する：
  目的変数  : 保健医療費（二人以上の世帯）[円/月]
              = 家計の医療・保健支出（健康意識・健康状態の代理指標）
              スポーツ活動が盛んな地域は予防医療に敏感であり、
              また医療アクセスも高い傾向があると仮説を立てる。
  説明変数1 : 消費支出_万 = 消費支出（二人以上の世帯）/ 10000
              （所得水準の代理：スポーツ参加には費用がかかる）
  説明変数2 : 教養娯楽費_千 = 教養娯楽費（二人以上の世帯）/ 1000
              （余暇・スポーツ活動への支出の代理指標）
  説明変数3 : 高齢化率 = 65歳以上人口 / 総人口 × 100 [%]
              （高齢化は体力低下・医療需要増加の直接要因）
  説明変数4 : 光熱水道費_千 = 光熱・水道費（二人以上の世帯）/ 1000
              （生活基盤・生活様式の指標）
"""

# ============================================================
# 【データの準備】実行前に以下のデータファイルを用意してください
#
#   必要ファイル:
#     ・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/2018_H4_katsuyo.py
# ============================================================


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
from scipy.stats import kruskal
from matplotlib.patches import Patch

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)

print("=== df_b.columns ===")
print(df_b.columns.tolist())
print()

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

def get_short_name(full_name):
    for suffix in ['都', '道', '府', '県']:
        if full_name.endswith(suffix):
            return full_name[:-1]
    return full_name

region_order = ['北海道・東北', '関東', '中部', '近畿', '中国・四国', '九州・沖縄']
region_colors = {
    '北海道・東北': '#4e9af1',
    '関東':         '#e05c5c',
    '中部':         '#f0a500',
    '近畿':         '#5cb85c',
    '中国・四国':   '#9b59b6',
    '九州・沖縄':   '#f39c12'
}

# ── 分析用データセット作成（2023年断面） ─────────────────────────
latest_yr = 2023
df_cross = df_b[df_b['年度'] == latest_yr].copy()

df_cross['short_name'] = df_cross['都道府県'].apply(get_short_name)
df_cross['region']     = df_cross['short_name'].map(region_map_raw)

# プロキシ変数の計算
df_cross['高齢化率']      = df_cross['65歳以上人口'] / df_cross['総人口'] * 100
df_cross['消費支出_万']   = df_cross['消費支出（二人以上の世帯）'] / 10000
df_cross['教養娯楽費_千'] = df_cross['教養娯楽費（二人以上の世帯）'] / 1000
df_cross['光熱水道費_千'] = df_cross['光熱・水道費（二人以上の世帯）'] / 1000

y_col = '保健医療費（二人以上の世帯）'
analysis_cols = [y_col, '消費支出_万', '教養娯楽費_千', '高齢化率', '光熱水道費_千', 'region']
df_cross = df_cross.dropna(subset=analysis_cols).copy()

print(f"分析サンプル数: {len(df_cross)} 都道府県（{latest_yr}年）")
print()
print("=== 基本統計量 ===")
print(df_cross[[y_col, '消費支出_万', '教養娯楽費_千', '高齢化率', '光熱水道費_千']].describe().round(2))
print()

national_avg = df_cross[y_col].mean()
print(f"保健医療費（全国平均）: {national_avg:.0f} 円/月")
print()

# Top/Bottom 5 都道府県
df_sorted_rank = df_cross.sort_values(y_col, ascending=False)
print("=== 保健医療費 上位5都道府県 ===")
print(df_sorted_rank[['short_name', y_col]].head(5).to_string(index=False))
print("\n=== 保健医療費 下位5都道府県 ===")
print(df_sorted_rank[['short_name', y_col]].tail(5).to_string(index=False))
print()

# ════════════════════════════════════════════════════════════════════
# Figure 1: 都道府県別保健医療費ランキング棒グラフ
# ════════════════════════════════════════════════════════════════════
df_sorted = df_cross.sort_values(y_col, ascending=False).copy()

fig1, ax1 = plt.subplots(figsize=(14, 7))
colors_bar = [region_colors[r] for r in df_sorted['region']]
ax1.bar(range(len(df_sorted)), df_sorted[y_col] / 1000,
        color=colors_bar, edgecolor='white', linewidth=0.5, alpha=0.88)

ax1.axhline(national_avg / 1000, color='#333333', linewidth=1.8,
            linestyle='--', label=f'全国平均: {national_avg/1000:.1f} 千円/月')

ax1.set_xticks(range(len(df_sorted)))
ax1.set_xticklabels(df_sorted['short_name'], rotation=55, ha='right', fontsize=8.5)
ax1.set_ylabel('保健医療費（千円/月）', fontsize=12)
ax1.set_title('図1：都道府県別 保健医療費（2023年）\n'
              '〈健康意識・健康状態の代理指標：二人以上の世帯の月次保健医療支出〉',
              fontsize=13, fontweight='bold', pad=14)
ax1.set_xlim(-0.8, len(df_sorted) - 0.2)
ax1.grid(axis='y', alpha=0.3, linewidth=0.8)
ax1.spines['top'].set_visible(False)
ax1.spines['right'].set_visible(False)

legend_elements = [Patch(facecolor=region_colors[r], label=r, alpha=0.88)
                   for r in region_order]
legend_elements.append(
    plt.Line2D([0], [0], color='#333333', linewidth=1.8, linestyle='--',
               label=f'全国平均: {national_avg/1000:.1f} 千円/月')
)
ax1.legend(handles=legend_elements, fontsize=9, loc='upper right', ncol=2)

plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2018_H4_fig1.png'), bbox_inches='tight')
plt.close(fig1)
print("fig1 saved.")

# ════════════════════════════════════════════════════════════════════
# Figure 2: 消費支出 vs 保健医療費の散布図（47都道府県ラベル付き）
# ════════════════════════════════════════════════════════════════════
x2 = df_cross['消費支出_万'].values
y2 = df_cross[y_col].values / 1000  # 千円
r2, p2 = stats.pearsonr(x2, y2)

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

for _, row in df_cross.iterrows():
    c = region_colors[row['region']]
    ax2.scatter(row['消費支出_万'], row[y_col] / 1000,
                color=c, s=72, alpha=0.85, zorder=3,
                edgecolors='white', linewidths=0.8)
    ax2.annotate(row['short_name'],
                 (row['消費支出_万'], row[y_col] / 1000),
                 xytext=(3, 3), textcoords='offset points',
                 fontsize=6.5, color='#333', zorder=4)

m, b_int = np.polyfit(x2, y2, 1)
xr = np.linspace(x2.min() - 0.5, x2.max() + 0.5, 200)
ax2.plot(xr, m * xr + b_int, 'k-', linewidth=1.8, alpha=0.7)

p2_str = f'{p2:.4f}' if p2 >= 0.0001 else f'{p2:.2e}'
ax2.set_xlabel('消費支出（万円/月）', fontsize=12)
ax2.set_ylabel('保健医療費（千円/月）', fontsize=12)
ax2.set_title(f'図2：消費支出 vs 保健医療費（2023年、47都道府県）\n'
              f'Pearson r = {r2:.3f},  p = {p2_str}',
              fontsize=13, fontweight='bold', pad=12)

legend_elem2 = [Patch(facecolor=region_colors[r], label=r, alpha=0.85)
                for r in region_order]
ax2.legend(handles=legend_elem2, fontsize=9, loc='upper left',
           framealpha=0.85, ncol=2)
ax2.grid(alpha=0.25, linewidth=0.8)
ax2.spines['top'].set_visible(False)
ax2.spines['right'].set_visible(False)

plt.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2018_H4_fig2.png'), bbox_inches='tight')
plt.close(fig2)
print(f"fig2 saved.  r={r2:.3f}, p={p2_str}")

# ════════════════════════════════════════════════════════════════════
# Figure 3: 地域ブロック別 保健医療費 箱ひげ図（Kruskal-Wallis検定）
# ════════════════════════════════════════════════════════════════════
groups = [df_cross[df_cross['region'] == r][y_col].dropna().values
          for r in region_order]
groups_nonempty = [g for g in groups if len(g) > 0]
kw_stat, kw_p = kruskal(*groups_nonempty)
kw_p_str = f'{kw_p:.4f}' if kw_p >= 0.0001 else f'{kw_p:.2e}'
print(f"\nKruskal-Wallis: H={kw_stat:.3f}, p={kw_p_str}")

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

bp_data = [g / 1000 for g in groups]
bp = ax3.boxplot(bp_data, patch_artist=True, notch=False,
                 medianprops=dict(color='black', linewidth=2),
                 whiskerprops=dict(linewidth=1.4),
                 capprops=dict(linewidth=1.4),
                 flierprops=dict(marker='o', markersize=5, alpha=0.6))

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

ax3.set_xticks(range(1, len(region_order) + 1))
ax3.set_xticklabels(region_order, rotation=20, ha='right', fontsize=10)
ax3.set_ylabel('保健医療費（千円/月）', fontsize=12)
ax3.set_title(f'図3：地域ブロック別 保健医療費の分布（2023年）\n'
              f'Kruskal-Wallis検定: H = {kw_stat:.2f},  p = {kw_p_str}',
              fontsize=13, fontweight='bold', pad=12)
ax3.grid(axis='y', alpha=0.3, linewidth=0.8)
ax3.spines['top'].set_visible(False)
ax3.spines['right'].set_visible(False)

# n数表示（ylimを先に取得）
y3_lim = ax3.get_ylim()
for i, (region, grp) in enumerate(zip(region_order, groups)):
    ax3.text(i + 1, y3_lim[0] + (y3_lim[1] - y3_lim[0]) * 0.01,
             f'n={len(grp)}', ha='center', va='bottom', fontsize=9, color='#555')

plt.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2018_H4_fig3.png'), bbox_inches='tight')
plt.close(fig3)
print("fig3 saved.")

# ════════════════════════════════════════════════════════════════════
# Figure 4: 重回帰分析 – 標準化偏回帰係数（横棒グラフ + 95%CI）
# ════════════════════════════════════════════════════════════════════
X_vars = ['消費支出_万', '教養娯楽費_千', '高齢化率', '光熱水道費_千']
var_labels = {
    '消費支出_万':   '消費支出（万円/月）\n〈所得水準の代理〉',
    '教養娯楽費_千': '教養娯楽費（千円/月）\n〈余暇・スポーツ支出の代理〉',
    '高齢化率':      '高齢化率（%）\n〈年齢構成・体力水準の代理〉',
    '光熱水道費_千': '光熱・水道費（千円/月）\n〈生活様式の指標〉',
}

df_reg = df_cross[X_vars + [y_col]].dropna().copy()

# 標準化
for col in X_vars + [y_col]:
    df_reg[col + '_z'] = (df_reg[col] - df_reg[col].mean()) / df_reg[col].std()

X_std = sm.add_constant(df_reg[[v + '_z' for v in X_vars]])
y_std = df_reg[y_col + '_z']
ols_result = sm.OLS(y_std, X_std).fit()

print("\n=== OLS重回帰結果（標準化変数） ===")
print(ols_result.summary())

coefs    = ols_result.params[1:]
ci       = ols_result.conf_int().iloc[1:]
pvals    = ols_result.pvalues[1:]
r2_val   = ols_result.rsquared
adj_r2   = ols_result.rsquared_adj

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

labels_plot = [var_labels[v] for v in X_vars]
coef_vals   = coefs.values
ci_lo       = ci.iloc[:, 0].values
ci_hi       = ci.iloc[:, 1].values
xerr_lo     = coef_vals - ci_lo
xerr_hi     = ci_hi - coef_vals

colors4 = ['#e05c5c' if c > 0 else '#4e9af1' for c in coef_vals]

fig4, ax4 = plt.subplots(figsize=(10, 5.5))
y4_pos = np.arange(len(X_vars))

ax4.barh(y4_pos, coef_vals, xerr=[xerr_lo, xerr_hi],
         color=colors4, alpha=0.82, edgecolor='white',
         error_kw=dict(ecolor='#333', linewidth=1.5, capsize=5), height=0.55)
ax4.axvline(0, color='black', linewidth=1.2, linestyle='-')

for i, (c, p) in enumerate(zip(coef_vals, pvals.values)):
    star = sig_star(p)
    offset = 0.03 if c >= 0 else -0.03
    ha_a = 'left' if c >= 0 else 'right'
    color_s = '#c0392b' if star != 'n.s.' else '#888'
    ax4.text(c + offset, i, star, va='center', ha=ha_a,
             fontsize=12, fontweight='bold', color=color_s)

ax4.set_yticks(y4_pos)
ax4.set_yticklabels(labels_plot, fontsize=10)
ax4.set_xlabel('標準化偏回帰係数（β）', fontsize=12)
ax4.set_title(f'図4：重回帰分析 — 標準化偏回帰係数と95%信頼区間（2023年）\n'
              f'R² = {r2_val:.3f},  Adj.R² = {adj_r2:.3f}　'
              f'* p<0.05  ** p<0.01  *** p<0.001',
              fontsize=12, fontweight='bold', pad=12)
ax4.grid(axis='x', alpha=0.3, linewidth=0.8)
ax4.spines['top'].set_visible(False)
ax4.spines['right'].set_visible(False)

plt.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2018_H4_fig4.png'), bbox_inches='tight')
plt.close(fig4)
print("fig4 saved.")

# ── 最終サマリー出力 ──────────────────────────────────────────────
print("\n=== 分析結果サマリー ===")
print(f"分析年度: {latest_yr}年  N={len(df_cross)}")
print(f"目的変数: 保健医療費（二人以上の世帯）[代理変数]")
print(f"消費支出との相関（Pearson r）: r={r2:.3f}, p={p2_str}")
print(f"Kruskal-Wallis（地域別）: H={kw_stat:.3f}, p={kw_p_str}")
print(f"重回帰 R²={r2_val:.3f}, Adj.R²={adj_r2:.3f}")
print("\n標準化偏回帰係数:")
for v, c, p in zip(X_vars, coef_vals, pvals.values):
    short_label = var_labels[v].split('\n')[0]
    print(f"  {short_label:25s}: β={c:+.4f}  {sig_star(p)}  (p={p:.4f})")

print("\nDONE: 2018_H4_katsuyo")