クラスタリングを用いた
男子中学生の持久力低下の要因特定

data_file = 'data/2025_U3/2025_U3_panel.csv'
df_raw = pd.read_csv(data_file)
# linearmodels 用：都道府県名を数値IDに変換（MultiIndex の entity）
df_raw['pref_id'] = pd.Categorical(df_raw['pref']).codes + 1

print("=" * 65)
print("■ パネルデータ（実データ：N=47都道府県 × T=3年）")
print("  出典: スポーツ庁体力調査 / SSDSE-B / 愛媛県オープンデータ")
print("=" * 65)
print(df_raw[['shuttle', 'bmi_rate', 'stu_per_school', 'sports_club', 'breakfast',
              'avg_temp', 'max_temp']].describe().round(3))

# Z-score は data_prep.py で既に計算済み（{変数名}_z 列として保存）
SCALE_VARS = ['avg_temp', 'max_temp', 'rain_days', 'rain_vol',
              'bmi_rate', 'stu_per_school', 'sports_club', 'breakfast']
Z_VARS = [v + '_z' for v in SCALE_VARS]

VAR_LABELS = {
    'avg_temp_z':        '年平均気温',
    'max_temp_z':        '最高気温',
    'rain_days_z':       '降水日数（年間）',
    'rain_vol_z':        '降水量（年間）',
    'bmi_rate_z':        '県肥満率',
    'stu_per_school_z':  '1学校当たり生徒数',
    'sports_club_z':     '運動部加入率',
    'breakfast_z':       '朝食を毎日食べる割合',
}

print("図1: シャトルラン時系列推移グラフを作成中...")

# スポーツ庁「全国体力・運動能力、運動習慣等調査」の公表集計値
# 出典: https://www.mext.go.jp/sports/b_menu/toukei/kodomo/zencyo/1368222_00002.htm
# 2020年はコロナにより調査実施なし（None）
hist_years   = [2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015,
                2016, 2017, 2018, 2019, 2020, 2021, 2022]
shuttle_vals = [85.8, 85.4, 85.7, 84.4, 84.7, 84.7, 84.5, 85.0,
                85.3, 85.5, 85.8, 83.7, None, 79.5, 77.8]

# 都道府県別・年度別 シャトルラン平均（実パネルデータから計算）
pref_year_avg = df_raw.groupby(['pref', 'year'])['shuttle'].mean().reset_index()

fig1, axes1 = plt.subplots(1, 2, figsize=(13, 5))
fig1.suptitle('男子中学生の持久力テスト結果の推移\n（出典：スポーツ庁 全国体力・運動能力、運動習慣等調査）',
              fontsize=12, fontweight='bold')

# (a) 全国平均の時系列（公表値）
years_v = [y for y, v in zip(hist_years, shuttle_vals) if v is not None]
vals_v  = [v for v in shuttle_vals if v is not None]
ax1a = axes1[0]
ax1a.plot(years_v, vals_v, 'o-', color='#1565C0', linewidth=2.2, markersize=6,
          markerfacecolor='white', markeredgewidth=2, label='20mシャトルラン（全国平均）')
ax1a.axvspan(2020, 2022, alpha=0.12, color='red', label='コロナ禍影響期')
ax1a.axhline(vals_v[0], color='gray', linestyle='--', linewidth=0.8, alpha=0.6,
             label=f'2008年水準 ({vals_v[0]:.1f}回)')
ax1a.set_xlabel('年度', fontsize=10)
ax1a.set_ylabel('平均回数（回）', fontsize=10)
ax1a.set_title('全国平均 時系列推移（公表値）', fontsize=11, fontweight='bold')
ax1a.legend(fontsize=8.5)
ax1a.grid(True, alpha=0.3)
ax1a.set_ylim(74, 90)
for y, v in zip(years_v[-3:], vals_v[-3:]):
    ax1a.annotate(f'{v:.1f}', (y, v), textcoords='offset points', xytext=(0, 8),
                  fontsize=8.5, ha='center', color='#C62828', fontweight='bold')

