""" 2019_H5_3_shorei.py 都道府県別ごみ排出量とリサイクル率の決定要因：EKC仮説の検証 審査員奨励賞（高校生部門）2019年度 使用データ: SSDSE-B-2026.csv (実データのみ・合成データ禁止) 【分析概要】 データ : SSDSE-B-2026.csv（都道府県パネル、47都道府県×12年度） 目的変数 : 1人1日あたりごみ排出量（g/人/日）& ごみのリサイクル率（%） 所得代理 : 消費支出（二人以上の世帯）（円/月） EKC仮説 : 所得（消費支出）の上昇に伴いごみ排出量が逆U字型に推移するか検証（2次項） 【変数】 1人1日排出量 ← 消費支出（1次・2次項）＋ 光熱水道費 ＋ 高齢化率 ＋ 年平均気温 リサイクル率 ← 消費支出（Pearson相関） 【データ出典】 SSDSE-B-2026.csv: 社会・人口統計体系（都道府県データ） 統計センターより公表の実データ 【DS学習ポイント】 1. EKC仮説（環境クズネッツ曲線）の概念と2次項回帰の解釈 2. 2次項を含む回帰モデルの転換点計算 3. 寒冷地効果（北海道・東北の光熱水道費とごみ排出） 4. 3R政策と日本のごみ削減トレンド（時系列） """ # ============================================================ # 【データの準備】実行前に以下のデータファイルを用意してください # # 必要ファイル: # ・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/2019_H5_3_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 matplotlib.patches import Patch 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) print("=== 利用可能な列 ===") print(df_b.columns.tolist()) print(f"\n年度: {sorted(df_b['年度'].unique())}") print(f"都道府県数: {df_b['都道府県'].nunique()}") # 列名エイリアス COL_PERCAP = '1人1日当たりの排出量' # g/人/日 COL_RECYCLE = 'ごみのリサイクル率' # % COL_CONSUME = '消費支出（二人以上の世帯）' # 円/月 COL_ENERGY = '光熱・水道費（二人以上の世帯）' # 円/月 COL_POP = '総人口' COL_AGED = '65歳以上人口' COL_TEMP = '年平均気温' COL_PREF = '都道府県' COL_YEAR = '年度' # 都道府県名の正規化（語尾除去） def clean_pref(s): for suffix in ['県', '府', '都', '道']: if s.endswith(suffix): return s[:-1] return s df_b['都道府県短'] = df_b[COL_PREF].apply(clean_pref) # ─── 地域区分 ───────────────────────────────────────────────────── region_map = { '北海道': '北海道・東北', '青森': '北海道・東北', '岩手': '北海道・東北', '宮城': '北海道・東北', '秋田': '北海道・東北', '山形': '北海道・東北', '福島': '北海道・東北', '茨城': '関東', '栃木': '関東', '群馬': '関東', '埼玉': '関東', '千葉': '関東', '東京': '関東', '神奈川': '関東', '新潟': '中部', '富山': '中部', '石川': '中部', '福井': '中部', '山梨': '中部', '長野': '中部', '岐阜': '中部', '静岡': '中部', '愛知': '中部', '三重': '近畿', '滋賀': '近畿', '京都': '近畿', '大阪': '近畿', '兵庫': '近畿', '奈良': '近畿', '和歌山': '近畿', '鳥取': '中国・四国', '島根': '中国・四国', '岡山': '中国・四国', '広島': '中国・四国', '山口': '中国・四国', '徳島': '中国・四国', '香川': '中国・四国', '愛媛': '中国・四国', '高知': '中国・四国', '福岡': '九州・沖縄', '佐賀': '九州・沖縄', '長崎': '九州・沖縄', '熊本': '九州・沖縄', '大分': '九州・沖縄', '宮崎': '九州・沖縄', '鹿児島': '九州・沖縄', '沖縄': '九州・沖縄' } region_colors = { '北海道・東北': '#4e9af1', '関東': '#e05c5c', '中部': '#f0a500', '近畿': '#5cb85c', '中国・四国': '#9b59b6', '九州・沖縄': '#f39c12' } region_order = ['北海道・東北', '関東', '中部', '近畿', '中国・四国', '九州・沖縄'] df_b['地域'] = df_b['都道府県短'].map(region_map) # ─── 断面データ（最新年：2023年）──────────────────────────────── LATEST = df_b[COL_YEAR].max() df_cross = df_b[df_b[COL_YEAR] == LATEST].copy() print(f"\n最新年: {LATEST}, N={len(df_cross)}") # 派生変数の計算 df_cross['高齢化率'] = pd.to_numeric(df_cross[COL_AGED], errors='coerce') / \ pd.to_numeric(df_cross[COL_POP], errors='coerce') * 100 df_cross['消費支出'] = pd.to_numeric(df_cross[COL_CONSUME], errors='coerce') df_cross['光熱水道費'] = pd.to_numeric(df_cross[COL_ENERGY], errors='coerce') df_cross['気温'] = pd.to_numeric(df_cross[COL_TEMP], errors='coerce') df_cross['1人排出量'] = pd.to_numeric(df_cross[COL_PERCAP], errors='coerce') df_cross['リサイクル率'] = pd.to_numeric(df_cross[COL_RECYCLE], errors='coerce') df_cross['消費支出_sq'] = df_cross['消費支出'] ** 2 # ─── 目的変数・説明変数の設定 ──────────────────────────────────── Y_COL = '1人排出量' Y_LABEL = '1人1日あたりごみ排出量（g/人/日）' X_COL = '消費支出' X_LABEL = '消費支出（円/月，二人以上世帯）—— 所得代理変数' CTRL_COLS = ['光熱水道費', '高齢化率', '気温'] CTRL_LABELS = ['光熱水道費', '高齢化率', '年平均気温'] need_cols = [Y_COL, X_COL, '消費支出_sq'] + CTRL_COLS df_reg = df_cross[need_cols + ['都道府県短', '地域', COL_PREF, 'リサイクル率']]\ .dropna(subset=need_cols).copy() print(f"回帰用 N={len(df_reg)}") print(df_reg[need_cols].describe().round(1)) # ─── Pearson相関 ───────────────────────────────────────────────── print("\n=== Pearson相関（1人排出量 vs 説明変数） ===") for col, label in [(X_COL, '消費支出'), ('光熱水道費', '光熱水道費'), ('高齢化率', '高齢化率'), ('気温', '年平均気温')]: r, p = stats.pearsonr(df_reg[col], df_reg[Y_COL]) sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else 'n.s.' print(f" {label}: r={r:.3f}, p={p:.4f} {sig}") df_rec = df_reg.dropna(subset=['リサイクル率']) print("\n=== Pearson相関（リサイクル率 vs 説明変数） ===") for col, label in [(X_COL, '消費支出'), ('高齢化率', '高齢化率'), ('気温', '年平均気温')]: r, p = stats.pearsonr(df_rec[col], df_rec['リサイクル率']) sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else 'n.s.' print(f" {label}: r={r:.3f}, p={p:.4f} {sig}") # ─── EKC 重回帰分析（2次項含む）───────────────────────────────── print("\n=== EKC 重回帰分析（1人1日排出量）===") exog_cols = [X_COL, '消費支出_sq'] + CTRL_COLS exog_labels = ['消費支出（β₁）', '消費支出²（β₂）'] + CTRL_LABELS # 標準化 df_std = (df_reg[exog_cols + [Y_COL]] - df_reg[exog_cols + [Y_COL]].mean()) / \ df_reg[exog_cols + [Y_COL]].std() X_ekc = sm.add_constant(df_std[exog_cols].astype(float)) res_ekc = sm.OLS(df_std[Y_COL].astype(float), X_ekc).fit() print(res_ekc.summary()) coef_std = res_ekc.params[1:] ci = res_ekc.conf_int().iloc[1:] pvals = res_ekc.pvalues[1:] print("\n=== 標準化偏回帰係数 ===") for name, c, p in zip(exog_labels, coef_std, pvals): sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else 'n.s.' print(f" {name}: β={c:.4f}, p={p:.4f} {sig}") print(f"\nR² = {res_ekc.rsquared:.4f}, Adj.R² = {res_ekc.rsquared_adj:.4f}") print(f"F({res_ekc.df_model:.0f},{res_ekc.df_resid:.0f}) = {res_ekc.fvalue:.3f}, p = {res_ekc.f_pvalue:.6f}") # EKC判定 beta_sq = coef_std.iloc[1] p_sq = pvals.iloc[1] if beta_sq < 0 and p_sq < 0.05: ekc_result = "EKC仮説を支持（逆U字型：β₂ < 0, 有意）" elif beta_sq > 0 and p_sq < 0.05: ekc_result = "EKC仮説を棄却（U字型：β₂ > 0, 有意）" else: ekc_result = "EKC仮説の証拠なし（β₂ 非有意）" print(f"\nEKC検定: {ekc_result}") # 転換点（非標準化） X_raw_reg = sm.add_constant(df_reg[exog_cols].astype(float)) res_raw = sm.OLS(df_reg[Y_COL].astype(float), X_raw_reg).fit() a1 = res_raw.params[X_COL] a2 = res_raw.params['消費支出_sq'] if a2 != 0: turning_point = -a1 / (2 * a2) print(f"転換点（消費支出）: {turning_point:,.0f} 円/月") print(f" → データ範囲: {df_reg[X_COL].min():,.0f}〜{df_reg[X_COL].max():,.0f} 円/月") # ─── 図1：1人1日ごみ排出量ランキング棒グラフ ───────────────────── fig1, ax1 = plt.subplots(figsize=(14, 7)) df_rank = df_reg[['都道府県短', '地域', Y_COL]].sort_values(Y_COL, ascending=True).copy() national_mean = df_reg[Y_COL].mean() colors_bar = [region_colors.get(r, '#999') for r in df_rank['地域']] ax1.barh(df_rank['都道府県短'], df_rank[Y_COL], color=colors_bar, edgecolor='white', linewidth=0.4) ax1.axvline(national_mean, color='#333', linestyle='--', linewidth=1.8, label=f'全国平均 {national_mean:.0f} g/人/日', zorder=5) legend_elements = [Patch(facecolor=region_colors[r], label=r) for r in region_order] legend_elements.append( plt.Line2D([0], [0], color='#333', linestyle='--', linewidth=1.8, label=f'全国平均 {national_mean:.0f} g')) ax1.legend(handles=legend_elements, loc='lower right', fontsize=9, framealpha=0.9) ax1.set_xlabel('1人1日あたりごみ排出量（g/人/日）', fontsize=12) ax1.set_title(f'図1：都道府県別 1人1日あたりごみ排出量ランキング（{LATEST}年度）\n' f'北海道・東北が上位に集中（寒冷地効果・ごみ焼却量の影響）', fontsize=13, fontweight='bold') ax1.tick_params(axis='y', labelsize=8) ax1.tick_params(axis='x', labelsize=10) ax1.set_xlim(0, df_rank[Y_COL].max() * 1.12) ax1.grid(axis='x', alpha=0.3) max_val = df_rank[Y_COL].max() min_val = df_rank[Y_COL].min() max_pref = df_rank.loc[df_rank[Y_COL].idxmax(), '都道府県短'] min_pref = df_rank.loc[df_rank[Y_COL].idxmin(), '都道府県短'] n_rank = len(df_rank) ax1.annotate(f'最大: {max_pref} {max_val:.0f} g', xy=(max_val, n_rank - 1), xytext=(max_val - 90, n_rank - 4), fontsize=8.5, color='#c0392b', fontweight='bold') ax1.annotate(f'最小: {min_pref} {min_val:.0f} g', xy=(min_val, 0), xytext=(min_val + 8, 3), fontsize=8.5, color='#2980b9', fontweight='bold') plt.tight_layout() fig1.savefig(os.path.join(FIG_DIR, '2019_H5_3_fig1.png'), bbox_inches='tight') plt.close(fig1) print("図1 保存完了") # ─── 図2：消費支出 vs 1人1日排出量（EKC散布図）────────────────── fig2, ax2 = plt.subplots(figsize=(9, 7)) for region in region_order: df_r = df_reg[df_reg['地域'] == region] ax2.scatter(df_r[X_COL], df_r[Y_COL], color=region_colors[region], label=region, s=65, alpha=0.88, edgecolors='white', linewidth=0.5, zorder=4) # 都道府県ラベル（上位・下位 + 特徴的な都市） label_prefs = set() for col_val, ascending in [(Y_COL, False), (Y_COL, True), (X_COL, False), (X_COL, True)]: for _, row in df_reg.sort_values(col_val, ascending=ascending).head(4).iterrows(): label_prefs.add(row['都道府県短']) for _, row in df_reg.iterrows(): if row['都道府県短'] in label_prefs: ax2.annotate(row['都道府県短'], (row[X_COL], row[Y_COL]), textcoords='offset points', xytext=(5, 2), fontsize=8, color='#333', fontweight='bold') x_sorted = np.linspace(df_reg[X_COL].min(), df_reg[X_COL].max(), 300) # 線形フィット slope_l, intercept_l, r_l, p_l, _ = stats.linregress(df_reg[X_COL], df_reg[Y_COL]) ax2.plot(x_sorted, intercept_l + slope_l * x_sorted, 'k--', linewidth=1.5, alpha=0.6, label='線形フィット') # 2次フィット（EKC） z2 = np.polyfit(df_reg[X_COL], df_reg[Y_COL], 2) p2 = np.poly1d(z2) ax2.plot(x_sorted, p2(x_sorted), 'r-', linewidth=2.2, alpha=0.85, label='2次フィット（EKC）') # 転換点マーク if z2[0] != 0: tp = -z2[1] / (2 * z2[0]) x_min, x_max = df_reg[X_COL].min(), df_reg[X_COL].max() if x_min <= tp <= x_max: ax2.axvline(tp, color='#c0392b', linestyle=':', linewidth=1.5, alpha=0.8) ax2.text(tp + 2000, df_reg[Y_COL].min() + 5, f'転換点\n{tp/10000:.1f}万円', fontsize=9, color='#c0392b', va='bottom', fontweight='bold') r_lin, p_lin = stats.pearsonr(df_reg[X_COL], df_reg[Y_COL]) sig_lin = '***' if p_lin < 0.001 else '**' if p_lin < 0.01 else '*' if p_lin < 0.05 else 'n.s.' ax2.text(0.04, 0.97, f'Pearson r = {r_lin:.3f} ({sig_lin})\np = {p_lin:.4f}', transform=ax2.transAxes, fontsize=10, va='top', bbox=dict(boxstyle='round', facecolor='white', alpha=0.85, edgecolor='#ccc')) ax2.set_xlabel(X_LABEL, fontsize=11) ax2.set_ylabel(Y_LABEL, fontsize=11) ax2.set_title(f'図2：消費支出（所得代理）とごみ排出量の関係（EKC検討）\n{LATEST}年度 47都道府県', fontsize=12, fontweight='bold') ax2.legend(loc='upper right', fontsize=9, framealpha=0.9) ax2.grid(alpha=0.3) plt.tight_layout() fig2.savefig(os.path.join(FIG_DIR, '2019_H5_3_fig2.png'), bbox_inches='tight') plt.close(fig2) print("図2 保存完了") # ─── 図3：地域別時系列推移（2012〜2023年）──────────────────────── fig3, ax3 = plt.subplots(figsize=(11, 6)) years_range = sorted(df_b[COL_YEAR].unique()) # COVID グレー帯（2020〜2021） ax3.axvspan(2019.5, 2021.5, alpha=0.12, color='gray', label='COVID-19期間（参考）') # 全国平均 national_ts = df_b.groupby(COL_YEAR)[COL_PERCAP].apply( lambda x: pd.to_numeric(x, errors='coerce').mean()) ax3.plot(years_range, national_ts[years_range].values, 'k-', linewidth=3, alpha=0.5, label='全国平均', zorder=3) for region in region_order: df_r = df_b[df_b['地域'] == region].copy() df_r[COL_PERCAP] = pd.to_numeric(df_r[COL_PERCAP], errors='coerce') ts = df_r.groupby(COL_YEAR)[COL_PERCAP].mean().reindex(years_range) ax3.plot(years_range, ts.values, marker='o', markersize=4, color=region_colors[region], label=region, linewidth=2, alpha=0.9) ax3.set_xlabel('年度', fontsize=12) ax3.set_ylabel('1人1日あたりごみ排出量（g/人/日）', fontsize=12) ax3.set_title('図3：地域別 1人1日ごみ排出量の時系列推移（2012〜2023年度）\n' '全地域で一貫した削減トレンド（3R政策の効果）', fontsize=13, fontweight='bold') ax3.legend(loc='upper right', fontsize=9, framealpha=0.9) ax3.set_xticks(years_range) ax3.tick_params(axis='x', rotation=45) ax3.grid(alpha=0.3) y_start = float(national_ts.iloc[0]) y_end = float(national_ts.iloc[-1]) change_pct = (y_end - y_start) / y_start * 100 ax3.annotate(f'全国平均 {change_pct:.1f}%\n（{y_start:.0f}→{y_end:.0f} g）', xy=(years_range[-1], y_end), xytext=(years_range[-4], y_end + 30), arrowprops=dict(arrowstyle='->', color='#333', lw=1.5), fontsize=9, color='#333', fontweight='bold') plt.tight_layout() fig3.savefig(os.path.join(FIG_DIR, '2019_H5_3_fig3.png'), bbox_inches='tight') plt.close(fig3) print("図3 保存完了") # ─── 図4：標準化偏回帰係数の横棒グラフ ───────────────────────────── fig4, ax4 = plt.subplots(figsize=(8, 5)) colors_coef = [] for c, p in zip(coef_std.values, pvals.values): if p < 0.05: colors_coef.append('#e05c5c' if c > 0 else '#4e9af1') else: colors_coef.append('#bbb') y_pos = np.arange(len(exog_labels)) ci_lo = ci.iloc[:, 0].values ci_hi = ci.iloc[:, 1].values xerr_lo = coef_std.values - ci_lo xerr_hi = ci_hi - coef_std.values ax4.barh(y_pos, coef_std.values, color=colors_coef, edgecolor='white', xerr=np.array([xerr_lo, xerr_hi]), error_kw={'ecolor': '#555', 'capsize': 4, 'lw': 1.5}, height=0.6) ax4.axvline(0, color='#333', linewidth=1.2) for i, (c, p) in enumerate(zip(coef_std.values, pvals.values)): sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else '' if sig: x_pos = ci_hi[i] + 0.02 if c >= 0 else ci_lo[i] - 0.02 ha = 'left' if c >= 0 else 'right' ax4.text(x_pos, i, sig, va='center', ha=ha, fontsize=13, color='#333', fontweight='bold') ax4.set_yticks(y_pos) ax4.set_yticklabels(exog_labels, fontsize=11) ax4.set_xlabel('標準化偏回帰係数（β）', fontsize=12) ax4.set_title(f'図4：ごみ排出量の決定要因 — EKC重回帰（標準化係数）{LATEST}年度\n' f'R²={res_ekc.rsquared:.3f}, Adj.R²={res_ekc.rsquared_adj:.3f} ' f'（N=47都道府県）', fontsize=12, fontweight='bold') ax4.grid(axis='x', alpha=0.3) legend_coef = [ Patch(facecolor='#e05c5c', label='正の効果（有意, p<0.05）'), Patch(facecolor='#4e9af1', label='負の効果（有意, p<0.05）'), Patch(facecolor='#bbb', label='非有意'), ] ax4.legend(handles=legend_coef, loc='lower right', fontsize=9) plt.tight_layout() fig4.savefig(os.path.join(FIG_DIR, '2019_H5_3_fig4.png'), bbox_inches='tight') plt.close(fig4) print("図4 保存完了") # ─── サマリー出力 ───────────────────────────────────────────────── print("\n=== 1人1日排出量 上位・下位5都道府県 ===") df_summary = df_reg[['都道府県短', '地域', Y_COL, 'リサイクル率', X_COL, '高齢化率']]\ .sort_values(Y_COL, ascending=False) print("上位5（排出量多）:") print(df_summary.head(5).to_string(index=False)) print("下位5（排出量少）:") print(df_summary.tail(5).to_string(index=False)) print("\nDONE: 2019_H5_3_shorei")