"""
2021_H2 分析スクリプト: 地域別健康格差と社会経済要因の関係
=============================================================
【論文】地域別健康格差と社会経済要因の関係
        2021年度（令和3年度）統計データ分析コンペティション 優秀賞（高校生の部）

【手法】
  1. 散布図: 高齢化率 vs 保健医療費（地域別色分け・都道府県ラベル・回帰線）
  2. 相関ヒートマップ: 健康・社会経済指標のPearson相関行列
  3. 重回帰分析（OLS）: 保健医療費を目的変数とした標準化偏回帰係数プロット
  4. Ward法デンドログラム: 47都道府県の健康特性クラスタリング

【入力】
  data/raw/SSDSE-B-2026.csv（2022年断面データ、47都道府県）

【出力】
  html/figures/2021_H2_fig1.png ... 散布図（高齢化率 vs 保健医療費）
  html/figures/2021_H2_fig2.png ... 相関ヒートマップ
  html/figures/2021_H2_fig3.png ... 標準化偏回帰係数プロット
  html/figures/2021_H2_fig4.png ... Ward法デンドログラム

【注意】
  合成データ（np.random.seed等）は一切使用しない。実データのみ使用。
"""

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


import os
import warnings
import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import statsmodels.api as sm
from scipy import stats
from scipy.cluster.hierarchy import dendrogram, linkage
from sklearn.preprocessing import StandardScaler

warnings.filterwarnings('ignore')

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)

# 2022年断面データ（47都道府県）
df = df_b[df_b['年度'] == 2022].copy().reset_index(drop=True)
assert len(df) == 47, f"Expected 47 prefectures, got {len(df)}"

print(f"=== データ読み込み完了 ===")
print(f"  都道府県数: {len(df)}")
print(f"  年度: {df['年度'].unique()}")

# ── 変数の構築 ─────────────────────────────────────────────────────
# 人口関連
pop = df['総人口'].values.astype(float)

# 健康指標
# (1) 保健医療費（二人以上の世帯）[円/月] → 家計支出ベース
health_spending = df['保健医療費（二人以上の世帯）'].values.astype(float)

# (2) 保健医療費の消費支出比率 [%]
total_spending = df['消費支出（二人以上の世帯）'].values.astype(float)
health_ratio = health_spending / total_spending * 100

# (3) 医療施設数（一般病院数 + 一般診療所数）/ 人口10万対
hospital_cnt = df['一般病院数'].values.astype(float)
clinic_cnt   = df['一般診療所数'].values.astype(float)
medical_per10k = (hospital_cnt + clinic_cnt) / pop * 100000

# 社会経済指標
# (4) 高齢化率 [%]
elderly_pop = df['65歳以上人口'].values.astype(float)
aging_rate = elderly_pop / pop * 100

# (5) 消費支出水準（所得代理変数）[円/月] → 所得の代理として消費支出を使用
# SSDSE-Bには現金給与総額がないため消費支出を使用
income_proxy = total_spending  # 消費支出（所得代理）

# (6) 人口密度: 一般病院数を面積代理として使わず、人口絶対数で代用
# SSDSE-Bに面積変数がないため、都道府県の総人口を用いて対数変換で密度感を表現
# 転入者数（移動者）を用いた都市化指標（流入多いほど都市）
inflow = df['転入者数（日本人移動者）'].values.astype(float)
urban_idx = inflow / pop * 1000   # 転入率（都市化代理）[‰]

# データフレームに格納
df['高齢化率'] = aging_rate
df['保健医療費'] = health_spending
df['保健医療費率'] = health_ratio
df['医療施設数10万対'] = medical_per10k
df['消費支出'] = income_proxy
df['転入率'] = urban_idx

pref_names = df['都道府県'].values

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

df['地域区分'] = df['都道府県'].map(region_map)

# ════════════════════════════════════════════════════════════════════
# ■ Figure 1: 散布図（高齢化率 vs 保健医療費、地域色分け・ラベル・回帰線）
# ════════════════════════════════════════════════════════════════════
fig, ax = plt.subplots(figsize=(12, 8))

x = aging_rate
y = health_spending