# (b) 都道府県別の分布（実パネルデータ）
ax1b = axes1[1]
pref_colors = plt.cm.tab20(np.linspace(0, 1, N))
nat_mean = pref_year_avg.groupby('year')['shuttle'].mean()
for i, pref in enumerate(df_raw['pref'].unique()):
    pref_d = pref_year_avg[pref_year_avg['pref'] == pref].sort_values('year')
    ax1b.plot(pref_d['year'], pref_d['shuttle'],
              color=pref_colors[i % 20], alpha=0.3, linewidth=0.9)
ax1b.plot(nat_mean.index, nat_mean.values,
          color='black', linewidth=2.5, linestyle='--', label='47都道府県平均', zorder=5)
ax1b.set_xlabel('年度', fontsize=10)
ax1b.set_ylabel('シャトルラン平均回数（回）', fontsize=10)
ax1b.set_title('都道府県別 推移（実データ）', fontsize=11, fontweight='bold')
ax1b.legend(fontsize=9)
ax1b.grid(True, alpha=0.3)
ax1b.set_xticks(YEARS)

plt.tight_layout()
save_fig('fig1_trend')

print("図2: 固定効果モデル選択と回帰係数グラフを作成中...")

fig2, axes2 = plt.subplots(1, 3, figsize=(16, 5))
fig2.suptitle('Step1〜2. パネルモデル選択（F検定・Hausman検定）と固定効果モデルの係数\n'
              '被説明変数：男子中学生 20mシャトルラン平均回数（実データ）',
              fontsize=11, fontweight='bold')

# (a) F統計量の棒グラフ
ax2a = axes2[0]
f_labels    = ['Pooled OLS\nvs 時間固定効果', 'Pooled OLS\nvs 個体固定効果', 'Pooled OLS\nvs two-way FE']
f_vals_ref  = [f_time, f_ind, f_twoway]
p_vals_ref  = [p_time, p_ind, p_twoway]
bar_colors  = ['#E53935' if p < 0.05 else '#9E9E9E' for p in p_vals_ref]
ax2a.barh(np.arange(3), f_vals_ref, color=bar_colors, alpha=0.85, edgecolor='white')
crit_f = stats.f.ppf(0.95, 2, 130)
ax2a.axvline(crit_f, color='orange', linestyle='--', linewidth=1.5, label=f'5%臨界値 ({crit_f:.2f})')
ax2a.set_yticks(np.arange(3))
ax2a.set_yticklabels(f_labels, fontsize=9)
ax2a.set_xlabel('F統計量', fontsize=10)
ax2a.set_title('F検定：固定効果の有意性', fontsize=10, fontweight='bold')
ax2a.legend(fontsize=8)
ax2a.grid(axis='x', alpha=0.3)
ax2a.invert_yaxis()
for i, (fv, pv) in enumerate(zip(f_vals_ref, p_vals_ref)):
    sig = '***' if pv < 0.001 else '**' if pv < 0.01 else '*' if pv < 0.05 else 'n.s.'
    ax2a.text(fv + max(f_vals_ref) * 0.02, i, f'F={fv:.1f} {sig}',
              va='center', fontsize=8.5, fontweight='bold')
red_patch  = mpatches.Patch(color='#E53935', label='有意(p<0.05)')
gray_patch = mpatches.Patch(color='#9E9E9E', label='非有意')
ax2a.legend(handles=[red_patch, gray_patch], fontsize=8, loc='lower right')

# (b) 決定係数の比較 + Hausman検定
ax2b = axes2[1]
model_names = ['Pooled OLS', '時間FE\n(論文:0.39)', '個体FE\n(論文:0.75)', 'Two-way FE']
r2_vals     = [ols_pooled.rsquared, ols_time_fe.rsquared, ols_ind_fe.rsquared, ols_twoway.rsquared]
r2_colors   = ['#9E9E9E', '#FF8F00', '#6A1B9A', '#757575']
ax2b.bar(np.arange(4), r2_vals, color=r2_colors, alpha=0.85, edgecolor='white')
ax2b.set_xticks(np.arange(4))
ax2b.set_xticklabels(model_names, fontsize=8.5)
ax2b.set_ylabel('決定係数 R²', fontsize=10)
ax2b.set_title('各モデルの決定係数', fontsize=10, fontweight='bold')
ax2b.grid(axis='y', alpha=0.3)
ax2b.set_ylim(0, 1.0)
for i, v in enumerate(r2_vals):
    ax2b.text(i, v + 0.02, f'{v:.3f}', ha='center', fontsize=9, fontweight='bold')
