""" 2021年度 統計データ分析コンペティション 審査員奨励賞（大学生の部） 「ソーシャルキャピタルと地域コミュニティ活動の関係分析」 再現・教育用スクリプト 使用データ: SSDSE-B-2026.csv（47都道府県断面, 2022年度） """ # ============================================================ # 【データの準備】実行前に以下のデータファイルを用意してください # # 必要ファイル: # ・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_U5_5_shorei.py # ============================================================ import os import numpy as np import pandas as pd import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt import statsmodels.api as sm from scipy import stats from sklearn.preprocessing import StandardScaler from scipy.cluster.hierarchy import dendrogram, linkage from scipy.spatial.distance import pdist plt.rcParams['font.family'] = 'Hiragino Sans' plt.rcParams['axes.unicode_minus'] = False plt.rcParams['figure.dpi'] = 150 FIG_DIR = 'html/figures' DATA_B = 'data/raw/SSDSE-B-2026.csv' os.makedirs(FIG_DIR, exist_ok=True) # --------------------------------------------------------------- # データ読み込み # --------------------------------------------------------------- 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_b['年度'] = df_b['年度'].astype(int) # 2022年度断面データ df = df_b[df_b['年度'] == 2022].copy().reset_index(drop=True) assert len(df) == 47, f"47都道府県が揃っていません: {len(df)}件" # --------------------------------------------------------------- # 地域区分 # --------------------------------------------------------------- region_map = { '北海道': '北海道・東北', '青森県': '北海道・東北', '岩手県': '北海道・東北', '宮城県': '北海道・東北', '秋田県': '北海道・東北', '山形県': '北海道・東北', '福島県': '北海道・東北', '茨城県': '関東', '栃木県': '関東', '群馬県': '関東', '埼玉県': '関東', '千葉県': '関東', '東京都': '関東', '神奈川県': '関東', '新潟県': '中部', '富山県': '中部', '石川県': '中部', '福井県': '中部', '山梨県': '中部', '長野県': '中部', '岐阜県': '中部', '静岡県': '中部', '愛知県': '中部', '三重県': '近畿', '滋賀県': '近畿', '京都府': '近畿', '大阪府': '近畿', '兵庫県': '近畿', '奈良県': '近畿', '和歌山県': '近畿', '鳥取県': '中国・四国', '島根県': '中国・四国', '岡山県': '中国・四国', '広島県': '中国・四国', '山口県': '中国・四国', '徳島県': '中国・四国', '香川県': '中国・四国', '愛媛県': '中国・四国', '高知県': '中国・四国', '福岡県': '九州・沖縄', '佐賀県': '九州・沖縄', '長崎県': '九州・沖縄', '熊本県': '九州・沖縄', '大分県': '九州・沖縄', '宮崎県': '九州・沖縄', '鹿児島県': '九州・沖縄', '沖縄県': '九州・沖縄' } region_colors = { '北海道・東北': '#4e9af1', '関東': '#e05c5c', '中部': '#f0a500', '近畿': '#5cb85c', '中国・四国': '#9b59b6', '九州・沖縄': '#f39c12' } df['地域'] = df['都道府県'].map(region_map) # --------------------------------------------------------------- # 派生変数の作成（人口1万人あたり等） # --------------------------------------------------------------- pop = df['総人口'] # 単位: 人 # 保育所密度（保育所等数 / 総人口 × 10,000） df['保育所密度'] = df['保育所等数'] / pop * 10000 # 高齢化率（65歳以上 / 総人口） df['高齢化率'] = df['65歳以上人口'] / pop * 100 # 年少人口率（15歳未満 / 総人口） df['年少人口率'] = df['15歳未満人口'] / pop * 100 # 医療施設密度（一般病院数 + 診療所数 / 万人） df['医療施設密度'] = (df['一般病院数'] + df['一般診療所数']) / pop * 10000 # リサイクル率（環境への意識 = ソーシャルキャピタルの代理） df['リサイクル率'] = df['ごみのリサイクル率'] # 教育進学率（高校卒業→進学率） df['進学率'] = df['高等学校卒業者のうち進学者数'] / df['高等学校卒業者数'] * 100 # 婚姻率（人口千人あたり） df['婚姻率'] = df['婚姻件数'] / pop * 1000 # 消費支出に占める教養娯楽費率（文化活動への関与） df['教養娯楽費率'] = df['教養娯楽費（二人以上の世帯）'] / df['消費支出（二人以上の世帯）'] * 100 # 消費支出に占める教育費率 df['教育費率'] = df['教育費（二人以上の世帯）'] / df['消費支出（二人以上の世帯）'] * 100 # 旅行・観光活動（延べ宿泊者 / 万人） df['観光密度'] = df['延べ宿泊者数'] / pop * 10000 # 保育所定員充足率（在所児 / 定員） df['保育充足率'] = df['保育所等在所児数'] / df['保育所等定員数'] * 100 # --------------------------------------------------------------- # 分析用データフレーム（欠損値を含む行を除外） # --------------------------------------------------------------- vars_all = [ '合計特殊出生率', '保育所密度', '高齢化率', '年少人口率', '医療施設密度', 'リサイクル率', '進学率', '婚姻率', '教養娯楽費率', '教育費率', '観光密度', '保育充足率' ] df_ana = df[['都道府県', '地域'] + vars_all].dropna().copy() print(f"分析対象: {len(df_ana)}都道府県") # --------------------------------------------------------------- # 相関行列の計算（ヒートマップ用） # --------------------------------------------------------------- sc_vars = [ '保育所密度', '高齢化率', 'リサイクル率', '進学率', '婚姻率', '教養娯楽費率', '教育費率', '医療施設密度' ] corr_labels = [ '保育所密度', '高齢化率', 'リサイクル率', '進学率', '婚姻率', '教養娯楽費率', '教育費率', '医療施設密度' ] corr_matrix = df_ana[sc_vars].corr() # --------------------------------------------------------------- # 重回帰分析 # --------------------------------------------------------------- y_var = '合計特殊出生率' x_vars = ['保育所密度', '高齢化率', 'リサイクル率', '進学率', '婚姻率', '教養娯楽費率', '医療施設密度'] y = df_ana[y_var] X_raw = df_ana[x_vars] # 標準化 scaler = StandardScaler() X_std = scaler.fit_transform(X_raw) X_std_df = pd.DataFrame(X_std, columns=x_vars, index=df_ana.index) X_std_df = sm.add_constant(X_std_df) ols_model = sm.OLS(y, X_std_df).fit() print(ols_model.summary()) # 標準化偏回帰係数と95%CI coef = ols_model.params[1:] # 定数項を除く ci = ols_model.conf_int().iloc[1:] pval = ols_model.pvalues[1:] r2 = ols_model.rsquared r2_adj = ols_model.rsquared_adj f_stat = ols_model.fvalue f_pval = ols_model.f_pvalue print(f"\nR² = {r2:.3f}, Adj R² = {r2_adj:.3f}") print(f"F統計量 = {f_stat:.2f}, p値 = {f_pval:.4f}") # --------------------------------------------------------------- # クラスタリング用変数 # --------------------------------------------------------------- cluster_vars = ['保育所密度', '高齢化率', 'リサイクル率', '婚姻率', '医療施設密度', '進学率'] X_clust = df_ana[cluster_vars].values X_clust_std = StandardScaler().fit_transform(X_clust) Z = linkage(X_clust_std, method='ward') # 個々の相関係数（fig1用） r_birth_hoiku, p_birth_hoiku = stats.pearsonr(df_ana['保育所密度'], df_ana['合計特殊出生率']) print(f"\n保育所密度 vs 合計特殊出生率: r={r_birth_hoiku:.3f}, p={p_birth_hoiku:.4f}") # ================================================================ # 図1: 散布図 — 保育所密度 vs 合計特殊出生率 # ================================================================ fig, ax = plt.subplots(figsize=(9, 6.5)) for region, grp in df_ana.groupby('地域'): color = region_colors[region] ax.scatter(grp['保育所密度'], grp['合計特殊出生率'], color=color, s=60, alpha=0.85, label=region, zorder=3) # ラベル for _, row in df_ana.iterrows(): ax.annotate(row['都道府県'], xy=(row['保育所密度'], row['合計特殊出生率']), fontsize=6.5, ha='left', va='bottom', xytext=(2, 2), textcoords='offset points', color='#333333') # 回帰線 x_fit = np.linspace(df_ana['保育所密度'].min(), df_ana['保育所密度'].max(), 200) slope, intercept, r_val, p_val, se = stats.linregress(df_ana['保育所密度'], df_ana['合計特殊出生率']) y_fit = slope * x_fit + intercept ax.plot(x_fit, y_fit, color='#333333', linewidth=1.5, linestyle='--', zorder=2) ax.set_xlabel('保育所密度（保育所等数 / 人口1万人）', fontsize=12) ax.set_ylabel('合計特殊出生率', fontsize=12) ax.set_title('図1: 保育所密度と合計特殊出生率の関係（2022年度, 47都道府県）', fontsize=13, pad=12) # 凡例 handles, labels = ax.get_legend_handles_labels() ax.legend(handles, labels, loc='upper right', fontsize=9, framealpha=0.9) # 相関係数注記 ax.text(0.03, 0.97, f'r = {r_val:.3f} (p = {p_val:.4f})', transform=ax.transAxes, fontsize=10, va='top', ha='left', bbox=dict(boxstyle='round,pad=0.4', fc='white', ec='#aaaaaa', alpha=0.9)) ax.grid(True, alpha=0.3) plt.tight_layout() fig.savefig(os.path.join(FIG_DIR, '2021_U5_5_fig1.png'), bbox_inches='tight') plt.close(fig) print("fig1 saved") # ================================================================ # 図2: 相関ヒートマップ # ================================================================ fig, ax = plt.subplots(figsize=(9, 7.5)) n = len(sc_vars) im = ax.imshow(corr_matrix.values, cmap='RdBu_r', vmin=-1, vmax=1, aspect='auto') plt.colorbar(im, ax=ax, shrink=0.8, label='Pearson r') ax.set_xticks(range(n)) ax.set_yticks(range(n)) ax.set_xticklabels(corr_labels, rotation=40, ha='right', fontsize=10) ax.set_yticklabels(corr_labels, fontsize=10) for i in range(n): for j in range(n): val = corr_matrix.values[i, j] text_color = 'white' if abs(val) > 0.6 else 'black' ax.text(j, i, f'{val:.2f}', ha='center', va='center', fontsize=9, color=text_color, fontweight='bold') ax.set_title('図2: ソーシャルキャピタル関連変数のPearson相関行列（2022年度）', fontsize=12, pad=12) plt.tight_layout() fig.savefig(os.path.join(FIG_DIR, '2021_U5_5_fig2.png'), bbox_inches='tight') plt.close(fig) print("fig2 saved") # ================================================================ # 図3: 標準化偏回帰係数プロット（横棒グラフ） # ================================================================ labels_jp = { '保育所密度': '保育所密度\n（保育所等数/万人）', '高齢化率': '高齢化率', 'リサイクル率': 'リサイクル率\n（環境意識）', '進学率': '進学率\n（教育志向）', '婚姻率': '婚姻率', '教養娯楽費率': '教養娯楽費率\n（文化活動）', '医療施設密度': '医療施設密度', } order = coef.abs().sort_values(ascending=True).index coef_sorted = coef[order] ci_sorted = ci.loc[order] pval_sorted = pval[order] ylabels = [labels_jp.get(v, v) for v in order] colors_bar = ['#e05c5c' if p < 0.05 else '#aec6e8' for p in pval_sorted] fig, ax = plt.subplots(figsize=(9, 5.5)) y_pos = np.arange(len(order)) bars = ax.barh(y_pos, coef_sorted.values, color=colors_bar, height=0.55, zorder=3) ax.errorbar(coef_sorted.values, y_pos, xerr=np.array([(coef_sorted.values - ci_sorted.iloc[:, 0].values), (ci_sorted.iloc[:, 1].values - coef_sorted.values)]), fmt='none', color='#333333', linewidth=1.3, capsize=4, zorder=4) ax.axvline(0, color='black', linewidth=1.0, zorder=2) ax.set_yticks(y_pos) ax.set_yticklabels(ylabels, fontsize=10) ax.set_xlabel('標準化偏回帰係数（β）', fontsize=11) ax.set_title(f'図3: 標準化偏回帰係数と95%信頼区間\n（目的変数: 合計特殊出生率, R²={r2:.3f}）', fontsize=12, pad=10) from matplotlib.patches import Patch legend_elements = [ Patch(facecolor='#e05c5c', label='有意 (p < 0.05)'), Patch(facecolor='#aec6e8', label='非有意 (p >= 0.05)') ] ax.legend(handles=legend_elements, loc='lower right', fontsize=10) ax.grid(True, axis='x', alpha=0.35) plt.tight_layout() fig.savefig(os.path.join(FIG_DIR, '2021_U5_5_fig3.png'), bbox_inches='tight') plt.close(fig) print("fig3 saved") # ================================================================ # 図4: Ward法デンドログラム # ================================================================ fig, ax = plt.subplots(figsize=(14, 6)) pref_names = df_ana['都道府県'].tolist() dn = dendrogram(Z, labels=pref_names, ax=ax, leaf_rotation=90, leaf_font_size=9, color_threshold=Z[-4, 2]) ax.set_title('図4: Ward法による47都道府県のコミュニティ特性クラスタリング（2022年度）', fontsize=13, pad=12) ax.set_ylabel('クラスター間距離（Ward距離）', fontsize=11) ax.set_xlabel('都道府県', fontsize=11) ax.grid(True, axis='y', alpha=0.3) plt.tight_layout() fig.savefig(os.path.join(FIG_DIR, '2021_U5_5_fig4.png'), bbox_inches='tight') plt.close(fig) print("fig4 saved") # ================================================================ # 統計サマリー出力（HTML埋め込み用） # ================================================================ print("\n========== HTML用 統計サマリー ==========") print(f"分析年度: 2022年度") print(f"分析対象: {len(df_ana)}都道府県") print(f"決定係数 R²: {r2:.3f}") print(f"自由度調整済み R²: {r2_adj:.3f}") print(f"F統計量: {f_stat:.2f} p値: {f_pval:.4f}") print(f"保育所密度 vs 合計特殊出生率: r = {r_birth_hoiku:.3f}, p = {p_birth_hoiku:.4f}") print("\n-- 標準化偏回帰係数 --") for v in x_vars: b = ols_model.params[v] p = ols_model.pvalues[v] ci_lo = ols_model.conf_int().loc[v, 0] ci_hi = ols_model.conf_int().loc[v, 1] sig = '*' if p < 0.05 else '' print(f" {v}: β={b:.3f} CI=[{ci_lo:.3f}, {ci_hi:.3f}] p={p:.4f} {sig}") print("\n-- 相関行列（ソーシャルキャピタル変数）--") for i, v1 in enumerate(sc_vars): for j, v2 in enumerate(sc_vars): if j > i: r, p = stats.pearsonr(df_ana[v1], df_ana[v2]) if abs(r) > 0.4: sig = '**' if p < 0.01 else ('*' if p < 0.05 else '') print(f" {v1} ↔ {v2}: r={r:.3f} p={p:.4f} {sig}") print("\n-- 都道府県別クラスター（閾値 = Z[-4,2]）--") from scipy.cluster.hierarchy import fcluster labels_clust = fcluster(Z, t=4, criterion='maxclust') for i, pref in enumerate(pref_names): print(f" {pref}: クラスター{labels_clust[i]}") print("\n全図の保存先:", FIG_DIR) print("スクリプト完了")