""" 教育用再現コード: 2025年 統計データ分析コンペティション 統計数理賞 ================================================================= 論文タイトル：クラスタリングを用いた男子中学生の持久力低下の要因特定 著者：新井陽登・陸家傑・宮下真翔（上智大学大学院理工学研究科） 岡田蒼未・柴田悠生（上智大学理工学部情報理工学科） 【分析概要】 データ：47都道府県パネルデータ（2019, 2021, 2022年度） 目的変数：男子中学生 20mシャトルラン平均回数 Step1. モデル選択（F検定 + Hausman検定） - 5モデル比較：Pooled OLS / 個体固定効果 / 時間固定効果 / Two-way FE / 変量効果 - F検定：時間FE有意(F≈10.3)、個体FE有意(F≈31.4)、two-way FE非有意 - Hausman検定：固定効果モデルが適切（p<0.05） Step2. 有意な説明変数の特定 - 共通有意：県肥満率(−***)、運動部加入率(+***)、朝食を毎日食べる割合(+***) - 時間FEのみ：最高気温(+*)、1学校あたり生徒数(−**) Step3. K-meansクラスタリング（k=3）で47都道府県を類型化 - Cluster1「標準型」: シャトルラン平均81.9（最高） - Cluster2「生活習慣課題型」: 77.8（最低）、高肥満率・低朝食・低運動部 - Cluster3「食文化・環境要因型」: 78.3（中）、東北・北関東に集中 【データサイエンス学習ポイント】 1. パネルデータの5モデル（Pooled OLS/個体FE/時間FE/Two-way FE/RE） 2. F検定によるモデル比較（固定効果の有意性検定） 3. Hausman検定（FEとREの選択基準） 4. K-meansクラスタリング（エルボー法・シルエットスコア） 5. 変数の標準化（Z-score）の意義 6. クラスタ解釈：重心プロファイルとレーダーチャート 【データ】 data/2025_U3/2025_U3_panel.csv（2025_U3_data_prep.py で生成） スポーツ庁・SSDSE-B・愛媛県オープンデータをもとに構築した実データ。 ※ 合成データ・乱数生成（np.random.seed等）は一切使用しない。 ================================================================= """ # ============================================================ # 【データの準備】実行前に以下のデータファイルを用意してください # # 必要ファイル: data/raw/ 以下に SSDSE CSVを配置 # ダウンロード先: https://www.nstac.go.jp/use/literacy/ssdse/ # # フォルダ配置（プロジェクトルートからの相対パス）: # code/ ← このスクリプトの場所 # data/raw/ ← CSV ファイルをここに配置 # html/figures/ ← 図の出力先（自動生成） # # 実行方法（ファイルを一切編集せず実行可能）: # python3 code/2025_U3_suri.py # ============================================================ 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) # ================================================================ # ■ 実データ読み込み # data/2025_U3/2025_U3_panel.csv（2025_U3_data_prep.py で生成） # ================================================================ 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': '朝食を毎日食べる割合', } # ================================================================ # ■ Step1. パネルモデル推定と F検定・Hausman検定 # ================================================================ 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}") # ================================================================ # ■ Step2. 有意な説明変数の回帰結果 # ================================================================ 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}") # ================================================================ # ■ Step3. K-meansクラスタリング（k=3） # ================================================================ 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())) # ================================================================ # ■ 図の生成（4枚） # fig1_trend : シャトルラン時系列推移 # fig2_fe : パネルモデル選択（F検定 + Hausman + 個体FE係数） # fig3_elbow : エルボー法 + シルエットスコア # fig4_cluster : クラスター散布図 + ボックスプロット # ================================================================ print("\n\n" + "=" * 65) print("■ 図の生成（4枚）") print("=" * 65) # ──────────────────────────────────────────────────────────────── # 図1：シャトルラン結果の時系列推移（スポーツ庁公表値） # ──────────────────────────────────────────────────────────────── 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') # ──────────────────────────────────────────────────────────────── # 図2：固定効果モデルの選択と係数（F検定 + Hausman + 個体FE係数） # ──────────────────────────────────────────────────────────────── 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') # ──────────────────────────────────────────────────────────────── # 図3：エルボー法 + シルエットスコア # ──────────────────────────────────────────────────────────────── 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') # ──────────────────────────────────────────────────────────────── # 図4：クラスター散布図 + ボックスプロット（クラスター解釈） # ──────────────────────────────────────────────────────────────── 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}回")