"""
教育用再現コード: 2025年 統計データ分析コンペティション 審査員奨励賞（大学生）
=================================================================================
論文タイトル：女性の社会進出を支援する環境要因分析と改善案の提言
受賞区分  ：審査員奨励賞 ［大学生・一般の部］

【分析概要】
  データ：SSDSE-B-2026.csv（47都道府県クロスセクション、2022年）
  目的   ：女性就業率を目的変数とし、LASSO→重回帰→PCA→階層クラスタリングの
            パイプラインで環境要因を特定し、都道府県の政策類型を明らかにする

  Step1. データ読み込みと変数構築
    - 女性就業率 = 女性就業者数 ÷ (総人口(女) - 15歳未満人口(女))
    - 幼稚園教員平均負担人数 = 幼稚園在園者数 ÷ 幼稚園教員数
    - 保育士平均負担人数 = 保育所等在所児数 ÷ 保育所等保育士数
    - 保育所入園倍率 = 保育所等在所児数 ÷ 保育所等定員数
  Step2. LASSO回帰（変数選択）
  Step3. 重回帰分析（選択変数で）
  Step4. 主成分分析（バリマックス回転）
  Step5. 階層クラスタリング（Ward法）

【データサイエンス学習ポイント】
  1. LASSO回帰（L1正則化による変数選択）
  2. 重回帰分析と調整済みR²
  3. 主成分分析（PCA）とバリマックス回転
  4. 階層クラスタリング（デンドログラム・Ward法）
  5. pandas によるデータ加工・変数構築

【入力データ】
  data/raw/SSDSE-B-2026.csv（2022年データを使用）

【出力図】
  html/figures/2025_U5_2_fig1_lasso.png  ... LASSO正則化パスとCV結果
  html/figures/2025_U5_2_fig2_reg.png    ... 重回帰係数プロット
  html/figures/2025_U5_2_fig3_pca.png    ... 主成分分析バイプロット
  html/figures/2025_U5_2_fig4_cluster.png ... 階層クラスタリングデンドログラム
=================================================================================
"""

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


import os
import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from matplotlib.gridspec import GridSpec
import warnings
warnings.filterwarnings('ignore')

from sklearn.linear_model import Lasso, LassoCV
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from sklearn.linear_model import LinearRegression
from scipy.cluster.hierarchy import dendrogram, linkage, fcluster
from scipy import stats

# ── パス設定 ─────────────────────────────────────────────────────────────────
DATA_RAW = 'data/raw/SSDSE-B-2026.csv'
DATA_D   = 'data/raw/SSDSE-D-2023.csv'
FIG_DIR  = 'html/figures'
os.makedirs(FIG_DIR, exist_ok=True)

plt.rcParams.update({
    'font.family':        'Hiragino Sans',
    'axes.unicode_minus': False,
    'figure.dpi':         150,
    'axes.spines.top':    False,
    'axes.spines.right':  False,
})

# ═══════════════════════════════════════════════════════════════════════════════
# ■ Step1. データ読み込みと変数構築
# ═══════════════════════════════════════════════════════════════════════════════

print("=" * 65)
print("■ Step1. SSDSE-B-2026.csv 読み込みと変数構築")
print("=" * 65)

df_raw = pd.read_csv(DATA_RAW, encoding='cp932', header=0, skiprows=[1])
df_raw = df_raw.rename(columns={'SSDSE-B-2026': 'year'})

# 2022年のみ抽出（最新の完全データ）
df = df_raw[df_raw['year'] == 2022].copy().reset_index(drop=True)
print(f"\n2022年データ: {len(df)} 都道府県")

# ── 変数構築 ─────────────────────────────────────────────────────────────────
# 総人口（女）: A110102 / 15歳未満人口（女）: A130102
# 15-64歳人口（女）: A130202 / 65歳以上人口（女）: A130302
# 幼稚園教員数: E1301 / 幼稚園在園者数: E1501
# 保育所等在所児数: J2506 / 保育所等保育士数: J2526 / 定員数: J2505
# 待機児童数: J250502

