"""
教育用再現コード: 2025年 統計データ分析コンペティション 審査員奨励賞（大学生）
=================================================================
論文タイトル：人口動態の地域差に着目した基礎的分析
著者：山形陽生（武蔵精密工業株式会社）

【分析概要】
  データ：47都道府県（2022年）
  目的変数：自然増減率・社会増減率（SSDSE-Bから実計算）

  Step1. ボックスプロットで外れ値検出（沖縄・秋田・東京）
  Step2. ガウス混合モデル（GMM）クラスタリング
         → 外れ値3都県を独立クラスタ、残り44都道府県を3クラスタ
  Step3. PLSregression（Partial Least Squares）回帰分析
  Step4. VIF（多重共線性確認）
  Step5. バイオリンプロット

  Key findings:
    - 沖縄（高自然増）・秋田（低自然増）・東京（高社会増）が外れ値
    - 44都道府県をGMMで3クラスタに類型化
    - クラスタごとにPLS回帰で寄与変数が異なる → 地域特性対応施策の必要性

【データ】
  data/raw/SSDSE-B-2026.csv（SSDSE-B: 47都道府県 2022年データ）
  data/raw/SSDSE-E-2026.csv（SSDSE-E: 1人当たり県民所得）
  ※ 合成データ・乱数は一切使用しない
=================================================================
"""

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


import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from sklearn.mixture import GaussianMixture
from sklearn.cross_decomposition import PLSRegression
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import cross_val_score
from scipy import stats
import warnings
warnings.filterwarnings('ignore')

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

import os
DATA_DIR = 'data/raw'
FIG_DIR  = 'html/figures'
os.makedirs(FIG_DIR, exist_ok=True)

# ================================================================
# ■ Step 0. データ構築（SSDSE-B 2022年 + SSDSE-E）
# ================================================================

print("=" * 65)
print("■ データ読み込み・構築（SSDSE-B & SSDSE-E 実データのみ）")
print("=" * 65)

# ── SSDSE-B-2026 読み込み ────────────────────────────────────
df_b_raw = pd.read_csv(
    os.path.join(DATA_DIR, 'SSDSE-B-2026.csv'),
    encoding='cp932', header=1
)

# 47都道府県のみ
df_b_raw = df_b_raw[df_b_raw['地域コード'].str.match(r'^R\d{5}$', na=False)].copy()
df_b_raw['年度'] = pd.to_numeric(df_b_raw['年度'], errors='coerce')

# 2022年データを使用
df_b = df_b_raw[df_b_raw['年度'] == 2022].copy().reset_index(drop=True)
print(f"SSDSE-B 2022年: {len(df_b)}都道府県")

# 必要列を数値変換
num_cols_b = [
    '総人口', '65歳以上人口', '15～64歳人口', '15歳未満人口',
    '出生数', '死亡数',
    '転入者数（日本人移動者）', '転出者数（日本人移動者）',
    '保育所等数', '保育所等定員数',
    '婚姻件数', '年平均気温',
    '高等学校卒業者数', '高等学校卒業者のうち進学者数',
    '合計特殊出生率',
    '一般病院数', '月間有効求人数（一般）',
    '標準価格（平均価格）（住宅地）',
    '消費支出（二人以上の世帯）',
    '着工新設持家数',
    '教育費（二人以上の世帯）',
]
for c in num_cols_b:
    df_b[c] = pd.to_numeric(df_b[c], errors='coerce')

# ── SSDSE-E-2026 読み込み ────────────────────────────────────
df_e_raw = pd.read_csv(
    os.path.join(DATA_DIR, 'SSDSE-E-2026.csv'),
    encoding='cp932', header=1
)
df_e = df_e_raw.iloc[1:].copy()
df_e.columns = df_e_raw.iloc[0].values
df_e = df_e[df_e['都道府県'] != '全国'].reset_index(drop=True)
df_e['県民所得_1人'] = pd.to_numeric(df_e['1人当たり県民所得（平成27年基準）'], errors='coerce')
df_e['医師数'] = pd.to_numeric(df_e['医師数'], errors='coerce')
df_e['面積_km2'] = pd.to_numeric(df_e['総面積（北方地域及び竹島を除く）'], errors='coerce') / 100.0
print(f"SSDSE-E: {len(df_e)}都道府県")