hausman_note = (f"Hausman検定\n"
                f"χ²={H_stat:.2f}, p={p_H:.3f}\n"
                f"→ {'個体FE採用' if p_H < 0.05 else 'RE採用'}")
ax2b.text(0.97, 0.04, hausman_note, transform=ax2b.transAxes, fontsize=8.5,
          va='bottom', ha='right',
          bbox=dict(boxstyle='round', facecolor='#FFF9C4', alpha=0.9, edgecolor='#F9A825'))

# (c) 個体固定効果モデルの回帰係数（95% CI付き）
ax2c = axes2[2]
y_pos2 = np.arange(len(Z_VARS))
coef_i_list = [coef_ind[zv]     for zv in Z_VARS]
pval_i_list = [pval_ind[zv]     for zv in Z_VARS]
se_i_list   = [ols_ind_fe.bse[zv] for zv in Z_VARS]
ci_lo_list  = [c - 1.96 * s for c, s in zip(coef_i_list, se_i_list)]
ci_hi_list  = [c + 1.96 * s for c, s in zip(coef_i_list, se_i_list)]
bar_c2      = ['#6A1B9A' if p < 0.05 else '#BDBDBD' for p in pval_i_list]

ax2c.barh(y_pos2, coef_i_list, color=bar_c2, alpha=0.8, height=0.6)
ax2c.errorbar(coef_i_list, y_pos2,
              xerr=[np.array(coef_i_list) - np.array(ci_lo_list),
                    np.array(ci_hi_list)  - np.array(coef_i_list)],
              fmt='none', color='black', capsize=3, linewidth=1)
ax2c.axvline(0, color='black', linewidth=0.8, linestyle='--')
ax2c.set_yticks(y_pos2)
ax2c.set_yticklabels([VAR_LABELS[zv] for zv in Z_VARS], fontsize=9)
ax2c.set_xlabel('標準化係数（95% CI）', fontsize=9)
ax2c.set_title(f'個体固定効果モデル\nR²={ols_ind_fe.rsquared:.3f}', fontsize=10, fontweight='bold')
ax2c.grid(axis='x', alpha=0.3)
ax2c.invert_yaxis()
for i, (c, s, p) in enumerate(zip(coef_i_list, se_i_list, pval_i_list)):
    if p < 0.05:
        sig = '***' if p < 0.001 else '**' if p < 0.01 else '*'
        ax2c.text(c + 1.96 * s + 0.05, i, sig, va='center', fontsize=9, color='#4A148C')

plt.tight_layout()
save_fig('fig2_fe')

import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
import matplotlib.patches as mpatches
import statsmodels.api as sm
from linearmodels.panel import PooledOLS, PanelOLS, RandomEffects
from sklearn.cluster import KMeans
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import silhouette_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
FIGURE_DIR = 'html/figures'
os.makedirs(FIGURE_DIR, exist_ok=True)

def save_fig(name):
    path = os.path.join(FIGURE_DIR, f'2025_U3_{name}.png')
    plt.savefig(path, bbox_inches='tight', dpi=150)
    plt.close()
    print(f"  → {path} 保存完了")

C_CLUSTER = {1: '#E53935', 2: '#43A047', 3: '#1E88E5'}
C_LABELS  = {1: 'Cluster1: 標準型', 2: 'Cluster2: 生活習慣課題型', 3: 'Cluster3: 食文化・環境要因型'}

YEARS = [2019, 2021, 2022]
N     = 47
T     = len(YEARS)

