"""
2021年度 統計データ分析コンペティション
審査員奨励賞（高校生の部）
「都道府県別自殺死亡率と社会経済要因の関係分析」

SSDSE-B-2026.csv（都道府県別データ）を使用し、
死亡率（千人当たり）を目的変数として社会経済要因との関係を分析する。
注: SSDSE-Bには自殺死亡率の直接データは含まれないため、
    粗死亡率と高齢化率・雇用環境・保健医療費等との関係を分析する。
"""

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

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年断面データ
df = df_b[df_b['年度'] == 2022].copy()
assert len(df) == 47, f"47都道府県のデータが必要です。実際: {len(df)}"

# ---- 変数の計算 ----
# 目的変数: 粗死亡率（千人当たり）
df['死亡率_千人'] = df['死亡数'] / df['総人口'] * 1000

# 説明変数
# 1. 求人倍率（月間有効求人数 / 月間有効求職者数）
df['求人倍率'] = df['月間有効求人数（一般）'] / df['月間有効求職者数（一般）']

# 2. 高齢化率（65歳以上人口 / 総人口 × 100）
df['高齢化率'] = df['65歳以上人口'] / df['総人口'] * 100

# 3. 保健医療費（二人以上の世帯、円）
df['保健医療費'] = df['保健医療費（二人以上の世帯）']

# 4. 消費支出（生活水準の代理）
df['消費支出'] = df['消費支出（二人以上の世帯）']

# 5. 離婚率（千人当たり、社会的孤立の代理）
df['離婚率_千人'] = df['離婚件数'] / df['総人口'] * 1000

# 6. 出生率（合計特殊出生率）の逆数（社会活力の喪失指標）
df['出生率'] = df['合計特殊出生率']

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

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

# ---- 相関分析（Pearson） ----
print("=" * 60)
print("相関分析（目的変数: 死亡率 千人当たり）")
print("=" * 60)
explan_vars = {
    '求人倍率': '求人倍率',
    '高齢化率': '高齢化率',
    '保健医療費': '保健医療費',
    '消費支出': '消費支出',
    '離婚率_千人': '離婚率_千人',
    '出生率': '出生率'
}
corr_results = {}
for varname, col in explan_vars.items():
    r, p = stats.pearsonr(df[col], df['死亡率_千人'])
    corr_results[varname] = {'r': r, 'p': p}
    sig = '*' if p < 0.05 else ' '
    print(f"  {varname:10s}: r={r:+.4f}, p={p:.4f} {sig}")

# ---- 重回帰分析（OLS） ----
print()
print("=" * 60)
print("重回帰分析（OLS）")
print("=" * 60)
X_vars = ['求人倍率', '高齢化率', '保健医療費', '消費支出', '離婚率_千人', '出生率']
X = df[X_vars].copy()

# 標準化（標準化偏回帰係数用）
X_std = (X - X.mean()) / X.std()
y = df['死亡率_千人']

# モデル推定（非標準化）
X_const = sm.add_constant(X)
model = sm.OLS(y, X_const).fit()
print(model.summary())

# 標準化モデル
X_std_const = sm.add_constant(X_std)
model_std = sm.OLS(y, X_std_const).fit()

print()
print("標準化偏回帰係数:")
for v, coef, p in zip(X_vars, model_std.params[1:], model_std.pvalues[1:]):
    sig = '**' if p < 0.01 else ('*' if p < 0.05 else '')
    print(f"  {v:12s}: β={coef:+.4f}, p={p:.4f} {sig}")
print(f"  R² = {model_std.rsquared:.4f}")
print(f"  Adj.R² = {model_std.rsquared_adj:.4f}")

# ==================================================================
# 図1: 都道府県別死亡率ランキング棒グラフ（地域色分け）
# ==================================================================
df_sorted = df.sort_values('死亡率_千人', ascending=True).reset_index(drop=True)

fig, ax = plt.subplots(figsize=(10, 14))
bars = ax.barh(
    range(len(df_sorted)),
    df_sorted['死亡率_千人'],
    color=df_sorted['地域色'],
    edgecolor='white', linewidth=0.4, height=0.75
)
ax.set_yticks(range(len(df_sorted)))
ax.set_yticklabels(df_sorted['都道府県'], fontsize=8.5)
ax.set_xlabel('死亡率（人口千人当たり、2022年）', fontsize=11)
ax.set_title('都道府県別 粗死亡率ランキング（2022年）\n目的変数：精神的健康の社会経済的代理指標', fontsize=12, fontweight='bold')

# 凡例
from matplotlib.patches import Patch
legend_elements = [Patch(facecolor=c, label=r) for r, c in region_colors.items()]
ax.legend(handles=legend_elements, loc='lower right', fontsize=9, framealpha=0.9)

