"""
教育用再現コード: 2025年 統計データ分析コンペティション 審査員奨励賞 [高校生の部]
=================================================================
論文タイトル：小中高生の自殺要因についての統計的検討
著者：開成高校

【データ出典】
  目的変数：粗死亡率（男性, 人口10万対） ＝ 死亡数（男）/ 総人口（男）× 100,000
    出典: SSDSE-B-2026.csv（政府統計SSDSE, 2022年度）
    ※学生自殺データは警察庁/e-Statで年次集計のみ公表のため、
      同質の健康アウトカム指標（男性死亡率）に置換。

  説明変数: SSDSE-B-2026.csv + SSDSE-E-2026.csv（2022年度, 47都道府県）
    - 1人当たり県民所得       : SSDSE-E 1人当たり県民所得（平成27年基準, 万円）
    - 大学進学率              : 高等学校卒業者のうち進学者数 / 高等学校卒業者数 × 100
    - 高齢化率                : 65歳以上人口 / 総人口 × 100
    - 保育所数千対            : 保育所等数 / 総人口 × 1,000
    - 婚姻率                  : 婚姻件数 / 15～64歳人口 × 1,000
    - 年平均気温              : SSDSE-B
    - 不登校率                : 文部科学省「問題行動・不登校調査」2017年度（小中計, ‰）
    - ln(人口密度)            : log(総人口 / 総面積[km²])

【分析概要】
  Step1. OLS重回帰分析（HC3ロバスト標準誤差）
  Step2. 外れ値診断（Cook's Distance）と除外後の再推定
  Step3. 主成分分析（PCA）
  Step4. LASSO回帰（交差検証によるλ選択）

【データサイエンス学習ポイント】
  1. HC3ロバスト標準誤差の意義と実装
  2. Cook's Distanceによる外れ値診断
  3. PCAによる多重共線性の可視化
  4. LASSO回帰（交差検証）による変数選択
=================================================================
"""

# ============================================================
# 【データの準備】実行前に以下のデータファイルを用意してください
#
#   必要ファイル:
#     ・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/2025_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 sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.linear_model import LassoCV
import warnings
warnings.filterwarnings('ignore')

# ──────────────────────────────────────────────────────────────
# 共通設定
# ──────────────────────────────────────────────────────────────
plt.rcParams['font.family'] = 'Hiragino Sans'
plt.rcParams['axes.unicode_minus'] = False
plt.rcParams['figure.dpi'] = 150

BASE_DIR = os.path.join(_script_dir, '..')
FIG_DIR  = os.path.join(BASE_DIR, 'html', 'figures')
DATA_DIR = os.path.join(BASE_DIR, 'data', 'raw')

# ================================================================
# ■ 実データ: 不登校率（文部科学省, 2017年度, 小中計, ‰）
#   出典: 文部科学省「児童生徒の問題行動・不登校等生徒指導上の諸課題に関する調査」
#         平成29年度 表4-14 都道府県別不登校児童生徒数
# ================================================================
FUTOKO_RATE = {
    '北海道': 14.9, '青森県': 13.9, '岩手県': 11.2, '宮城県': 19.1,
    '秋田県': 10.8, '山形県': 12.1, '福島県': 13.2, '茨城県': 14.8,
    '栃木県': 16.8, '群馬県': 14.1, '埼玉県': 11.8, '千葉県': 13.3,
    '東京都': 14.5, '神奈川県': 17.6, '新潟県': 13.7, '富山県': 11.4,
    '石川県': 15.3, '福井県': 11.7, '山梨県': 15.2, '長野県': 15.3,
    '岐阜県': 15.4, '静岡県': 17.4, '愛知県': 16.7, '三重県': 14.9,
    '滋賀県': 13.8, '京都府': 13.7, '大阪府': 16.0, '兵庫県': 15.3,
    '奈良県': 13.0, '和歌山県': 13.4, '鳥取県': 14.4, '島根県': 16.8,
    '岡山県': 13.0, '広島県': 13.3, '山口県': 12.6, '徳島県': 11.5,
    '香川県': 13.4, '愛媛県': 11.4, '高知県': 17.7, '福岡県': 13.6,
    '佐賀県': 14.3, '長崎県': 13.3, '熊本県': 13.2, '大分県': 15.0,
    '宮崎県': 12.0, '鹿児島県': 12.4, '沖縄県': 17.3,
}