# ── 2つのデータを結合 ─────────────────────────────────────────
df_pref = df_b[['都道府県', '総人口', '65歳以上人口', '15～64歳人口', '15歳未満人口',
                '出生数', '死亡数', '転入者数（日本人移動者）', '転出者数（日本人移動者）',
                '保育所等数', '保育所等定員数', '婚姻件数', '年平均気温',
                '高等学校卒業者数', '高等学校卒業者のうち進学者数', '合計特殊出生率',
                '一般病院数', '月間有効求人数（一般）',
                '標準価格（平均価格）（住宅地）',
                '消費支出（二人以上の世帯）', '着工新設持家数',
                '教育費（二人以上の世帯）',
                ]].copy()
df_pref = df_pref.rename(columns={'都道府県': 'pref'})
df_pref = df_pref.merge(
    df_e[['都道府県', '県民所得_1人', '医師数', '面積_km2']].rename(columns={'都道府県': 'pref'}),
    on='pref', how='left'
)

# ── 主要指標の計算（実データのみ） ───────────────────────────
# 自然増減率: (出生数 - 死亡数) / 総人口 × 1000  (‰)
df_pref['自然増減率'] = (df_pref['出生数'] - df_pref['死亡数']) / df_pref['総人口'] * 1000

# 社会増減率: (転入者数 - 転出者数) / 総人口 × 1000  (‰)
df_pref['社会増減率'] = (df_pref['転入者数（日本人移動者）'] - df_pref['転出者数（日本人移動者）']) / df_pref['総人口'] * 1000

# 高齢化率: 65歳以上人口 / 総人口 × 100
df_pref['高齢化率'] = df_pref['65歳以上人口'] / df_pref['総人口'] * 100

# 保育所充実度: 保育所等数 / 総人口 × 10000
df_pref['保育所充実度'] = df_pref['保育所等数'] / df_pref['総人口'] * 10000

# 保育所定員率: 保育所等定員数 / 総人口 × 1000
df_pref['保育所定員率'] = df_pref['保育所等定員数'] / df_pref['総人口'] * 1000

# 大学進学率: 高校卒業者うち進学者 / 高校卒業者 × 100
df_pref['大学進学率'] = df_pref['高等学校卒業者のうち進学者数'] / df_pref['高等学校卒業者数'] * 100

# 婚姻率: 婚姻件数 / 15〜64歳人口 × 1000
df_pref['婚姻率'] = df_pref['婚姻件数'] / df_pref['15～64歳人口'] * 1000

# 医師密度: 医師数 / 総人口 × 10000
df_pref['医師密度'] = df_pref['医師数'] / df_pref['総人口'] * 10000

# 人口密度: 総人口 / 面積km²
df_pref['人口密度'] = df_pref['総人口'] / df_pref['面積_km2']

print(f"\n構築データ: {len(df_pref)}都道府県")
print(f"変数数: {len(df_pref.columns)}列")
print("\n記述統計（主要変数）:")
print(df_pref[['自然増減率','社会増減率','高齢化率','婚姻率']].describe().round(3))

# ================================================================
# ■ Step 1. 外れ値の特定（ボックスプロット基準）
# ================================================================

print("\n" + "=" * 65)
print("■ Step1. 外れ値の特定")
print("=" * 65)

df_pref = df_pref.dropna(subset=['自然増減率', '社会増減率']).reset_index(drop=True)

def iqr_outliers(series, factor=1.5):
    q1 = series.quantile(0.25)
    q3 = series.quantile(0.75)
    iqr = q3 - q1
    lower = q1 - factor * iqr
    upper = q3 + factor * iqr
    return series[(series < lower) | (series > upper)].index.tolist()

nat_outlier_idx = iqr_outliers(df_pref['自然増減率'])
soc_outlier_idx = iqr_outliers(df_pref['社会増減率'])
nat_outlier_prefs = df_pref.loc[nat_outlier_idx, 'pref'].tolist()
soc_outlier_prefs = df_pref.loc[soc_outlier_idx, 'pref'].tolist()

print(f"自然増減率の外れ値: {nat_outlier_prefs}")
print(f"社会増減率の外れ値: {soc_outlier_prefs}")