# 女性就業率（SSDSE-D 社会生活基本調査2021 — 都道府県別女性「仕事」平均時間を代理変数として使用）
df['pop_f_15plus']    = df['A110102'] - df['A130102']          # 15歳以上女性人口
df_d = pd.read_csv(DATA_D, encoding='cp932', header=1)
df_d_f = df_d[(df_d['男女の別'] == '2_女') & (df_d['地域コード'] != 'R00000')].copy()
df_d_f['仕事_min'] = pd.to_numeric(df_d_f['仕事'], errors='coerce')
d_map = df_d_f.set_index('地域コード')['仕事_min']
df['fem_employment_rate'] = df['地域コード'].map(d_map)
df['fem_employment_rate'] = pd.to_numeric(df['fem_employment_rate'], errors='coerce')
df['fem_employment_rate'].fillna(df['fem_employment_rate'].median(), inplace=True)

# 幼稚園教員平均負担人数（在園者÷教員数）
df['kindergarten_burden'] = (df['E1501'] / df['E1301'].replace(0, np.nan)).fillna(0)

# 保育士平均負担人数
df['nursery_burden']      = (df['J2506'] / df['J2526'].replace(0, np.nan)).fillna(0)

# 保育所入園倍率（在所÷定員、1超=超過充足・待機あり）
df['nursery_fill_rate']   = (df['J2506'] / df['J2505'].replace(0, np.nan)).fillna(0)

# 保育所等落選率（待機÷申込者≈待機÷在所）代理変数
df['nursery_wait_rate']   = (df['J250502'] / df['J2506'].replace(0, np.nan)).clip(0, 0.3).fillna(0)

# 15歳未満人口割合（女）
df['child_ratio_f']       = df['A130102'] / df['A110102']

# 第3次産業就業率の代理：直接変数が無いため、労働力データから近似
# F3105=就職件数(一般), F3104=充足数(一般), F3102=有効求職者数
# → 有効求人数/有効求職者数 = 有効求人倍率（高いほど第3次産業盛ん）
df['job_offer_ratio']     = (df['F3103'] / df['F3102'].replace(0, np.nan)).fillna(1.0)

# 大学卒業者割合（高学歴化 = 第3次産業就業と相関）
df['univ_grad_share']     = (df['E6502'] / df['A110102'] * 1000)  # 人口千人あたり

print(f"\n【構築変数の記述統計】")
key_vars = ['fem_employment_rate', 'child_ratio_f', 'kindergarten_burden',
            'nursery_burden', 'nursery_fill_rate', 'nursery_wait_rate',
            'job_offer_ratio', 'univ_grad_share']
print(df[key_vars].describe().round(4))

# 欠損値確認と除去
df_clean = df[key_vars + ['Prefecture']].dropna()
print(f"\n有効都道府県数: {len(df_clean)}")

# ═══════════════════════════════════════════════════════════════════════════════
# ■ Step2. LASSO回帰（変数選択）
# ═══════════════════════════════════════════════════════════════════════════════

print("\n" + "=" * 65)
print("■ Step2. LASSO回帰（LassoCV で最適 α を選択）")
print("=" * 65)

y = df_clean['fem_employment_rate'].values
X_names = ['child_ratio_f', 'kindergarten_burden', 'nursery_burden',
           'nursery_fill_rate', 'nursery_wait_rate', 'job_offer_ratio', 'univ_grad_share']
X = df_clean[X_names].values

scaler = StandardScaler()
X_std  = scaler.fit_transform(X)

# LassoCV: 5-fold CV で最適なαを選択
alphas_grid = np.logspace(-4, 0, 100)
lasso_cv    = LassoCV(alphas=alphas_grid, cv=5, random_state=2025, max_iter=5000)
lasso_cv.fit(X_std, y)
alpha_opt   = lasso_cv.alpha_

print(f"\n最適 α（LassoCV）: {alpha_opt:.6f}")

# 最適αでフィット
lasso_final = Lasso(alpha=alpha_opt, max_iter=5000)
lasso_final.fit(X_std, y)
lasso_coef  = lasso_final.coef_