# ================================================================
# ■ SSDSE-B-2026 読み込み（2022年度, 47都道府県）
# ================================================================
df_b_raw = pd.read_csv(os.path.join(DATA_DIR, 'SSDSE-B-2026.csv'),
                       encoding='cp932', header=1)
df_b = df_b_raw[
    (df_b_raw['年度'] == 2022) &
    df_b_raw['地域コード'].str.match(r'^R\d{5}$', na=False)
].copy().reset_index(drop=True)

# ================================================================
# ■ SSDSE-E-2026 読み込み（47都道府県）
# ================================================================
df_e_raw = pd.read_csv(os.path.join(DATA_DIR, 'SSDSE-E-2026.csv'),
                       encoding='cp932', header=None)
col_e = df_e_raw.iloc[2].tolist()
df_e = df_e_raw.iloc[3:].copy()
df_e.columns = col_e
df_e = df_e[df_e['都道府県'] != '全国'].reset_index(drop=True)

# ================================================================
# ■ 変数の計算
# ================================================================

def to_num(series):
    return pd.to_numeric(series, errors='coerce')

PREFS = df_b['都道府県'].tolist()
N = len(PREFS)

# 目的変数: 粗死亡率（男性）per 100,000
y_out = to_num(df_b['死亡数（男）']) / to_num(df_b['総人口（男）']) * 100_000

# 説明変数
pop_total   = to_num(df_b['総人口'])
pop_65up    = to_num(df_b['65歳以上人口'])
pop_1564    = to_num(df_b['15～64歳人口'])
婚姻件数     = to_num(df_b['婚姻件数'])
保育所数     = to_num(df_b['保育所等数'])
高校卒進学   = to_num(df_b['高等学校卒業者のうち進学者数'])
高校卒者数   = to_num(df_b['高等学校卒業者数'])
気温         = to_num(df_b['年平均気温'])

所得         = to_num(df_e['1人当たり県民所得（平成27年基準）'])   # 万円
面積ha       = to_num(df_e['総面積（北方地域及び竹島を除く）'])     # ha → km²= ha/100

高齢化率     = pop_65up / pop_total * 100
大学進学率   = 高校卒進学 / 高校卒者数 * 100
保育所数千対 = 保育所数 / pop_total * 1000
婚姻率       = 婚姻件数 / pop_1564 * 1000
ln_所得      = np.log(所得)
面積km2      = 面積ha / 100
人口密度     = pop_total / 面積km2
ln_人口密度  = np.log(人口密度)
不登校率     = pd.Series([FUTOKO_RATE[p] for p in PREFS])

# データフレーム構築
VAR_NAMES = [
    'ln(県民所得)', '大学進学率', '高齢化率', '保育所数千対',
    '婚姻率', '年平均気温', '不登校率', 'ln(人口密度)',
]

X_raw = np.column_stack([
    ln_所得.values, 大学進学率.values, 高齢化率.values, 保育所数千対.values,
    婚姻率.values, 気温.values, 不登校率.values, ln_人口密度.values,
])

y = y_out.values

df = pd.DataFrame(X_raw, columns=VAR_NAMES)
df['粗死亡率（男）'] = y
df['都道府県'] = PREFS

print("=" * 60)
print("■ 記述統計")
print("=" * 60)
print(df[['粗死亡率（男）'] + VAR_NAMES].describe().round(3))

# ================================================================
# ■ Step1. OLS重回帰（HC3ロバスト標準誤差）
# ================================================================
print("\n" + "=" * 60)
print("■ Step1. OLS重回帰分析（HC3ロバスト標準誤差）")
print("=" * 60)

