都道府県別CO2排出量の決定要因と
地球温暖化対策の分析

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)

# 2022年度断面データ
df = df_b[df_b['年度'] == 2022].copy().reset_index(drop=True)
print(f"[データ確認] 2022年度 都道府県数: {len(df)}")

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

df['地域区分'] = df['都道府県'].map(region_map)
df['地域色']   = df['地域区分'].map(region_colors)

df['y_gomi_pp']    = df['1人1日当たりの排出量']

# 説明変数
df['x_konetsu']    = df['光熱・水道費（二人以上の世帯）']      # エネルギー消費代理
df['x_kotsu']      = df['交通・通信費（二人以上の世帯）']       # 自動車利用代理
df['x_recycle']    = df['ごみのリサイクル率']                   # 廃棄物対策
df['x_temp']       = df['年平均気温']                           # 気候（冷暖房需要）
df['x_rain']       = df['降水日数（年間）']                     # 降水日数

# 分析データセット
ana_cols = ['都道府県', '地域区分', '地域色',
            'y_gomi_pp', 'x_konetsu', 'x_kotsu', 'x_recycle', 'x_temp', 'x_rain']
dfa = df[ana_cols].dropna().copy()
print(f"[分析対象] {len(dfa)}都道府県（欠損除外後）")

# 基本統計量の表示
print("\n=== 基本統計量 ===")
desc_cols = {
    '1人1日ごみ排出量 (g/人・日)': 'y_gomi_pp',
    '光熱・水道費 (円/月)': 'x_konetsu',
    '交通・通信費 (円/月)': 'x_kotsu',
    'ごみリサイクル率 (%)': 'x_recycle',
    '年平均気温 (℃)': 'x_temp',
    '降水日数 (日/年)': 'x_rain',
}
for label, col in desc_cols.items():
    s = dfa[col]
    print(f"  {label}: 平均={s.mean():.1f}, 標準偏差={s.std():.1f}, "
          f"最小={s.min():.1f}, 最大={s.max():.1f}")

print("\n=== 地域別 1人1日ごみ排出量 平均 ===")
region_stats = dfa.groupby('地域区分')['y_gomi_pp'].agg(['mean', 'std', 'count'])
region_stats.columns = ['平均', '標準偏差', '都道府県数']
print(region_stats.sort_values('平均', ascending=False).to_string())

# 上位・下位都道府県
print("\n=== 上位5都道府県（ごみ排出量多） ===")
top5 = dfa_sorted.head(5)[['都道府県', 'y_gomi_pp', '地域区分']]
top5.columns = ['都道府県', '1人1日排出量(g)', '地域区分']
print(top5.to_string(index=False))

print("\n=== 下位5都道府県（ごみ排出量少） ===")
bot5 = dfa_sorted.tail(5)[['都道府県', 'y_gomi_pp', '地域区分']]
bot5.columns = ['都道府県', '1人1日排出量(g)', '地域区分']
print(bot5.to_string(index=False))

# 特化係数（LQ: Location Quotient）の計算例
# LQ = (都道府県iの光熱水道費シェア) / (全国の光熱水道費シェア)
# ここでは光熱水道費の全国消費支出に占める割合の特化係数を計算
df_lq = df.copy()
national_share = df_lq['光熱・水道費（二人以上の世帯）'].mean() / df_lq['消費支出（二人以上の世帯）'].mean()
df_lq['LQ_konetsu'] = (df_lq['光熱・水道費（二人以上の世帯）'] / df_lq['消費支出（二人以上の世帯）']) / national_share
print("\n=== 光熱・水道費特化係数（LQ）上位10都道府県 ===")
lq_top = df_lq[['都道府県', 'LQ_konetsu']].sort_values('LQ_konetsu', ascending=False).head(10)
print(lq_top.to_string(index=False))

print("\n[完了] 全4枚の図を保存しました。")
print(f"  保存先: {FIG_DIR}")

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

for region, color in region_colors.items():
    sub = dfa[dfa['地域区分'] == region]
    ax1.scatter(sub['x_konetsu'], sub['y_gomi_pp'],
                color=color, label=region, s=60, zorder=3, alpha=0.9)
    for _, row in sub.iterrows():
        pref = row['都道府県'].replace('県', '').replace('都', '').replace('府', '').replace('道', '')
        ax1.annotate(pref,
                     xy=(row['x_konetsu'], row['y_gomi_pp']),
                     xytext=(3, 3), textcoords='offset points',
                     fontsize=6.5, color='#444', zorder=4)

