都道府県別ごみ排出量・環境負荷と
経済発展の関係

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

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)

# 地域マップ
region_map = {
    '北海道': '北海道・東北', '青森県': '北海道・東北', '岩手県': '北海道・東北',
    '宮城県': '北海道・東北', '秋田県': '北海道・東北', '山形県': '北海道・東北',
    '福島県': '北海道・東北', '茨城県': '関東', '栃木県': '関東', '群馬県': '関東',
    '埼玉県': '関東', '千葉県': '関東', '東京都': '関東', '神奈川県': '関東',
    '新潟県': '中部', '富山県': '中部', '石川県': '中部', '福井県': '中部',
    '山梨県': '中部', '長野県': '中部', '岐阜県': '中部', '静岡県': '中部', '愛知県': '中部',
    '三重県': '近畿', '滋賀県': '近畿', '京都府': '近畿', '大阪府': '近畿',
    '兵庫県': '近畿', '奈良県': '近畿', '和歌山県': '近畿',
    '鳥取県': '中国・四国', '島根県': '中国・四国', '岡山県': '中国・四国',
    '広島県': '中国・四国', '山口県': '中国・四国', '徳島県': '中国・四国',
    '香川県': '中国・四国', '愛媛県': '中国・四国', '高知県': '中国・四国',
    '福岡県': '九州・沖縄', '佐賀県': '九州・沖縄', '長崎県': '九州・沖縄',
    '熊本県': '九州・沖縄', '大分県': '九州・沖縄', '宮崎県': '九州・沖縄',
    '鹿児島県': '九州・沖縄', '沖縄県': '九州・沖縄'
}
region_colors = {
    '北海道・東北': '#4e9af1', '関東': '#e05c5c', '中部': '#f0a500',
    '近畿': '#5cb85c', '中国・四国': '#9b59b6', '九州・沖縄': '#f39c12'
}

# 全データに地域列を追加
df_b['地域'] = df_b['都道府県'].map(region_map)

# ── 2022年断面データ ────────────────────────────────────────────────────────────
df22 = df_b[df_b['年度'] == 2022].copy()
df22 = df22.dropna(subset=['1人1日当たりの排出量', 'ごみのリサイクル率',
                            '消費支出（二人以上の世帯）', '光熱・水道費（二人以上の世帯）'])

# ── 基本統計の出力 ───────────────────────────────────────────────────────────────
print("=" * 60)
print("2022年断面データ 基本統計（47都道府県）")
print("=" * 60)
cols_show = ['1人1日当たりの排出量', 'ごみのリサイクル率',
             '消費支出（二人以上の世帯）', '光熱・水道費（二人以上の世帯）']
print(df22[cols_show].describe().round(2))

# ── EKC仮説: 消費支出 × ごみ排出量 相関 ────────────────────────────────────────
x_ekc = df22['消費支出（二人以上の世帯）'].values
y_ekc = df22['1人1日当たりの排出量'].values
r_ekc, p_ekc = stats.pearsonr(x_ekc, y_ekc)
print(f"\n消費支出 vs ごみ排出量: r = {r_ekc:.3f}, p = {p_ekc:.4f}")

# ── 重回帰分析 ──────────────────────────────────────────────────────────────────
# 目的変数: 1人1日当たりごみ排出量
# 説明変数: 消費支出, 光熱・水道費（温暖さの代理）, 食料費, 総人口（都市規模）
reg_df = df22[['1人1日当たりの排出量',
               '消費支出（二人以上の世帯）',
               '光熱・水道費（二人以上の世帯）',
               '食料費（二人以上の世帯）',
               '総人口']].dropna().copy()

# 標準化（Zスコア）
reg_z = (reg_df - reg_df.mean()) / reg_df.std()

y_reg = reg_z['1人1日当たりの排出量'].values
X_reg = reg_z[['消費支出（二人以上の世帯）',
               '光熱・水道費（二人以上の世帯）',
               '食料費（二人以上の世帯）',
               '総人口']].values
X_reg_sm = sm.add_constant(X_reg)
ols_model = sm.OLS(y_reg, X_reg_sm).fit()
print("\n" + "=" * 60)
print("重回帰分析 OLS結果（標準化係数）")
print("=" * 60)
print(ols_model.summary())

coef_names = ['消費支出', '光熱・水道費', '食料費', '総人口']
coefs = ols_model.params[1:]
pvals = ols_model.pvalues[1:]
conf_arr = ols_model.conf_int()
conf = pd.DataFrame(conf_arr, columns=['lo', 'hi']).iloc[1:].reset_index(drop=True)