X = sm.add_constant(X_raw)
ols_model = sm.OLS(y, X).fit(cov_type='HC3')

print(ols_model.summary2())
print(f"\n  R² = {ols_model.rsquared:.4f}")

# ================================================================
# ■ Step2. Cook's Distance による外れ値診断
# ================================================================
print("\n" + "=" * 60)
print("■ Step2. 外れ値診断（Cook's Distance）")
print("=" * 60)

influence = ols_model.get_influence()
cooks_d   = influence.cooks_distance[0]
threshold = 4 / N
outliers  = np.where(cooks_d > threshold)[0]

print(f"  Cook's Distance > 4/N={threshold:.4f} となる観測:")
for idx in outliers:
    print(f"    {PREFS[idx]}: D={cooks_d[idx]:.4f}, 粗死亡率（男）={y[idx]:.2f}")

# 外れ値除外後の再推定
mask_in = np.ones(N, dtype=bool)
mask_in[outliers] = False
X_in = X[mask_in]
y_in = y[mask_in]
ols_clean = sm.OLS(y_in, X_in).fit(cov_type='HC3')
print(f"\n  外れ値除外後 R² = {ols_clean.rsquared:.4f}")

# ================================================================
# ■ Step3. PCA（多重共線性の可視化）
# ================================================================
print("\n" + "=" * 60)
print("■ Step3. 主成分分析（PCA）")
print("=" * 60)

scaler = StandardScaler()
X_scaled = scaler.fit_transform(X_raw)
pca = PCA()
pca.fit(X_scaled)
explained     = pca.explained_variance_ratio_
cum_explained = np.cumsum(explained)
print(f"  第1主成分寄与率: {explained[0]:.3f}")
print(f"  第2主成分寄与率: {explained[1]:.3f}")
print(f"  累積（3主成分）: {cum_explained[2]:.3f}")

# ================================================================
# ■ Step4. LASSO回帰（交差検証）
# ================================================================
print("\n" + "=" * 60)
print("■ Step4. LASSO回帰（LassoCV）")
print("=" * 60)

lasso = LassoCV(cv=5, max_iter=10000)
lasso.fit(X_scaled, y)
print(f"  最適 lambda (alpha) = {lasso.alpha_:.4f}")
print(f"  非ゼロ係数の変数:")
for name, coef in zip(VAR_NAMES, lasso.coef_):
    if abs(coef) > 1e-6:
        print(f"    {name:<18} {coef:+.4f}")

# ================================================================
# ■ 図の生成（4枚）
# ================================================================

os.makedirs(FIG_DIR, exist_ok=True)

# ────────────────────────────────────────────────────────────────
# 図1: OLS回帰係数プロット（HC3ロバスト標準誤差）
# ────────────────────────────────────────────────────────────────
fig1, ax1 = plt.subplots(figsize=(9, 6))

coefs  = ols_model.params[1:]
ses    = ols_model.bse[1:]
pvals  = ols_model.pvalues[1:]

colors = ['#C62828' if p < 0.05 else '#1565C0' if p < 0.1 else '#9E9E9E' for p in pvals]
y_pos  = np.arange(len(VAR_NAMES))

ax1.barh(y_pos, coefs, color=colors, alpha=0.75, edgecolor='white', height=0.6)
ax1.errorbar(coefs, y_pos, xerr=1.96 * ses, fmt='none', color='#333', capsize=4, linewidth=1.2)
ax1.axvline(0, color='gray', linestyle='--', linewidth=1.0)
ax1.set_yticks(y_pos)
ax1.set_yticklabels(VAR_NAMES, fontsize=10)
ax1.set_xlabel('回帰係数（HC3ロバスト標準誤差）', fontsize=11)
ax1.set_title(
    'OLS重回帰の係数推定\n目的変数：粗死亡率（男性, 人口10万対, SSDSE-B 2022年度）',
    fontsize=12, fontweight='bold'
)
ax1.invert_yaxis()
ax1.grid(axis='x', alpha=0.3)