# 地域ごとに色付きでプロット
for region, color in region_colors.items():
    mask = df['地域区分'] == region
    ax.scatter(x[mask], y[mask], color=color, s=70, alpha=0.85,
               label=region, zorder=3, edgecolors='white', linewidths=0.5)

# 都道府県ラベル
for i, pref in enumerate(pref_names):
    short = pref.replace('県', '').replace('府', '').replace('都', '').replace('道', '')
    ax.annotate(short, (x[i], y[i]),
                xytext=(3, 3), textcoords='offset points',
                fontsize=6.5, alpha=0.85, zorder=4)

# 回帰線
slope, intercept, r_val, p_val, se = stats.linregress(x, y)
x_line = np.linspace(x.min(), x.max(), 200)
y_line = slope * x_line + intercept
ax.plot(x_line, y_line, color='#333333', lw=2, ls='--', alpha=0.8,
        label=f'回帰線  r={r_val:.3f}  (p={p_val:.4f})')

ax.set_xlabel('高齢化率（65歳以上人口比率）[%]', fontsize=12)
ax.set_ylabel('保健医療費（二人以上世帯）[円/月]', fontsize=12)
ax.set_title('高齢化率と保健医療費の関係\n（2022年・47都道府県）',
             fontsize=14, fontweight='bold')
ax.legend(fontsize=9, loc='upper left', framealpha=0.9)
ax.grid(True, alpha=0.3, lw=0.7)

plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2021_H2_fig1.png'), bbox_inches='tight')
plt.close()
print("Figure 1 saved: 散布図（高齢化率 vs 保健医療費）")

print(f"\n  高齢化率 vs 保健医療費: r={r_val:.3f}, p={p_val:.4f}")
print(f"  回帰式: y = {slope:.2f}x + {intercept:.2f}")

# ════════════════════════════════════════════════════════════════════
# ■ 相関行列の計算（Figure 2用）
# ════════════════════════════════════════════════════════════════════
var_labels = [
    '保健医療費', '保健医療費率',
    '医療施設数\n10万対', '高齢化率',
    '消費支出', '転入率'
]
var_data = np.column_stack([
    health_spending, health_ratio,
    medical_per10k, aging_rate,
    income_proxy, urban_idx
])

corr_matrix = np.corrcoef(var_data, rowvar=False)

# p値行列の計算
n = len(health_spending)
pval_matrix = np.ones_like(corr_matrix)
for i in range(len(var_labels)):
    for j in range(len(var_labels)):
        if i != j:
            _, pval_matrix[i, j] = stats.pearsonr(var_data[:, i], var_data[:, j])

print("\n=== 相関行列 ===")
for i, li in enumerate(var_labels):
    for j, lj in enumerate(var_labels):
        if i < j:
            label_i = li.replace('\n', '')
            label_j = lj.replace('\n', '')
            sig = '***' if pval_matrix[i,j]<0.001 else ('**' if pval_matrix[i,j]<0.01 else ('*' if pval_matrix[i,j]<0.05 else ''))
            print(f"  {label_i} vs {label_j}: r={corr_matrix[i,j]:.3f} {sig}")

# ════════════════════════════════════════════════════════════════════
# ■ Figure 2: 相関ヒートマップ
# ════════════════════════════════════════════════════════════════════
fig, ax = plt.subplots(figsize=(9, 7))

im = ax.imshow(corr_matrix, cmap='RdBu_r', vmin=-1, vmax=1, aspect='auto')

ax.set_xticks(range(len(var_labels)))
ax.set_xticklabels(var_labels, fontsize=10)
ax.set_yticks(range(len(var_labels)))
ax.set_yticklabels(var_labels, fontsize=10)

for i in range(len(var_labels)):
    for j in range(len(var_labels)):
        val = corr_matrix[i, j]
        p   = pval_matrix[i, j]
        sig_mark = '***' if p < 0.001 else ('**' if p < 0.01 else ('*' if p < 0.05 else ''))
        text_color = 'white' if abs(val) > 0.5 else 'black'
        ax.text(j, i, f'{val:.2f}{sig_mark}',
                ha='center', va='center', fontsize=9, color=text_color, fontweight='bold')

