"""
2018_H1_daijin.py
都道府県別死亡率の規定要因分析：高齢化・医療アクセス・経済格差の複合効果
総務大臣賞（高校生部門） 2018年度

使用データ：SSDSE-B-2026.csv（実データのみ、合成データ不使用）
"""

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

# 図の保存先フォルダを用意（プロジェクトルートから実行する前提）
os.makedirs('html/figures', exist_ok=True)

# ── データ読み込み ─────────────────────────────────────────────────────────────
df_b = pd.read_csv('data/raw/SSDSE-B-2026.csv', 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:", df_b.columns.tolist())

# 最新年（2023年）断面データ
latest_year = df_b['年度'].max()
df_cross = df_b[df_b['年度'] == latest_year].copy().reset_index(drop=True)
print(f"最新年度: {latest_year}, n={len(df_cross)}")

# ── 変数計算 ──────────────────────────────────────────────────────────────────
# 目的変数：死亡率（人口千人当たり）
df_cross['死亡率'] = df_cross['死亡数'] / df_cross['総人口'] * 1000

# 高齢化率（65歳以上人口割合 %）
df_cross['高齢化率'] = df_cross['65歳以上人口'] / df_cross['総人口'] * 100

# 医療アクセス指標：人口10万人当たり一般病院数
df_cross['病院数_10万対'] = df_cross['一般病院数'] / df_cross['総人口'] * 100000

# 保健医療費（世帯）
df_cross['保健医療費'] = df_cross['保健医療費（二人以上の世帯）']

# 有効求人倍率（月間有効求人数 / 月間有効求職者数）
df_cross['有効求人倍率'] = df_cross['月間有効求人数（一般）'] / df_cross['月間有効求職者数（一般）']

# 消費支出
df_cross['消費支出'] = df_cross['消費支出（二人以上の世帯）']

# 都道府県名の「県」「都」「道」「府」を省略（表示用）
def shorten_pref(name):
    # 北海道 → 北海道、京都府 → 京都 のように末尾の都道府県を除く
    for suffix in ['道', '都', '府', '県']:
        if name.endswith(suffix) and len(name) > len(suffix):
            return name[:-len(suffix)]
    return name

df_cross['pref_short'] = df_cross['都道府県'].apply(shorten_pref)

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

df_cross['地域'] = df_cross['pref_short'].map(region_map)
# 北海道は pref_short='北海' になるので都道府県名でも補完
region_map_full = {
    '北海道': '北海道・東北', '青森県': '北海道・東北', '岩手県': '北海道・東北', '宮城県': '北海道・東北',
    '秋田県': '北海道・東北', '山形県': '北海道・東北', '福島県': '北海道・東北',
    '茨城県': '関東', '栃木県': '関東', '群馬県': '関東', '埼玉県': '関東', '千葉県': '関東',
    '東京都': '関東', '神奈川県': '関東',
    '新潟県': '中部', '富山県': '中部', '石川県': '中部', '福井県': '中部', '山梨県': '中部',
    '長野県': '中部', '岐阜県': '中部', '静岡県': '中部', '愛知県': '中部',
    '三重県': '近畿', '滋賀県': '近畿', '京都府': '近畿', '大阪府': '近畿',
    '兵庫県': '近畿', '奈良県': '近畿', '和歌山県': '近畿',
    '鳥取県': '中国・四国', '島根県': '中国・四国', '岡山県': '中国・四国', '広島県': '中国・四国',
    '山口県': '中国・四国', '徳島県': '中国・四国', '香川県': '中国・四国',
    '愛媛県': '中国・四国', '高知県': '中国・四国',
    '福岡県': '九州・沖縄', '佐賀県': '九州・沖縄', '長崎県': '九州・沖縄', '熊本県': '九州・沖縄',
    '大分県': '九州・沖縄', '宮崎県': '九州・沖縄', '鹿児島県': '九州・沖縄', '沖縄県': '九州・沖縄'
}
df_cross['地域'] = df_cross['都道府県'].map(region_map_full)
df_cross['color'] = df_cross['地域'].map(region_colors)

# 確認
print("\n変数サマリー:")
print(df_cross[['都道府県','死亡率','高齢化率','病院数_10万対','有効求人倍率','保健医療費']].describe())

# ══════════════════════════════════════════════════════════════════════════════
# 図1：死亡率ランキング棒グラフ（47都道府県、地域色分け、全国平均線）
# ══════════════════════════════════════════════════════════════════════════════
print("\n図1 作成中...")
df_sorted = df_cross.sort_values('死亡率', ascending=False).reset_index(drop=True)
national_avg = df_cross['死亡率'].mean()

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

bars = ax.bar(
    range(len(df_sorted)),
    df_sorted['死亡率'],
    color=df_sorted['color'],
    edgecolor='white', linewidth=0.5
)

# 全国平均線
ax.axhline(national_avg, color='black', linestyle='--', linewidth=1.5,
           label=f'全国平均 {national_avg:.2f}‰')

# 軸設定
ax.set_xticks(range(len(df_sorted)))
ax.set_xticklabels(df_sorted['pref_short'], rotation=90, fontsize=8)
ax.set_ylabel('死亡率（人口千人当たり）', fontsize=12)
ax.set_title(f'図1：都道府県別死亡率ランキング（{latest_year}年）\n人口千人当たり死亡数（地域別色分け）', fontsize=13, fontweight='bold')
ax.set_ylim(0, df_sorted['死亡率'].max() * 1.15)

# 凡例（地域）
legend_handles = [
    plt.Rectangle((0, 0), 1, 1, color=c, label=r)
    for r, c in region_colors.items()
]
legend_handles.append(
    plt.Line2D([0], [0], color='black', linestyle='--', linewidth=1.5, label=f'全国平均 {national_avg:.2f}‰')
)
ax.legend(handles=legend_handles, loc='upper right', fontsize=9, ncol=2)

ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
plt.tight_layout()
fig.savefig('html/figures/2018_H1_fig1.png', dpi=150, bbox_inches='tight')
plt.close()
print("  → 2018_H1_fig1.png 保存完了")

# ══════════════════════════════════════════════════════════════════════════════
# 図2：Lorenz曲線（保健医療費の都道府県格差、Gini係数付き）
# ══════════════════════════════════════════════════════════════════════════════
print("図2 作成中...")

def gini(arr):
    v = np.sort(arr)
    n = len(v)
    cumv = np.cumsum(v)
    return (n + 1 - 2 * np.sum(cumv) / cumv[-1]) / n

# 保健医療費のGini係数・Lorenz曲線
health_vals = df_cross['保健医療費'].dropna().values
g_health = gini(health_vals)

# 病院数のGini係数・Lorenz曲線
hosp_vals = df_cross['病院数_10万対'].dropna().values
g_hosp = gini(hosp_vals)

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

plot_configs = [
    (axes[0], health_vals, '保健医療費（世帯月額・円）', g_health, '#1565C0'),
    (axes[1], hosp_vals,   '一般病院数（人口10万人当たり）', g_hosp,   '#00695C'),
]
for ax, vals, label, g, color in plot_configs:
    sorted_vals = np.sort(vals)
    cum_pop = np.linspace(0, 1, len(sorted_vals) + 1)
    cum_val = np.concatenate([[0], np.cumsum(sorted_vals) / sorted_vals.sum()])

    ax.fill_between(cum_pop, cum_pop, cum_val, alpha=0.3, color=color, label='格差面積')
    ax.plot(cum_pop, cum_val, color=color, linewidth=2.5, label=f'Lorenz曲線 (Gini={g:.3f})')
    ax.plot([0, 1], [0, 1], 'k--', linewidth=1.5, label='完全平等線')
    ax.set_xlabel('都道府県の累積人口割合', fontsize=11)
    ax.set_ylabel('累積シェア', fontsize=11)
    ax.set_title(f'{label}\nGini係数 = {g:.3f}', fontsize=11, fontweight='bold')
    ax.legend(fontsize=9)
    ax.set_xlim(0, 1)
    ax.set_ylim(0, 1)
    ax.grid(alpha=0.3)
    ax.spines['top'].set_visible(False)
    ax.spines['right'].set_visible(False)

fig.suptitle('図2：医療関連指標のLorenz曲線と格差（Gini係数）', fontsize=13, fontweight='bold', y=1.02)
plt.tight_layout()
fig.savefig('html/figures/2018_H1_fig2.png', dpi=150, bbox_inches='tight')
plt.close()
print(f"  → 2018_H1_fig2.png 保存完了 (保健医療費Gini={g_health:.3f}, 病院数Gini={g_hosp:.3f})")

# ══════════════════════════════════════════════════════════════════════════════
# 図3：有効求人倍率 vs 死亡率（47都道府県、地域色分け、都道府県ラベル、回帰直線）
# ══════════════════════════════════════════════════════════════════════════════
print("図3 作成中...")

x_var = df_cross['有効求人倍率'].values
y_var = df_cross['死亡率'].values

r_val, p_val = stats.pearsonr(x_var, y_var)
slope, intercept, _, _, _ = stats.linregress(x_var, y_var)
x_line = np.linspace(x_var.min(), x_var.max(), 100)
y_line = slope * x_line + intercept

fig, ax = plt.subplots(figsize=(11, 8))

for _, row in df_cross.iterrows():
    ax.scatter(row['有効求人倍率'], row['死亡率'],
               color=row['color'], s=60, zorder=3, alpha=0.85, edgecolors='white', linewidth=0.5)
    ax.annotate(row['pref_short'],
                xy=(row['有効求人倍率'], row['死亡率']),
                xytext=(3, 3), textcoords='offset points',
                fontsize=7, color='#333333')

ax.plot(x_line, y_line, 'k-', linewidth=1.8, alpha=0.7,
        label=f'回帰直線 (r={r_val:.3f}, p={p_val:.4f})')

# 地域凡例
legend_handles = [
    plt.scatter([], [], color=c, s=60, label=r)
    for r, c in region_colors.items()
]
legend_handles.append(
    plt.Line2D([0], [0], color='black', linewidth=1.8, label=f'回帰直線 (r={r_val:.3f}, p={p_val:.4f})')
)
ax.legend(handles=legend_handles, fontsize=9, loc='upper right')

ax.set_xlabel('有効求人倍率（月間有効求人数 / 有効求職者数）', fontsize=12)
ax.set_ylabel('死亡率（人口千人当たり）', fontsize=12)
ax.set_title(f'図3：有効求人倍率と死亡率の関係（{latest_year}年、47都道府県）\nr={r_val:.3f}, p={p_val:.4f}',
             fontsize=13, fontweight='bold')
ax.grid(alpha=0.25)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
plt.tight_layout()
fig.savefig('html/figures/2018_H1_fig3.png', dpi=150, bbox_inches='tight')
plt.close()
print(f"  → 2018_H1_fig3.png 保存完了 (r={r_val:.3f}, p={p_val:.4f})")

# ══════════════════════════════════════════════════════════════════════════════
# 重回帰分析（OLS）：標準化偏回帰係数
# ══════════════════════════════════════════════════════════════════════════════
print("\n重回帰分析...")
pred_cols = {
    '高齢化率 (%)': '高齢化率',
    '病院数\n(10万対)': '病院数_10万対',
    '有効求人倍率': '有効求人倍率',
    '保健医療費\n(円)': '保健医療費',
    '消費支出\n(円)': '消費支出',
}

df_reg = df_cross[['死亡率', '高齢化率', '病院数_10万対', '有効求人倍率', '保健医療費', '消費支出']].dropna()

# 標準化
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
X_cols = ['高齢化率', '病院数_10万対', '有効求人倍率', '保健医療費', '消費支出']
y_col  = '死亡率'
X_scaled = scaler.fit_transform(df_reg[X_cols])
y_scaled  = (df_reg[y_col].values - df_reg[y_col].mean()) / df_reg[y_col].std()

X_sm = sm.add_constant(X_scaled)
model = sm.OLS(y_scaled, X_sm).fit()
print(model.summary())

betas    = model.params[1:]   # 標準化偏回帰係数
pvalues  = model.pvalues[1:]
conf_int = model.conf_int()[1:]
var_labels = ['高齢化率 (%)', '病院数\n(10万対)', '有効求人倍率', '保健医療費\n(円)', '消費支出\n(円)']

# ══════════════════════════════════════════════════════════════════════════════
# 図4：標準化偏回帰係数（横棒グラフ + 95%CI + 有意マーク）
# ══════════════════════════════════════════════════════════════════════════════
print("図4 作成中...")

fig, ax = plt.subplots(figsize=(9, 5))

# 色分け（正=青、負=赤）
bar_colors = ['#1565C0' if b > 0 else '#C62828' for b in betas]

y_pos = np.arange(len(betas))
bars = ax.barh(y_pos, betas, color=bar_colors, alpha=0.8, height=0.55)

# 95% CI
for i, (lo, hi) in enumerate(conf_int):
    ax.plot([lo, hi], [i, i], 'k-', linewidth=2, zorder=4)
    ax.plot([lo, lo], [i - 0.12, i + 0.12], 'k-', linewidth=2)
    ax.plot([hi, hi], [i - 0.12, i + 0.12], 'k-', linewidth=2)

# 有意マーク
for i, p in enumerate(pvalues):
    mark = ''
    if p < 0.001:
        mark = '***'
    elif p < 0.01:
        mark = '**'
    elif p < 0.05:
        mark = '*'
    else:
        mark = 'n.s.'
    x_off = betas[i] + (0.03 if betas[i] >= 0 else -0.03)
    ha = 'left' if betas[i] >= 0 else 'right'
    ax.text(x_off, i, mark, va='center', ha=ha, fontsize=11,
            color='#C62828' if mark != 'n.s.' else '#888888', fontweight='bold')

ax.axvline(0, color='black', linewidth=1.0, linestyle='-')
ax.set_yticks(y_pos)
ax.set_yticklabels(var_labels, fontsize=11)
ax.set_xlabel('標準化偏回帰係数（β）', fontsize=12)
ax.set_title(f'図4：重回帰分析の標準化偏回帰係数（n=47都道府県、{latest_year}年）\n目的変数：死亡率（人口千人当たり）',
             fontsize=12, fontweight='bold')
ax.grid(axis='x', alpha=0.3)
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)