print(f"\n【LASSO 係数（標準化）】")
var_labels = {
    'child_ratio_f':        '15歳未満人口割合（女）',
    'kindergarten_burden':  '幼稚園教員平均負担人数',
    'nursery_burden':       '保育士平均負担人数',
    'nursery_fill_rate':    '保育所入園倍率',
    'nursery_wait_rate':    '保育所落選率（代理）',
    'job_offer_ratio':      '有効求人倍率（第3次産業代理）',
    'univ_grad_share':      '大学卒業者数（千人あたり）',
}
selected_vars = []
for vn, coef in zip(X_names, lasso_coef):
    status = '✓ 選択' if abs(coef) > 1e-6 else '  ゼロ'
    print(f"  {status}  {var_labels[vn]:<20} β={coef:+.6f}")
    if abs(coef) > 1e-6:
        selected_vars.append(vn)

print(f"\nLASSOにより選択された変数数: {len(selected_vars)}")

# 正則化パス（複数αでの係数変化）
alphas_path = np.logspace(-4, 0, 50)
coef_path   = []
for a in alphas_path:
    l = Lasso(alpha=a, max_iter=5000)
    l.fit(X_std, y)
    coef_path.append(l.coef_.copy())
coef_path = np.array(coef_path)

# ═══════════════════════════════════════════════════════════════════════════════
# ■ Step3. 重回帰分析（LASSO選択変数を使用）
# ═══════════════════════════════════════════════════════════════════════════════

print("\n" + "=" * 65)
print("■ Step3. 重回帰分析（LASSO選択変数）")
print("=" * 65)

# LASSOで選択された変数が少ない場合は上位変数を使用
if len(selected_vars) < 3:
    # 係数の絶対値上位4変数を使用
    abs_coefs = np.abs(lasso_coef)
    top_idx   = np.argsort(abs_coefs)[::-1][:4]
    selected_vars = [X_names[i] for i in top_idx]
    print(f"  (LASSOで選択が少ないため上位係数変数{len(selected_vars)}個を使用)")

X_sel     = df_clean[selected_vars].values
X_sel_std = scaler.fit_transform(X_sel)

# statsmodels で詳細な検定結果を得る
import statsmodels.api as sm
X_ols = sm.add_constant(X_sel_std)
ols   = sm.OLS(y, X_ols).fit()

print(f"\n  調整済みR²: {ols.rsquared_adj:.4f}  (論文参考値: 0.647)")
print(f"  F統計量: {ols.fvalue:.4f}  (p={ols.f_pvalue:.4f})")
print(f"\n  {'変数名':<22} {'係数':>8} {'SE':>7} {'t値':>7} {'p値':>8} {'有意'}")
print("  " + "-" * 60)
print(f"  {'定数項':<22} {ols.params[0]:>8.4f} {ols.bse[0]:>7.4f} {ols.tvalues[0]:>7.3f} {ols.pvalues[0]:>8.4f}")
for i, vn in enumerate(selected_vars):
    p   = ols.pvalues[i + 1]
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else '†' if p < 0.1 else ''
    print(f"  {var_labels[vn]:<22} {ols.params[i+1]:>8.4f} {ols.bse[i+1]:>7.4f}"
          f" {ols.tvalues[i+1]:>7.3f} {p:>8.4f} {sig}")

# ═══════════════════════════════════════════════════════════════════════════════
# ■ Step4. 主成分分析（PCA）+ バリマックス回転
# ═══════════════════════════════════════════════════════════════════════════════

print("\n" + "=" * 65)
print("■ Step4. 主成分分析（PCA）+ バリマックス回転")
print("=" * 65)

X_all_std = StandardScaler().fit_transform(df_clean[X_names].values)
pca       = PCA(n_components=3)
scores    = pca.fit_transform(X_all_std)
loadings  = pca.components_.T  # shape: (n_vars, n_components)

print(f"\n【寄与率】")
for i, (ev, r, cr) in enumerate(zip(
    pca.explained_variance_,
    pca.explained_variance_ratio_,
    np.cumsum(pca.explained_variance_ratio_)
)):
    print(f"  PC{i+1}: 固有値={ev:.4f}, 寄与率={r*100:.1f}%, 累積={cr*100:.1f}%")