plt.colorbar(im, ax=ax, label='Pearson相関係数', shrink=0.8)
ax.set_title('健康・社会経済指標のPearson相関行列\n（2022年・47都道府県）'
             '\n※ *p<0.05, **p<0.01, ***p<0.001',
             fontsize=12, fontweight='bold')
plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2021_H2_fig2.png'), bbox_inches='tight')
plt.close()
print("Figure 2 saved: 相関ヒートマップ")

# ════════════════════════════════════════════════════════════════════
# ■ 重回帰分析（OLS）
# ════════════════════════════════════════════════════════════════════
# 目的変数: 保健医療費
# 説明変数: 高齢化率, 消費支出, 転入率, 医療施設数10万対
# 標準化して標準化偏回帰係数を推定

y_raw = health_spending.copy()

X_raw = np.column_stack([
    aging_rate,       # 高齢化率
    income_proxy,     # 消費支出（所得代理）
    urban_idx,        # 転入率（都市化指標）
    medical_per10k,   # 医療施設数10万対
])

feature_names_reg = ['高齢化率', '消費支出\n（所得代理）', '転入率\n（都市化）', '医療施設数\n10万対']

# 標準化
scaler = StandardScaler()
X_std = scaler.fit_transform(X_raw)
y_std = (y_raw - y_raw.mean()) / y_raw.std()

# OLS 推定
X_ols = sm.add_constant(X_std)
res = sm.OLS(y_std, X_ols).fit()

print("\n=== 重回帰分析結果（標準化） ===")
print(res.summary2().tables[1][['Coef.', 'Std.Err.', 't', 'P>|t|']].to_string())
print(f"\n  R² = {res.rsquared:.4f}, 調整済みR² = {res.rsquared_adj:.4f}")
print(f"  F統計量 = {res.fvalue:.3f}, p = {res.f_pvalue:.6f}")

coefs = np.array(res.params[1:])   # 定数項を除く
pvals = np.array(res.pvalues[1:])
cis   = np.array(res.conf_int())[1:]

# ════════════════════════════════════════════════════════════════════
# ■ Figure 3: 標準化偏回帰係数プロット
# ════════════════════════════════════════════════════════════════════
fig, ax = plt.subplots(figsize=(9, 5.5))

n_feat = len(feature_names_reg)
y_pos = range(n_feat)
bar_colors = ['#e05c5c' if p < 0.05 else '#aaaaaa' for p in pvals]

bars = ax.barh(y_pos, coefs,
               color=bar_colors, alpha=0.85,
               height=0.5,
               xerr=[coefs - cis[:, 0], cis[:, 1] - coefs],
               capsize=4, error_kw={'ecolor': '#555', 'lw': 1.5})

ax.axvline(0, color='black', lw=1.2)
ax.set_yticks(y_pos)
ax.set_yticklabels(feature_names_reg, fontsize=11)
ax.set_xlabel('標準化偏回帰係数（β）', fontsize=11)
ax.set_title(f'重回帰分析: 保健医療費の規定要因\n（目的変数：保健医療費, R²={res.rsquared:.3f}, n=47）',
             fontsize=12, fontweight='bold')

# 値ラベル
for i, (coef, pval) in enumerate(zip(coefs, pvals)):
    sig = '**' if pval < 0.01 else ('*' if pval < 0.05 else 'n.s.')
    offset = 0.01
    ha = 'left' if coef >= 0 else 'right'
    xpos = coef + offset if coef >= 0 else coef - offset
    ax.text(xpos, i, f'β={coef:.3f}\n({sig})', va='center', ha=ha,
            fontsize=9, color='#333333')

sig_patch = mpatches.Patch(color='#e05c5c', alpha=0.85, label='p < 0.05（有意）')
ns_patch  = mpatches.Patch(color='#aaaaaa', alpha=0.85, label='p ≥ 0.05（非有意）')
ax.legend(handles=[sig_patch, ns_patch], loc='lower right', fontsize=9)
ax.grid(True, axis='x', alpha=0.3, lw=0.7)

plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2021_H2_fig3.png'), bbox_inches='tight')
plt.close()
print("Figure 3 saved: 標準化偏回帰係数プロット")