print("\n" + "=" * 65)
print("■ Step1. モデル選択（F検定 + Hausman検定）")
print("=" * 65)

# ── Pooled OLS（statsmodels） ────────────────────────────────
df_flat = df_raw.copy()
X_ols_df = sm.add_constant(df_flat[Z_VARS])
ols_pooled = sm.OLS(df_flat['shuttle'], X_ols_df).fit()

# ── 個体固定効果モデル（明示的ダミー変数）───────────────────
pref_dummies = pd.get_dummies(df_flat['pref_id'], prefix='pref', drop_first=True).astype(float)
X_ind_fe_df = pd.concat([df_flat[Z_VARS].reset_index(drop=True),
                          pref_dummies.reset_index(drop=True)], axis=1)
X_ind_fe_df = sm.add_constant(X_ind_fe_df)
ols_ind_fe = sm.OLS(df_flat['shuttle'], X_ind_fe_df).fit()

# ── 時間固定効果モデル（明示的ダミー変数）───────────────────
year_dummies = pd.get_dummies(df_flat['year'], prefix='yr', drop_first=True).astype(float)
X_time_fe_df = pd.concat([df_flat[Z_VARS].reset_index(drop=True),
                           year_dummies.reset_index(drop=True)], axis=1)
X_time_fe_df = sm.add_constant(X_time_fe_df)
ols_time_fe = sm.OLS(df_flat['shuttle'], X_time_fe_df).fit()

# ── Two-way 固定効果モデル ────────────────────────────────────
X_twoway_df = pd.concat([df_flat[Z_VARS].reset_index(drop=True),
                          pref_dummies.reset_index(drop=True),
                          year_dummies.reset_index(drop=True)], axis=1)
X_twoway_df = sm.add_constant(X_twoway_df)
ols_twoway = sm.OLS(df_flat['shuttle'], X_twoway_df).fit()

# ── F検定 ────────────────────────────────────────────────────
def panel_f_test(ols_restricted, ols_unrestricted, df_added):
    """
    F検定: 固定効果ダミーの結合有意性
    H₀: 追加ダミー変数の係数がすべてゼロ（固定効果不要）
    """
    rss_r = ols_restricted.ssr
    rss_u = ols_unrestricted.ssr
    df_num = df_added
    df_den = ols_unrestricted.df_resid
    if rss_r <= rss_u or df_num <= 0:
        return 0.0, 1.0
    f_stat = ((rss_r - rss_u) / df_num) / (rss_u / df_den)
    p_val  = 1 - stats.f.cdf(f_stat, df_num, df_den)
    return float(f_stat), float(p_val)

f_time,  p_time  = panel_f_test(ols_pooled, ols_time_fe, df_added=T-1)
f_ind,   p_ind   = panel_f_test(ols_pooled, ols_ind_fe,  df_added=N-1)
f_twoway,p_twoway= panel_f_test(ols_pooled, ols_twoway,  df_added=N+T-2)

print(f"\n【F検定結果（Pooled OLS との比較）】")
print(f"  {'比較':<30} {'F統計量':>9} {'p値':>12} {'判定':>8}")
print("  " + "-" * 65)
print(f"  {'Pooled OLS vs 時間固定効果':<30} {f_time:>9.3f} {p_time:>12.2e}"
      f"  {'有意*' if p_time < 0.05 else '非有意':>8}")
print(f"  {'Pooled OLS vs 個体固定効果':<30} {f_ind:>9.3f} {p_ind:>12.2e}"
      f"  {'有意*' if p_ind < 0.05 else '非有意':>8}")
print(f"  {'Pooled OLS vs two-way 固定効果':<30} {f_twoway:>9.3f} {p_twoway:>12.4f}"
      f"  {'有意*' if p_twoway < 0.05 else '非有意':>8}")
print(f"\n  論文の参考値: 時間FE F=10.269, 個体FE F=31.431, Two-way F=1.011")