OUTLIERS = list(set(nat_outlier_prefs + soc_outlier_prefs))
print(f"\n外れ値として独立クラスタ扱い: {OUTLIERS}")

# ================================================================
# ■ Step 2. GMM クラスタリング（外れ値除外後の都道府県）
# ================================================================

print("\n" + "=" * 65)
print("■ Step2. ガウス混合モデル（GMM）クラスタリング")
print("=" * 65)

df_main = df_pref[~df_pref['pref'].isin(OUTLIERS)].copy().reset_index(drop=True)
print(f"クラスタリング対象: {len(df_main)}都道府県")

CLUSTER_FEATURES = ['自然増減率', '社会増減率', '高齢化率']
X_gmm = df_main[CLUSTER_FEATURES].values
scaler_gmm = StandardScaler()
X_gmm_scaled = scaler_gmm.fit_transform(X_gmm)

# BIC でコンポーネント数を選択
bics = []
for k in range(1, 6):
    gmm_k = GaussianMixture(n_components=k, covariance_type='full', random_state=2025)
    gmm_k.fit(X_gmm_scaled)
    bics.append(gmm_k.bic(X_gmm_scaled))
    print(f"  k={k}: BIC={gmm_k.bic(X_gmm_scaled):.1f}")

best_k = np.argmin(bics) + 1
print(f"\n最適コンポーネント数（BIC最小）: k={best_k}")

# k=3 で実行（論文に合わせる）
gmm = GaussianMixture(n_components=3, covariance_type='full', random_state=2025)
gmm.fit(X_gmm_scaled)
labels_gmm = gmm.predict(X_gmm_scaled)
probs  = gmm.predict_proba(X_gmm_scaled)

# クラスタを高齢化率順に整列（クラスタ1=低高齢化、クラスタ3=高高齢化）
cluster_aging = {k: X_gmm[labels_gmm == k, 2].mean() for k in range(3)}
order = sorted(cluster_aging, key=lambda x: cluster_aging[x])
remap = {order[0]: 0, order[1]: 1, order[2]: 2}
df_main['cluster'] = [remap[l] + 1 for l in labels_gmm]

print("\nクラスター別プロファイル:")
summary = df_main.groupby('cluster')[['自然増減率','社会増減率','高齢化率']].mean().round(3)
print(summary)
for cl in [1, 2, 3]:
    prefs_cl = df_main[df_main['cluster'] == cl]['pref'].tolist()
    print(f"\nクラスタ{cl}({len(prefs_cl)}都道府県): {', '.join(prefs_cl[:8])}{'...' if len(prefs_cl) > 8 else ''}")

# ================================================================
# ■ Step 3. PLS 回帰分析（実データの変数のみ）
# ================================================================

print("\n" + "=" * 65)
print("■ Step3. PLS回帰分析（クラスタ別）")
print("=" * 65)

# 自然増減率モデルの説明変数: 実データのみ
PLS_FEATURES_NAT = ['高齢化率', '保育所定員率', '医師密度', '婚姻率', '合計特殊出生率']

# 社会増減率モデルの説明変数: 実データのみ
PLS_FEATURES_SOC = ['月間有効求人数（一般）', '大学進学率',
                    '標準価格（平均価格）（住宅地）', '着工新設持家数',
                    '県民所得_1人']

def run_pls(df_sub, X_features, y_col, n_components=2, label=''):
    df_clean = df_sub[X_features + [y_col]].dropna()
    if len(df_clean) < n_components + 2:
        print(f"  {label}: サンプル数不足 (n={len(df_clean)})")
        return 0.0, np.zeros(len(X_features)), None
    X = df_clean[X_features].values
    y = df_clean[y_col].values
    scaler = StandardScaler()
    X_scaled = scaler.fit_transform(X)
    nc = min(n_components, X.shape[1], len(df_clean) - 1)
    pls = PLSRegression(n_components=nc)
    pls.fit(X_scaled, y)
    y_pred = pls.predict(X_scaled).ravel()
    ss_res = np.sum((y - y_pred)**2)
    ss_tot = np.sum((y - y.mean())**2)
    r2 = 1 - ss_res/ss_tot if ss_tot > 0 else 0.0
    coefs = pls.coef_.ravel() if hasattr(pls, 'coef_') else np.zeros(len(X_features))
    print(f"  {label}: R²={r2:.3f}, n={len(df_clean)}")
    return r2, coefs, pls

