"""
2019_H5_1_shorei.py
審査員奨励賞（高校生部門）
「都道府県別交通事故・安全指標の規定要因：経済・都市化・医療アクセスの複合分析」

代理変数について:
  SSDSE-B-2026 には交通事故件数の直接データがないため、
  安全・健康指標の代理として「死亡率（人口千人あたり死亡数）」を目的変数として使用。
  死亡率は高齢化・都市化・医療アクセス・経済水準を反映する複合指標であり、
  元論文の「経済・都市化・医療アクセス」の複合分析の枠組みに適合する。

説明変数（元論文の枠組みを踏襲）:
  - 有効求人倍率：月間有効求人数 / 月間有効求職者数（経済指標）
  - 高齢化率：65歳以上人口 / 総人口 × 100（%）
  - 人口密度：総人口 / 可住地面積（人/km²）  ← SSDSE-E より
  - 消費支出：二人以上の世帯の月額消費支出（万円・経済水準の代理）
  - 人口10万対病院数：一般病院数 / 総人口 × 100,000（医療アクセス）

使用データ:
  SSDSE-B-2026.csv（2023年度 47都道府県・最新横断面）
  SSDSE-E-2026.csv（47都道府県・横断面）
"""

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

# ============================================================
# 1. データ読み込み
# ============================================================
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("SSDSE-B columns:", df_b.columns.tolist())
print("Years available:", sorted(df_b['年度'].unique()))

# 最新年（2023年）データ抽出
latest_year = df_b['年度'].max()
d = df_b[df_b['年度'] == latest_year].copy().reset_index(drop=True)
print(f"\n{latest_year}年 都道府県数: {len(d)}")

# SSDSE-E（横断面・都道府県別）
df_e_raw = pd.read_csv(DATA_E, encoding='cp932', header=0)
col_names_e = df_e_raw.iloc[1].tolist()
df_e = df_e_raw.iloc[2:].copy()
df_e.columns = col_names_e
df_e = df_e[df_e['地域コード'].str.match(r'R\d{5}', na=False)].copy()
df_e = df_e[df_e['地域コード'] != 'R00000'].copy().reset_index(drop=True)

# ============================================================
# 2. 都道府県名の正規化
# ============================================================
def normalize_pref(name):
    name = str(name).strip()
    if name == '北海道':
        return '北海道'
    for s in ['都', '道', '府', '県']:
        if name.endswith(s):
            return name[:-1]
    return name

d['pref_key'] = d['都道府県'].apply(normalize_pref)
df_e['pref_key'] = df_e['都道府県'].apply(normalize_pref)

print("\nSSDSE-B 都道府県 sample:", d['pref_key'].tolist()[:5])
print("SSDSE-E 都道府県 sample:", df_e['pref_key'].tolist()[:5])

# ============================================================
# 3. 変数の作成（SSDSE-B）
# ============================================================
d['総人口_num']      = pd.to_numeric(d['総人口'], errors='coerce')
d['65歳以上_num']    = pd.to_numeric(d['65歳以上人口'], errors='coerce')
d['死亡数_num']      = pd.to_numeric(d['死亡数'], errors='coerce')
d['有効求人数_num']  = pd.to_numeric(d['月間有効求人数（一般）'], errors='coerce')
d['有効求職者数_num'] = pd.to_numeric(d['月間有効求職者数（一般）'], errors='coerce')
d['病院数_num']      = pd.to_numeric(d['一般病院数'], errors='coerce')
d['消費支出_num']    = pd.to_numeric(d['消費支出（二人以上の世帯）'], errors='coerce')

# 目的変数: 死亡率（人口千人あたり）
d['死亡率'] = d['死亡数_num'] / d['総人口_num'] * 1000

# 説明変数1: 高齢化率（%）
d['高齢化率'] = d['65歳以上_num'] / d['総人口_num'] * 100

# 説明変数2: 有効求人倍率
d['有効求人倍率'] = d['有効求人数_num'] / d['有効求職者数_num']

# 説明変数3: 人口10万対病院数（医療アクセス）
d['人口10万対病院数'] = d['病院数_num'] / d['総人口_num'] * 100000

# 説明変数4: 消費支出（万円換算）
d['消費支出_万円'] = d['消費支出_num'] / 10000

# ============================================================
# 4. 変数の作成（SSDSE-E）
# ============================================================
df_e['可住地面積_num'] = pd.to_numeric(df_e['可住地面積'], errors='coerce')
df_e['総人口_e_num']  = pd.to_numeric(df_e['総人口'], errors='coerce')

