""" 2024_U5_1_shorei.py ごみの削減とリサイクルを推進する要因――環境ボランティアは、ごみ削減の効果を持つのか―― 審査員奨励賞 [大学生・一般の部] 市村遼ほか（中央大学商学部・文学部） 実データ（SSDSE-B-2026.csv）による教育用再現コード ※合成データ・乱数生成は一切使用しない """ # ============================================================ # 【データの準備】実行前に以下のデータファイルを用意してください # # 必要ファイル: # ・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/2024_U5_1_shorei.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 scipy import stats from numpy.linalg import lstsq plt.rcParams['font.family'] = 'Hiragino Sans' plt.rcParams['axes.unicode_minus'] = False import os FIGDIR = os.path.normpath('html/figures') + os.sep DATA_PATH = 'data/raw/SSDSE-B-2026.csv' os.makedirs(FIGDIR, exist_ok=True) # ---------------------------------------------------------------- # 実データ読み込み: SSDSE-B-2026（2022年度、都道府県レベル） # ---------------------------------------------------------------- df_b = pd.read_csv(DATA_PATH, encoding='cp932', header=1) # 都道府県レベルのみ（地域コードが R+5桁数字） df_b = df_b[df_b['地域コード'].str.match(r'^R\d{5}$', na=False)] df = df_b[df_b['年度'] == 2022].copy().reset_index(drop=True) # 数値変換 for col in ['総人口', '65歳以上人口', '小学校数', '一般病院数', '保育所等数', '1人1日当たりの排出量', 'ごみのリサイクル率', '消費支出（二人以上の世帯）']: df[col] = pd.to_numeric(df[col], errors='coerce') # ---------------------------------------------------------------- # 派生変数の作成 # ---------------------------------------------------------------- # 高齢化率: 65歳以上人口 / 総人口 df['高齢化率'] = df['65歳以上人口'] / df['総人口'] # 小学校数割合（人口1万人あたり）: 地域コミュニティ密度・環境ボランティア代理変数 df['小学校数割合\n（1万人あたり）'] = df['小学校数'] / df['総人口'] * 10000 # 病院数割合（人口10万人あたり）: 都市化率の代理変数 df['病院数割合\n（10万人あたり）'] = df['一般病院数'] / df['総人口'] * 100000 # 保育所数割合（人口1万人あたり）: 都市インフラ整備度の代理変数 df['保育所数割合\n（1万人あたり）'] = df['保育所等数'] / df['総人口'] * 10000 # 消費支出（円）: 所得水準の代理変数 df['消費支出\n（円/月）'] = df['消費支出（二人以上の世帯）'] # 目的変数 df['1人1日排出量\n（g/人・日）'] = df['1人1日当たりの排出量'] df['リサイクル率\n（%）'] = df['ごみのリサイクル率'] # 欠損除去 feature_cols = [ '高齢化率', '小学校数割合\n（1万人あたり）', '病院数割合\n（10万人あたり）', '保育所数割合\n（1万人あたり）', '消費支出\n（円/月）', ] target_cols = ['1人1日排出量\n（g/人・日）', 'リサイクル率\n（%）'] all_cols = feature_cols + target_cols df = df.dropna(subset=all_cols).reset_index(drop=True) # 短い表示名（軸ラベル用） feature_labels = [ '高齢化率', '小学校数割合\n(1万人あたり)', '病院数割合\n(10万人あたり)', '保育所数割合\n(1万人あたり)', '消費支出\n(円/月)', ] prefectures = df['都道府県'].tolist() n = len(df) print(f"分析対象: {n} 都道府県（2022年度）") X = df[feature_cols].values y_waste = df['1人1日排出量\n（g/人・日）'].values y_recycle = df['リサイクル率\n（%）'].values # 標準化 X_mean = X.mean(axis=0) X_std_v = X.std(axis=0, ddof=0) X_std = (X - X_mean) / X_std_v # ---------------------------------------------------------------- # OLS 回帰関数 # ---------------------------------------------------------------- def ols_full(X_mat, y): X_c = np.column_stack([np.ones(len(y)), X_mat]) b, _, _, _ = lstsq(X_c, y, rcond=None) res = y - X_c @ b n_r, k = X_c.shape sigma2 = res @ res / (n_r - k) cov = sigma2 * np.linalg.inv(X_c.T @ X_c) se = np.sqrt(np.diag(cov))[1:] t_s = b[1:] / se p_v = 2 * (1 - stats.t.cdf(np.abs(t_s), df=n_r - k)) rss = res @ res ss_tot = np.sum((y - y.mean()) ** 2) r2 = 1 - rss / ss_tot intercept = b[0] return b[1:], se, t_s, p_v, r2, intercept coef_w, se_w, t_w, p_w, r2_w, int_w = ols_full(X_std, y_waste) coef_r, se_r, t_r, p_r, r2_r, int_r = ols_full(X_std, y_recycle) # ---------------------------------------------------------------- # 回帰結果の表示 # ---------------------------------------------------------------- def sig_star(p): if p < 0.001: return "***" if p < 0.01: return "**" if p < 0.05: return "*" return "n.s." print("\n=== OLS 回帰結果: 1人1日排出量 ===") print(f"R² = {r2_w:.4f}") print(f"{'変数':<28} {'標準化係数':>10} {'SE':>8} {'t値':>8} {'p値':>8} {'有意':>5}") short_names = ['高齢化率', '小学校数割合', '病院数割合', '保育所数割合', '消費支出'] for nm, c, s, t, p in zip(short_names, coef_w, se_w, t_w, p_w): print(f" {nm:<26} {c:>10.4f} {s:>8.4f} {t:>8.4f} {p:>8.4f} {sig_star(p):>5}") print("\n=== OLS 回帰結果: リサイクル率 ===") print(f"R² = {r2_r:.4f}") print(f"{'変数':<28} {'標準化係数':>10} {'SE':>8} {'t値':>8} {'p値':>8} {'有意':>5}") for nm, c, s, t, p in zip(short_names, coef_r, se_r, t_r, p_r): print(f" {nm:<26} {c:>10.4f} {s:>8.4f} {t:>8.4f} {p:>8.4f} {sig_star(p):>5}") print("\n=== 相関係数（各説明変数 × 目的変数）===") print(f"{'変数':<20} {'r（排出量）':>12} {'p':>8} {'r（リサイクル率）':>16} {'p':>8}") for nm, fc in zip(short_names, feature_cols): r1, p1 = stats.pearsonr(df[fc], y_waste) r2, p2 = stats.pearsonr(df[fc], y_recycle) print(f" {nm:<18} {r1:>12.4f} {p1:>8.4f} {r2:>16.4f} {p2:>8.4f}") # ================================================================ # 図1: 相関行列ヒートマップ（全変数） # ================================================================ corr_cols_labels = feature_labels + ['1人1日排出量\n（g/人・日）', 'リサイクル率\n（%）'] corr_data = df[feature_cols + target_cols].copy() corr_data.columns = corr_cols_labels corr_mat = corr_data.corr() fig, ax = plt.subplots(figsize=(9, 7)) im = ax.imshow(corr_mat.values, cmap='RdBu_r', vmin=-1, vmax=1, aspect='auto') plt.colorbar(im, ax=ax, shrink=0.8, label='Pearson r') n_vars = len(corr_cols_labels) ax.set_xticks(range(n_vars)) ax.set_yticks(range(n_vars)) ax.set_xticklabels(corr_cols_labels, fontsize=9, rotation=45, ha='right') ax.set_yticklabels(corr_cols_labels, fontsize=9) for i in range(n_vars): for j in range(n_vars): val = corr_mat.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=8, color=text_color) ax.set_title("図1: 相関行列ヒートマップ（47都道府県, 2022年度）", fontsize=13, pad=12) plt.tight_layout() plt.savefig(FIGDIR + "2024_U5_1_fig1_scatter.png", dpi=150) plt.close() print("\nfig1 saved") # ================================================================ # 図2: 説明変数 vs 1人1日排出量（回帰直線付き散布図） # ================================================================ fig, axes = plt.subplots(2, 3, figsize=(14, 9)) axes_flat = axes.flatten() highlight = {'北海道', '東京都', '大阪府', '秋田県', '沖縄県'} for idx, (fc, label) in enumerate(zip(feature_cols, feature_labels)): ax = axes_flat[idx] x_vals = df[fc].values y_vals = y_waste ax.scatter(x_vals, y_vals, color='#1565C0', alpha=0.65, s=50, zorder=3) slope, intercept_lr, r_val, p_val, _ = stats.linregress(x_vals, y_vals) x_line = np.linspace(x_vals.min(), x_vals.max(), 200) ax.plot(x_line, intercept_lr + slope * x_line, color='#C62828', lw=2, label=f'r={r_val:.2f}, p={p_val:.3f}') for pref, xi, yi in zip(prefectures, x_vals, y_vals): if pref in highlight: ax.annotate(pref, (xi, yi), xytext=(xi, yi + 8), fontsize=8, ha='center', arrowprops=dict(arrowstyle='->', color='gray', lw=0.8)) xlabel = label.replace('\n', ' ') ax.set_xlabel(xlabel, fontsize=9) ax.set_ylabel("1人1日排出量（g/人・日）", fontsize=9) ax.set_title(f"{xlabel}\nvs 排出量", fontsize=9) ax.legend(fontsize=8) ax.grid(True, alpha=0.3) # 6番目のパネルは空白にして説明テキスト ax6 = axes_flat[5] ax6.axis('off') ax6.text(0.5, 0.7, "【データ出典】\nSSDSE-B-2026\n（2022年度, 47都道府県）", ha='center', va='center', fontsize=11, transform=ax6.transAxes, bbox=dict(boxstyle='round', facecolor='#E3F2FD', alpha=0.8)) ax6.text(0.5, 0.3, "【変数の意味】\n高齢化率 → 人口高齢化\n小学校数割合 → コミュニティ密度\n" "病院数割合 → 都市化度\n保育所数割合 → インフラ整備\n消費支出 → 所得水準", ha='center', va='center', fontsize=9, transform=ax6.transAxes) fig.suptitle("図2: 各説明変数 vs 1人1日ごみ排出量（47都道府県, 2022年度）", fontsize=13) plt.tight_layout() plt.savefig(FIGDIR + "2024_U5_1_fig2_vif.png", dpi=150) plt.close() print("fig2 saved") # ================================================================ # 図3: 重回帰係数の比較（ごみ排出量 vs リサイクル率） # ================================================================ fig, axes = plt.subplots(1, 2, figsize=(13, 5)) label_short = ['高齢化率', '小学校数割合\n(1万人あたり)', '病院数割合\n(10万人あたり)', '保育所数割合\n(1万人あたり)', '消費支出\n(円/月)'] for ax, coefs, ses, p_vals_m, title, r2_v in [ (axes[0], coef_w, se_w, p_w, f"目的変数①: 1人1日排出量（g/人・日）\n(R²={r2_w:.3f})", r2_w), (axes[1], coef_r, se_r, p_r, f"目的変数②: リサイクル率（%）\n(R²={r2_r:.3f})", r2_r), ]: colors_c = ['#C62828' if c > 0 else '#1565C0' for c in coefs] bars = ax.barh(label_short, coefs, color=colors_c, alpha=0.8, xerr=1.96 * ses, capsize=4) ax.axvline(0, color='black', lw=1.2) for i, (c, se_i, p_v) in enumerate(zip(coefs, ses, p_vals_m)): star = sig_star(p_v) offset = max(abs(c) * 0.05, 0.05) ha_pos = 'left' if c >= 0 else 'right' x_pos = c + (offset if c >= 0 else -offset) color_s = '#C62828' if star != 'n.s.' else 'gray' ax.text(x_pos, i, star, va='center', ha=ha_pos, fontsize=11, color=color_s) ax.set_xlabel("標準化回帰係数（95% CI）", fontsize=10) ax.set_title(title, fontsize=10) ax.grid(True, axis='x', alpha=0.3) red_patch = mpatches.Patch(color='#C62828', alpha=0.8, label='正の係数（排出増・リサイクル増）') blue_patch = mpatches.Patch(color='#1565C0', alpha=0.8, label='負の係数（排出減・リサイクル減）') ax.legend(handles=[red_patch, blue_patch], fontsize=8, loc='lower right') fig.suptitle("図3: 重回帰分析の標準化係数（47都道府県, 2022年度）\n" "*** p<0.001 ** p<0.01 * p<0.05 n.s. 有意でない", fontsize=12) plt.tight_layout() plt.savefig(FIGDIR + "2024_U5_1_fig3_coef.png", dpi=150) plt.close() print("fig3 saved") # ================================================================ # 図4: リサイクル率 vs 消費支出（散布図）+ 残差ドットプロット # ================================================================ fig, axes = plt.subplots(1, 2, figsize=(13, 6)) # 左: リサイクル率 vs 消費支出 ax = axes[0] x_spend = df['消費支出\n（円/月）'].values / 10000 # 万円単位 slope_s, intercept_s, r_s, p_s, _ = stats.linregress(x_spend, y_recycle) x_line_s = np.linspace(x_spend.min(), x_spend.max(), 200) ax.scatter(x_spend, y_recycle, color='#2E7D32', alpha=0.7, s=55, zorder=3) ax.plot(x_line_s, intercept_s + slope_s * x_line_s, color='#C62828', lw=2.5, label=f'回帰直線 (r={r_s:.2f}, p={p_s:.3f})') for pref, xi, yi in zip(prefectures, x_spend, y_recycle): if pref in {'東京都', '大阪府', '秋田県', '富山県', '沖縄県', '愛知県'}: ax.annotate(pref, (xi, yi), xytext=(xi + 0.02, yi + 0.4), fontsize=8, arrowprops=dict(arrowstyle='->', color='gray', lw=0.8)) ax.set_xlabel("消費支出（万円/月）", fontsize=11) ax.set_ylabel("リサイクル率（%）", fontsize=11) ax.set_title("リサイクル率 vs 消費支出\n（所得水準の代理変数）", fontsize=11) ax.legend(fontsize=10) ax.grid(True, alpha=0.3) # 右: ごみ排出量モデルの都道府県別残差 X_c_full = np.column_stack([np.ones(n), X_std]) bw_full, _, _, _ = lstsq(X_c_full, y_waste, rcond=None) residuals_w = y_waste - X_c_full @ bw_full ax2 = axes[1] sorted_idx = np.argsort(residuals_w) sorted_prefs = [prefectures[i] for i in sorted_idx] sorted_res = residuals_w[sorted_idx] colors_res = ['#C62828' if r > 0 else '#1565C0' for r in sorted_res] ax2.barh(sorted_prefs, sorted_res, color=colors_res, alpha=0.8) ax2.axvline(0, color='black', lw=1.5) ax2.set_xlabel("残差（g/人・日）", fontsize=10) ax2.set_title("ごみ排出量モデルの都道府県別残差\n（赤=過大推定, 青=過小推定）", fontsize=10) ax2.grid(True, axis='x', alpha=0.3) ax2.tick_params(axis='y', labelsize=7.5) fig.suptitle("図4: リサイクル率と消費支出の関係 / ごみ排出量モデルの残差（47都道府県, 2022年度）", fontsize=12) plt.tight_layout() plt.savefig(FIGDIR + "2024_U5_1_fig4_map.png", dpi=150) plt.close() print("fig4 saved") print("\nAll figures saved. 実データ（SSDSE-B-2026）による分析完了。")