print("\n【自然増減率へのPLS回帰】")
r2_nat_all, coef_nat_all, _ = run_pls(df_main, PLS_FEATURES_NAT, '自然増減率', label='全体（外れ値除外）')
r2_results = {}
for cl in [1, 2, 3]:
    df_cl = df_main[df_main['cluster'] == cl]
    r2, coefs, _ = run_pls(df_cl, PLS_FEATURES_NAT, '自然増減率',
                            n_components=min(2, len(df_cl)-1), label=f'クラスタ{cl}')
    r2_results[f'nat_cl{cl}'] = r2

print("\n【社会増減率へのPLS回帰】")
r2_soc_all, coef_soc_all, _ = run_pls(df_main, PLS_FEATURES_SOC, '社会増減率', label='全体（外れ値除外）')
for cl in [1, 2, 3]:
    df_cl = df_main[df_main['cluster'] == cl]
    r2, coefs, _ = run_pls(df_cl, PLS_FEATURES_SOC, '社会増減率',
                            n_components=min(2, len(df_cl)-1), label=f'クラスタ{cl}')
    r2_results[f'soc_cl{cl}'] = r2

# ================================================================
# ■ VIF 計算（実データのみ）
# ================================================================

print("\n" + "=" * 65)
print("■ VIF（分散拡大係数）の計算")
print("=" * 65)

from statsmodels.stats.outliers_influence import variance_inflation_factor
import statsmodels.api as sm

def calc_vif(df_sub, features):
    df_clean = df_sub[features].dropna()
    X_vif = (df_clean - df_clean.mean()) / df_clean.std()
    X_vif = sm.add_constant(X_vif)
    vifs = {}
    for i, feat in enumerate(features):
        col_idx = i + 1  # +1 for const
        vif = variance_inflation_factor(X_vif.values, col_idx)
        vifs[feat] = vif
    return vifs

vif_nat = calc_vif(df_main, PLS_FEATURES_NAT)
vif_soc = calc_vif(df_main, PLS_FEATURES_SOC)
print("\n自然増減率 説明変数のVIF:")
for feat, vif in vif_nat.items():
    flag = '★多重共線性の疑い' if vif > 5 else ''
    print(f"  {feat:<25} VIF={vif:.2f} {flag}")
print("\n社会増減率 説明変数のVIF:")
for feat, vif in vif_soc.items():
    flag = '★多重共線性の疑い' if vif > 5 else ''
    print(f"  {feat:<25} VIF={vif:.2f} {flag}")

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

print("\n" + "=" * 65)
print("■ 図の生成（4枚）")
print("=" * 65)

# クラスタ色マップ
C_GMM = {1: '#E53935', 2: '#43A047', 3: '#1E88E5', 'outlier': '#9C27B0'}

# ─── 図1: ボックスプロット（外れ値の可視化）──────────────
print("図1: ボックスプロットによる外れ値検出を作成中...")

fig1, axes1 = plt.subplots(1, 2, figsize=(12, 6))
fig1.suptitle('自然増減率・社会増減率の分布と外れ値（2022年, SSDSE-B実データ）',
              fontsize=13, fontweight='bold')