# 説明変数5: 人口密度（人/km²）　可住地面積は ha 単位 → / 100 で km²
df_e['人口密度'] = df_e['総人口_e_num'] / (df_e['可住地面積_num'] / 100)

# ============================================================
# 5. データ結合・整形
# ============================================================
df_merge = d[['pref_key', '都道府県', '死亡率', '高齢化率', '有効求人倍率',
              '人口10万対病院数', '消費支出_万円']].merge(
    df_e[['pref_key', '人口密度']],
    on='pref_key', how='inner'
)

print(f"\nマージ後データ件数: {len(df_merge)}")
print(df_merge[['pref_key', '死亡率', '高齢化率', '有効求人倍率',
                '人口密度', '消費支出_万円', '人口10万対病院数']].head(10).to_string())

# 欠損値確認・除去
print("\n欠損値:")
print(df_merge.isnull().sum())
df_merge = df_merge.dropna().reset_index(drop=True)
print(f"欠損除去後: {len(df_merge)} 都道府県")

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

df_merge['地域'] = df_merge['pref_key'].map(region_map)
df_merge['地域色'] = df_merge['地域'].map(region_colors)

# ============================================================
# 7. 記述統計・相関分析
# ============================================================
vars_analysis = ['死亡率', '高齢化率', '有効求人倍率', '人口密度', '消費支出_万円', '人口10万対病院数']
var_labels = {
    '死亡率':         '死亡率\n(‰)',
    '高齢化率':       '高齢化率\n(%)',
    '有効求人倍率':    '有効求人\n倍率',
    '人口密度':       '人口密度\n(人/km²)',
    '消費支出_万円':  '消費支出\n(万円/月)',
    '人口10万対病院数': '病院数\n(10万対)',
}

print("\n=== 記述統計 ===")
print(df_merge[vars_analysis].describe().round(3).to_string())

# Pearson相関行列
print("\n=== Pearson相関行列 ===")
corr_mat = df_merge[vars_analysis].corr()
print(corr_mat.round(3).to_string())

# ============================================================
# 8. 重回帰分析（OLS）
# ============================================================
y = df_merge['死亡率'].values
X_vars = ['高齢化率', '有効求人倍率', '人口密度', '消費支出_万円', '人口10万対病院数']
X = df_merge[X_vars].values

# 標準化（標準化偏回帰係数のため）
y_std = (y - y.mean()) / y.std()
X_std = (X - X.mean(axis=0)) / X.std(axis=0)

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

print("\n=== 重回帰分析（標準化） ===")
print(model_std.summary())

# 非標準化モデル（Cook's distance 等）
X_const = sm.add_constant(X)
model = sm.OLS(y, X_const).fit()

print("\n=== 重回帰分析（非標準化） ===")
print(model.summary())

# 標準化偏回帰係数の抽出
betas      = model_std.params[1:]
beta_se    = model_std.bse[1:]
beta_pvals = model_std.pvalues[1:]

print("\n=== 標準化偏回帰係数 ===")
for i, v in enumerate(X_vars):
    sig = ('***' if beta_pvals[i] < 0.001 else
           '**'  if beta_pvals[i] < 0.01  else
           '*'   if beta_pvals[i] < 0.05  else 'n.s.')
    print(f"  {v}: β={betas[i]:.4f}, SE={beta_se[i]:.4f}, p={beta_pvals[i]:.4f} {sig}")

print(f"\nR² = {model_std.rsquared:.4f}")
print(f"Adj.R² = {model_std.rsquared_adj:.4f}")

# Cook's distance
influence = model.get_influence()
cooks_d = influence.cooks_distance[0]
print("\n=== Cook's Distance（上位5） ===")
top5_cooks = pd.DataFrame({'都道府県': df_merge['pref_key'].values, 'Cooks_D': cooks_d})
print(top5_cooks.nlargest(5, 'Cooks_D').to_string(index=False))

# 上位・下位 5 都道府県（死亡率）
print("\n=== 死亡率 上位5都道府県 ===")
print(df_merge[['pref_key', '死亡率']].nlargest(5, '死亡率').to_string(index=False))
print("\n=== 死亡率 下位5都道府県 ===")
print(df_merge[['pref_key', '死亡率']].nsmallest(5, '死亡率').to_string(index=False))

# ============================================================
# 9. 図1: 都道府県別死亡率ランキング棒グラフ（地域色分け・全国平均線）
# ============================================================
df_sorted = df_merge.sort_values('死亡率', ascending=False).reset_index(drop=True)
colors_bar = [region_colors[r] for r in df_sorted['地域']]
mean_val = df_merge['死亡率'].mean()