# ── Hausman検定（linearmodels の FE vs RE）───────────────────
df_panel = df_raw.set_index(['pref_id', 'year'])
y_p  = df_panel['shuttle']
X_p  = sm.add_constant(df_panel[Z_VARS])

fe_res = PanelOLS(y_p, X_p, entity_effects=True).fit(cov_type='unadjusted')
re_res = RandomEffects(y_p, X_p).fit(cov_type='unadjusted')

def hausman_test(fe_result, re_result):
    """
    Hausman検定: FE と RE の係数の差を利用した検定
    H₀: 固有効果と説明変数が無相関（RE が一致・有効）
    H₁: 相関あり（RE はバイアスを持つ → FE を使うべき）
    """
    fe_coef = fe_result.params
    re_coef = re_result.params
    shared  = fe_coef.index.intersection(re_coef.index)
    diff    = (fe_coef - re_coef)[shared].values

    V_fe = pd.DataFrame(fe_result.cov.values, index=fe_result.params.index,
                         columns=fe_result.params.index).loc[shared, shared].values
    V_re = pd.DataFrame(re_result.cov.values, index=re_result.params.index,
                         columns=re_result.params.index).loc[shared, shared].values
    V_diff = V_fe - V_re

    eigvals, eigvecs = np.linalg.eigh(V_diff)
    pos_mask = eigvals > 1e-10
    if pos_mask.sum() == 0:
        V_inv = np.linalg.pinv(V_diff)
    else:
        V_inv = (eigvecs[:, pos_mask]
                 @ np.diag(1.0 / eigvals[pos_mask])
                 @ eigvecs[:, pos_mask].T)

    H_stat = float(diff @ V_inv @ diff)
    df_H   = int(pos_mask.sum())
    p_H    = 1 - stats.chi2.cdf(H_stat, df_H)
    return H_stat, df_H, p_H

H_stat, df_H, p_H = hausman_test(fe_res, re_res)
sig_h = '***' if p_H < 0.01 else '**' if p_H < 0.05 else '*' if p_H < 0.1 else 'n.s.'
print(f"\n【Hausman検定（個体固定効果 vs 変量効果）】")
print(f"  χ²={H_stat:.3f}, df={df_H}, p={p_H:.4f} {sig_h}")
if p_H < 0.05:
    print("  → 固定効果モデルを採用（個体効果と説明変数に相関がある）")
else:
    print("  → 変量効果モデルが一致・有効（固有効果と説明変数に相関なし）")

print(f"  論文の参考値: Hausman χ²≈17.1, p=0.047")

# ── モデル決定係数のまとめ ───────────────────────────────────
print(f"\n【決定係数のまとめ】")
print(f"  Pooled OLS   : R²={ols_pooled.rsquared:.4f}")
print(f"  時間固定効果 : R²={ols_time_fe.rsquared:.4f}  (論文 0.3872)")
print(f"  個体固定効果 : R²={ols_ind_fe.rsquared:.4f}  (論文 0.7451)")
print(f"  Two-way FE   : R²={ols_twoway.rsquared:.4f}")

print("\n" + "=" * 65)
print("■ Step2. 回帰分析結果（時間FE vs 個体FE）")
print("=" * 65)

print(f"\n  {'説明変数':<22} {'時間FE 係数':>10} {'時間FE p':>9}"
      f" {'個体FE 係数':>12} {'個体FE p':>9} {'有意':>6}")
print("  " + "-" * 75)

coef_time = dict(zip(Z_VARS, ols_time_fe.params[Z_VARS]))
pval_time = dict(zip(Z_VARS, ols_time_fe.pvalues[Z_VARS]))
coef_ind  = dict(zip(Z_VARS, ols_ind_fe.params[Z_VARS]))
pval_ind  = dict(zip(Z_VARS, ols_ind_fe.pvalues[Z_VARS]))