for ax, col, title in zip(axes1, ['自然増減率','社会増減率'], ['自然増減率（‰）','社会増減率（‰）']):
    bp = ax.boxplot(df_pref[col].dropna().values, patch_artist=True, widths=0.4)
    bp['boxes'][0].set_facecolor('#BBDEFB')
    bp['boxes'][0].set_alpha(0.8)
    for flier in bp['fliers']:
        flier.set(marker='D', color='#E53935', markersize=8, alpha=0.8)

    # 都道府県名をプロット（jitterはなし、x=1に揃える）
    x_jitter = np.linspace(-0.12, 0.12, len(df_pref))  # 決定論的オフセット
    for i, (pref, val) in enumerate(zip(df_pref['pref'], df_pref[col])):
        if pd.isna(val):
            continue
        color = '#E53935' if pref in OUTLIERS else '#1565C0'
        ax.scatter(1 + x_jitter[i], val, color=color, alpha=0.7, s=40, zorder=3)
        if pref in OUTLIERS:
            ax.annotate(pref.replace('県','').replace('都','').replace('道',''),
                       (1 + x_jitter[i], val), fontsize=9, fontweight='bold',
                       color='#E53935',
                       xytext=(8, 0), textcoords='offset points')

    ax.set_ylabel(title, fontsize=11)
    ax.set_title(title, fontsize=11, fontweight='bold')
    ax.set_xticks([])
    ax.grid(axis='y', alpha=0.3)
    q1 = df_pref[col].quantile(0.25)
    q3 = df_pref[col].quantile(0.75)
    iqr_ = q3 - q1
    ax.axhline(q3 + 1.5*iqr_, color='orange', linestyle='--', linewidth=1.2, alpha=0.7, label='IQR×1.5')
    ax.axhline(q1 - 1.5*iqr_, color='orange', linestyle='--', linewidth=1.2, alpha=0.7)
    ax.legend(fontsize=9)

blue_p = mpatches.Patch(color='#1565C0', alpha=0.7, label='通常都道府県')
red_p  = mpatches.Patch(color='#E53935', alpha=0.7, label='外れ値')
fig1.legend(handles=[blue_p, red_p], loc='lower center', ncol=2, fontsize=10, bbox_to_anchor=(0.5, -0.02))

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

# ─── 図2: GMMクラスタリング結果（散布図 + BIC）────────────
print("図2: GMMクラスタリング結果を作成中...")

fig2, axes2 = plt.subplots(1, 2, figsize=(13, 5.5))
fig2.suptitle('ガウス混合モデル（GMM）クラスタリング（外れ値除外後の都道府県）',
              fontsize=13, fontweight='bold')

# (a) BIC曲線
ax2a = axes2[0]
ks_plot = range(1, 6)
ax2a.plot(list(ks_plot), bics, 'o-', color='#1565C0', linewidth=2.2,
          markersize=8, markerfacecolor='white', markeredgewidth=2)
ax2a.axvline(3, color='#E53935', linestyle='--', linewidth=1.5, label='k=3 選択')
ax2a.set_xlabel('コンポーネント数 k', fontsize=11)
ax2a.set_ylabel('BIC（情報量基準）', fontsize=11)
ax2a.set_title('BIC によるコンポーネント数の選択', fontsize=11, fontweight='bold')
ax2a.legend(fontsize=10)
ax2a.grid(True, alpha=0.3)
for k, bic in zip(ks_plot, bics):
    ax2a.annotate(f'{bic:.0f}', (k, bic), textcoords='offset points',
                  xytext=(0, 8), ha='center', fontsize=8.5)

# (b) 散布図（自然増減率 vs 社会増減率）
ax2b = axes2[1]
for cl in [1, 2, 3]:
    mask = df_main['cluster'] == cl
    ax2b.scatter(df_main.loc[mask, '自然増減率'],
                 df_main.loc[mask, '社会増減率'],
                 color=C_GMM[cl], s=80, alpha=0.85,
                 edgecolors='white', linewidth=0.5,
                 label=f'クラスタ{cl}', zorder=3)

# 外れ値を表示
for pref in OUTLIERS:
    row = df_pref[df_pref['pref'] == pref]
    if len(row) == 0:
        continue
    ax2b.scatter(row['自然増減率'].values, row['社会増減率'].values,
                 color=C_GMM['outlier'], s=120, marker='*',
                 edgecolors='white', linewidth=0.5, zorder=4)
    ax2b.annotate(pref.replace('県','').replace('都','').replace('道',''),
                  (row['自然増減率'].values[0], row['社会増減率'].values[0]),
                  fontsize=9, fontweight='bold', color=C_GMM['outlier'],
                  xytext=(6, 4), textcoords='offset points')

