""" 教育用再現コード: 2022年度 統計データ分析コンペティション 審査員奨励賞（高校生の部） ================================================================= 論文タイトル：SDGs指標の地域比較：都道府県別の持続可能性スコア分析 【分析概要】 データ：SSDSE-B-2026.csv（都道府県データ、2022年断面） 目的：SDGs関連指標を用いた都道府県の持続可能性スコア分析 SDGs関連7指標の構築： 経済指標： 消費支出_log = log(消費支出（二人以上の世帯）) → 生活水準の代理 求人倍率 = 月間有効求人数 / 月間有効求職者数 → 雇用充実度 環境指標： ごみリサイクル率 = ごみのリサイクル率（%） ごみ排出量_inv = 1 / (1人1日当たりの排出量) ← 小さいほど良い → 逆数で高スコア化 社会指標： 大学進学率 = 高校卒業者のうち進学者数 / 高校卒業者数 × 100 保育所密度 = 保育所等数 / 総人口 × 10000 （人口1万対） 人口持続： 合計特殊出生率 （出生率が高いほど持続可能性が高い） 分析の流れ： Step1. 主成分分析（PCA）→ 都道府県の総合持続可能性スコア Step2. PCAスコアプロット（PC1 vs PC2, 地域色分け） Step3. 変数の主成分負荷量ヒートマップ（どの変数がPCを構成するか） Step4. Ward法クラスタリング デンドログラム（類似都道府県のグループ化） Step5. 地域別レーダーチャート（6地域ブロック × 4指標の特徴比較） 【変数定義（SSDSE-B-2026列名）】 SSDSE-B-2026 : 年度 Prefecture : 都道府県 A1101 : 総人口 A4103 : 合計特殊出生率 E4601 : 高校卒業者数 E4602 : 高校卒業者のうち進学者数 F3102 : 月間有効求職者数（一般） F3103 : 月間有効求人数（一般） H5610 : 1人1日当たりの排出量（g） H5614 : ごみのリサイクル率（%） J2503 : 保育所等数 L3221 : 消費支出（二人以上の世帯）（円/月） 【派生指標の計算式】 消費支出_log = np.log(L3221) 求人倍率 = F3103 / F3102 ごみリサイクル率 = H5614 ごみ排出量_inv = 1 / H5610 ← 排出量が少ないほど逆数が大きい 大学進学率 = E4602 / E4601 × 100 保育所密度 = J2503 / A1101 × 10000 合計特殊出生率 = A4103 【データ出典】 SSDSE-B-2026.csv: 社会・人口統計体系（都道府県データ） 統計数理研究所・統計センターより公表の実データ 【データサイエンス学習ポイント】 1. 複合指標の構築（PCAによる合成スコアの作り方） ─ 相関のある複数変数を直交する主成分に変換する ─ 第1主成分スコアを「総合持続可能性指数」として解釈 2. 因子負荷量の解釈（どの変数がPC1・PC2を構成するか） ─ 負荷量の大きい変数がその主成分の「意味」を規定する ─ ヒートマップで変数×成分の関係を可視化 3. レーダーチャートの作り方（matplotlib.projections.polar の使い方） ─ ax = fig.add_subplot(projection='polar') で極座標系を取得 ─ 閉じたポリゴンを描くには末尾に先頭の値を追加する 4. SDGs達成度の統計的評価（ベンチマークと格差の可視化） ─ 地域ブロック別の平均スコアを比較してSDGs達成の「地域格差」を測定 ─ Ward法クラスタリングで「類似した持続可能性プロファイル」の都道府県を識別 ================================================================= """ # ============================================================ # 【データの準備】実行前に以下のデータファイルを用意してください # # 必要ファイル: # ・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_H5_15_shorei.py # ============================================================ import numpy as np import pandas as pd import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt from sklearn.preprocessing import StandardScaler from sklearn.decomposition import PCA from scipy.cluster.hierarchy import dendrogram, linkage from scipy.spatial.distance import pdist import warnings warnings.filterwarnings('ignore') # ────────────────────────────────────────────────────────────── # 共通設定 # ────────────────────────────────────────────────────────────── plt.rcParams['font.family'] = 'Hiragino Sans' plt.rcParams['axes.unicode_minus'] = False plt.rcParams['figure.dpi'] = 150 import os FIG_DIR = 'html/figures' DATA_B = 'data/raw/SSDSE-B-2026.csv' os.makedirs(FIG_DIR, exist_ok=True) # ================================================================ # ■ 実データ読み込み（SSDSE-B-2026: 都道府県データ） # ================================================================ df_raw = pd.read_csv(DATA_B, encoding='shift-jis', header=0) # 行0がラベル行（日本語名）、行1以降が実データ df_all = df_raw.iloc[1:].copy() df_all.columns = df_raw.columns # 必要な数値列に変換 NUM_COLS = ['A1101', 'A4103', 'E4601', 'E4602', 'F3102', 'F3103', 'H5610', 'H5614', 'J2503', 'L3221'] for c in NUM_COLS: df_all[c] = pd.to_numeric(df_all[c], errors='coerce') df_all['年度'] = df_all['SSDSE-B-2026'].astype(int) df_all['都道府県'] = df_all['Prefecture'].astype(str) print("=" * 65) print("■ SSDSE-B-2026 読み込み完了") print(f" 年度範囲: {df_all['年度'].min()}〜{df_all['年度'].max()}") print(f" 都道府県数（各年）: {df_all.groupby('年度')['都道府県'].count().iloc[0]}件") print("=" * 65) # ================================================================ # ■ 2022年データ抽出・SDGs指標構築 # ================================================================ df = df_all[df_all['年度'] == 2022].copy().reset_index(drop=True) # --- SDGs関連7指標の計算 --- df['消費支出_log'] = np.log(df['L3221']) # 経済: 生活水準（対数変換） df['求人倍率'] = df['F3103'] / df['F3102'] # 経済: 雇用充実度 df['ごみリサイクル率'] = df['H5614'] # 環境: リサイクル率(%) df['ごみ排出量_inv'] = 1.0 / df['H5610'] # 環境: 排出少＝良い → 逆数 df['大学進学率'] = df['E4602'] / df['E4601'] * 100 # 社会: 教育到達度(%) df['保育所密度'] = df['J2503'] / df['A1101'] * 10000 # 社会: 人口1万対保育所数 df['合計特殊出生率'] = df['A4103'] # 人口持続性 INDICATOR_COLS = [ '消費支出_log', '求人倍率', 'ごみリサイクル率', 'ごみ排出量_inv', '大学進学率', '保育所密度', '合計特殊出生率' ] INDICATOR_LABELS = [ '消費支出\n(log)', '求人倍率', 'ごみ\nリサイクル率', 'ごみ排出量\n(逆数)', '大学進学率', '保育所密度', '合計特殊\n出生率' ] # 欠損値がある行を除外 df_sdgs = df[['都道府県'] + INDICATOR_COLS].dropna().reset_index(drop=True) prefectures = df_sdgs['都道府県'].values X_raw = df_sdgs[INDICATOR_COLS].values print(f"\n■ 2022年 SDGs指標データ: {len(df_sdgs)}都道府県") print(f" 7指標の基本統計:") for col in INDICATOR_COLS: print(f" {col}: 平均={df_sdgs[col].mean():.3f}, SD={df_sdgs[col].std():.3f}") # ================================================================ # ■ 標準化 + 主成分分析（PCA） # ================================================================ scaler = StandardScaler() X_std = scaler.fit_transform(X_raw) pca = PCA(n_components=7) pca.fit(X_std) scores = pca.transform(X_std) # 主成分スコア (47 × 7) loadings = pca.components_.T # 因子負荷行列 (7変数 × 7成分) ev_ratio = pca.explained_variance_ratio_ print(f"\n■ PCA 寄与率:") for i, ev in enumerate(ev_ratio[:5]): cum = ev_ratio[:i+1].sum() print(f" PC{i+1}: {ev*100:.1f}% (累積 {cum*100:.1f}%)") # ================================================================ # ■ 地方区分マップ（6地域ブロック） # ================================================================ REGION_MAP = { '北海道': '北海道・東北', '青森県': '北海道・東北', '岩手県': '北海道・東北', '宮城県': '北海道・東北', '秋田県': '北海道・東北', '山形県': '北海道・東北', '福島県': '北海道・東北', '茨城県': '関東', '栃木県': '関東', '群馬県': '関東', '埼玉県': '関東', '千葉県': '関東', '東京都': '関東', '神奈川県': '関東', '新潟県': '中部', '富山県': '中部', '石川県': '中部', '福井県': '中部', '山梨県': '中部', '長野県': '中部', '岐阜県': '中部', '静岡県': '中部', '愛知県': '中部', '三重県': '近畿', '滋賀県': '近畿', '京都府': '近畿', '大阪府': '近畿', '兵庫県': '近畿', '奈良県': '近畿', '和歌山県': '近畿', '鳥取県': '中国・四国', '島根県': '中国・四国', '岡山県': '中国・四国', '広島県': '中国・四国', '山口県': '中国・四国', '徳島県': '中国・四国', '香川県': '中国・四国', '愛媛県': '中国・四国', '高知県': '中国・四国', '福岡県': '九州・沖縄', '佐賀県': '九州・沖縄', '長崎県': '九州・沖縄', '熊本県': '九州・沖縄', '大分県': '九州・沖縄', '宮崎県': '九州・沖縄', '鹿児島県': '九州・沖縄', '沖縄県': '九州・沖縄', } REGION_ORDER = ['北海道・東北', '関東', '中部', '近畿', '中国・四国', '九州・沖縄'] REGION_COLORS = { '北海道・東北': '#1565C0', '関東': '#E65100', '中部': '#2E7D32', '近畿': '#6A1B9A', '中国・四国': '#795548', '九州・沖縄': '#00838F', } regions = np.array([REGION_MAP.get(p, '不明') for p in prefectures]) df_sdgs['地方'] = regions # ================================================================ # ■ 図1: PCAスコアプロット（PC1 vs PC2, 地域色分け） # ================================================================ print("\n" + "=" * 65) print("■ 図1: PCAスコアプロット（PC1 vs PC2）") print("=" * 65) fig1, ax1 = plt.subplots(figsize=(12, 8)) for region in REGION_ORDER: mask = regions == region if mask.sum() == 0: continue ax1.scatter(scores[mask, 0], scores[mask, 1], color=REGION_COLORS[region], s=80, alpha=0.85, edgecolors='white', linewidth=0.7, label=f'{region} (n={mask.sum()})', zorder=3) # 都道府県ラベル（全都道府県に短縮名を表示） for i, pref in enumerate(prefectures): short = pref.replace('県', '').replace('都', '').replace('道', '').replace('府', '') ax1.annotate(short, (scores[i, 0], scores[i, 1]), xytext=(3, 3), textcoords='offset points', fontsize=6.5, color='#333', alpha=0.9) ax1.axhline(0, color='gray', linestyle='--', linewidth=0.8, alpha=0.5) ax1.axvline(0, color='gray', linestyle='--', linewidth=0.8, alpha=0.5) ax1.set_xlabel(f'PC1（寄与率 {ev_ratio[0]*100:.1f}%）── 総合持続可能性スコア', fontsize=12) ax1.set_ylabel(f'PC2（寄与率 {ev_ratio[1]*100:.1f}%）── 経済 vs 環境・人口', fontsize=12) ax1.set_title( 'SDGs指標の主成分スコアプロット（2022年, 都道府県別）\n' '〜PC1: 総合持続可能性スコア。右ほど持続可能性が高い傾向〜', fontsize=13, fontweight='bold' ) ax1.legend(loc='lower left', fontsize=9, framealpha=0.88, bbox_to_anchor=(0.0, 0.0)) ax1.grid(True, alpha=0.25) plt.tight_layout() fig1.savefig(os.path.join(FIG_DIR, '2022_H5_15_fig1_pca_score.png'), bbox_inches='tight', dpi=150) plt.close(fig1) print("図1保存: 2022_H5_15_fig1_pca_score.png") # PC1スコア上位・下位を表示 df_pca = df_sdgs.copy() df_pca['PC1'] = scores[:, 0] df_pca['PC2'] = scores[:, 1] top5 = df_pca.nlargest(5, 'PC1')[['都道府県', '地方', 'PC1', 'PC2']] bot5 = df_pca.nsmallest(5, 'PC1')[['都道府県', '地方', 'PC1', 'PC2']] print(f"\n PC1スコア TOP5（持続可能性が高い都道府県）:") print(top5.to_string(index=False)) print(f"\n PC1スコア BOTTOM5（持続可能性が低い都道府県）:") print(bot5.to_string(index=False)) # ================================================================ # ■ 図2: 変数の主成分負荷量ヒートマップ（因子負荷行列） # ================================================================ print("\n" + "=" * 65) print("■ 図2: 主成分負荷量ヒートマップ") print("=" * 65) N_PC_SHOW = 4 # PC1〜PC4を表示 load_show = loadings[:, :N_PC_SHOW] fig2, ax2 = plt.subplots(figsize=(9, 6)) im = ax2.imshow(load_show, cmap='RdBu_r', vmin=-1, vmax=1, aspect='auto') plt.colorbar(im, ax=ax2, label='因子負荷量', shrink=0.85) ax2.set_xticks(range(N_PC_SHOW)) ax2.set_xticklabels( [f'PC{i+1}\n({ev_ratio[i]*100:.1f}%)' for i in range(N_PC_SHOW)], fontsize=11 ) ax2.set_yticks(range(len(INDICATOR_LABELS))) ax2.set_yticklabels(INDICATOR_LABELS, fontsize=11) # セル内に数値を表示 for i in range(len(INDICATOR_COLS)): for j in range(N_PC_SHOW): val = load_show[i, j] color = 'white' if abs(val) > 0.55 else '#222' ax2.text(j, i, f'{val:.2f}', ha='center', va='center', fontsize=10, fontweight='bold', color=color) ax2.set_title( 'SDGs指標の主成分負荷量ヒートマップ（PC1〜PC4）\n' '〜赤：正の負荷（その指標がPCと同方向）、青：負の負荷〜', fontsize=13, fontweight='bold' ) ax2.set_xlabel('主成分（説明した分散の割合）', fontsize=11) ax2.set_ylabel('SDGs指標', fontsize=11) plt.tight_layout() fig2.savefig(os.path.join(FIG_DIR, '2022_H5_15_fig2_loadings.png'), bbox_inches='tight', dpi=150) plt.close(fig2) print("図2保存: 2022_H5_15_fig2_loadings.png") print(f"\n 因子負荷量（PC1〜PC4）:") df_load = pd.DataFrame(load_show, index=INDICATOR_COLS, columns=[f'PC{i+1}' for i in range(N_PC_SHOW)]) print(df_load.round(3).to_string()) # ================================================================ # ■ 図3: Ward法クラスタリング デンドログラム # ================================================================ print("\n" + "=" * 65) print("■ 図3: Ward法クラスタリング デンドログラム") print("=" * 65) # Ward法でクラスタリング（PCA後のスコア全7成分を使用） Z = linkage(X_std, method='ward', metric='euclidean') # 都道府県の短縮名リスト pref_short = [p.replace('県', '').replace('都', '').replace('道', '').replace('府', '') for p in prefectures] fig3, ax3 = plt.subplots(figsize=(16, 7)) dn = dendrogram(Z, labels=pref_short, ax=ax3, color_threshold=0.5 * Z[-5:, 2].mean(), leaf_font_size=8.5, leaf_rotation=75) ax3.set_title( 'Ward法クラスタリング デンドログラム（SDGs7指標, 47都道府県）\n' '〜類似した持続可能性プロファイルを持つ都道府県がグループ化される〜', fontsize=13, fontweight='bold' ) ax3.set_xlabel('都道府県', fontsize=11) ax3.set_ylabel('Ward距離（ユークリッド距離）', fontsize=11) ax3.yaxis.grid(True, alpha=0.3) plt.tight_layout() fig3.savefig(os.path.join(FIG_DIR, '2022_H5_15_fig3_dendrogram.png'), bbox_inches='tight', dpi=150) plt.close(fig3) print("図3保存: 2022_H5_15_fig3_dendrogram.png") # ================================================================ # ■ 図4: 地域別レーダーチャート（6地域ブロック × 4指標） # ================================================================ print("\n" + "=" * 65) print("■ 図4: 地域別レーダーチャート（6地域 × 4指標）") print("=" * 65) # レーダーチャートに使う4指標（解釈しやすいものを選択） RADAR_COLS = ['求人倍率', 'ごみリサイクル率', '大学進学率', '合計特殊出生率'] RADAR_LABELS = ['求人倍率\n（雇用）', 'ごみリサイクル率\n（環境）', '大学進学率\n（教育）', '合計特殊出生率\n（人口）'] # 地域別平均を計算（正規化用に全国平均・標準偏差を使用） df_radar = df_sdgs[['地方'] + RADAR_COLS].copy() region_means = df_radar.groupby('地方')[RADAR_COLS].mean() # 全国平均・SD で0〜1正規化（各列を min-max スケール） col_min = df_sdgs[RADAR_COLS].min() col_max = df_sdgs[RADAR_COLS].max() region_norm = (region_means - col_min) / (col_max - col_min) print("\n 地域別 正規化スコア（0〜1スケール）:") print(region_norm.round(3).to_string()) # レーダーチャートの角度計算 N_VARS = len(RADAR_COLS) angles = np.linspace(0, 2 * np.pi, N_VARS, endpoint=False).tolist() angles += angles[:1] # 閉じるために先頭を末尾に追加 fig4, ax4 = plt.subplots(figsize=(10, 8), subplot_kw=dict(projection='polar')) # 全体の平均をグレーの参照線として描画 national_mean = np.ones(N_VARS) * 0.5 # 正規化後の中央値（参考） vals_nat = national_mean.tolist() + national_mean[:1].tolist() ax4.plot(angles, vals_nat, 'k--', linewidth=1.0, alpha=0.4, label='_nolegend_') ax4.fill(angles, vals_nat, color='gray', alpha=0.08) # 各地域をプロット for region in REGION_ORDER: if region not in region_norm.index: continue vals = region_norm.loc[region, RADAR_COLS].values.tolist() vals += vals[:1] # 閉じる ax4.plot(angles, vals, color=REGION_COLORS[region], linewidth=2.2, label=region, alpha=0.9) ax4.fill(angles, vals, color=REGION_COLORS[region], alpha=0.12) # 軸ラベル設定 ax4.set_xticks(angles[:-1]) ax4.set_xticklabels(RADAR_LABELS, fontsize=10.5) ax4.set_ylim(0, 1) ax4.set_yticks([0.25, 0.5, 0.75]) ax4.set_yticklabels(['25%', '50%', '75%'], fontsize=8, color='#666') ax4.set_title( '地域別 SDGs4指標 レーダーチャート（2022年）\n' '〜各指標を全国min-max正規化。グレー点線は中央値（0.5）〜', fontsize=13, fontweight='bold', pad=20 ) ax4.legend(loc='lower left', bbox_to_anchor=(-0.18, -0.08), fontsize=9.5, framealpha=0.88) ax4.grid(True, alpha=0.3) plt.tight_layout() fig4.savefig(os.path.join(FIG_DIR, '2022_H5_15_fig4_radar.png'), bbox_inches='tight', dpi=150) plt.close(fig4) print("図4保存: 2022_H5_15_fig4_radar.png") # ================================================================ # ■ 完了メッセージ # ================================================================ print("\n" + "=" * 65) print("■ 全図生成完了（4枚）") print(f" fig1_pca_score.png : PCAスコアプロット（PC1 vs PC2, 地域色分け）") print(f" fig2_loadings.png : 変数の主成分負荷量ヒートマップ") print(f" fig3_dendrogram.png : Ward法クラスタリング デンドログラム") print(f" fig4_radar.png : 地域別レーダーチャート（6地域 × 4指標）") print("=" * 65) print(f"\n■ PCA寄与率サマリ:") for i in range(4): print(f" PC{i+1}: {ev_ratio[i]*100:.1f}% (累積: {ev_ratio[:i+1].sum()*100:.1f}%)") print(f"\n■ PC1スコア上位5都道府県（持続可能性が高い）:") print(top5[['都道府県', '地方', 'PC1']].to_string(index=False)) print(f"\n■ PC1スコア下位5都道府県（持続可能性が低い）:") print(bot5[['都道府県', '地方', 'PC1']].to_string(index=False))