fig1, ax1 = plt.subplots(figsize=(15, 7))
ax1.bar(range(len(df_sorted)), df_sorted['死亡率'],
        color=colors_bar, width=0.7, edgecolor='white', linewidth=0.4)
ax1.axhline(mean_val, color='#333', linewidth=2.0, linestyle='--', zorder=5)

ax1.set_xticks(range(len(df_sorted)))
ax1.set_xticklabels(df_sorted['pref_key'], rotation=90, fontsize=8.5)
ax1.set_ylabel('死亡率（人口千人あたり）', fontsize=11)
ax1.set_title(f'都道府県別 死亡率ランキング（{latest_year}年・47都道府県）\n'
              '高齢化・医療アクセス・経済水準の複合指標', fontsize=13, fontweight='bold', pad=12)

ax1.set_xlim(-0.7, len(df_sorted) - 0.3)
ax1.set_ylim(0, df_sorted['死亡率'].max() * 1.15)
ax1.spines['top'].set_visible(False)
ax1.spines['right'].set_visible(False)
ax1.yaxis.grid(True, alpha=0.3)
ax1.set_axisbelow(True)

# 全国平均ラベル
ax1.text(len(df_sorted) - 1, mean_val + 0.2, f'全国平均\n{mean_val:.2f}‰',
         ha='right', va='bottom', fontsize=9, color='#333')

# 地域凡例
legend_els = [Patch(facecolor=c, label=r) for r, c in region_colors.items()]
legend_els.append(plt.Line2D([0], [0], color='#333', linewidth=2.0,
                              linestyle='--', label=f'全国平均 {mean_val:.2f}‰'))
ax1.legend(handles=legend_els, fontsize=9, loc='upper right', framealpha=0.9)

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

# ============================================================
# 10. 図2: 有効求人倍率 vs 死亡率 散布図（都道府県ラベル付き・回帰直線）
# ============================================================
fig2, ax2 = plt.subplots(figsize=(11, 7))

for region, color in region_colors.items():
    mask = df_merge['地域'] == region
    ax2.scatter(df_merge.loc[mask, '有効求人倍率'],
                df_merge.loc[mask, '死亡率'],
                color=color, s=65, alpha=0.85,
                edgecolors='white', linewidths=0.5, label=region, zorder=3)

# 都道府県ラベル
for _, row in df_merge.iterrows():
    ax2.annotate(row['pref_key'],
                 (row['有効求人倍率'], row['死亡率']),
                 fontsize=6.5, ha='left', va='bottom',
                 xytext=(2, 2), textcoords='offset points', color='#333')

# 回帰直線
x_sc = df_merge['有効求人倍率'].values
y_sc = df_merge['死亡率'].values
r_val, p_val = stats.pearsonr(x_sc, y_sc)
slope, intercept, _, _, _ = stats.linregress(x_sc, y_sc)
x_line = np.linspace(x_sc.min() - 0.05, x_sc.max() + 0.05, 200)
ax2.plot(x_line, slope * x_line + intercept,
         color='#333', linewidth=1.8, linestyle='--', zorder=2)

p_str = f'p = {p_val:.4f}' if p_val >= 0.001 else 'p < 0.001'
ax2.text(0.05, 0.95, f'r = {r_val:.3f}\n{p_str}',
         transform=ax2.transAxes, fontsize=11, va='top', ha='left',
         bbox=dict(boxstyle='round,pad=0.4', facecolor='white',
                   edgecolor='#aaa', alpha=0.9))

ax2.set_xlabel('有効求人倍率（月間有効求人数 / 月間有効求職者数）', fontsize=11)
ax2.set_ylabel('死亡率（人口千人あたり）', fontsize=11)
ax2.set_title(f'有効求人倍率 vs 死亡率（{latest_year}年・47都道府県）\n'
              '経済的活力と健康・安全指標の地域差', fontsize=13, fontweight='bold', pad=12)

legend_els2 = [Patch(facecolor=c, label=r) for r, c in region_colors.items()]
ax2.legend(handles=legend_els2, fontsize=9, loc='lower right', framealpha=0.9)
ax2.spines['top'].set_visible(False)
ax2.spines['right'].set_visible(False)

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

# ============================================================
# 11. 図3: 相関ヒートマップ（Pearson相関行列・有意マーク付き）
# ============================================================
from matplotlib.colors import TwoSlopeNorm

