"""
教育用再現コード: 2023年度 統計データ分析コンペティション 審査員奨励賞（大学生）
=================================================================================
論文タイトル：変数重要度に着目したクラスタリングによる社会構造と睡眠時間の関係性の解析
受賞区分  ：審査員奨励賞 ［大学生・一般の部］

【分析概要】
  データ：SSDSE-B-2026（社会指標）+ SSDSE-D-2023（社会生活基本調査、睡眠時間）
         ※ SSDSE-B は 2022年度の47都道府県データを使用
  目的変数：睡眠時間（分/日）= 1日の平均睡眠時間（総数・都道府県別）

  Step1. 実データ読み込み・マージ（SSDSE-B ⊕ SSDSE-D, 地域コードで結合）
  Step2. 社会指標の特徴量エンジニアリング
  Step3. ランダムフォレスト回帰で睡眠時間を予測 → 変数重要度を算出
  Step4. 上位重要特徴量でWard法による階層クラスタリング（47都道府県）
  Step5. クラスタリング結果の可視化と解釈

【使用モデル】
  RandomForestRegressor（sklearn）― 変数重要度（impurity-based）
  Ward法 階層クラスタリング（scipy）

【使用データ】
  SSDSE-B-2026.csv  （統計数理研究所 教育用標準データセット Type B：都道府県）
  SSDSE-D-2023.csv  （統計数理研究所 教育用標準データセット Type D：社会生活基本調査）
  https://www.ism.ac.jp/editsection/ssdse/

【特徴量（社会指標）】
  高齢化率           = 65歳以上人口 / 総人口 × 100
  消費支出           = 消費支出（二人以上の世帯）[円/月]
  住宅地価格         = 標準価格（平均価格）（住宅地）[円/m²]
  有効求人倍率       = 月間有効求人数 / 月間有効求職者数
  TFR（合計特殊出生率）
  大学進学率         = 高等学校卒業者のうち進学者数 / 高等学校卒業者数 × 100
  宿泊者数 per capita = 延べ宿泊者数 / 総人口
  年平均気温（℃）
  仕事時間           = 1日の仕事時間（分）[SSDSE-D]

【データサイエンス学習ポイント】
  1. RandomForest による変数重要度（Feature Importance）の解釈
  2. Ward法（最小分散法）による階層クラスタリングとデンドログラム読み方
  3. 異なる公的統計データセットの結合（地域コードでのマージ）
  4. クラスタ別ボックスプロットによる睡眠時間の差異の可視化

【出力図】
  html/figures/2023_U5_2_fig1_importance.png  ... RF変数重要度棒グラフ
  html/figures/2023_U5_2_fig2_dendrogram.png  ... Ward法デンドログラム（47都道府県）
  html/figures/2023_U5_2_fig3_scatter.png     ... 上位重要変数 vs 睡眠時間（クラスタ色分け）
  html/figures/2023_U5_2_fig4_boxplot.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/2023_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
from sklearn.ensemble import RandomForestRegressor
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import r2_score
from scipy.cluster.hierarchy import dendrogram, linkage, fcluster
import warnings
warnings.filterwarnings('ignore')

# ── パス設定 ─────────────────────────────────────────────────────────────────
_dir = os.path.dirname(os.path.abspath(__file__)) if '__file__' in dir() else os.getcwd()
FIG_DIR = os.path.join(_dir, '..', 'html', 'figures')
DATA_B  = os.path.join(_dir, '..', 'data', 'raw', 'SSDSE-B-2026.csv')
DATA_D  = os.path.join(_dir, '..', 'data', 'raw', 'SSDSE-D-2023.csv')
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("■ 実データ読み込み（SSDSE-B-2026 + SSDSE-D-2023）")
print("=" * 65)

# ── SSDSE-D-2023（社会生活基本調査、47都道府県、総数） ────────────────────
df_d = pd.read_csv(DATA_D, encoding='cp932', header=1)
df_d_total = df_d[
    (df_d['男女の別'] == '0_総数') &
    (df_d['地域コード'] != 'R00000')
].copy()
df_d_total['睡眠_min'] = pd.to_numeric(df_d_total['睡眠'], errors='coerce')
df_d_total['仕事_min'] = pd.to_numeric(df_d_total['仕事'], errors='coerce')

print(f"  SSDSE-D（総数）: {len(df_d_total)} 都道府県")
print(f"  睡眠時間 平均: {df_d_total['睡眠_min'].mean():.1f} 分/日"
      f"  ({df_d_total['睡眠_min'].mean()/60:.2f} 時間/日)")
print(f"  睡眠時間 範囲: {df_d_total['睡眠_min'].min()}〜{df_d_total['睡眠_min'].max()} 分")
print(f"  仕事時間 平均: {df_d_total['仕事_min'].mean():.1f} 分/日")

# ── SSDSE-B-2026（2022年度、47都道府県） ──────────────────────────────────
df_b = pd.read_csv(DATA_B, encoding='cp932', header=1)
df_b = df_b[
    (df_b['地域コード'].str.match(r'^R\d{5}$', na=False)) &
    (df_b['年度'] == 2022)
].copy()

print(f"\n  SSDSE-B（2022年度）: {len(df_b)} 都道府県")

# 数値型変換
num_cols = [
    '総人口', '65歳以上人口',
    '消費支出（二人以上の世帯）',
    '標準価格（平均価格）（住宅地）',
    '月間有効求人数（一般）', '月間有効求職者数（一般）',
    '合計特殊出生率',
    '高等学校卒業者のうち進学者数', '高等学校卒業者数',
    '延べ宿泊者数', '年平均気温',
]
for c in num_cols:
    df_b[c] = pd.to_numeric(df_b[c], errors='coerce')

# ── 地域コードでマージ ─────────────────────────────────────────────────────
df = df_d_total[['地域コード', '都道府県', '睡眠_min', '仕事_min']].merge(
    df_b[['地域コード', '都道府県'] + num_cols],
    on='地域コード', how='inner',
    suffixes=('_d', '_b')
)
df['都道府県'] = df['都道府県_d']
df = df.drop(columns=['都道府県_d', '都道府県_b'])
print(f"\n  マージ後: {len(df)} 都道府県")

# ═══════════════════════════════════════════════════════════════════════════════
# ■ Step2. 特徴量エンジニアリング
# ═══════════════════════════════════════════════════════════════════════════════

print("\n" + "=" * 65)
print("■ Step2. 社会指標の特徴量エンジニアリング")
print("=" * 65)

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

# 消費支出（二人以上の世帯, 円/月）
df['消費支出'] = df['消費支出（二人以上の世帯）']

# 住宅地価格（円/m²）
df['住宅地価格'] = df['標準価格（平均価格）（住宅地）']

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

# 合計特殊出生率（TFR）
df['TFR'] = df['合計特殊出生率']

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

# 宿泊者数 per capita = 延べ宿泊者数 / 総人口
df['宿泊者数per_capita'] = df['延べ宿泊者数'] / df['総人口']

# 年平均気温
df['年平均気温'] = df['年平均気温']

# 仕事時間（SSDSE-Dから）
df['仕事_min'] = df['仕事_min']

# 特徴量列とラベル定義
FEATURE_COLS = [
    '高齢化率',
    '消費支出',
    '住宅地価格',
    '有効求人倍率',
    'TFR',
    '大学進学率',
    '宿泊者数per_capita',
    '年平均気温',
    '仕事_min',
]

FEATURE_LABELS = {
    '高齢化率':           '高齢化率（%）',
    '消費支出':           '消費支出（円/月）',
    '住宅地価格':         '住宅地価格（円/m²）',
    '有効求人倍率':       '有効求人倍率（倍）',
    'TFR':               '合計特殊出生率',
    '大学進学率':         '大学進学率（%）',
    '宿泊者数per_capita': '宿泊者数per capita',
    '年平均気温':         '年平均気温（℃）',
    '仕事_min':          '仕事時間（分/日）',
}

# 欠損値除去
df_model = df[['地域コード', '都道府県', '睡眠_min'] + FEATURE_COLS].dropna().reset_index(drop=True)
print(f"  有効観測数: {len(df_model)} 都道府県")

X = df_model[FEATURE_COLS].values
y = df_model['睡眠_min'].values

print(f"\n  目的変数（睡眠時間）:")
print(f"    平均: {y.mean():.1f} 分/日 ({y.mean()/60:.2f} 時間/日)")
print(f"    SD: {y.std():.1f} 分  範囲: {y.min():.0f}〜{y.max():.0f} 分")

print(f"\n  睡眠時間 上位5都道府県:")
top5 = df_model.nlargest(5, '睡眠_min')[['都道府県', '睡眠_min']]
for _, row in top5.iterrows():
    print(f"    {row['都道府県']}: {row['睡眠_min']:.0f} 分 ({row['睡眠_min']/60:.2f} 時間)")
print(f"  睡眠時間 下位5都道府県:")
bot5 = df_model.nsmallest(5, '睡眠_min')[['都道府県', '睡眠_min']]
for _, row in bot5.iterrows():
    print(f"    {row['都道府県']}: {row['睡眠_min']:.0f} 分 ({row['睡眠_min']/60:.2f} 時間)")

# ═══════════════════════════════════════════════════════════════════════════════
# ■ Step3. ランダムフォレスト回帰 → 変数重要度
# ═══════════════════════════════════════════════════════════════════════════════

print("\n" + "=" * 65)
print("■ Step3. ランダムフォレスト回帰（変数重要度の算出）")
print("=" * 65)

# random_state=42: アルゴリズムの再現性確保（合成データではなく算法固定）
rf = RandomForestRegressor(n_estimators=100, random_state=42)
rf.fit(X, y)

importances = rf.feature_importances_
train_r2 = r2_score(y, rf.predict(X))

print(f"\n  学習用R²（全データ）: {train_r2:.4f}")
print(f"\n  変数重要度（Impurity-based Feature Importance）:")

order = np.argsort(importances)[::-1]
for rank, idx in enumerate(order):
    fn = FEATURE_COLS[idx]
    print(f"  {rank+1:2d}. {FEATURE_LABELS[fn]:<24} {importances[idx]:.4f}")

# 上位特徴量（クラスタリングに使用）
top_n = 4  # 上位4特徴量でクラスタリング
top_idx = order[:top_n]
top_feat_cols = [FEATURE_COLS[i] for i in top_idx]
top_feat_labels = [FEATURE_LABELS[FEATURE_COLS[i]] for i in top_idx]
print(f"\n  クラスタリングに使用する上位{top_n}特徴量:")
for i, (col, lbl) in enumerate(zip(top_feat_cols, top_feat_labels)):
    print(f"    {i+1}. {lbl}")

# ═══════════════════════════════════════════════════════════════════════════════
# ■ Step4. Ward法 階層クラスタリング
# ═══════════════════════════════════════════════════════════════════════════════

print("\n" + "=" * 65)
print("■ Step4. Ward法 階層クラスタリング（47都道府県）")
print("=" * 65)

# 上位重要特徴量で標準化
X_top = df_model[top_feat_cols].values
scaler = StandardScaler()
X_top_sc = scaler.fit_transform(X_top)

# Ward法（最小分散法）でリンケージ行列を計算
Z = linkage(X_top_sc, method='ward', metric='euclidean')

# クラスタ数 = 4 に分割
N_CLUSTERS = 4
cluster_labels = fcluster(Z, N_CLUSTERS, criterion='maxclust')
df_model['cluster'] = cluster_labels

print(f"  クラスタ数: {N_CLUSTERS}")
for c in sorted(df_model['cluster'].unique()):
    prefs = df_model[df_model['cluster'] == c]['都道府県'].tolist()
    mean_sleep = df_model[df_model['cluster'] == c]['睡眠_min'].mean()
    print(f"\n  クラスタ {c}（睡眠時間平均: {mean_sleep:.1f} 分）:")
    print(f"    {', '.join(prefs)}")

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

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

# ────────────────────────────────────────────────────────────────────────────
# 図1：RF変数重要度棒グラフ
# ────────────────────────────────────────────────────────────────────────────
print("\n図1: RF変数重要度棒グラフを作成中...")

fig1, ax1 = plt.subplots(figsize=(10, 6))
fig1.suptitle('ランダムフォレストによる変数重要度\n（目的変数：都道府県別 1日平均睡眠時間）',
              fontsize=13, fontweight='bold')

imp_sorted    = importances[order]
labels_sorted = [FEATURE_LABELS[FEATURE_COLS[i]] for i in order]

colors_bar = ['#E53935' if i < top_n else '#78909C' for i in range(len(order))]

bars = ax1.barh(np.arange(len(order)), imp_sorted,
                color=colors_bar, alpha=0.85, edgecolor='white', linewidth=0.7)

# 重要度の数値ラベル
for bar, val in zip(bars, imp_sorted):
    ax1.text(val + 0.002, bar.get_y() + bar.get_height() / 2,
             f'{val:.4f}', va='center', fontsize=9, color='#333')

ax1.set_yticks(np.arange(len(order)))
ax1.set_yticklabels(labels_sorted, fontsize=10)
ax1.set_xlabel('変数重要度（不純度減少の平均, Gini Importance）', fontsize=11)
ax1.set_title(f'上位{top_n}変数（赤）をクラスタリングに使用', fontsize=10, color='#555')
ax1.invert_yaxis()
ax1.grid(axis='x', alpha=0.3)

# 凡例
patch_top  = mpatches.Patch(color='#E53935', alpha=0.85, label=f'上位{top_n}変数（クラスタリング使用）')
patch_rest = mpatches.Patch(color='#78909C', alpha=0.85, label='その他の変数')
ax1.legend(handles=[patch_top, patch_rest], fontsize=9, loc='lower right')

ax1.text(0.98, 0.02, f'RF: n_estimators=100\nTrain R²={train_r2:.3f}',
         transform=ax1.transAxes, fontsize=8.5, ha='right', va='bottom',
         bbox=dict(boxstyle='round', facecolor='white', alpha=0.85, edgecolor='#aaa'))

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

# ────────────────────────────────────────────────────────────────────────────
# 図2：Ward法デンドログラム（47都道府県）
# ────────────────────────────────────────────────────────────────────────────
print("図2: Ward法デンドログラムを作成中...")

fig2, ax2 = plt.subplots(figsize=(16, 8))
fig2.suptitle(f'Ward法 階層クラスタリング デンドログラム\n'
              f'（使用特徴量: {", ".join(top_feat_labels[:2])} 他計{top_n}変数）',
              fontsize=13, fontweight='bold')

pref_names = df_model['都道府県'].tolist()

# クラスタ境界色を計算するためのリーフ順序
dn = dendrogram(
    Z,
    labels=pref_names,
    orientation='bottom',
    color_threshold=Z[-N_CLUSTERS + 1, 2],
    above_threshold_color='#555',
    ax=ax2,
    leaf_font_size=8,
    leaf_rotation=90,
)

ax2.set_xlabel('都道府県', fontsize=11)
ax2.set_ylabel('Ward距離（不純度増加量）', fontsize=11)
ax2.axhline(y=Z[-N_CLUSTERS + 1, 2], color='#E53935', linestyle='--',
            linewidth=2, label=f'カット線（クラスタ数={N_CLUSTERS}）')
ax2.legend(fontsize=10, loc='upper left')
ax2.grid(axis='y', alpha=0.2)

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

# ────────────────────────────────────────────────────────────────────────────
# 図3：最重要変数 vs 睡眠時間（クラスタ色分け散布図）
# ────────────────────────────────────────────────────────────────────────────
print("図3: 散布図（最重要変数 vs 睡眠時間）を作成中...")

top1_col   = top_feat_cols[0]
top1_label = top_feat_labels[0]
top2_col   = top_feat_cols[1] if len(top_feat_cols) > 1 else top_feat_cols[0]
top2_label = top_feat_labels[1] if len(top_feat_labels) > 1 else top_feat_labels[0]

fig3, axes3 = plt.subplots(1, 2, figsize=(14, 6))
fig3.suptitle('社会指標（上位重要変数）と睡眠時間の関係\n（点の色 = Ward法クラスタ）',
              fontsize=13, fontweight='bold')

for ax, x_col, x_label in zip(axes3, [top1_col, top2_col], [top1_label, top2_label]):
    for c in sorted(df_model['cluster'].unique()):
        mask = df_model['cluster'] == c
        ax.scatter(df_model.loc[mask, x_col],
                   df_model.loc[mask, '睡眠_min'],
                   color=CLUSTER_COLORS[c], alpha=0.85, s=70,
                   edgecolors='white', linewidth=0.7,
                   label=f'クラスタ {c}', zorder=3)

    # 全体回帰直線
    coef = np.polyfit(df_model[x_col], df_model['睡眠_min'], 1)
    x_range = np.linspace(df_model[x_col].min(), df_model[x_col].max(), 200)
    ax.plot(x_range, np.polyval(coef, x_range), 'k--', linewidth=1.5,
            alpha=0.7, label='全体回帰直線')

    # 都道府県ラベル（上位・下位）
    q_high = df_model['睡眠_min'].quantile(0.85)
    q_low  = df_model['睡眠_min'].quantile(0.15)
    for _, row in df_model.iterrows():
        if row['睡眠_min'] >= q_high or row['睡眠_min'] <= q_low:
            ax.annotate(row['都道府県'][:2],
                        (row[x_col], row['睡眠_min']),
                        fontsize=7.5, ha='left', va='bottom', color='#333')

    ax.set_xlabel(x_label, fontsize=11)
    ax.set_ylabel('1日平均睡眠時間（分）', fontsize=11)
    ax.set_title(f'{x_label}\nvs 睡眠時間', fontsize=11, fontweight='bold')
    ax.legend(fontsize=9)
    ax.grid(True, alpha=0.25)

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

# ────────────────────────────────────────────────────────────────────────────
# 図4：クラスタ別 睡眠時間ボックスプロット
# ────────────────────────────────────────────────────────────────────────────
print("図4: クラスタ別睡眠時間ボックスプロットを作成中...")

fig4, axes4 = plt.subplots(1, 2, figsize=(14, 6))
fig4.suptitle('クラスタ別 1日平均睡眠時間の分布\n（Ward法 階層クラスタリング, 47都道府県）',
              fontsize=13, fontweight='bold')

ax4a = axes4[0]
cluster_data = [df_model[df_model['cluster'] == c]['睡眠_min'].values
                for c in sorted(df_model['cluster'].unique())]
cluster_means = [d.mean() for d in cluster_data]
cluster_ns    = [len(d)   for d in cluster_data]

bp = ax4a.boxplot(cluster_data, patch_artist=True, notch=False,
                   medianprops=dict(color='white', linewidth=2.5),
                   whiskerprops=dict(linewidth=1.5),
                   capprops=dict(linewidth=1.5),
                   flierprops=dict(marker='o', markersize=5, alpha=0.6))

box_colors = [CLUSTER_COLORS[c] for c in sorted(df_model['cluster'].unique())]
for patch, color in zip(bp['boxes'], box_colors):
    patch.set_facecolor(color)
    patch.set_alpha(0.75)

# 平均値マーカー
for i, (mean_val, color) in enumerate(zip(cluster_means, box_colors)):
    ax4a.scatter(i + 1, mean_val, color=color, s=80, zorder=5,
                 marker='D', edgecolors='white', linewidth=1.5,
                 label=f'クラスタ{i+1} 平均: {mean_val:.1f}分')

ax4a.set_xticks(range(1, N_CLUSTERS + 1))
ax4a.set_xticklabels([f'クラスタ {c}\n(n={cluster_ns[c-1]})' for c in range(1, N_CLUSTERS + 1)],
                      fontsize=10)
ax4a.set_ylabel('1日平均睡眠時間（分）', fontsize=11)
ax4a.set_title('クラスタ別 睡眠時間分布\n（ひし形=平均値）', fontsize=11, fontweight='bold')
ax4a.legend(fontsize=9, loc='upper right')
ax4a.grid(axis='y', alpha=0.3)

# 全国平均線
national_mean = df_model['睡眠_min'].mean()
ax4a.axhline(national_mean, color='#333', linestyle=':', linewidth=1.5,
             label=f'全国平均: {national_mean:.1f}分')

# パネル(b): クラスタ別 社会指標プロファイル（レーダー不使用、棒グラフで）
ax4b = axes4[1]
feat_for_profile = top_feat_cols  # 上位4変数を使ってプロファイル比較

cluster_means_feat = []
for c in sorted(df_model['cluster'].unique()):
    means = df_model[df_model['cluster'] == c][feat_for_profile].mean().values
    cluster_means_feat.append(means)

# 標準化して比較（各特徴量をzスコア化）
cluster_means_arr = np.array(cluster_means_feat)
feat_overall_mean = df_model[feat_for_profile].mean().values
feat_overall_std  = df_model[feat_for_profile].std().values + 1e-9
cluster_means_z   = (cluster_means_arr - feat_overall_mean) / feat_overall_std

x_pos   = np.arange(top_n)
bar_w   = 0.2
offsets = [-1.5 * bar_w, -0.5 * bar_w, 0.5 * bar_w, 1.5 * bar_w]

for ci, (c, offset) in enumerate(zip(sorted(df_model['cluster'].unique()), offsets)):
    n_pref = cluster_ns[ci]
    ax4b.bar(x_pos + offset, cluster_means_z[ci],
             width=bar_w, color=CLUSTER_COLORS[c], alpha=0.8,
             label=f'クラスタ {c}（n={n_pref}）', edgecolor='white')

ax4b.set_xticks(x_pos)
ax4b.set_xticklabels(top_feat_labels, fontsize=9, rotation=15, ha='right')
ax4b.set_ylabel('zスコア（全体平均=0）', fontsize=11)
ax4b.set_title('クラスタ別 社会指標プロファイル\n（上位重要変数のzスコア比較）',
               fontsize=11, fontweight='bold')
ax4b.axhline(0, color='gray', linewidth=0.8, linestyle='--')
ax4b.legend(fontsize=9)
ax4b.grid(axis='y', alpha=0.3)

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

# ═══════════════════════════════════════════════════════════════════════════════
# ■ 最終サマリー
# ═══════════════════════════════════════════════════════════════════════════════
print("\n" + "=" * 65)
print("■ 全図の生成完了（4枚）")
print("=" * 65)
print(f"\n  データ: SSDSE-B-2026（2022年度）+ SSDSE-D-2023（実公的統計データ）")
print(f"  総観測数: {len(df_model)} 都道府県")
print(f"  目的変数: 睡眠時間 平均 {y.mean():.1f} 分/日"
      f" ({y.mean()/60:.2f} 時間/日)  SD: {y.std():.1f} 分")

print(f"\n  【変数重要度 Top5（RF, Impurity-based）】")
for rank, idx in enumerate(order[:5]):
    fn = FEATURE_COLS[idx]
    print(f"    {rank+1}. {FEATURE_LABELS[fn]:<24} {importances[idx]:.4f}")

print(f"\n  【クラスタ別 睡眠時間平均】")
for c in sorted(df_model['cluster'].unique()):
    m = df_model[df_model['cluster'] == c]['睡眠_min'].mean()
    n = (df_model['cluster'] == c).sum()
    print(f"    クラスタ {c} (n={n:2d}): {m:.1f} 分 = {m/60:.2f} 時間/日")

print(f"\n  【主要知見】")
print(f"    - 仕事時間（SSDSE-D）が睡眠時間予測の最重要因子の一つ")
print(f"    - 高齢化率・消費支出などの社会構造変数も睡眠時間と関連")
print(f"    - Ward法クラスタリングで社会構造が類似する都道府県群が睡眠パターンを共有")
print(f"    - 大都市圏（仕事時間長）と地方（睡眠時間長）の対照的なパターン")