# バリマックス回転（手動実装）
def varimax(Phi, gamma=1.0, q=20, tol=1e-6):
    """
    バリマックス回転（直交回転）
    Phi: (n_vars, n_components) の因子負荷行列
    """
    p, k = Phi.shape
    R    = np.eye(k)
    d    = 0
    for _ in range(q):
        d_old = d
        for i in range(k):
            for j in range(i + 1, k):
                x = Phi @ R
                u = x[:, i] ** 2 - x[:, j] ** 2
                v = 2 * x[:, i] * x[:, j]
                A = np.sum(u)
                B = np.sum(v)
                C = np.sum(u ** 2 - v ** 2)
                D = np.sum(2 * u * v)
                num = D - 2 * gamma / p * A * B
                den = C - gamma / p * (A ** 2 - B ** 2)
                theta = np.arctan2(num, den) / 4
                Rot   = np.eye(k)
                c, s  = np.cos(theta), np.sin(theta)
                Rot[i, i] =  c; Rot[j, j] =  c
                Rot[i, j] = -s; Rot[j, i] =  s
                R = R @ Rot
        d = np.sum(pca.explained_variance_ratio_)
        if abs(d - d_old) < tol:
            break
    return Phi @ R, R

L_rotated, R_mat = varimax(loadings)
print(f"\n【バリマックス回転後の因子負荷（上位変数のみ）】")
print(f"  {'変数名':<22} {'F1':>8} {'F2':>8} {'F3':>8}")
print("  " + "-" * 50)
for i, vn in enumerate(X_names):
    print(f"  {var_labels[vn]:<22} {L_rotated[i,0]:>8.4f} {L_rotated[i,1]:>8.4f} {L_rotated[i,2]:>8.4f}")

# 回転後のスコア
scores_rot = scores @ R_mat

# ═══════════════════════════════════════════════════════════════════════════════
# ■ Step5. 階層クラスタリング（Ward法）
# ═══════════════════════════════════════════════════════════════════════════════

print("\n" + "=" * 65)
print("■ Step5. 階層クラスタリング（Ward法）")
print("=" * 65)

Z    = linkage(scores_rot, method='ward')
labels_cluster = fcluster(Z, t=4, criterion='maxclust')

df_clean = df_clean.copy()
df_clean['cluster'] = labels_cluster
df_clean['PC1']     = scores_rot[:, 0]
df_clean['PC2']     = scores_rot[:, 1]
df_clean['fem_emp'] = y

print(f"\n【クラスター別の都道府県】")
for cl in range(1, 5):
    prefs = df_clean[df_clean['cluster'] == cl]['Prefecture'].tolist()
    mean_emp = df_clean[df_clean['cluster'] == cl]['fem_emp'].mean()
    print(f"  Cluster{cl} (n={len(prefs)}, 女性就業率平均={mean_emp:.3f}): {', '.join(prefs)}")

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

CLUSTER_COLORS = {1: '#1565C0', 2: '#2E7D32', 3: '#E65100', 4: '#6A1B9A'}

# ────────────────────────────────────────────────────────────────────────────
# 図1：LASSO正則化パスとCVエラー
# ────────────────────────────────────────────────────────────────────────────
print("\n図1: LASSO正則化パスを作成中...")

fig1, axes1 = plt.subplots(1, 2, figsize=(13, 5))
fig1.suptitle('LASSO回帰による変数選択', fontsize=13, fontweight='bold')

# (a) 正則化パス
ax = axes1[0]
colors_path = plt.cm.tab10(np.linspace(0, 1, len(X_names)))
for i, (vn, col) in enumerate(zip(X_names, colors_path)):
    ax.plot(np.log10(alphas_path), coef_path[:, i], color=col,
            linewidth=1.8, label=var_labels[vn])
ax.axvline(np.log10(alpha_opt), color='red', linestyle='--', linewidth=1.5,
           label=f'最適 α={alpha_opt:.5f}')
ax.set_xlabel('log₁₀(α)', fontsize=11)
ax.set_ylabel('標準化係数', fontsize=11)
ax.set_title('LASSO正則化パス\n（αが大きいほど係数がゼロに縮小）', fontsize=11, fontweight='bold')
ax.legend(fontsize=7.5, loc='upper left', ncol=1)
ax.grid(True, alpha=0.3)
ax.axhline(0, color='gray', linewidth=0.8)

