"""
2019_H4_katsuyo.py
都道府県別社会保障給付費と生活保護受給率の規定要因分析
統計活用奨励賞（高校生部門）2019年度

【プロキシ変数の説明】
SSDSE-B-2026には社会保障給付費や生活保護受給率の直接データが含まれないため、
以下の代理変数を使用する：
  目的変数  : 保健医療費（二人以上の世帯）[円/月] = 家計の医療・保健支出（社会保障需要の代理）
  説明変数1 : 高齢化率 = 65歳以上人口 / 総人口 × 100 [%]（高齢化の代理）
  説明変数2 : 有効求人倍率 = 月間有効求人数 / 月間有効求職者数（労働市場の逼迫度）
  説明変数3 : 病院数（人口10万対）= 一般病院数 / 総人口 × 100000（医療インフラ）
  説明変数4 : 消費支出（万円/月）= 消費支出（二人以上の世帯）/ 10000（所得水準の代理）
"""

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

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()
print("=== df_b.head(3) ===")
print(df_b.head(3))
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['月間有効求人数（一般）'] / df_cross['月間有効求職者数（一般）']
df_cross['病院数_per10万'] = df_cross['一般病院数'] / df_cross['総人口'] * 100000
df_cross['消費支出_万'] = df_cross['消費支出（二人以上の世帯）'] / 10000

# 目的変数：保健医療費（代理変数）
y_col = '保健医療費（二人以上の世帯）'
df_cross = df_cross.dropna(subset=[y_col, '高齢化率', '有効求人倍率', '病院数_per10万', '消費支出_万', 'region'])

print(f"分析サンプル数: {len(df_cross)} 都道府県（{latest_yr}年）")
print()
print("=== 基本統計量 ===")
stats_cols = [y_col, '高齢化率', '有効求人倍率', '病院数_per10万', '消費支出_万']
print(df_cross[stats_cols].describe().round(2))
print()

# ── 全国平均 ──────────────────────────────────────────────────────
national_avg = df_cross[y_col].mean()
print(f"保健医療費（全国平均）: {national_avg:.0f} 円/月")
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']]
bars = 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.legend(fontsize=11, loc='upper right')
ax1.grid(axis='y', alpha=0.3, linewidth=0.8)
ax1.spines['top'].set_visible(False)
ax1.spines['right'].set_visible(False)

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

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

# ════════════════════════════════════════════════════════════════════
# Figure 2: 高齢化率 vs 保健医療費の散布図（回帰直線つき）
# ════════════════════════════════════════════════════════════════════
x2 = df_cross['高齢化率'].values
y2 = df_cross[y_col].values / 1000  # 千円
r2, p2 = stats.pearsonr(x2, y2)

fig2, ax2 = plt.subplots(figsize=(9, 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 = np.polyfit(x2, y2, 1)
xr = np.linspace(x2.min() - 0.5, x2.max() + 0.5, 200)
ax2.plot(xr, m * xr + b, 'k-', linewidth=1.8, alpha=0.7, label=f'回帰直線')

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'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, '2019_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)
print(f"\nKruskal-Wallis: H={kw_stat:.3f}, p={kw_p:.4f}")

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)
kw_p_str = f'{kw_p:.4f}' if kw_p >= 0.0001 else f'{kw_p:.2e}'
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数表示
for i, (region, grp) in enumerate(zip(region_order, groups)):
    ax3.text(i + 1, ax3.get_ylim()[0] + 0.2,
             f'n={len(grp)}', ha='center', va='bottom', fontsize=9, color='#555')

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

# ════════════════════════════════════════════════════════════════════
# Figure 4: 重回帰分析 – 標準化偏回帰係数（横棒グラフ）
# ════════════════════════════════════════════════════════════════════
# 変数の標準化
X_vars = ['高齢化率', '有効求人倍率', '病院数_per10万', '消費支出_万']
var_labels = {
    '高齢化率': '高齢化率（%）',
    '有効求人倍率': '有効求人倍率',
    '病院数_per10万': '病院数（人口10万対）',
    '消費支出_万': '消費支出（万円/月）'
}

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())

# 標準化係数・信頼区間・p値
coefs    = ols_result.params[1:]      # 定数項除く
ci       = ols_result.conf_int().iloc[1:]  # 95%CI
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=(9, 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.025 if c >= 0 else -0.025
    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=11)
ax4.set_xlabel('標準化偏回帰係数（β）', fontsize=12)
ax4.set_title(f'図4：重回帰分析 — 標準化偏回帰係数（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, '2019_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"高齢化率との相関: 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):
    print(f"  {var_labels[v]:20s}: β={c:+.4f}  {sig_star(p)}  (p={p:.4f})")

print("\nDONE: 2019_H4_katsuyo")