# 出力保存用変数
r2 = ols_model.rsquared
adj_r2 = ols_model.rsquared_adj
f_stat = ols_model.fvalue
f_pval = ols_model.f_pvalue
n_obs = int(ols_model.nobs)

print(f"\nR² = {r2:.3f}, Adj.R² = {adj_r2:.3f}")
print(f"F統計量 = {f_stat:.2f}, p = {f_pval:.4f}")
print(f"サンプル数 = {n_obs}")

# 係数テーブル
coef_table = pd.DataFrame({
    '変数': coef_names,
    '標準化係数': coefs.round(3),
    'p値': pvals.round(4),
    '下限95%CI': conf['lo'].round(3),
    '上限95%CI': conf['hi'].round(3)
})
print(coef_table.to_string(index=False))

# ── リサイクル率: トップ/ボトム確認 ────────────────────────────────────────────
print("\n" + "=" * 60)
print("リサイクル率ランキング（2022年）")
print("=" * 60)
rec_rank = df22[['都道府県', 'ごみのリサイクル率', '地域']].sort_values(
    'ごみのリサイクル率', ascending=False)
print(rec_rank.to_string(index=False))
print(f"\n平均リサイクル率: {df22['ごみのリサイクル率'].mean():.1f}%")
print(f"最高: {df22['ごみのリサイクル率'].max():.1f}% ({df22.loc[df22['ごみのリサイクル率'].idxmax(), '都道府県']})")
print(f"最低: {df22['ごみのリサイクル率'].min():.1f}% ({df22.loc[df22['ごみのリサイクル率'].idxmin(), '都道府県']})")

# ── 時系列データ ────────────────────────────────────────────────────────────────
ts_years = list(range(2012, 2024))
df_ts = df_b[df_b['年度'].isin(ts_years)].copy()
ts_region = df_ts.groupby(['年度', '地域'])['1人1日当たりの排出量'].mean().reset_index()
ts_all = df_ts.groupby('年度')['1人1日当たりの排出量'].mean()

print("\n" + "=" * 60)
print("時系列: 全国平均ごみ排出量（2012〜2023年）")
print("=" * 60)
print(ts_all.round(1))
print(f"\n変化量: {ts_all[2023] - ts_all[2012]:.1f}g/人/日")
print(f"変化率: {(ts_all[2023] / ts_all[2012] - 1) * 100:.1f}%")

# ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━

fig1, ax1 = plt.subplots(figsize=(10, 6))

region_order = ['北海道・東北', '関東', '中部', '近畿', '中国・四国', '九州・沖縄']
for reg in region_order:
    sub = ts_region[ts_region['地域'] == reg].sort_values('年度')
    ax1.plot(sub['年度'], sub['1人1日当たりの排出量'],
             color=region_colors[reg], linewidth=2.2, marker='o', markersize=5,
             label=reg, zorder=3)

# 全国平均
ts_all_plot = ts_all.reset_index()
ax1.plot(ts_all_plot['年度'], ts_all_plot['1人1日当たりの排出量'],
         color='#333333', linewidth=2.5, linestyle='--', marker='s', markersize=4,
         label='全国平均', zorder=4)

# 2020年（コロナ禍）の縦線
ax1.axvline(x=2020, color='#aaaaaa', linestyle=':', linewidth=1.5, alpha=0.8)
ax1.text(2020.1, ax1.get_ylim()[0] + 5, 'COVID-19', fontsize=9, color='#888888')

ax1.set_xlabel('年度', fontsize=12)
ax1.set_ylabel('1人1日当たりごみ排出量（g）', fontsize=12)
ax1.set_title('図1: 地域別ごみ排出量の時系列推移（2012〜2023年）', fontsize=13, fontweight='bold')
ax1.set_xlim(2011.5, 2023.5)
ax1.set_xticks(ts_years)
ax1.set_xticklabels([str(y) for y in ts_years], fontsize=9, rotation=45)
ax1.legend(loc='upper right', fontsize=10, ncol=2)
ax1.grid(axis='y', alpha=0.4, linewidth=0.8)
ax1.grid(axis='x', alpha=0.2, linewidth=0.5)
ax1.spines['top'].set_visible(False)
ax1.spines['right'].set_visible(False)

fig1.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2020_H5_4_fig1.png'), dpi=150, bbox_inches='tight')
plt.close(fig1)
print("\n[FIG1] 保存完了: 2020_H5_4_fig1.png")

# ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━

fig2, ax2 = plt.subplots(figsize=(10, 8))