# (b) CV エラー曲線
ax2 = axes1[1]
mse_path = lasso_cv.mse_path_.mean(axis=1)
mse_std  = lasso_cv.mse_path_.std(axis=1)
ax2.semilogx(lasso_cv.alphas_, mse_path, color='#1565C0', linewidth=2)
ax2.fill_between(lasso_cv.alphas_, mse_path - mse_std, mse_path + mse_std,
                 alpha=0.2, color='#1565C0')
ax2.axvline(alpha_opt, color='red', linestyle='--', linewidth=1.5,
            label=f'最適 α={alpha_opt:.5f}')
ax2.set_xlabel('α（正則化強度）', fontsize=11)
ax2.set_ylabel('平均二乗誤差（MSE, CV）', fontsize=11)
ax2.set_title('5-fold CV によるα選択\n（MSE最小点を採用）', fontsize=11, fontweight='bold')
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)

plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2025_U5_2_fig1_lasso.png'), bbox_inches='tight', dpi=150)
plt.close(fig1)
print("  → 2025_U5_2_fig1_lasso.png 保存完了")

# ────────────────────────────────────────────────────────────────────────────
# 図2：重回帰係数プロット
# ────────────────────────────────────────────────────────────────────────────
print("図2: 重回帰係数プロットを作成中...")

fig2, axes2 = plt.subplots(1, 2, figsize=(13, 5))
fig2.suptitle('重回帰分析の結果（女性就業率を目的変数）', fontsize=13, fontweight='bold')

# (a) 係数と95%信頼区間
ax3 = axes2[0]
params   = np.asarray(ols.params)[1:]   # 定数項除く
ci       = np.asarray(ols.conf_int())
ci_lower = ci[1:, 0]
ci_upper = ci[1:, 1]
pvals    = np.asarray(ols.pvalues)[1:]
labels_sel = [var_labels[vn] for vn in selected_vars]

y_pos = np.arange(len(params))
colors_reg = ['#E53935' if p < 0.05 else '#90A4AE' for p in pvals]
ax3.barh(y_pos, params, color=colors_reg, alpha=0.8, edgecolor='white')
xerr_lo = np.maximum(0, params - ci_lower)
xerr_hi = np.maximum(0, ci_upper - params)
ax3.errorbar(params, y_pos, xerr=[xerr_lo, xerr_hi],
             fmt='none', color='black', capsize=4, linewidth=1.5)
ax3.axvline(0, color='gray', linestyle='--', linewidth=1.0)
ax3.set_yticks(y_pos)
ax3.set_yticklabels(labels_sel, fontsize=9)
ax3.set_xlabel('標準化偏回帰係数（±95%CI）', fontsize=11)
ax3.set_title(f'標準化偏回帰係数\n調整済みR²={ols.rsquared_adj:.3f}（論文参考値: 0.647）',
              fontsize=11, fontweight='bold')
red_patch  = mpatches.Patch(color='#E53935', label='有意（p<0.05）')
gray_patch = mpatches.Patch(color='#90A4AE', label='非有意')
ax3.legend(handles=[red_patch, gray_patch], fontsize=9)
ax3.grid(axis='x', alpha=0.3)

# 有意性マーカー
for i, (c, p) in enumerate(zip(params, pvals)):
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else ''
    if sig:
        offset = 0.003 if c >= 0 else -0.003
        ax3.text(ci_upper[i] + abs(offset), i, sig, va='center', fontsize=10,
                 color='#E53935', fontweight='bold')

# (b) 実績値 vs 予測値
ax4 = axes2[1]
y_pred  = ols.fittedvalues
residuals = y - y_pred
ax4.scatter(y_pred, y, color='#1565C0', alpha=0.6, s=50, edgecolors='white', linewidth=0.5)
lim = [min(y_pred.min(), y.min()) - 0.01, max(y_pred.max(), y.max()) + 0.01]
ax4.plot(lim, lim, 'r--', linewidth=1.5, label='完全予測線')
# 都道府県名ラベル（外れ値）
residual_std = np.std(residuals)
for i, (pref, yp, ya) in enumerate(zip(df_clean['Prefecture'], y_pred, y)):
    if abs(ya - yp) > residual_std * 1.5:
        ax4.annotate(pref, (yp, ya), fontsize=7.5, color='#555',
                     xytext=(5, 3), textcoords='offset points')