ax2b.axhline(0, color='gray', linewidth=0.8, alpha=0.5)
ax2b.axvline(0, color='gray', linewidth=0.8, alpha=0.5)
ax2b.set_xlabel('自然増減率（‰）', fontsize=11)
ax2b.set_ylabel('社会増減率（‰）', fontsize=11)
ax2b.set_title('GMM クラスタリング結果\n（★：外れ値として除外）', fontsize=11, fontweight='bold')
outlier_names = [p.replace('県','').replace('都','') for p in OUTLIERS]
outlier_p = mpatches.Patch(color=C_GMM['outlier'], label=f'外れ値（{", ".join(outlier_names)}）')
handles, labels_leg = ax2b.get_legend_handles_labels()
ax2b.legend(handles=handles+[outlier_p], fontsize=9)
ax2b.grid(True, alpha=0.3)

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

# ─── 図3: PLS回帰係数のクラスタ別比較（バイオリンプロット含む）
print("図3: PLS係数とバイオリンプロットを作成中...")

fig3, axes3 = plt.subplots(1, 2, figsize=(14, 6))
fig3.suptitle('PLS 回帰係数の比較とクラスタ別分布（バイオリンプロット）',
              fontsize=13, fontweight='bold')

# (a) PLS係数（全体、自然増減率）
ax3a = axes3[0]
feat_labels_nat = ['高齢化率', '保育所定員率', '医師密度\n(医師数/人口万)', '婚姻率', '合計特殊出生率']
coef_vals = coef_nat_all[:len(feat_labels_nat)]
bar_colors_nat = ['#E53935' if c < 0 else '#43A047' for c in coef_vals]
ax3a.barh(np.arange(len(feat_labels_nat)), coef_vals,
          color=bar_colors_nat, alpha=0.8, edgecolor='white')
ax3a.axvline(0, color='black', linewidth=1.0)
ax3a.set_yticks(np.arange(len(feat_labels_nat)))
ax3a.set_yticklabels(feat_labels_nat, fontsize=10)
ax3a.set_xlabel('PLS 回帰係数', fontsize=11)
ax3a.set_title(f'PLS 回帰係数（自然増減率）\n全体 R²={r2_nat_all:.3f}', fontsize=11, fontweight='bold')
ax3a.grid(axis='x', alpha=0.3)
ax3a.invert_yaxis()

# (b) バイオリンプロット（クラスタ別 自然増減率分布）
ax3b = axes3[1]
violin_data = [df_main[df_main['cluster'] == cl]['自然増減率'].dropna().values for cl in [1, 2, 3]]
violin_data = [d for d in violin_data if len(d) > 1]
if len(violin_data) >= 2:
    vp = ax3b.violinplot(violin_data, positions=range(1, len(violin_data)+1),
                         showmedians=True, showextrema=True)
    for i, (body, cl) in enumerate(zip(vp['bodies'], [1, 2, 3])):
        body.set_facecolor(C_GMM[cl])
        body.set_alpha(0.7)
    vp['cmedians'].set_color('white')
    vp['cmedians'].set_linewidth(2.5)

ax3b.set_xticks([1, 2, 3])
ax3b.set_xticklabels([f'クラスタ{cl}\n(n={len(df_main[df_main["cluster"]==cl])})' for cl in [1,2,3]],
                     fontsize=10)
ax3b.set_ylabel('自然増減率（‰）', fontsize=11)
ax3b.set_title('クラスタ別 自然増減率の分布\n（バイオリンプロット）', fontsize=11, fontweight='bold')
ax3b.axhline(0, color='gray', linestyle='--', linewidth=0.8, alpha=0.7)
ax3b.grid(axis='y', alpha=0.3)

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

# ─── 図4: クラスタ別PLS R²の比較とVIF ──────────────────────
print("図4: クラスタ別PLS R²とVIF比較を作成中...")

fig4, axes4 = plt.subplots(1, 2, figsize=(13, 5))
fig4.suptitle('クラスタ別 PLS 回帰性能（R²）と VIF（多重共線性確認）',
              fontsize=13, fontweight='bold')