for _, row in df22.iterrows():
    reg = row['地域']
    col = region_colors.get(reg, '#999999')
    ax2.scatter(row['消費支出（二人以上の世帯）'] / 10000,
                row['1人1日当たりの排出量'],
                color=col, s=80, alpha=0.85, zorder=3, edgecolors='white', linewidth=0.5)
    # 都道府県名ラベル（小さめ）
    ax2.annotate(row['都道府県'].replace('県', '').replace('府', '').replace('都', '').replace('道', ''),
                 (row['消費支出（二人以上の世帯）'] / 10000,
                  row['1人1日当たりの排出量']),
                 fontsize=7, ha='center', va='bottom',
                 xytext=(0, 4), textcoords='offset points', color='#444444')

# 線形回帰フィット
x_fit = df22['消費支出（二人以上の世帯）'].values / 10000
y_fit = df22['1人1日当たりの排出量'].values
slope, intercept, r_val, p_val, se = stats.linregress(x_fit, y_fit)
x_line = np.linspace(x_fit.min(), x_fit.max(), 100)
ax2.plot(x_line, intercept + slope * x_line,
         color='#333333', linewidth=1.8, linestyle='--', alpha=0.7, label=f'回帰直線 (r={r_val:.3f})')

# 凡例（地域）
from matplotlib.lines import Line2D
handles = [Line2D([0], [0], marker='o', color='w', markerfacecolor=region_colors[r],
                  markersize=9, label=r) for r in region_order]
handles.append(Line2D([0], [0], color='#333333', linewidth=1.8, linestyle='--',
                      label=f'回帰直線 (r={r_val:.3f}, p={p_val:.3f})'))
ax2.legend(handles=handles, loc='upper left', fontsize=9, framealpha=0.9)

ax2.set_xlabel('消費支出（万円/月・世帯）', fontsize=12)
ax2.set_ylabel('1人1日当たりごみ排出量（g）', fontsize=12)
ax2.set_title('図2: 消費支出とごみ排出量の関係（2022年・47都道府県）\n'
              '─ EKC仮説（環境クズネッツ曲線）の検証 ─', fontsize=12, fontweight='bold')
ax2.grid(alpha=0.3, linewidth=0.8)
ax2.spines['top'].set_visible(False)
ax2.spines['right'].set_visible(False)

# 注釈: EKC
ax2.text(0.97, 0.05,
         f'Pearson r = {r_val:.3f}\np = {p_val:.3f}\nn = 47都道府県',
         transform=ax2.transAxes, fontsize=10, ha='right', va='bottom',
         bbox=dict(boxstyle='round,pad=0.5', facecolor='#FFF9C4', alpha=0.9))

fig2.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2020_H5_4_fig2.png'), dpi=150, bbox_inches='tight')
plt.close(fig2)
print("[FIG2] 保存完了: 2020_H5_4_fig2.png")

# ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━

rec_sorted = df22.sort_values('ごみのリサイクル率', ascending=True).reset_index(drop=True)
colors_bar = [region_colors.get(r, '#999999') for r in rec_sorted['地域']]
labels_bar = [p.replace('県', '').replace('府', '').replace('都', '').replace('道', '')
              for p in rec_sorted['都道府県']]

fig3, ax3 = plt.subplots(figsize=(9, 12))

bars = ax3.barh(range(len(rec_sorted)), rec_sorted['ごみのリサイクル率'],
                color=colors_bar, edgecolor='white', linewidth=0.5, height=0.8)
ax3.set_yticks(range(len(rec_sorted)))
ax3.set_yticklabels(labels_bar, fontsize=9)

# 値ラベル
for i, (val, bar) in enumerate(zip(rec_sorted['ごみのリサイクル率'], bars)):
    ax3.text(val + 0.2, i, f'{val:.1f}%', va='center', fontsize=8, color='#333333')

# 全国平均線
mean_rec = df22['ごみのリサイクル率'].mean()
ax3.axvline(x=mean_rec, color='#333333', linewidth=1.5, linestyle='--', alpha=0.8)
ax3.text(mean_rec + 0.3, len(rec_sorted) - 0.5, f'全国平均\n{mean_rec:.1f}%',
         fontsize=9, color='#333333', va='top')

# 地域凡例
from matplotlib.patches import Patch
legend_patches = [Patch(color=region_colors[r], label=r) for r in region_order]
ax3.legend(handles=legend_patches, loc='lower right', fontsize=9, framealpha=0.9)