# 回帰直線
x_line = np.linspace(dfa['x_konetsu'].min(), dfa['x_konetsu'].max(), 100)
r_val, p_val = stats.pearsonr(dfa['x_konetsu'], dfa['y_gomi_pp'])
slope, intercept, *_ = stats.linregress(dfa['x_konetsu'], dfa['y_gomi_pp'])
y_line = intercept + slope * x_line
ax1.plot(x_line, y_line, color='#333', linewidth=1.5, linestyle='--',
         label=f'回帰直線 (r={r_val:.2f})')

ax1.set_xlabel('光熱・水道費（円/月）\n（エネルギー消費の代理変数）', fontsize=11)
ax1.set_ylabel('1人1日当たりごみ排出量（g/人・日）\n（CO2排出の代理変数）', fontsize=11)
ax1.set_title('光熱・水道費と1人1日ごみ排出量の散布図\n（2022年度 都道府県別）', fontsize=12, fontweight='bold')
ax1.legend(loc='upper left', fontsize=9, framealpha=0.8)
ax1.grid(True, alpha=0.3)

fig1.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2020_H5_5_fig1.png'), dpi=150, bbox_inches='tight')
plt.close(fig1)
print("\n[図1] 散布図 → 保存完了")

dfa_sorted = dfa.sort_values('y_gomi_pp', ascending=False).reset_index(drop=True)

fig2, ax2 = plt.subplots(figsize=(12, 7))
bars = ax2.bar(range(len(dfa_sorted)),
               dfa_sorted['y_gomi_pp'],
               color=dfa_sorted['地域色'].values,
               edgecolor='white', linewidth=0.5, zorder=3)

# 都道府県ラベル（X軸）
pref_labels = [p.replace('県', '').replace('都', '').replace('府', '').replace('道', '')
               for p in dfa_sorted['都道府県']]
ax2.set_xticks(range(len(dfa_sorted)))
ax2.set_xticklabels(pref_labels, rotation=90, fontsize=8)

# 凡例（地域区分）
from matplotlib.patches import Patch
legend_patches = [Patch(color=c, label=r) for r, c in region_colors.items()]
ax2.legend(handles=legend_patches, loc='upper right', fontsize=9, framealpha=0.8)

# 全国平均線
mean_val = dfa['y_gomi_pp'].mean()
ax2.axhline(mean_val, color='#333', linestyle='--', linewidth=1.5, zorder=4)
ax2.text(len(dfa_sorted) - 0.5, mean_val + 5,
         f'全国平均\n{mean_val:.0f}g', ha='right', fontsize=8.5, color='#333')

ax2.set_xlabel('都道府県', fontsize=11)
ax2.set_ylabel('1人1日当たりごみ排出量（g/人・日）', fontsize=11)
ax2.set_title('都道府県別 1人1日当たりごみ排出量ランキング\n（2022年度・CO2排出の代理変数）', fontsize=12, fontweight='bold')
ax2.set_xlim(-0.7, len(dfa_sorted) - 0.3)
ax2.set_ylim(700, max(dfa_sorted['y_gomi_pp']) * 1.12)
ax2.grid(True, axis='y', alpha=0.3, zorder=0)

fig2.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2020_H5_5_fig2.png'), dpi=150, bbox_inches='tight')
plt.close(fig2)
print("[図2] ランキング棒グラフ → 保存完了")

heat_vars = ['y_gomi_pp', 'x_konetsu', 'x_kotsu', 'x_recycle', 'x_temp', 'x_rain']
heat_labels = ['1人1日\nごみ排出量', '光熱・\n水道費', '交通・\n通信費', 'リサイ\nクル率', '年平均\n気温', '降水\n日数']

corr_matrix = dfa[heat_vars].corr()

fig3, ax3 = plt.subplots(figsize=(7, 6))
n = len(heat_vars)
im = ax3.imshow(corr_matrix.values, cmap='RdBu_r', vmin=-1, vmax=1)
plt.colorbar(im, ax=ax3, shrink=0.85, label='Pearson r')

ax3.set_xticks(range(n))
ax3.set_yticks(range(n))
ax3.set_xticklabels(heat_labels, fontsize=9)
ax3.set_yticklabels(heat_labels, fontsize=9)

for i in range(n):
    for j in range(n):
        val = corr_matrix.values[i, j]
        r_ij, p_ij = stats.pearsonr(dfa[heat_vars[i]], dfa[heat_vars[j]])
        sig_mark = '**' if p_ij < 0.01 else ('*' if p_ij < 0.05 else '')
        txt_color = 'white' if abs(val) > 0.65 else 'black'
        ax3.text(j, i, f'{val:.2f}{sig_mark}', ha='center', va='center',
                 fontsize=8, color=txt_color, fontweight='bold' if sig_mark else 'normal')

