""" 教育用再現コード: 2021年度 統計数理賞（大学生・一般の部） ================================================================= 論文タイトル：地域産業構造と経済成長：特化係数・主成分分析による産業集積の分析 受賞：統計数理賞（大学生・一般の部） 【分析概要】 データ：SSDSE-B-2026（47都道府県・2022年度） 本論文は地域産業構造を特化係数（LQ）で定量化し、主成分分析（PCA）で多次元的な 産業集積パターンを抽出、Ward法クラスタリングで地域類型を導出する。 SSDSE-Bには直接的な産業別就業者数がないため、入手可能な統計変数を代理指標として活用。 【変数設計（代理指標）】 - 観光特化 = 旅館営業施設数（ホテルを含む）/ 総人口 × 10000 観光産業の集積を表す代理変数（施設密度） - 農業特化代理 = 1000 / (1人1日当たりごみ排出量) 農業地域は都市化が低くごみ排出量が少ない傾向を利用 - 教育特化 = 大学学生数 / 総人口 × 1000 高等教育産業（大学）の集積密度 - 医療特化 = (一般病院数 + 一般診療所数) / 総人口 × 10000 医療・福祉産業の集積密度 - 建設特化代理 = 着工新設住宅戸数 / 総人口 × 10000 建設需要・活動の代理変数 - 消費支出_log = log(消費支出（二人以上の世帯）) 地域の消費活動水準（経済成長の代理） 【学習ポイント 4項目】 1. 特化係数（Location Quotient）の定義：地域の産業集積を全国平均との比率で測定 2. レーダーチャートの作成：matplotlib の polar axes を用いた多変数比較 3. PCAによる多次元産業構造の次元削減：因子負荷量の解釈と主成分スコア 4. Ward法クラスタリングによる地域類型化：政策グループとしての活用 【使用データ】 SSDSE-B-2026.csv（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/2021_U3_suri.py # ============================================================ import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt import matplotlib.patches as mpatches import numpy as np import pandas as pd from sklearn.preprocessing import StandardScaler from sklearn.decomposition import PCA from scipy.cluster.hierarchy import dendrogram, linkage from scipy import stats import warnings warnings.filterwarnings('ignore') # ────────────────────────────────────────────────────────────── # パス設定 # ────────────────────────────────────────────────────────────── import os FIG_DIR = 'html/figures' DATA_B = 'data/raw/SSDSE-B-2026.csv' os.makedirs(FIG_DIR, exist_ok=True) # ────────────────────────────────────────────────────────────── # 共通設定 # ────────────────────────────────────────────────────────────── plt.rcParams['font.family'] = 'Hiragino Sans' plt.rcParams['axes.unicode_minus'] = False plt.rcParams['figure.dpi'] = 150 PREFIX = '2021_U3' def save_fig(label): path = os.path.join(FIG_DIR, f'{PREFIX}_{label}.png') plt.savefig(path, bbox_inches='tight', dpi=150) plt.close() print(f' → {path} 保存完了') # ────────────────────────────────────────────────────────────── # 地域ブロック定義 # ────────────────────────────────────────────────────────────── REGION_MAP = { '北海道': '北海道', '青森県': '東北', '岩手県': '東北', '宮城県': '東北', '秋田県': '東北', '山形県': '東北', '福島県': '東北', '茨城県': '関東', '栃木県': '関東', '群馬県': '関東', '埼玉県': '関東', '千葉県': '関東', '東京都': '関東', '神奈川県': '関東', '新潟県': '中部', '富山県': '中部', '石川県': '中部', '福井県': '中部', '山梨県': '中部', '長野県': '中部', '岐阜県': '中部', '静岡県': '中部', '愛知県': '中部', '三重県': '近畿', '滋賀県': '近畿', '京都府': '近畿', '大阪府': '近畿', '兵庫県': '近畿', '奈良県': '近畿', '和歌山県': '近畿', '鳥取県': '中国・四国', '島根県': '中国・四国', '岡山県': '中国・四国', '広島県': '中国・四国', '山口県': '中国・四国', '徳島県': '中国・四国', '香川県': '中国・四国', '愛媛県': '中国・四国', '高知県': '中国・四国', '福岡県': '九州・沖縄', '佐賀県': '九州・沖縄', '長崎県': '九州・沖縄', '熊本県': '九州・沖縄', '大分県': '九州・沖縄', '宮崎県': '九州・沖縄', '鹿児島県': '九州・沖縄', '沖縄県': '九州・沖縄', } REGION_COLORS = { '北海道': '#E53935', '東北': '#FB8C00', '関東': '#43A047', '中部': '#1E88E5', '近畿': '#8E24AA', '中国・四国': '#00ACC1', '九州・沖縄': '#6D4C41', } REGION_ORDER = ['北海道', '東北', '関東', '中部', '近畿', '中国・四国', '九州・沖縄'] # ────────────────────────────────────────────────────────────── # データ読み込み・前処理 # ────────────────────────────────────────────────────────────── print('=== データ読み込み ===') df_raw = pd.read_csv(DATA_B, encoding='cp932', header=1) # 都道府県行のみ抽出（地域コードが R+5桁、かつ全国コードR00000を除外） mask = (df_raw['地域コード'].str.match(r'^R\d{5}$', na=False) & (df_raw['地域コード'] != 'R00000')) df_all = df_raw[mask].copy() # 2022年データ（消費支出が揃っている最新年） df = df_all[df_all['年度'] == 2022].copy().reset_index(drop=True) print(f'2022年データ: {len(df)}都道府県') # 数値変換 num_cols = ['総人口', '旅館営業施設数（ホテルを含む）', '大学学生数', '一般病院数', '一般診療所数', '着工新設住宅戸数', '1人1日当たりの排出量', '消費支出（二人以上の世帯）'] for c in num_cols: df[c] = pd.to_numeric(df[c], errors='coerce') # 欠損値確認 print('欠損値確認:') for c in num_cols: n_na = df[c].isna().sum() if n_na > 0: print(f' {c}: {n_na}件') # ────────────────────────────────────────────────────────────── # 特化係数（Location Quotient）代理指標の計算 # ────────────────────────────────────────────────────────────── print('\n=== 特化係数代理指標の計算 ===') pop = df['総人口'] # 観光特化: 旅館施設数 / 総人口 × 10000 df['観光特化'] = df['旅館営業施設数（ホテルを含む）'] / pop * 10000 # 農業特化代理: ごみ排出量の逆数（農村部はごみが少ない） df['農業特化代理'] = 1000 / df['1人1日当たりの排出量'] # 教育特化: 大学学生数 / 総人口 × 1000 df['教育特化'] = df['大学学生数'] / pop * 1000 # 医療特化: (病院数+診療所数) / 総人口 × 10000 df['医療特化'] = (df['一般病院数'] + df['一般診療所数']) / pop * 10000 # 建設特化代理: 新設住宅戸数 / 総人口 × 10000 df['建設特化代理'] = df['着工新設住宅戸数'] / pop * 10000 # 消費支出（log変換） df['消費支出_log'] = np.log(df['消費支出（二人以上の世帯）']) # 地域ブロック df['地域'] = df['都道府県'].map(REGION_MAP) # 変数リスト VARS = ['観光特化', '農業特化代理', '教育特化', '医療特化', '建設特化代理', '消費支出_log'] VAR_LABELS = ['観光\n特化', '農業\n特化代理', '教育\n特化', '医療\n特化', '建設\n特化代理', '消費支出\n(log)'] # 欠損除去 df_clean = df[['都道府県', '地域'] + VARS + ['消費支出（二人以上の世帯）']].dropna().reset_index(drop=True) print(f'分析対象: {len(df_clean)}都道府県') # 統計サマリー print('\n--- 代理指標サマリー ---') print(df_clean[VARS].describe().round(4).to_string()) # ────────────────────────────────────────────────────────────── # LQ正規化（全国平均=1.0 として比較） # ────────────────────────────────────────────────────────────── # 各指標を全国平均で割って相対特化度（LQ的）を算出 lq_df = df_clean[VARS].copy() for v in VARS: nat_mean = lq_df[v].mean() lq_df[v] = lq_df[v] / nat_mean # LQ: 地域値 / 全国平均 lq_df['地域'] = df_clean['地域'].values lq_df['都道府県'] = df_clean['都道府県'].values print('\n--- LQ値サマリー（全国平均=1.0）---') print(lq_df[VARS].describe().round(3).to_string()) # ────────────────────────────────────────────────────────────── # Fig1: レーダーチャート（6地域ブロック別平均 LQ） # ────────────────────────────────────────────────────────────── print('\n=== Fig1: レーダーチャート ===') region_lq = lq_df.groupby('地域')[VARS].mean() N_vars = len(VARS) angles = np.linspace(0, 2 * np.pi, N_vars, endpoint=False).tolist() angles += angles[:1] # close the polygon fig, axes = plt.subplots(2, 3, figsize=(14, 9), subplot_kw=dict(polar=True)) fig.suptitle('地域ブロック別 産業特化指標（LQ代理）レーダーチャート\n（2022年, 全国平均=1.0）', fontsize=14, fontweight='bold', y=1.01) axes_flat = axes.flatten() for i, region in enumerate(REGION_ORDER[:6]): ax = axes_flat[i] if region not in region_lq.index: ax.set_visible(False) continue vals = region_lq.loc[region, VARS].tolist() vals += vals[:1] color = REGION_COLORS[region] ax.plot(angles, vals, color=color, linewidth=2.0) ax.fill(angles, vals, color=color, alpha=0.25) # 参照線（LQ=1.0） ax.plot(angles, [1.0] * len(angles), color='gray', linewidth=0.8, linestyle='--', alpha=0.6) ax.set_xticks(angles[:-1]) ax.set_xticklabels(VAR_LABELS, fontsize=8) ax.set_ylim(0, None) ax.set_title(region, fontsize=11, fontweight='bold', color=color, pad=12) ax.tick_params(axis='y', labelsize=7) plt.tight_layout() save_fig('fig1_radar') # ────────────────────────────────────────────────────────────── # PCA 準備（標準化） # ────────────────────────────────────────────────────────────── print('\n=== PCA ===') X_raw = df_clean[VARS].values scaler = StandardScaler() X_scaled = scaler.fit_transform(X_raw) pca = PCA(n_components=min(len(VARS), len(df_clean))) scores = pca.fit_transform(X_scaled) loadings = pca.components_ # shape: (n_components, n_vars) ev_ratio = pca.explained_variance_ratio_ print('寄与率:') for i, ev in enumerate(ev_ratio[:4]): print(f' PC{i+1}: {ev*100:.1f}% (累積: {ev_ratio[:i+1].sum()*100:.1f}%)') # 因子負荷量 print('\n--- 因子負荷量（PC1, PC2）---') for j, v in enumerate(VARS): print(f' {v:12s}: PC1={loadings[0,j]:+.3f}, PC2={loadings[1,j]:+.3f}') # ────────────────────────────────────────────────────────────── # Fig2: PCAスコアプロット（PC1 vs PC2） # ────────────────────────────────────────────────────────────── print('\n=== Fig2: PCAスコアプロット ===') fig, ax = plt.subplots(figsize=(11, 8)) for region in REGION_ORDER: mask_r = df_clean['地域'].values == region if mask_r.sum() == 0: continue ax.scatter(scores[mask_r, 0], scores[mask_r, 1], color=REGION_COLORS[region], s=80, alpha=0.85, label=region, zorder=3) for idx in np.where(mask_r)[0]: pref = df_clean['都道府県'].iloc[idx] ax.annotate(pref, (scores[idx, 0], scores[idx, 1]), fontsize=7, xytext=(4, 3), textcoords='offset points', color=REGION_COLORS[region]) # 原点に十字線 ax.axhline(0, color='gray', linewidth=0.6, linestyle='--', alpha=0.5) ax.axvline(0, color='gray', linewidth=0.6, linestyle='--', alpha=0.5) ax.set_xlabel(f'PC1（寄与率 {ev_ratio[0]*100:.1f}%）', fontsize=12) ax.set_ylabel(f'PC2（寄与率 {ev_ratio[1]*100:.1f}%）', fontsize=12) ax.set_title('主成分分析（PCA）スコアプロット\n産業特化指標の都道府県別分布（2022年）', fontsize=13, fontweight='bold') ax.legend(title='地域ブロック', bbox_to_anchor=(1.01, 1), loc='upper left', fontsize=9, title_fontsize=9) # 因子負荷量の矢印（biplot） scale_factor = 2.5 for j, v in enumerate(VARS): lx, ly = loadings[0, j] * scale_factor, loadings[1, j] * scale_factor ax.annotate('', xy=(lx, ly), xytext=(0, 0), arrowprops=dict(arrowstyle='->', color='#E53935', lw=1.5)) ax.text(lx * 1.08, ly * 1.08, v.replace('代理', '代理\n').replace('特化', '特化\n'), fontsize=8, color='#E53935', ha='center') ax.grid(True, alpha=0.3) plt.tight_layout() save_fig('fig2_pca') # ────────────────────────────────────────────────────────────── # Fig3: Ward法クラスタリング デンドログラム # ────────────────────────────────────────────────────────────── print('\n=== Fig3: Ward法デンドログラム ===') linked = linkage(X_scaled, method='ward') fig, ax = plt.subplots(figsize=(14, 6)) dend = dendrogram( linked, labels=df_clean['都道府県'].values, ax=ax, leaf_rotation=90, leaf_font_size=8, color_threshold=0.7 * max(linked[:, 2]), ) ax.set_title('Ward法階層クラスタリング デンドログラム\n産業特化指標による都道府県類型（2022年）', fontsize=13, fontweight='bold') ax.set_xlabel('都道府県', fontsize=11) ax.set_ylabel('非類似度（Ward距離）', fontsize=11) # k=4 でカットする水準を示す補助線 # Ward距離の上位ギャップを探してカット点を示す sorted_dist = np.sort(linked[:, 2]) gaps = np.diff(sorted_dist) # 上位5ギャップのうち最大を採用（k=4カット） n_clusters_target = 4 cutoff_height = (linked[-n_clusters_target + 1, 2] + linked[-n_clusters_target, 2]) / 2 ax.axhline(y=cutoff_height, color='red', linestyle='--', linewidth=1.5, label=f'k={n_clusters_target} クラスタ分割線') ax.legend(fontsize=9) plt.tight_layout() save_fig('fig3_dendro') # ────────────────────────────────────────────────────────────── # Fig4: 観光特化 vs 消費支出（2022年）散布図 # ────────────────────────────────────────────────────────────── print('\n=== Fig4: 観光特化 vs 消費支出 散布図 ===') x_plot = df_clean['観光特化'].values y_plot = df_clean['消費支出（二人以上の世帯）'].values / 1000 # 千円単位 # 相関係数 r_val, p_val = stats.pearsonr(x_plot, y_plot) print(f'相関係数: r={r_val:.4f}, p={p_val:.4f}') # 回帰直線 slope, intercept, _, _, stderr = stats.linregress(x_plot, y_plot) x_fit = np.linspace(x_plot.min(), x_plot.max(), 100) y_fit = slope * x_fit + intercept fig, ax = plt.subplots(figsize=(10, 7)) for region in REGION_ORDER: mask_r = df_clean['地域'].values == region if mask_r.sum() == 0: continue ax.scatter(x_plot[mask_r], y_plot[mask_r], color=REGION_COLORS[region], s=80, alpha=0.85, label=region, zorder=3) for idx in np.where(mask_r)[0]: pref = df_clean['都道府県'].iloc[idx] ax.annotate(pref, (x_plot[idx], y_plot[idx]), fontsize=7, xytext=(4, 3), textcoords='offset points', color=REGION_COLORS[region]) # 回帰直線 line_style = '-' if p_val < 0.05 else '--' ax.plot(x_fit, y_fit, color='black', linewidth=1.8, linestyle=line_style, label=f'回帰直線 (r={r_val:.3f}, p={p_val:.3f})') ax.set_xlabel('観光特化指標（旅館施設数/万人）', fontsize=12) ax.set_ylabel('消費支出（千円, 二人以上世帯）', fontsize=12) ax.set_title('観光産業集積と消費支出水準の関係（2022年）\n都道府県別クロスセクション分析', fontsize=13, fontweight='bold') ax.legend(title='地域ブロック', bbox_to_anchor=(1.01, 1), loc='upper left', fontsize=9, title_fontsize=9) ax.grid(True, alpha=0.3) # 相関注釈 sig_str = '**' if p_val < 0.01 else ('*' if p_val < 0.05 else 'n.s.') ax.text(0.05, 0.95, f'r = {r_val:.3f} {sig_str}\np = {p_val:.3f}\nN = {len(df_clean)}', transform=ax.transAxes, fontsize=10, verticalalignment='top', bbox=dict(boxstyle='round', facecolor='lightyellow', alpha=0.8)) plt.tight_layout() save_fig('fig4_scatter') print('\n=== 全図の保存完了 ===') print('Fig1: 産業特化指標レーダーチャート（地域ブロック別）') print('Fig2: PCAスコアプロット（都道府県ラベル付き biplot）') print('Fig3: Ward法クラスタリング デンドログラム') print('Fig4: 観光特化 vs 消費支出 散布図（2022年）')