# ════════════════════════════════════════════════════════════════════
# ■ Figure 4: Ward法デンドログラム
# ════════════════════════════════════════════════════════════════════
# クラスタリング用変数（健康・医療系指標）
cluster_vars = np.column_stack([
    health_spending / health_spending.std(),    # 保健医療費（標準化）
    health_ratio / health_ratio.std(),          # 保健医療費率（標準化）
    aging_rate / aging_rate.std(),              # 高齢化率（標準化）
    medical_per10k / medical_per10k.std(),      # 医療施設数（標準化）
])

Z = linkage(cluster_vars, method='ward')

fig, ax = plt.subplots(figsize=(14, 6))

# 地域色のリスト（都道府県順）
leaf_colors_dict = {}
for i, pref in enumerate(pref_names):
    region = region_map.get(pref, '九州・沖縄')
    leaf_colors_dict[i] = region_colors[region]

dend = dendrogram(
    Z,
    labels=pref_names,
    leaf_rotation=90,
    leaf_font_size=8,
    color_threshold=Z[-4, 2],
    ax=ax,
    above_threshold_color='#aaaaaa',
)

# 切断線（4クラスター）
threshold = Z[-3, 2]
ax.axhline(threshold, color='#e05c5c', ls='--', lw=2, alpha=0.9,
           label=f'切断点（4クラスター, 距離={threshold:.2f}）')

ax.set_title('Ward法デンドログラム: 都道府県の健康特性クラスタリング\n（2022年・47都道府県, 健康指標4変数）',
             fontsize=13, fontweight='bold')
ax.set_ylabel('距離（Ward法）', fontsize=11)
ax.legend(fontsize=10)

# 地域凡例
legend_handles = [
    mpatches.Patch(color=c, label=r) for r, c in region_colors.items()
]
ax.legend(handles=legend_handles + [
    mpatches.Patch(color='#e05c5c', label=f'切断点（距離={threshold:.2f}）',
                   linestyle='--', fill=False, edgecolor='#e05c5c')
], fontsize=8, loc='upper right', ncol=2)

plt.tight_layout()
plt.savefig(os.path.join(FIG_DIR, '2021_H2_fig4.png'), bbox_inches='tight')
plt.close()
print("Figure 4 saved: Ward法デンドログラム")

# ════════════════════════════════════════════════════════════════════
# ■ 統計サマリー出力
# ════════════════════════════════════════════════════════════════════
print("\n" + "="*60)
print("■ 統計サマリー（2022年・47都道府県）")
print("="*60)
print(f"  保健医療費 平均: {health_spending.mean():.0f} 円/月")
print(f"  保健医療費 最大: {health_spending.max():.0f} ({pref_names[np.argmax(health_spending)]})")
print(f"  保健医療費 最小: {health_spending.min():.0f} ({pref_names[np.argmin(health_spending)]})")
print(f"  保健医療費率 平均: {health_ratio.mean():.2f} %")
print(f"  高齢化率 平均: {aging_rate.mean():.2f} %")
print(f"  高齢化率 最大: {aging_rate.max():.2f} % ({pref_names[np.argmax(aging_rate)]})")
print(f"  高齢化率 最小: {aging_rate.min():.2f} % ({pref_names[np.argmin(aging_rate)]})")
print(f"  医療施設数 平均 (10万対): {medical_per10k.mean():.2f}")
print(f"  消費支出 平均: {income_proxy.mean():.0f} 円/月")
print(f"\n  重回帰分析:")
print(f"    R²={res.rsquared:.4f}, 調整済みR²={res.rsquared_adj:.4f}")
print(f"    F={res.fvalue:.3f}, p={res.f_pvalue:.6f}")
for fname, coef, pval in zip(feature_names_reg, coefs, pvals):
    fname_clean = fname.replace('\n', '')
    sig = '***' if pval<0.001 else ('**' if pval<0.01 else ('*' if pval<0.05 else 'n.s.'))
    print(f"    {fname_clean}: β={coef:.4f}, p={pval:.4f} {sig}")

print("\n分析完了。html/figures/ に4枚の図を保存しました。")