変数重要度に着目したクラスタリングによる
社会構造と睡眠時間の関係性の解析

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

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"    - 大都市圏（仕事時間長）と地方（睡眠時間長）の対照的なパターン")

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 保存完了")

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 保存完了")

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 保存完了")

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

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)} 都道府県")

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

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

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

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 保存完了")