ax3.set_title('環境・生活指標の相関ヒートマップ\n（Pearson r、*p<0.05, **p<0.01）', fontsize=11, fontweight='bold')
fig3.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2020_H5_5_fig3.png'), dpi=150, bbox_inches='tight')
plt.close(fig3)
print("[図3] 相関ヒートマップ → 保存完了")

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
import warnings
warnings.filterwarnings('ignore')

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)

print("\n=== 相関分析（目的変数: 1人1日ごみ排出量） ===")
corr_vars = ['x_konetsu', 'x_kotsu', 'x_recycle', 'x_temp', 'x_rain']
corr_labels = {
    'x_konetsu': '光熱・水道費',
    'x_kotsu':   '交通・通信費',
    'x_recycle': 'リサイクル率',
    'x_temp':    '年平均気温',
    'x_rain':    '降水日数',
}
corr_results = {}
for col in corr_vars:
    r, p = stats.pearsonr(dfa['y_gomi_pp'], dfa[col])
    corr_results[col] = (r, p)
    sig = '**' if p < 0.01 else ('*' if p < 0.05 else 'n.s.')
    print(f"  {corr_labels[col]:12s}: r={r:.3f}, p={p:.4f} {sig}")

X_cols = corr_vars
X = dfa[X_cols].values
y = dfa['y_gomi_pp'].values

# 標準化
X_mean = X.mean(axis=0)
X_std  = X.std(axis=0, ddof=1)
y_mean = y.mean()
y_std  = y.std(ddof=1)

X_std_arr = (X - X_mean) / X_std
y_std_arr = (y - y_mean) / y_std

X_with_const = sm.add_constant(X_std_arr)
model = sm.OLS(y_std_arr, X_with_const).fit()

print("\n=== 重回帰分析結果 (標準化係数) ===")
print(f"  R² = {model.rsquared:.4f}, 調整済みR² = {model.rsquared_adj:.4f}")
print(f"  F値 = {model.fvalue:.3f}, p値 = {model.f_pvalue:.6f}")
print()
coef_labels = ['定数項'] + [corr_labels[c] for c in X_cols]
for i, (label, coef, pval) in enumerate(zip(coef_labels,
                                             model.params,
                                             model.pvalues)):
    sig = '**' if pval < 0.01 else ('*' if pval < 0.05 else 'n.s.')
    print(f"  {label:14s}: β={coef:+.4f}, p={pval:.4f} {sig}")

coef_names = [corr_labels[c] for c in X_cols]
coef_vals   = model.params[1:]   # 定数項除く
pvals_vals  = model.pvalues[1:]
conf_int    = np.array(model.conf_int())[1:]  # 信頼区間

# 係数の絶対値で降順ソート
order = np.argsort(np.abs(coef_vals))[::-1]
coef_vals_s  = coef_vals[order]
pvals_vals_s = pvals_vals[order]
names_s      = [coef_names[i] for i in order]
ci_s         = conf_int[order]

colors_bar = ['#e05c5c' if p < 0.05 else '#aaa' for p in pvals_vals_s]

fig4, ax4 = plt.subplots(figsize=(7, 5))
y_pos = range(len(names_s))
ax4.barh(y_pos, coef_vals_s, color=colors_bar, edgecolor='white', height=0.6, zorder=3)

# 信頼区間（95%CI）
for i, (ci_l, ci_u) in enumerate(ci_s):
    ax4.plot([ci_l, ci_u], [i, i], color='#333', linewidth=1.5, zorder=4)

ax4.set_yticks(list(y_pos))
ax4.set_yticklabels(names_s, fontsize=11)
ax4.axvline(0, color='black', linewidth=0.8)
ax4.set_xlabel('標準化偏回帰係数 β', fontsize=11)
ax4.set_title('重回帰分析：標準化偏回帰係数\n（赤：p<0.05有意、灰色：非有意）', fontsize=11, fontweight='bold')
ax4.grid(True, axis='x', alpha=0.3, zorder=0)

# p値ラベル
for i, (coef, pval) in enumerate(zip(coef_vals_s, pvals_vals_s)):
    sig = '**' if pval < 0.01 else ('*' if pval < 0.05 else 'n.s.')
    xpos = coef + (0.02 if coef >= 0 else -0.02)
    ha = 'left' if coef >= 0 else 'right'
    ax4.text(xpos, i, sig, va='center', ha=ha, fontsize=9, color='#333')

# 凡例
from matplotlib.patches import Patch
ax4.legend(handles=[
    Patch(color='#e05c5c', label='有意 (p<0.05)'),
    Patch(color='#aaa',    label='非有意')
], loc='lower right', fontsize=9)

fig4.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2020_H5_5_fig4.png'), dpi=150, bbox_inches='tight')
plt.close(fig4)
print("[図4] 標準化偏回帰係数プロット → 保存完了")