# 全国平均ライン
mean_val = df['死亡率_千人'].mean()
ax.axvline(mean_val, color='#c0392b', linestyle='--', linewidth=1.2, alpha=0.8)
ax.text(mean_val + 0.05, 1.5, f'全国平均\n{mean_val:.1f}‰', color='#c0392b', fontsize=9, va='bottom')

ax.set_xlim(0, df_sorted['死亡率_千人'].max() * 1.12)
ax.grid(axis='x', alpha=0.3)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
plt.tight_layout()
fig1_path = os.path.join(FIG_DIR, '2021_H5_4_fig1.png')
fig.savefig(fig1_path, bbox_inches='tight')
plt.close()
print(f"\n図1保存: {fig1_path}")

# ==================================================================
# 図2: 散布図 求人倍率 vs 死亡率（地域色分け・47都道府県ラベル・回帰線）
# ==================================================================
fig, ax = plt.subplots(figsize=(11, 8))

for region, color in region_colors.items():
    mask = df['地域区分'] == region
    ax.scatter(df.loc[mask, '求人倍率'], df.loc[mask, '死亡率_千人'],
               color=color, s=60, alpha=0.85, zorder=3, label=region)

# 都道府県名ラベル
for _, row in df.iterrows():
    pref = row['都道府県'].replace('県', '').replace('都', '').replace('道', '').replace('府', '')
    ax.annotate(pref, (row['求人倍率'], row['死亡率_千人']),
                fontsize=6.5, ha='center', va='bottom', xytext=(0, 3),
                textcoords='offset points', alpha=0.85)

# 回帰線
x_reg = df['求人倍率'].values
y_reg = df['死亡率_千人'].values
slope, intercept, r_val, p_val, se = stats.linregress(x_reg, y_reg)
x_line = np.linspace(x_reg.min(), x_reg.max(), 100)
ax.plot(x_line, slope * x_line + intercept, color='#2c3e50',
        linewidth=1.8, linestyle='-', zorder=2, label=f'回帰線 r={r_val:+.3f}')

ax.set_xlabel('求人倍率（月間有効求人数/有効求職者数、2022年）', fontsize=11)
ax.set_ylabel('粗死亡率（千人当たり、2022年）', fontsize=11)
ax.set_title('求人倍率と死亡率の関係（都道府県別、2022年）\n求人倍率が高い地域ほど死亡率が高い傾向（高齢化率と相関）', fontsize=11, fontweight='bold')
ax.legend(loc='upper left', fontsize=9, framealpha=0.9)
ax.text(0.98, 0.04, f'r = {r_val:+.3f}, p = {p_val:.4f}', transform=ax.transAxes,
        ha='right', va='bottom', fontsize=10,
        bbox=dict(boxstyle='round,pad=0.4', facecolor='#fff9c4', edgecolor='#f9a825', alpha=0.9))
ax.grid(alpha=0.3)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
plt.tight_layout()
fig2_path = os.path.join(FIG_DIR, '2021_H5_4_fig2.png')
fig.savefig(fig2_path, bbox_inches='tight')
plt.close()
print(f"図2保存: {fig2_path}")

# ==================================================================
# 図3: 相関行列ヒートマップ
# ==================================================================
var_labels = {
    '死亡率_千人': '死亡率\n(千人当たり)',
    '求人倍率': '求人倍率',
    '高齢化率': '高齢化率\n(%)',
    '保健医療費': '保健医療費\n(円/月)',
    '消費支出': '消費支出\n(円/月)',
    '離婚率_千人': '離婚率\n(千人当たり)',
    '出生率': '合計特殊\n出生率'
}
corr_cols = list(var_labels.keys())
corr_matrix = df[corr_cols].corr()

fig, ax = plt.subplots(figsize=(9, 7))
n = len(corr_cols)
labels = [var_labels[c] for c in corr_cols]

import matplotlib.colors as mcolors
cmap = plt.cm.RdBu_r

im = ax.imshow(corr_matrix.values, cmap=cmap, vmin=-1, vmax=1, aspect='auto')
plt.colorbar(im, ax=ax, shrink=0.8, label='Pearson 相関係数')

ax.set_xticks(range(n))
ax.set_yticks(range(n))
ax.set_xticklabels(labels, fontsize=9, rotation=0, ha='center')
ax.set_yticklabels(labels, fontsize=9)

for i in range(n):
    for j in range(n):
        val = corr_matrix.values[i, j]
        txt_color = 'white' if abs(val) > 0.6 else 'black'
        ax.text(j, i, f'{val:.2f}', ha='center', va='center',
                fontsize=9, color=txt_color, fontweight='bold' if abs(val) > 0.5 else 'normal')

ax.set_title('社会経済変数の相関行列（Pearson相関係数、2022年、47都道府県）',
             fontsize=11, fontweight='bold', pad=15)
