都道府県別ごみ排出量とリサイクル率の決定要因
— EKC仮説の検証 —

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)

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

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 保存完了")

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 保存完了")

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 保存完了")

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

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} 円/月")

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