for zv in Z_VARS:
    ct  = coef_time[zv]
    pt  = pval_time[zv]
    ci  = coef_ind[zv]
    pi  = pval_ind[zv]
    sig_t = '***' if pt < 0.01 else '**' if pt < 0.05 else '*' if pt < 0.1 else ''
    sig_i = '***' if pi < 0.01 else '**' if pi < 0.05 else '*' if pi < 0.1 else ''
    sig   = sig_t or sig_i
    vn    = VAR_LABELS[zv]
    print(f"  {vn:<22} {ct:>10.4f} {pt:>9.4f}"
          f" {ci:>12.4f} {pi:>9.4f} {sig:>6}")

print("\n" + "=" * 65)
print("■ Step3. K-meansクラスタリング（47都道府県を3類型に分類）")
print("=" * 65)

# 2022年データを使用（論文: 最新年度のデータでクラスタリング）
df_2022 = df_raw[df_raw['year'] == 2022].copy().reset_index(drop=True)

# 有意変数（論文の分析で特定された変数）でクラスタリング
CLUSTER_VARS = ['bmi_rate', 'sports_club', 'breakfast']
X_cluster    = df_2022[CLUSTER_VARS].values

# 標準化（KMeans 用）
scaler   = StandardScaler()
X_scaled = scaler.fit_transform(X_cluster)

# ── エルボー法（k=2〜6）────────────────────────────────────
print(f"\n【エルボー法とシルエットスコア（k=2〜6）】")
print(f"  {'k':>4} {'Inertia':>12} {'Silhouette Score':>18}")
print("  " + "-" * 36)
inertias, silhouettes = [], []
for k in range(2, 7):
    km = KMeans(n_clusters=k, random_state=2025, n_init=20)
    km.fit(X_scaled)
    inertias.append(km.inertia_)
    sil = silhouette_score(X_scaled, km.labels_)
    silhouettes.append(sil)
    print(f"  {k:>4} {km.inertia_:>12.2f} {sil:>18.4f}")

# ── k=3でクラスタリング ─────────────────────────────────────
kmeans = KMeans(n_clusters=3, random_state=2025, n_init=30)
kmeans.fit(X_scaled)
df_2022['cluster_raw'] = kmeans.labels_

# クラスター番号を論文と対応：シャトルラン平均で整列
shuttle_by_cluster = df_2022.groupby('cluster_raw')['shuttle'].mean().sort_values(ascending=False)
# 最高=Cluster1、最低=Cluster2、中間=Cluster3
rank_map = {shuttle_by_cluster.index[0]: 1,
            shuttle_by_cluster.index[2]: 2,
            shuttle_by_cluster.index[1]: 3}
df_2022['cluster'] = df_2022['cluster_raw'].map(rank_map)

print(f"\n【クラスター別の基本統計（論文 Table 4 に対応）】")
summary = df_2022.groupby('cluster').agg(
    都道府県数=('pref', 'count'),
    シャトルラン平均=('shuttle', 'mean'),
    肥満率平均=('bmi_rate', 'mean'),
    運動部加入率平均=('sports_club', 'mean'),
    朝食摂取率平均=('breakfast', 'mean'),
).round(2)
print(summary)
print(f"\n  論文の参考値: Cluster1=81.9, Cluster2=77.8, Cluster3=78.3")

print(f"\n【クラスター1（標準型）の都道府県】")
print("  " + ', '.join(df_2022[df_2022['cluster'] == 1]['pref'].tolist()))
print(f"【クラスター2（生活習慣課題型）の都道府県】")
print("  " + ', '.join(df_2022[df_2022['cluster'] == 2]['pref'].tolist()))
print(f"【クラスター3（食文化・環境要因型）の都道府県】")
print("  " + ', '.join(df_2022[df_2022['cluster'] == 3]['pref'].tolist()))

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

print("図3: エルボー法とシルエットスコアグラフを作成中...")

fig3, axes3 = plt.subplots(1, 2, figsize=(12, 5))
fig3.suptitle('Step3. K-meansクラスタリング: 最適クラスター数の選択（k=3）\n'
              '変数: 県肥満率・運動部加入率・朝食摂取率（2022年度、47都道府県）',
              fontsize=11, fontweight='bold')