ax3.set_xlabel('ごみリサイクル率（%）', fontsize=12)
ax3.set_title('図3: 都道府県別ごみリサイクル率ランキング（2022年）', fontsize=13, fontweight='bold')
ax3.set_xlim(0, rec_sorted['ごみのリサイクル率'].max() + 4)
ax3.grid(axis='x', alpha=0.4, linewidth=0.8)
ax3.spines['top'].set_visible(False)
ax3.spines['right'].set_visible(False)

fig3.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2020_H5_4_fig3.png'), dpi=150, bbox_inches='tight')
plt.close(fig3)
print("[FIG3] 保存完了: 2020_H5_4_fig3.png")

# ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━

fig4, ax4 = plt.subplots(figsize=(8, 5))

y_pos = range(len(coef_names))
bar_colors = ['#e05c5c' if c > 0 else '#4e9af1' for c in coefs]
ax4.barh(y_pos, coefs, color=bar_colors, edgecolor='white', height=0.55, alpha=0.85)

# 信頼区間
for i, (lo, hi) in enumerate(zip(conf['lo'], conf['hi'])):
    ax4.plot([lo, hi], [i, i], color='#333333', linewidth=2.5, solid_capstyle='round')

# 有意性マーク
for i, (c, p) in enumerate(zip(coefs, pvals)):
    mark = '**' if p < 0.01 else ('*' if p < 0.05 else '')
    sign = '†' if p < 0.10 else ''
    label = mark if mark else sign
    if label:
        ax4.text(c + (0.02 if c >= 0 else -0.02), i, label,
                 ha='left' if c >= 0 else 'right', va='center',
                 fontsize=12, fontweight='bold', color='#C62828')

ax4.axvline(x=0, color='#333333', linewidth=1.2, linestyle='-')
ax4.set_yticks(y_pos)
ax4.set_yticklabels(coef_names, fontsize=11)
ax4.set_xlabel('標準化偏回帰係数（β）', fontsize=12)
ax4.set_title('図4: ごみ排出量の重回帰分析 — 標準化偏回帰係数\n'
              f'(R² = {r2:.3f}, n = {n_obs}都道府県, 2022年)',
              fontsize=12, fontweight='bold')

ax4.text(0.97, 0.05,
         '** p<0.01  * p<0.05  † p<0.10\n横棒: 95%信頼区間',
         transform=ax4.transAxes, fontsize=9, ha='right', va='bottom',
         bbox=dict(boxstyle='round,pad=0.4', facecolor='#F8F9FA', alpha=0.9))

ax4.grid(axis='x', alpha=0.4, linewidth=0.8)
ax4.spines['top'].set_visible(False)
ax4.spines['right'].set_visible(False)

fig4.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2020_H5_4_fig4.png'), dpi=150, bbox_inches='tight')
plt.close(fig4)
print("[FIG4] 保存完了: 2020_H5_4_fig4.png")

print("\n" + "=" * 60)
print("全図の生成が完了しました。")
print(f"  図の保存先: {os.path.abspath(FIG_DIR)}")
print("=" * 60)

# ── HTML用の統計値まとめ出力 ─────────────────────────────────────────────────
print("\n[HTML用統計値]")
print(f"全国平均ごみ排出量（2022年）: {df22['1人1日当たりの排出量'].mean():.1f} g/人/日")
print(f"全国平均ごみ排出量（2012年）: {ts_all[2012]:.1f} g/人/日")
print(f"全国平均ごみ排出量（2023年）: {ts_all[2023]:.1f} g/人/日")
print(f"10年間減少量: {ts_all[2012] - ts_all[2023]:.1f} g")
print(f"減少率: {(1 - ts_all[2023]/ts_all[2012]) * 100:.1f}%")
print(f"全国平均リサイクル率（2022年）: {df22['ごみのリサイクル率'].mean():.1f}%")
print(f"リサイクル率最高: {df22['ごみのリサイクル率'].max():.1f}% - {df22.loc[df22['ごみのリサイクル率'].idxmax(), '都道府県']}")
print(f"リサイクル率最低: {df22['ごみのリサイクル率'].min():.1f}% - {df22.loc[df22['ごみのリサイクル率'].idxmin(), '都道府県']}")
print(f"消費支出 vs ごみ排出量: r = {r_ekc:.3f}, p = {p_ekc:.4f}")
print(f"重回帰 R² = {r2:.3f}, F = {f_stat:.2f}, p = {f_pval:.4f}")
for name, c, p in zip(coef_names, coefs, pvals):
    sig = '**' if p < 0.01 else ('*' if p < 0.05 else ('†' if p < 0.10 else 'ns'))
    print(f"  {name}: β = {c:.3f}, p = {p:.4f} [{sig}]")