ごみの削減とリサイクルを推進する要因
―環境ボランティアは、ごみ削減の効果を持つのか―

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)

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')

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

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}")

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")

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")

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")

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）による分析完了。")