ax.spines[:].set_visible(False)
plt.tight_layout()
fig3_path = os.path.join(FIG_DIR, '2021_H5_4_fig3.png')
fig.savefig(fig3_path, bbox_inches='tight')
plt.close()
print(f"図3保存: {fig3_path}")

# ==================================================================
# 図4: 標準化偏回帰係数プロット（横棒グラフ）
# ==================================================================
std_coefs = model_std.params[1:].values
std_pvals = model_std.pvalues[1:].values
std_conf  = model_std.conf_int().iloc[1:].values  # [lower, upper]

# エラーバー（95%CI）
err_lower = std_coefs - std_conf[:, 0]
err_upper = std_conf[:, 1] - std_coefs

# 変数名（表示用）
var_display = ['求人倍率', '高齢化率', '保健医療費\n(円/月)', '消費支出\n(円/月)', '離婚率\n(千人当たり)', '合計特殊\n出生率']

# 係数の絶対値でソート（大きい順）
sort_idx = np.argsort(np.abs(std_coefs))
coefs_s = std_coefs[sort_idx]
pvals_s = std_pvals[sort_idx]
labels_s = [var_display[i] for i in sort_idx]
err_lo_s = err_lower[sort_idx]
err_hi_s = err_upper[sort_idx]

colors_bar = ['#c0392b' if p < 0.05 else '#95a5a6' for p in pvals_s]

fig, ax = plt.subplots(figsize=(9, 5))
bars = ax.barh(range(len(coefs_s)), coefs_s, color=colors_bar,
               height=0.6, edgecolor='white', linewidth=0.5)
ax.errorbar(coefs_s, range(len(coefs_s)),
            xerr=[err_lo_s, err_hi_s],
            fmt='none', color='#2c3e50', capsize=4, linewidth=1.2)

ax.axvline(0, color='#2c3e50', linewidth=1.0, linestyle='--', alpha=0.6)
ax.set_yticks(range(len(labels_s)))
ax.set_yticklabels(labels_s, fontsize=10)
ax.set_xlabel('標準化偏回帰係数（β）', fontsize=11)
ax.set_title(f'重回帰分析：標準化偏回帰係数（95%信頼区間）\n目的変数: 粗死亡率（千人当たり）、R²={model_std.rsquared:.3f}',
             fontsize=11, fontweight='bold')

from matplotlib.patches import Patch
legend_elems = [
    Patch(facecolor='#c0392b', label='有意（p<0.05）'),
    Patch(facecolor='#95a5a6', label='非有意（p≥0.05）')
]
ax.legend(handles=legend_elems, loc='lower right', fontsize=9)

# p値をバーの横に表示
for i, (c, p) in enumerate(zip(coefs_s, pvals_s)):
    sign = '**' if p < 0.01 else ('*' if p < 0.05 else 'ns')
    offset = 0.02 if c >= 0 else -0.02
    ha = 'left' if c >= 0 else 'right'
    ax.text(c + offset, i, sign, ha=ha, va='center', fontsize=9,
            color='#2c3e50', fontweight='bold')

ax.grid(axis='x', alpha=0.3)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
plt.tight_layout()
fig4_path = os.path.join(FIG_DIR, '2021_H5_4_fig4.png')
fig.savefig(fig4_path, bbox_inches='tight')
plt.close()
print(f"図4保存: {fig4_path}")

# ---- 統計サマリーの出力 ----
print()
print("=" * 60)
print("統計サマリー（HTML作成用）")
print("=" * 60)
print(f"分析対象: 47都道府県（2022年度）")
print(f"目的変数: 粗死亡率（千人当たり）")
print(f"  平均: {df['死亡率_千人'].mean():.2f}‰")
print(f"  標準偏差: {df['死亡率_千人'].std():.2f}‰")
print(f"  最大: {df.loc[df['死亡率_千人'].idxmax(), '都道府県']} {df['死亡率_千人'].max():.2f}‰")
print(f"  最小: {df.loc[df['死亡率_千人'].idxmin(), '都道府県']} {df['死亡率_千人'].min():.2f}‰")
print()
print("相関分析:")
for k, v in corr_results.items():
    print(f"  {k}: r={v['r']:+.4f}, p={v['p']:.4f}")
print()
print(f"重回帰分析: R²={model_std.rsquared:.4f}, Adj.R²={model_std.rsquared_adj:.4f}")
print("標準化偏回帰係数:")
for v, coef, pval in zip(X_vars, model_std.params[1:], model_std.pvalues[1:]):
    sig = '**' if pval < 0.01 else ('*' if pval < 0.05 else 'ns')
    print(f"  {v}: β={coef:+.4f}, p={pval:.4f} {sig}")
print()
print("完了")