# R² 記載
ax.text(0.98, 0.04, f'R² = {model.rsquared:.3f}  adj.R² = {model.rsquared_adj:.3f}',
        transform=ax.transAxes, fontsize=10, ha='right',
        bbox=dict(boxstyle='round,pad=0.3', facecolor='#EFF3FF', alpha=0.8))

# 凡例（有意水準）
from matplotlib.patches import Patch
from matplotlib.lines import Line2D
leg_handles = [
    Patch(facecolor='#1565C0', alpha=0.8, label='正の効果'),
    Patch(facecolor='#C62828', alpha=0.8, label='負の効果'),
    Line2D([0], [0], color='k', linewidth=2, label='95% CI'),
]
ax.legend(handles=leg_handles, fontsize=9, loc='lower right')

plt.tight_layout()
fig.savefig('html/figures/2018_H1_fig4.png', dpi=150, bbox_inches='tight')
plt.close()
print("  → 2018_H1_fig4.png 保存完了")

# ── 相関係数まとめ ─────────────────────────────────────────────────────────────
print("\n=== Pearson相関係数（vs 死亡率）===")
for col in X_cols:
    r, p = stats.pearsonr(df_cross[col].dropna(), df_cross['死亡率'][df_cross[col].notna()])
    sig = '***' if p < 0.001 else ('**' if p < 0.01 else ('*' if p < 0.05 else 'n.s.'))
    print(f"  {col:20s}: r={r:+.4f}, p={p:.4f} {sig}")

print(f"\nGini係数:")
print(f"  保健医療費: {g_health:.3f}")
print(f"  病院数(10万対): {g_hosp:.3f}")

print("\nDONE: 2018_H1_daijin")