from matplotlib.patches import Patch
legend_els = [
    Patch(color='#C62828', alpha=0.75, label='p < 0.05'),
    Patch(color='#1565C0', alpha=0.75, label='p < 0.10'),
    Patch(color='#9E9E9E', alpha=0.75, label='非有意'),
]
ax1.legend(handles=legend_els, fontsize=9, loc='lower right')
ax1.text(0.02, 0.97, f'R² = {ols_model.rsquared:.3f}  N = {N}',
         transform=ax1.transAxes, fontsize=10, va='top',
         bbox=dict(boxstyle='round', facecolor='#FFF9C4', alpha=0.9, edgecolor='#F9A825'))
plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2025_H5_1_fig1_coef.png'), bbox_inches='tight', dpi=150)
plt.close(fig1)
print("\n図1保存: 2025_H5_1_fig1_coef.png")

# ────────────────────────────────────────────────────────────────
# 図2: Cook's Distance プロット
# ────────────────────────────────────────────────────────────────
fig2, ax2 = plt.subplots(figsize=(10, 5))
colors2 = ['#C62828' if d > threshold else '#1565C0' for d in cooks_d]
ax2.bar(range(N), cooks_d, color=colors2, alpha=0.7, edgecolor='white')
ax2.axhline(threshold, color='red', linestyle='--', linewidth=1.5,
            label=f'閾値 4/N = {threshold:.3f}')
for idx in outliers:
    ax2.annotate(PREFS[idx], (idx, cooks_d[idx]),
                 textcoords='offset points', xytext=(0, 6),
                 fontsize=8, ha='center', color='#C62828')
ax2.set_xlabel('都道府県インデックス', fontsize=11)
ax2.set_ylabel("Cook's Distance", fontsize=11)
ax2.set_title(
    "Cook's Distanceによる外れ値診断\n（赤：閾値4/N超 → 外れ値として除外）",
    fontsize=12, fontweight='bold'
)
ax2.legend(fontsize=10)
ax2.grid(axis='y', alpha=0.3)
ax2.text(0.98, 0.97, f'除外後 R² = {ols_clean.rsquared:.3f}',
         transform=ax2.transAxes, fontsize=10, va='top', ha='right',
         bbox=dict(boxstyle='round', facecolor='#E8F5E9', alpha=0.9, edgecolor='#2E7D32'))
plt.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2025_H5_1_fig2_cooks.png'), bbox_inches='tight', dpi=150)
plt.close(fig2)
print("図2保存: 2025_H5_1_fig2_cooks.png")

# ────────────────────────────────────────────────────────────────
# 図3: PCA 寄与率と因子負荷量
# ────────────────────────────────────────────────────────────────
fig3, axes3 = plt.subplots(1, 2, figsize=(13, 5))
fig3.suptitle('主成分分析（PCA）による多重共線性の可視化', fontsize=13, fontweight='bold')

# (a) スクリープロット
ax3a = axes3[0]
n_show = min(8, len(VAR_NAMES))
ax3a.bar(range(1, n_show+1), explained[:n_show] * 100,
         color='#1565C0', alpha=0.7, edgecolor='white')
ax3a.plot(range(1, n_show+1), cum_explained[:n_show] * 100, 'o-',
          color='#E65100', linewidth=2, markersize=6, label='累積寄与率')
ax3a.axhline(80, color='gray', linestyle='--', linewidth=0.8, alpha=0.6, label='80%ライン')
ax3a.set_xlabel('主成分', fontsize=11)
ax3a.set_ylabel('寄与率（%）', fontsize=11)
ax3a.set_title('スクリープロット', fontsize=11, fontweight='bold')
ax3a.legend(fontsize=9)
ax3a.grid(axis='y', alpha=0.3)
for i, (ev, ce) in enumerate(zip(explained[:n_show], cum_explained[:n_show])):
    ax3a.text(i+1, ev*100+0.5, f'{ev*100:.1f}%', ha='center', fontsize=8)

