""" 教育用再現コード: 2022年 統計データ分析コンペティション 統計数理賞（大学生・一般の部） ================================================================================= 論文タイトル：リサイクル活動に対する地域・政策要因の研究 -主成分分析・階層的クラスタリングを用いた市町村別分析- 【分析概要】 データ：47都道府県 2022年度データ（SSDSE-B-2026） 目的変数：ごみのリサイクル率 Step1. 相関分析（相関ヒートマップ） - リサイクル率とごみ・消費・人口・気温関連変数の相関を可視化 Step2. 主成分分析（PCA） - StandardScaler で標準化後に PCA を実施 - PC1・PC2 バイプロット（47都道府県ラベル付き） - 累積寄与率・因子負荷量の解釈 Step3. Ward法 階層的クラスタリング - 標準化データを用いてWard法デンドログラムを作成（47都道府県） - デンドログラムのカット高さからクラスター数を決定（k=4） Step4. クラスター別プロファイル比較 - クラスター別の主要変数平均値を棒グラフで比較 【データサイエンス学習ポイント】 1. 標準化（StandardScaler）の必要性：変数間のスケール差を除去 2. PCA の解釈：固有値・寄与率・因子負荷量（loadings） 3. バイプロット：サンプルと変数の関係を同一図で表示 4. Ward法の特徴：クラスター内分散を最小化する凝集型クラスタリング 5. デンドログラムの読み方：カット高さとクラスター数の選択 【データ】 data/raw/SSDSE-B-2026.csv（47都道府県レベル 都道府県別統計） 出典：統計数理研究所 SSDSE（社会・人口統計体系） ※ 合成データ・乱数生成（np.random.seed等）は一切使用しない。 ================================================================================= """ # ============================================================ # 【データの準備】実行前に以下のデータファイルを用意してください # # 必要ファイル: # ・SSDSE-B-2026.csv # → data/raw/SSDSE-B-2026.csv に配置 # # ダウンロード先: # https://www.nstac.go.jp/use/literacy/ssdse/ # （SSDSE-B（社会・人口統計体系 都道府県データ） の CSV をダウンロード） # # フォルダ配置（プロジェクトルートからの相対パス）: # code/ ← このスクリプトの場所 # data/raw/ ← CSV ファイルをここに配置 # html/figures/ ← 図の出力先（自動生成） # # 実行方法（ファイルを一切編集せず実行可能）: # python3 code/2022_U3_suri.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.preprocessing import StandardScaler from sklearn.decomposition import PCA from scipy.cluster.hierarchy import dendrogram, linkage, fcluster 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'2022_U3_{name}.png') plt.savefig(path, bbox_inches='tight', dpi=150) plt.close() print(f" → {path} 保存完了") # クラスターカラー（k=4） C_CLUSTER = {1: '#E53935', 2: '#1E88E5', 3: '#43A047', 4: '#FB8C00'} # ================================================================ # ■ 実データ読み込み（SSDSE-B-2026） # ================================================================ DATA_B = 'data/raw/SSDSE-B-2026.csv' 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)].copy() df_2022 = df_b[df_b['年度'] == 2022].copy() print("=" * 70) print("■ SSDSE-B-2026 実データ（2022年度、N=47都道府県）") print(" 出典: 統計数理研究所 SSDSE（社会・人口統計体系）") print("=" * 70) print(f" 行数: {len(df_2022)}, 列数: {df_2022.shape[1]}") # 分析変数の定義 TARGET = 'ごみのリサイクル率' FEATURES = { 'ごみ総排出量（総量）': 'ごみ総排出量', '1人1日当たりの排出量': '1人1日排出量', '消費支出（二人以上の世帯）': '消費支出', '食料費（二人以上の世帯）': '食料費', '光熱・水道費（二人以上の世帯）': '光熱・水道費', '教養娯楽費（二人以上の世帯）': '教養娯楽費', '年平均気温': '年平均気温', } # 高齢化率を計算（65歳以上人口 / 総人口） df_2022 = df_2022.copy() df_2022['高齢化率'] = df_2022['65歳以上人口'] / df_2022['総人口'] * 100 FEATURES['高齢化率'] = '高齢化率' # 分析用データフレームを作成 feat_cols = list(FEATURES.keys()) all_cols = [TARGET] + feat_cols df_ana = df_2022[['都道府県'] + all_cols].copy().reset_index(drop=True) df_ana.columns = ['都道府県', TARGET] + [FEATURES[c] for c in feat_cols] FEAT_SHORT = list(FEATURES.values()) # 短縮名リスト print(f"\n【分析変数一覧（N=47都道府県、2022年度）】") print(f" 目的変数: {TARGET}") print(f" 説明変数: {', '.join(FEAT_SHORT)}") print(f"\n【基本統計量】") print(df_ana[FEAT_SHORT + [TARGET]].describe().round(2)) # 標準化 scaler = StandardScaler() X_raw = df_ana[FEAT_SHORT].values X_scaled = scaler.fit_transform(X_raw) prefs = df_ana['都道府県'].tolist() # ================================================================ # ■ Step1. 相関分析 # ================================================================ print("\n" + "=" * 70) print("■ Step1. 相関分析（リサイクル率 vs 説明変数）") print("=" * 70) corr_df = df_ana[[TARGET] + FEAT_SHORT].corr() print(f"\n【リサイクル率との相関係数】") recycling_corr = corr_df[TARGET].drop(TARGET).sort_values(ascending=False) for var, r in recycling_corr.items(): direction = '+' if r > 0 else '-' strength = '強' if abs(r) >= 0.5 else '中' if abs(r) >= 0.3 else '弱' print(f" {var:<14}: r = {r:+.4f} ({direction} {strength})") # ================================================================ # ■ Step2. 主成分分析（PCA） # ================================================================ print("\n" + "=" * 70) print("■ Step2. 主成分分析（PCA）") print("=" * 70) pca = PCA() pca.fit(X_scaled) X_pca = pca.transform(X_scaled) explained_var = pca.explained_variance_ratio_ cumulative_var = np.cumsum(explained_var) loadings = pca.components_ # shape: (n_components, n_features) print(f"\n【主成分の寄与率】") print(f" {'PC':<6} {'固有値':>8} {'寄与率':>9} {'累積寄与率':>11}") print(" " + "-" * 38) eigenvalues = pca.explained_variance_ for i in range(min(6, len(explained_var))): print(f" PC{i+1:<4} {eigenvalues[i]:>8.4f} {explained_var[i]*100:>8.2f}% {cumulative_var[i]*100:>10.2f}%") print(f"\n【PC1・PC2 因子負荷量（loadings）】") print(f" {'変数':<14} {'PC1 負荷量':>11} {'PC2 負荷量':>11}") print(" " + "-" * 38) for j, feat in enumerate(FEAT_SHORT): print(f" {feat:<14} {loadings[0, j]:>+11.4f} {loadings[1, j]:>+11.4f}") print(f"\n PC1 寄与率: {explained_var[0]*100:.1f}% PC2 寄与率: {explained_var[1]*100:.1f}%") print(f" PC1+PC2 累積: {cumulative_var[1]*100:.1f}%") # ================================================================ # ■ Step3. Ward法 階層的クラスタリング # ================================================================ print("\n" + "=" * 70) print("■ Step3. Ward法 階層的クラスタリング（47都道府県）") print("=" * 70) Z = linkage(X_scaled, method='ward') # k=4 でクラスター割り当て N_CLUSTERS = 4 labels_ward = fcluster(Z, t=N_CLUSTERS, criterion='maxclust') df_ana['cluster'] = labels_ward print(f"\n【クラスター別の都道府県（k={N_CLUSTERS}）】") for cl in sorted(df_ana['cluster'].unique()): prefs_cl = df_ana[df_ana['cluster'] == cl]['都道府県'].tolist() print(f" クラスター{cl} (n={len(prefs_cl)}): {', '.join(prefs_cl)}") print(f"\n【クラスター別 リサイクル率の平均】") cl_recycling = df_ana.groupby('cluster')[TARGET].agg(['mean', 'std', 'min', 'max']) for cl, row in cl_recycling.iterrows(): print(f" クラスター{cl}: 平均={row['mean']:.2f}% SD={row['std']:.2f}% " f"[{row['min']:.1f}%~{row['max']:.1f}%]") # ================================================================ # ■ Step4. クラスター別プロファイル # ================================================================ print("\n" + "=" * 70) print("■ Step4. クラスター別プロファイル比較（標準化平均）") print("=" * 70) df_ana_z = df_ana.copy() for j, feat in enumerate(FEAT_SHORT): df_ana_z[feat + '_z'] = X_scaled[:, j] profile = df_ana_z.groupby('cluster')[[f + '_z' for f in FEAT_SHORT]].mean() print(f"\n【クラスター別 標準化変数平均値（高いほど全国平均より上）】") print(f" {'変数':<14}", end='') for cl in sorted(df_ana['cluster'].unique()): print(f" {'CL' + str(cl):>8}", end='') print() print(" " + "-" * (14 + 9 * N_CLUSTERS)) for feat in FEAT_SHORT: print(f" {feat:<14}", end='') for cl in sorted(df_ana['cluster'].unique()): v = profile.loc[cl, feat + '_z'] print(f" {v:>+8.3f}", end='') print() # ================================================================ # ■ 図の生成（4枚） # fig1_corr : 相関ヒートマップ（リサイクル率と説明変数） # fig2_pca : PCA バイプロット（PC1-PC2・都道府県ラベル） # fig3_dendrogram: Ward法デンドログラム（47都道府県） # fig4_cluster : クラスター別プロファイル比較（棒グラフ） # ================================================================ print("\n\n" + "=" * 70) print("■ 図の生成（4枚）") print("=" * 70) # ──────────────────────────────────────────────────────────────── # 図1: 相関ヒートマップ（リサイクル率と説明変数） # ──────────────────────────────────────────────────────────────── print("図1: 相関ヒートマップを作成中...") fig1, ax1 = plt.subplots(figsize=(9, 7)) # 相関行列（リサイクル率を含む全変数） all_vars_short = [TARGET] + FEAT_SHORT corr_mat = df_ana[[TARGET] + FEAT_SHORT].corr() im = ax1.imshow(corr_mat.values, cmap='RdBu_r', vmin=-1, vmax=1, aspect='auto') plt.colorbar(im, ax=ax1, fraction=0.046, pad=0.04, label='相関係数') tick_labels = ['リサイクル率', 'ごみ総排出量', '1人1日排出量', '消費支出', '食料費', '光熱・水道費', '教養娯楽費', '年平均気温', '高齢化率'] ax1.set_xticks(range(len(tick_labels))) ax1.set_yticks(range(len(tick_labels))) ax1.set_xticklabels(tick_labels, rotation=40, ha='right', fontsize=9) ax1.set_yticklabels(tick_labels, fontsize=9) # 数値を各セルに表示 for i in range(len(tick_labels)): for j in range(len(tick_labels)): val = corr_mat.values[i, j] color = 'white' if abs(val) > 0.6 else 'black' ax1.text(j, i, f'{val:.2f}', ha='center', va='center', fontsize=7.5, color=color, fontweight='bold') # リサイクル率の行・列をハイライト for pos in range(len(tick_labels)): ax1.add_patch(plt.Rectangle((pos - 0.5, -0.5), 1, 1, fill=False, edgecolor='#FFA000', linewidth=0, alpha=0)) ax1.axhline(y=0.5, color='#FFA000', linewidth=2) ax1.axhline(y=-0.5, color='#FFA000', linewidth=2) ax1.axvline(x=0.5, color='#FFA000', linewidth=2) ax1.axvline(x=-0.5, color='#FFA000', linewidth=2) ax1.set_title('相関係数行列ヒートマップ\n（2022年度 47都道府県・SSDSE-B）', fontsize=12, fontweight='bold', pad=14) fig1.text(0.5, 0.01, '出典: 統計数理研究所 SSDSE-B-2026（実データ）', ha='center', fontsize=8, color='#666') plt.tight_layout(rect=[0, 0.03, 1, 1]) save_fig('fig1_corr') # ──────────────────────────────────────────────────────────────── # 図2: PCA バイプロット（PC1-PC2・都道府県ラベル） # ──────────────────────────────────────────────────────────────── print("図2: PCA バイプロットを作成中...") fig2, axes2 = plt.subplots(1, 2, figsize=(16, 6.5)) fig2.suptitle('Step2. 主成分分析（PCA）バイプロット\n' '2022年度 47都道府県 ごみ・消費・人口・気温変数（SSDSE-B）', fontsize=12, fontweight='bold') # --- (a) バイプロット: スコア + ローディング --- ax2a = axes2[0] # クラスター色でサンプルを散布 for cl in sorted(df_ana['cluster'].unique()): mask = df_ana['cluster'] == cl ax2a.scatter(X_pca[mask, 0], X_pca[mask, 1], color=C_CLUSTER[cl], s=60, alpha=0.85, edgecolors='white', linewidth=0.5, zorder=3, label=f'クラスター{cl}') # 都道府県名ラベル for i, pref in enumerate(prefs): ax2a.annotate(pref, (X_pca[i, 0], X_pca[i, 1]), fontsize=5.5, ha='center', va='bottom', xytext=(0, 3), textcoords='offset points', color='#333333', zorder=4) # ローディングベクトル（矢印） scale = 2.5 for j, feat in enumerate(FEAT_SHORT): ax2a.annotate('', xy=(loadings[0, j] * scale, loadings[1, j] * scale), xytext=(0, 0), arrowprops=dict(arrowstyle='->', color='#C62828', lw=1.5, shrinkA=0, shrinkB=0)) ax2a.text(loadings[0, j] * scale * 1.12, loadings[1, j] * scale * 1.12, feat, fontsize=8, color='#C62828', fontweight='bold', ha='center', va='center') ax2a.axhline(0, color='gray', linewidth=0.6, linestyle='--', alpha=0.5) ax2a.axvline(0, color='gray', linewidth=0.6, linestyle='--', alpha=0.5) ax2a.set_xlabel(f'PC1（寄与率: {explained_var[0]*100:.1f}%）', fontsize=10) ax2a.set_ylabel(f'PC2（寄与率: {explained_var[1]*100:.1f}%）', fontsize=10) ax2a.set_title(f'バイプロット（PC1 vs PC2）\n累積寄与率: {cumulative_var[1]*100:.1f}%', fontsize=10, fontweight='bold') ax2a.legend(fontsize=8, loc='upper right') ax2a.grid(True, alpha=0.25) # --- (b) 累積寄与率の棒グラフ --- ax2b = axes2[1] n_pcs = min(8, len(explained_var)) pc_labels = [f'PC{i+1}' for i in range(n_pcs)] bars = ax2b.bar(pc_labels, explained_var[:n_pcs] * 100, color='#1565C0', alpha=0.8, edgecolor='white', label='寄与率') ax2b_twin = ax2b.twinx() ax2b_twin.plot(pc_labels, cumulative_var[:n_pcs] * 100, 'o-', color='#E65100', linewidth=2.2, markersize=7, markerfacecolor='white', markeredgewidth=2, label='累積寄与率') ax2b_twin.axhline(80, color='#E65100', linestyle='--', linewidth=1, alpha=0.5) ax2b_twin.set_ylabel('累積寄与率（%）', fontsize=10, color='#E65100') ax2b_twin.set_ylim(0, 105) ax2b_twin.tick_params(axis='y', labelcolor='#E65100') for i, (bar, ev) in enumerate(zip(bars, explained_var[:n_pcs])): ax2b.text(bar.get_x() + bar.get_width() / 2, bar.get_height() + 0.3, f'{ev*100:.1f}%', ha='center', va='bottom', fontsize=8, color='#1565C0') for i, cv in enumerate(cumulative_var[:n_pcs]): ax2b_twin.annotate(f'{cv*100:.1f}%', (i, cv * 100), textcoords='offset points', xytext=(6, 4), fontsize=7.5, color='#E65100') ax2b.set_xlabel('主成分', fontsize=10) ax2b.set_ylabel('寄与率（%）', fontsize=10, color='#1565C0') ax2b.set_title('主成分の寄与率・累積寄与率', fontsize=10, fontweight='bold') ax2b.grid(axis='y', alpha=0.3) ax2b.tick_params(axis='y', labelcolor='#1565C0') # 因子負荷量の凡例表として追加 load_text = '【因子負荷量 PC1 / PC2】\n' for j, feat in enumerate(FEAT_SHORT): load_text += f'{feat}: {loadings[0,j]:+.2f} / {loadings[1,j]:+.2f}\n' ax2b.text(0.01, 0.01, load_text.strip(), transform=ax2b.transAxes, fontsize=7, va='bottom', ha='left', fontfamily='monospace', bbox=dict(boxstyle='round', facecolor='#EFF3FF', alpha=0.9, edgecolor='#1565C0')) plt.tight_layout() save_fig('fig2_pca') # ──────────────────────────────────────────────────────────────── # 図3: Ward法デンドログラム（47都道府県） # ──────────────────────────────────────────────────────────────── print("図3: Ward法デンドログラムを作成中...") fig3, ax3 = plt.subplots(figsize=(16, 7)) # クラスター色対応のデンドログラム cl_colors_ward = {} for i, pref in enumerate(prefs): cl = df_ana.loc[i, 'cluster'] cl_colors_ward[pref] = C_CLUSTER[cl] # デンドログラムを描画 ddata = dendrogram( Z, labels=prefs, ax=ax3, leaf_rotation=90, leaf_font_size=8, color_threshold=Z[-N_CLUSTERS + 1, 2], # k=4 のカット高さ above_threshold_color='#9E9E9E', ) # カット線 cut_height = (Z[-N_CLUSTERS, 2] + Z[-N_CLUSTERS + 1, 2]) / 2 ax3.axhline(y=cut_height, color='red', linestyle='--', linewidth=2, label=f'カット高さ（k={N_CLUSTERS}クラスター）') ax3.set_title(f'Ward法 階層的クラスタリング デンドログラム\n' f'47都道府県 ごみ・消費・人口・気温変数（2022年度・SSDSE-B）', fontsize=12, fontweight='bold') ax3.set_xlabel('都道府県', fontsize=10) ax3.set_ylabel('Ward距離（クラスター内分散の増分）', fontsize=10) ax3.legend(fontsize=9, loc='upper right') ax3.grid(axis='y', alpha=0.25) # クラスター凡例 patch_list = [mpatches.Patch(color=C_CLUSTER[cl], label=f'クラスター{cl}') for cl in sorted(df_ana['cluster'].unique())] ax3.legend(handles=patch_list + [ mpatches.Patch(color='red', label=f'カット高さ（k={N_CLUSTERS}）') ], fontsize=9, loc='upper left') fig3.text(0.5, 0.01, '出典: 統計数理研究所 SSDSE-B-2026（実データ）| Ward法: scipy.cluster.hierarchy', ha='center', fontsize=8, color='#666') plt.tight_layout(rect=[0, 0.03, 1, 1]) save_fig('fig3_dendrogram') # ──────────────────────────────────────────────────────────────── # 図4: クラスター別プロファイル比較（棒グラフ） # ──────────────────────────────────────────────────────────────── print("図4: クラスター別プロファイル棒グラフを作成中...") fig4, axes4 = plt.subplots(2, 4, figsize=(16, 9)) fig4.suptitle(f'Step4. Ward法クラスタリング（k={N_CLUSTERS}）クラスター別プロファイル比較\n' '2022年度 47都道府県（SSDSE-B-2026 実データ）', fontsize=12, fontweight='bold') axes4_flat = axes4.flatten() # リサイクル率をプロファイル変数に追加 profile_vars = [TARGET] + FEAT_SHORT n_vars_p = len(profile_vars) cl_order = sorted(df_ana['cluster'].unique()) x_pos = np.arange(len(cl_order)) for ax_idx, var in enumerate(profile_vars): if ax_idx >= len(axes4_flat): break ax = axes4_flat[ax_idx] means = [df_ana[df_ana['cluster'] == cl][var].mean() for cl in cl_order] stds = [df_ana[df_ana['cluster'] == cl][var].std() for cl in cl_order] bars = ax.bar(x_pos, means, yerr=stds, capsize=4, color=[C_CLUSTER[cl] for cl in cl_order], alpha=0.85, edgecolor='white', error_kw={'elinewidth': 1.2, 'capthick': 1.2}) ax.set_xticks(x_pos) ax.set_xticklabels([f'CL{cl}' for cl in cl_order], fontsize=9) ax.set_title(var, fontsize=9, fontweight='bold') ax.grid(axis='y', alpha=0.3) # 全国平均との差を示す水平線 overall_mean = df_ana[var].mean() ax.axhline(overall_mean, color='black', linestyle='--', linewidth=1.2, alpha=0.7, label=f'全国平均 ({overall_mean:.1f})') ax.legend(fontsize=6.5, loc='upper right') # 値ラベル for bar, m, s in zip(bars, means, stds): ax.text(bar.get_x() + bar.get_width() / 2, m + s + (max(means) - min(means)) * 0.03, f'{m:.1f}', ha='center', va='bottom', fontsize=7, fontweight='bold') # 余った軸を非表示 for ax_idx in range(n_vars_p, len(axes4_flat)): axes4_flat[ax_idx].set_visible(False) # クラスター凡例（図全体） patch_list4 = [mpatches.Patch(color=C_CLUSTER[cl], label=f'クラスター{cl}') for cl in cl_order] fig4.legend(handles=patch_list4, fontsize=10, loc='lower right', bbox_to_anchor=(0.98, 0.02), title='クラスター', title_fontsize=10) fig4.text(0.5, 0.005, '出典: 統計数理研究所 SSDSE-B-2026 | Ward法（scipy）| エラーバー: ±1 SD', ha='center', fontsize=8, color='#666') plt.tight_layout(rect=[0, 0.03, 1, 0.96]) save_fig('fig4_cluster') # ================================================================ # ■ 最終サマリー # ================================================================ print("\n" + "=" * 70) print("✓ 全図の生成完了（4枚）") print(" fig1_corr.png : 相関ヒートマップ（リサイクル率と説明変数）") print(" fig2_pca.png : PCA バイプロット（PC1-PC2・都道府県ラベル）") print(" fig3_dendrogram.png : Ward法デンドログラム（47都道府県）") print(" fig4_cluster.png : クラスター別プロファイル比較（棒グラフ）") print("=" * 70) print("\n【最終サマリー】") print(f"\n ■ リサイクル率との相関（上位）") for var, r in recycling_corr.items(): if abs(r) >= 0.3: print(f" {var:<14}: r = {r:+.4f}") print(f"\n ■ PCA 結果") print(f" PC1 寄与率: {explained_var[0]*100:.1f}%（主にごみ量・消費支出系）") print(f" PC2 寄与率: {explained_var[1]*100:.1f}%") print(f" PC1+PC2 累積: {cumulative_var[1]*100:.1f}%") print(f"\n ■ Ward法クラスタリング（k={N_CLUSTERS}）") for cl in sorted(df_ana['cluster'].unique()): n_cl = (df_ana['cluster'] == cl).sum() m_rec = df_ana[df_ana['cluster'] == cl][TARGET].mean() prefs_cl = df_ana[df_ana['cluster'] == cl]['都道府県'].tolist() print(f" クラスター{cl}: n={n_cl}, リサイクル率平均={m_rec:.1f}%") print(f" → {', '.join(prefs_cl)}")