# (a) クラスタ別R²の比較
ax4a = axes4[0]
labels_r2 = ['全体\n(外れ値除外)', 'クラスタ1', 'クラスタ2', 'クラスタ3']
r2_nat_vals = [r2_nat_all, r2_results.get('nat_cl1', 0), r2_results.get('nat_cl2', 0), r2_results.get('nat_cl3', 0)]
r2_soc_vals = [r2_soc_all, r2_results.get('soc_cl1', 0), r2_results.get('soc_cl2', 0), r2_results.get('soc_cl3', 0)]
x_pos = np.arange(len(labels_r2))
width = 0.35
bars1 = ax4a.bar(x_pos - width/2, r2_nat_vals, width, color='#1565C0', alpha=0.8, label='自然増減率モデル')
bars2 = ax4a.bar(x_pos + width/2, r2_soc_vals, width, color='#E65100', alpha=0.8, label='社会増減率モデル')
ax4a.set_xticks(x_pos)
ax4a.set_xticklabels(labels_r2, fontsize=10)
ax4a.set_ylabel('R²（決定係数）', fontsize=11)
ax4a.set_title('PLS 回帰の R²（クラスタ別 vs 全体）\n値が大きいほど良い説明力', fontsize=11, fontweight='bold')
ax4a.legend(fontsize=9)
ax4a.grid(axis='y', alpha=0.3)
ax4a.set_ylim(0, 1.05)
for bars in [bars1, bars2]:
    for bar in bars:
        h = bar.get_height()
        ax4a.text(bar.get_x() + bar.get_width()/2, h + 0.02, f'{h:.2f}',
                  ha='center', fontsize=8.5, fontweight='bold')

# (b) VIF バーチャート
ax4b = axes4[1]
vif_combined = {}
feat_display = {
    '高齢化率': '高齢化率',
    '保育所定員率': '保育所定員率',
    '医師密度': '医師密度',
    '婚姻率': '婚姻率',
    '合計特殊出生率': '合計特殊出生率',
    '月間有効求人数（一般）': '有効求人数',
    '大学進学率': '大学進学率',
    '標準価格（平均価格）（住宅地）': '住宅地価格',
    '着工新設持家数': '持ち家着工数',
    '県民所得_1人': '1人当たり県民所得',
}
for feat in PLS_FEATURES_NAT:
    disp = feat_display.get(feat, feat)
    vif_combined[disp] = vif_nat.get(feat, 0)
for feat in PLS_FEATURES_SOC:
    disp = feat_display.get(feat, feat)
    vif_combined[disp] = vif_soc.get(feat, 0)

vif_features = list(vif_combined.keys())
vif_values   = list(vif_combined.values())
bar_col_vif  = ['#E53935' if v > 5 else '#FF8F00' if v > 3 else '#43A047' for v in vif_values]
ax4b.barh(np.arange(len(vif_features)), vif_values, color=bar_col_vif, alpha=0.8, edgecolor='white')
ax4b.axvline(5, color='#E53935', linestyle='--', linewidth=1.5, label='VIF=5（警戒）')
ax4b.axvline(3, color='#FF8F00', linestyle='--', linewidth=1.0, label='VIF=3（注意）')
ax4b.set_yticks(np.arange(len(vif_features)))
ax4b.set_yticklabels(vif_features, fontsize=9)
ax4b.set_xlabel('VIF（分散拡大係数）', fontsize=11)
ax4b.set_title('VIF による多重共線性確認\n（PLS回帰の必要性の根拠）', fontsize=11, fontweight='bold')
ax4b.legend(fontsize=9)
ax4b.grid(axis='x', alpha=0.3)
ax4b.invert_yaxis()
for i, v in enumerate(vif_values):
    ax4b.text(v + 0.05, i, f'{v:.1f}', va='center', fontsize=8.5)

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

print("\n" + "=" * 65)
print("全図の生成完了（4枚）")
print("=" * 65)
print(f"\n保存先: {FIG_DIR}")
print("  2025_U5_5_box.png  - ボックスプロット（外れ値検出）")
print("  2025_U5_5_gmm.png  - GMMクラスタリング結果")
print("  2025_U5_5_pls.png  - PLS係数とバイオリンプロット")
print("  2025_U5_5_vif.png  - クラスタ別R²とVIF")
print("\n※ 使用データ: SSDSE-B-2026.csv + SSDSE-E-2026.csv（合成データなし）")
print("  自然増減率 = (出生数 - 死亡数) / 総人口 × 1000  [‰]")
print("  社会増減率 = (転入者数 - 転出者数) / 総人口 × 1000  [‰]")