ax4.set_xlabel('予測値（女性就業率）', fontsize=11)
ax4.set_ylabel('実績値（女性就業率）', fontsize=11)
ax4.set_title(f'実績値 vs 予測値\nR²={ols.rsquared:.3f}', fontsize=11, fontweight='bold')
ax4.legend(fontsize=9)
ax4.grid(True, alpha=0.3)

plt.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2025_U5_2_fig2_reg.png'), bbox_inches='tight', dpi=150)
plt.close(fig2)
print("  → 2025_U5_2_fig2_reg.png 保存完了")

# ────────────────────────────────────────────────────────────────────────────
# 図3：主成分分析 バイプロット
# ────────────────────────────────────────────────────────────────────────────
print("図3: PCAバイプロットを作成中...")

fig3, axes3 = plt.subplots(1, 2, figsize=(13, 5))
fig3.suptitle('主成分分析（PCA）+ バリマックス回転', fontsize=13, fontweight='bold')

# (a) バイプロット（PC1 vs PC2）
ax5 = axes3[0]
scatter_colors = [CLUSTER_COLORS[c] for c in df_clean['cluster']]
sc = ax5.scatter(df_clean['PC1'], df_clean['PC2'], c=scatter_colors,
                 s=70, alpha=0.85, edgecolors='white', linewidth=0.5, zorder=3)
for i, (pref, x, y_val) in enumerate(zip(df_clean['Prefecture'],
                                          df_clean['PC1'], df_clean['PC2'])):
    ax5.annotate(pref[:2], (x, y_val), fontsize=6.5, ha='center', va='bottom',
                 xytext=(0, 4), textcoords='offset points', color='#333')

# 因子負荷ベクトル
scale = 3.0
for i, vn in enumerate(X_names):
    ax5.arrow(0, 0, L_rotated[i, 0] * scale, L_rotated[i, 1] * scale,
              head_width=0.08, head_length=0.05, fc='#E53935', ec='#E53935', alpha=0.7)
    ax5.text(L_rotated[i, 0] * scale * 1.12, L_rotated[i, 1] * scale * 1.12,
             var_labels[vn], fontsize=7, color='#C62828', ha='center')

ax5.axhline(0, color='gray', linewidth=0.6, alpha=0.5)
ax5.axvline(0, color='gray', linewidth=0.6, alpha=0.5)
ax5.set_xlabel(f'PC1（バリマックス回転後）', fontsize=11)
ax5.set_ylabel(f'PC2（バリマックス回転後）', fontsize=11)
ax5.set_title('バイプロット（PC1 × PC2）\nクラスターカラーで着色', fontsize=11, fontweight='bold')

for cl, col in CLUSTER_COLORS.items():
    ax5.scatter([], [], color=col, s=60, label=f'Cluster{cl}')
ax5.legend(fontsize=8, loc='upper left')
ax5.grid(True, alpha=0.2)

# (b) 寄与率（スクリープロット）
ax6 = axes3[1]
ev_ratio  = pca.explained_variance_ratio_ * 100
cum_ratio = np.cumsum(ev_ratio)
pc_labels = [f'PC{i+1}' for i in range(len(ev_ratio))]
ax6.bar(pc_labels, ev_ratio, color='#1565C0', alpha=0.8, edgecolor='white', label='各成分寄与率')
ax6.plot(pc_labels, cum_ratio, 'o-', color='#E53935', linewidth=2, markersize=7,
         markerfacecolor='white', markeredgewidth=2, label='累積寄与率')
ax6.axhline(70, color='gray', linestyle='--', linewidth=1.0, label='70%基準')
ax6.set_ylabel('寄与率（%）', fontsize=11)
ax6.set_title('固有値スクリープロット\n（3成分で累積寄与率の確認）', fontsize=11, fontweight='bold')
ax6.legend(fontsize=9)
ax6.grid(axis='y', alpha=0.3)
for i, (r, c) in enumerate(zip(ev_ratio, cum_ratio)):
    ax6.text(i, r + 1.0, f'{r:.1f}%', ha='center', fontsize=9, fontweight='bold', color='#1565C0')
    ax6.text(i, c + 1.5, f'{c:.1f}%', ha='center', fontsize=8, color='#E53935')