# (a) エルボー法
ks = range(2, 7)
ax3a = axes3[0]
ax3a.plot(list(ks), inertias, 'o-', color='#1565C0', linewidth=2.2, markersize=8,
          markerfacecolor='white', markeredgewidth=2.5)
ax3a.axvline(3, color='red', linestyle='--', linewidth=1.8, alpha=0.8, label='k=3 選択')
for k, ine in zip(ks, inertias):
    ax3a.annotate(f'{ine:.1f}', (k, ine), textcoords='offset points',
                  xytext=(0, 8), fontsize=8, ha='center')
ax3a.set_xlabel('クラスター数 k', fontsize=11)
ax3a.set_ylabel('Inertia（Within-cluster 平方和）', fontsize=10)
ax3a.set_title('エルボー法', fontsize=11, fontweight='bold')
ax3a.legend(fontsize=10)
ax3a.grid(True, alpha=0.3)
ax3a.set_xticks(list(ks))

# (b) シルエットスコア
ax3b = axes3[1]
ax3b.plot(list(ks), silhouettes, 's-', color='#E65100', linewidth=2.2, markersize=8,
          markerfacecolor='white', markeredgewidth=2.5)
ax3b.axvline(3, color='red', linestyle='--', linewidth=1.8, alpha=0.8, label='k=3 選択')
best_k = list(ks)[silhouettes.index(max(silhouettes))]
ax3b.axvline(best_k, color='green', linestyle=':', linewidth=1.5, alpha=0.7,
             label=f'Silhouette最大 k={best_k}')
for k, sil in zip(ks, silhouettes):
    ax3b.annotate(f'{sil:.3f}', (k, sil), textcoords='offset points',
                  xytext=(0, 8), fontsize=8, ha='center')
ax3b.set_xlabel('クラスター数 k', fontsize=11)
ax3b.set_ylabel('Silhouette Score', fontsize=10)
ax3b.set_title('シルエットスコア\n（高いほどクラスターの分離が明確）', fontsize=10, fontweight='bold')
ax3b.legend(fontsize=9)
ax3b.grid(True, alpha=0.3)
ax3b.set_xticks(list(ks))

plt.tight_layout()
save_fig('fig3_elbow')

print("図4: クラスター散布図とボックスプロットを作成中...")

fig4, axes4 = plt.subplots(1, 3, figsize=(15, 5))
fig4.suptitle('Step3. K-meansクラスタリング（k=3）の結果 ― 47都道府県の3類型\n'
              '出典: スポーツ庁体力調査・愛媛県オープンデータ（2022年度実データ）',
              fontsize=11, fontweight='bold')

# (a) 散布図: 肥満率 vs 運動部加入率
ax4a = axes4[0]
for cl in [1, 2, 3]:
    mask = df_2022['cluster'] == cl
    ax4a.scatter(df_2022.loc[mask, 'bmi_rate'], df_2022.loc[mask, 'sports_club'],
                 color=C_CLUSTER[cl], s=70, alpha=0.85, edgecolors='white', linewidth=0.5,
                 label=C_LABELS[cl], zorder=3)
ax4a.set_xlabel('県肥満率（%）', fontsize=10)
ax4a.set_ylabel('運動部加入率（%）', fontsize=10)
ax4a.set_title('肥満率 vs 運動部加入率', fontsize=10, fontweight='bold')
ax4a.legend(fontsize=8, loc='upper right')
ax4a.grid(True, alpha=0.3)

# (b) 散布図: 肥満率 vs 朝食摂取率
ax4b = axes4[1]
for cl in [1, 2, 3]:
    mask = df_2022['cluster'] == cl
    ax4b.scatter(df_2022.loc[mask, 'bmi_rate'], df_2022.loc[mask, 'breakfast'],
                 color=C_CLUSTER[cl], s=70, alpha=0.85, edgecolors='white', linewidth=0.5,
                 label=C_LABELS[cl], zorder=3)
