小中高生の自殺要因についての
統計的検討

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)

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))

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")

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")

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")

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')

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,
}

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}")

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}")

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}")

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}")

os.makedirs(FIG_DIR, exist_ok=True)

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枚）")