plt.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2025_U5_2_fig3_pca.png'), bbox_inches='tight', dpi=150)
plt.close(fig3)
print("  → 2025_U5_2_fig3_pca.png 保存完了")

# ────────────────────────────────────────────────────────────────────────────
# 図4：階層クラスタリング デンドログラム + クラスター別プロファイル
# ────────────────────────────────────────────────────────────────────────────
print("図4: 階層クラスタリングデンドログラムを作成中...")

fig4, axes4 = plt.subplots(1, 2, figsize=(15, 6))
fig4.suptitle('階層クラスタリング（Ward法）による都道府県の類型化', fontsize=13, fontweight='bold')

# (a) デンドログラム
ax7 = axes4[0]
pref_labels = df_clean['Prefecture'].tolist()
dend = dendrogram(Z, labels=pref_labels, orientation='left', ax=ax7,
                  color_threshold=Z[-3, 2],
                  leaf_font_size=7.5)
ax7.set_title('デンドログラム（Ward法）\n4クラスター', fontsize=11, fontweight='bold')
ax7.set_xlabel('距離（Ward）', fontsize=10)
ax7.grid(axis='x', alpha=0.3)

# (b) クラスター別の女性就業率と保育士負担
ax8 = axes4[1]
cl_means = df_clean.groupby('cluster')[['fem_emp', 'nursery_burden', 'child_ratio_f']].mean()

x_cl = np.arange(1, 5)
width = 0.28
ax8.bar(x_cl - width, cl_means['fem_emp'] * 100, width, label='女性就業率（%）',
        color='#1565C0', alpha=0.85, edgecolor='white')
ax8.bar(x_cl, cl_means['nursery_burden'], width, label='保育士1人あたり在所児数（人）',
        color='#E65100', alpha=0.85, edgecolor='white')
ax8.bar(x_cl + width, cl_means['child_ratio_f'] * 100, width, label='15歳未満人口割合（女, %）',
        color='#2E7D32', alpha=0.85, edgecolor='white')

ax8.set_xticks(x_cl)
ax8.set_xticklabels([f'Cluster{c}' for c in range(1, 5)], fontsize=10)
ax8.set_ylabel('各指標の値', fontsize=11)
ax8.set_title('クラスター別の主要指標（平均値）', fontsize=11, fontweight='bold')
ax8.legend(fontsize=8.5, loc='upper right')
ax8.grid(axis='y', alpha=0.3)

# クラスター別の都道府県数
for i, cl in enumerate(range(1, 5)):
    n = (df_clean['cluster'] == cl).sum()
    ax8.text(cl, max(cl_means.loc[cl]) + 0.5, f'n={n}', ha='center', fontsize=8)

plt.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2025_U5_2_fig4_cluster.png'), bbox_inches='tight', dpi=150)
plt.close(fig4)
print("  → 2025_U5_2_fig4_cluster.png 保存完了")

# ═══════════════════════════════════════════════════════════════════════════════
# ■ 最終サマリー
# ═══════════════════════════════════════════════════════════════════════════════
print("\n" + "=" * 65)
print("✓ 全図の生成完了（4枚）")
print("=" * 65)
print(f"\n  調整済みR²: {ols.rsquared_adj:.4f}（論文参考値: 0.647）")
print(f"  主成分3つの累積寄与率: {np.cumsum(pca.explained_variance_ratio_)[-1]*100:.1f}%")
print(f"  クラスター数: 4（Ward法）")
print(f"\n  【論文の主要知見（参考）】")
print(f"    - 第3次産業就業率が高いほど女性就業率が低下（仕事の構造的問題）")
print(f"    - 15歳未満人口割合が高いと女性就業率が高い（若年層が多い地域の活性化）")
print(f"    - 保育充足度（保育士負担・入園倍率）が女性就業に重要")