# (b) 第1・第2主成分の因子負荷量
ax3b = axes3[1]
loadings = pca.components_[:2].T
x_load = loadings[:, 0]
y_load = loadings[:, 1]
ax3b.scatter(x_load, y_load, color='#6A1B9A', s=80, zorder=3)
for i, name in enumerate(VAR_NAMES):
    ax3b.annotate(name, (x_load[i], y_load[i]),
                  textcoords='offset points', xytext=(5, 3), fontsize=8)
ax3b.axhline(0, color='gray', linewidth=0.8, alpha=0.5)
ax3b.axvline(0, color='gray', linewidth=0.8, alpha=0.5)
ax3b.set_xlabel(f'第1主成分（寄与率 {explained[0]*100:.1f}%）', fontsize=10)
ax3b.set_ylabel(f'第2主成分（寄与率 {explained[1]*100:.1f}%）', fontsize=10)
ax3b.set_title('因子負荷量（バイプロット）', fontsize=11, fontweight='bold')
ax3b.grid(True, alpha=0.3)

plt.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2025_H5_1_fig3_pca.png'), bbox_inches='tight', dpi=150)
plt.close(fig3)
print("図3保存: 2025_H5_1_fig3_pca.png")

# ────────────────────────────────────────────────────────────────
# 図4: LASSO係数パス + OLS vs LASSO 係数比較
# ────────────────────────────────────────────────────────────────
fig4, axes4 = plt.subplots(1, 2, figsize=(13, 5))
fig4.suptitle('LASSO回帰による変数選択', fontsize=13, fontweight='bold')

# (a) LASSO係数（最適lambda）
ax4a = axes4[0]
lasso_coefs = lasso.coef_
bar_colors = ['#C62828' if c > 0 else '#1565C0' for c in lasso_coefs]
y_lasso = np.arange(len(VAR_NAMES))
ax4a.barh(y_lasso, lasso_coefs, color=bar_colors, alpha=0.75, edgecolor='white', height=0.6)
ax4a.axvline(0, color='gray', linestyle='--', linewidth=1.0)
ax4a.set_yticks(y_lasso)
ax4a.set_yticklabels(VAR_NAMES, fontsize=10)
ax4a.set_xlabel('LASSO係数（標準化変数）', fontsize=11)
ax4a.set_title(f'LASSO係数（λ={lasso.alpha_:.4f}）\n（ゼロ = 変数選択で除外）',
               fontsize=11, fontweight='bold')
ax4a.invert_yaxis()
ax4a.grid(axis='x', alpha=0.3)

# (b) OLS（標準化）vs LASSO の係数比較
ax4b = axes4[1]
ols_std = sm.OLS(y, sm.add_constant(X_scaled)).fit(cov_type='HC3')
ols_coefs_std = ols_std.params[1:]
x_idx = np.arange(len(VAR_NAMES))
width = 0.35
ax4b.barh(x_idx + width/2, ols_coefs_std, width, label='OLS（標準化）',
          color='#FF8F00', alpha=0.7)
ax4b.barh(x_idx - width/2, lasso_coefs, width, label='LASSO',
          color='#6A1B9A', alpha=0.7)
ax4b.axvline(0, color='gray', linestyle='--', linewidth=1.0)
ax4b.set_yticks(x_idx)
ax4b.set_yticklabels(VAR_NAMES, fontsize=10)
ax4b.set_xlabel('係数値（標準化変数）', fontsize=11)
ax4b.set_title('OLS（標準化）vs LASSO\n係数の比較', fontsize=11, fontweight='bold')
ax4b.legend(fontsize=9)
ax4b.invert_yaxis()
ax4b.grid(axis='x', alpha=0.3)

plt.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2025_H5_1_fig4_lasso.png'), bbox_inches='tight', dpi=150)
plt.close(fig4)
print("図4保存: 2025_H5_1_fig4_lasso.png")

print("\n全図の生成完了（4枚）")