ax4b.set_xlabel('県肥満率（%）', fontsize=10)
ax4b.set_ylabel('朝食を毎日食べる割合（%）', fontsize=10)
ax4b.set_title('肥満率 vs 朝食摂取率', fontsize=10, fontweight='bold')
ax4b.legend(fontsize=8, loc='upper right')
ax4b.grid(True, alpha=0.3)

# (c) クラスター別シャトルラン平均（棒グラフ）
ax4c = axes4[2]
cl_order = [1, 2, 3]
sl_means = [df_2022[df_2022['cluster'] == cl]['shuttle'].mean() for cl in cl_order]
sl_stds  = [df_2022[df_2022['cluster'] == cl]['shuttle'].std()  for cl in cl_order]
cl_names = ['Cluster1\n標準型', 'Cluster2\n生活習慣\n課題型', 'Cluster3\n食文化・\n環境要因型']
ax4c.bar(np.arange(3), sl_means, yerr=sl_stds, capsize=5,
         color=[C_CLUSTER[c] for c in cl_order], alpha=0.85,
         edgecolor='white', error_kw={'elinewidth': 1.5})
ax4c.set_xticks(np.arange(3))
ax4c.set_xticklabels(cl_names, fontsize=9)
ax4c.set_ylabel('シャトルラン平均回数（回）', fontsize=10)
ax4c.set_title('クラスター別\nシャトルラン平均', fontsize=10, fontweight='bold')
ax4c.grid(axis='y', alpha=0.3)
y_min_c = min(sl_means) - max(sl_stds) * 1.2
y_max_c = max(sl_means) + max(sl_stds) * 1.8
ax4c.set_ylim(y_min_c, y_max_c)
for i, (m, s) in enumerate(zip(sl_means, sl_stds)):
    ax4c.text(i, m + s + 0.3, f'{m:.1f}', ha='center', fontsize=10, fontweight='bold')
ax4c.text(0.5, 0.03,
          '論文参考値: C1=81.9, C2=77.8, C3=78.3',
          transform=ax4c.transAxes, fontsize=8, ha='center',
          bbox=dict(boxstyle='round', facecolor='#FFF9C4', alpha=0.9, edgecolor='#F9A825'))

plt.tight_layout()
save_fig('fig4_cluster')


print("\n" + "=" * 65)
print("✓ 全図の生成完了（4枚）")
print("  fig1_trend.png   : シャトルラン時系列推移（公表値＋都道府県別実データ）")
print("  fig2_fe.png      : パネルモデル選択（F検定・Hausman）と固定効果係数")
print("  fig3_elbow.png   : エルボー法とシルエットスコア")
print("  fig4_cluster.png : クラスター散布図と類型別シャトルラン平均")
print("=" * 65)

print("\n【最終サマリー】")
print(f"\n  ■ モデル選択")
print(f"    時間固定効果 F={f_time:.2f} ({'有意' if p_time < 0.05 else '非有意'})")
print(f"    個体固定効果 F={f_ind:.2f} ({'有意' if p_ind < 0.05 else '非有意'})")
print(f"    Two-way FE   F={f_twoway:.2f} ({'有意' if p_twoway < 0.05 else '非有意'})")
print(f"    Hausman χ²={H_stat:.2f} p={p_H:.4f} → {'個体FE採用' if p_H < 0.05 else 'RE採用'}")
print(f"\n  ■ 個体固定効果モデルで有意な変数")
for zv in Z_VARS:
    pi = pval_ind[zv]
    if pi < 0.1:
        sig = '***' if pi < 0.001 else '**' if pi < 0.01 else '*' if pi < 0.05 else '†'
        print(f"    {VAR_LABELS[zv]:<22} β={coef_ind[zv]:+.4f}  p={pi:.4f} {sig}")
print(f"\n  ■ K-meansクラスタリング（k=3）")
for cl in [1, 2, 3]:
    n_cl = (df_2022['cluster'] == cl).sum()
    m_cl = df_2022[df_2022['cluster'] == cl]['shuttle'].mean()
    print(f"    Cluster{cl}: n={n_cl}都道府県, シャトルラン平均={m_cl:.1f}回")