n_v = len(vars_analysis)
pval_mat = np.ones((n_v, n_v))
for i in range(n_v):
    for j in range(n_v):
        if i != j:
            _, p_ij = stats.pearsonr(df_merge[vars_analysis[i]],
                                     df_merge[vars_analysis[j]])
            pval_mat[i, j] = p_ij

corr_arr = corr_mat.values
norm = TwoSlopeNorm(vmin=-1, vcenter=0, vmax=1)

fig3, ax3 = plt.subplots(figsize=(9, 7))
im = ax3.imshow(corr_arr, cmap='RdBu_r', norm=norm, aspect='auto')

tick_labels = [var_labels[v] for v in vars_analysis]
ax3.set_xticks(range(n_v))
ax3.set_yticks(range(n_v))
ax3.set_xticklabels(tick_labels, fontsize=9)
ax3.set_yticklabels(tick_labels, fontsize=9)

for i in range(n_v):
    for j in range(n_v):
        r_ij  = corr_arr[i, j]
        p_ij  = pval_mat[i, j]
        sig_mark = ('***' if p_ij < 0.001 else
                    '**'  if p_ij < 0.01  else
                    '*'   if p_ij < 0.05  else '')
        text_color = 'white' if abs(r_ij) > 0.6 else 'black'
        ax3.text(j, i, f'{r_ij:.2f}{sig_mark}',
                 ha='center', va='center', fontsize=9,
                 color=text_color,
                 fontweight='bold' if sig_mark else 'normal')

plt.colorbar(im, ax=ax3, shrink=0.8, label='Pearson r')
ax3.set_title(f'変数間Pearson相関行列（{latest_year}年・47都道府県）\n'
              '*p<0.05  **p<0.01  ***p<0.001',
              fontsize=12, fontweight='bold', pad=14)
plt.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2019_H5_1_fig3.png'), bbox_inches='tight')
plt.close(fig3)
print("Fig3 saved.")

# ============================================================
# 12. 図4: 重回帰の標準化偏回帰係数（横棒グラフ + 95%CI + 有意マーク）
# ============================================================
x_labels = ['高齢化率', '有効求人倍率', '人口密度\n(都市化)', '消費支出\n(万円/月)', '病院数\n(10万対)']
ci_95 = 1.96 * beta_se

fig4, ax4 = plt.subplots(figsize=(9, 6))
y_pos = np.arange(len(X_vars))[::-1]

bar_colors = []
for i, pv in enumerate(beta_pvals):
    if pv < 0.05:
        bar_colors.append('#e05c5c' if betas[i] > 0 else '#4e9af1')
    else:
        bar_colors.append('#cccccc')

ax4.barh(y_pos, betas, xerr=ci_95,
         color=bar_colors, height=0.55,
         edgecolor='white', linewidth=0.5,
         error_kw=dict(ecolor='#555', capsize=4, lw=1.5))
ax4.axvline(0, color='#333', linewidth=1.2)

ax4.set_yticks(y_pos)
ax4.set_yticklabels(x_labels, fontsize=10)
ax4.set_xlabel('標準化偏回帰係数（β）', fontsize=11)
ax4.set_title('重回帰分析：標準化偏回帰係数と95%信頼区間\n'
              '目的変数：死亡率（安全・健康指標）',
              fontsize=13, fontweight='bold', pad=12)

# 有意水準マーク
for i, (b, pv, ci) in enumerate(zip(betas, beta_pvals, ci_95)):
    sig_mark = ('***' if pv < 0.001 else
                '**'  if pv < 0.01  else
                '*'   if pv < 0.05  else 'n.s.')
    offset = 0.025
    x_text = b + ci + offset if b >= 0 else b - ci - offset
    ha = 'left' if b >= 0 else 'right'
    ax4.text(x_text, y_pos[i], sig_mark,
             ha=ha, va='center', fontsize=9,
             color='#c00' if pv < 0.05 else '#888')

legend_els4 = [
    Patch(facecolor='#e05c5c', label='正の有意効果 (p<0.05)'),
    Patch(facecolor='#4e9af1', label='負の有意効果 (p<0.05)'),
    Patch(facecolor='#cccccc', label='非有意 (n.s.)'),
]
ax4.legend(handles=legend_els4, fontsize=9, loc='lower right')
ax4.text(0.02, 0.02,
         f'R²={model_std.rsquared:.3f}  Adj.R²={model_std.rsquared_adj:.3f}  N={len(df_merge)}',
         transform=ax4.transAxes, fontsize=9, color='#555')
ax4.spines['top'].set_visible(False)
ax4.spines['right'].set_visible(False)

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

print("\nDONE: 2019_H5_1_shorei")