ごみの排出量と人口構成比に関する定量分析

import os
import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
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'
YEAR    = 2018            # 原論文が用いた年度（SSDSE-B-2021 収録の 2018 年度データ）
os.makedirs(FIG_DIR, exist_ok=True)

# ── データ読み込み（cp932・2行目に日本語列名 → skiprows=1 で日本語見出しを採用）──
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)

# ── 2018年断面（47都道府県）を抽出 ──────────────────────────────
d = df_b[df_b['年度'] == YEAR].copy()

d['男性割合']     = d['総人口（男）']                 / d['総人口']
d['女性割合']     = d['総人口（女）']                 / d['総人口']
d['外国人割合']   = (d['総人口'] - d['日本人人口'])    / d['総人口']
d['15歳未満割合'] = d['15歳未満人口']                 / d['総人口']
d['1564歳割合']   = d['15～64歳人口']                 / d['総人口']
d['65歳以上割合'] = d['65歳以上人口']                 / d['総人口']

RATIOS = ['男性割合', '外国人割合', '15歳未満割合', '1564歳割合', '65歳以上割合']

fig1, axes = plt.subplots(2, 3, figsize=(12, 7))
titles = {'男性割合': '総人口に占める男性の割合',
          '外国人割合': '総人口に占める外国人人口の割合',
          '15歳未満割合': '15歳未満人口の割合',
          '1564歳割合': '15〜64歳人口の割合',
          '65歳以上割合': '65歳以上人口の割合'}
for ax, col in zip(axes.ravel(), RATIOS):
    ax.hist(d[col], bins=12, color='#4e79c5', edgecolor='white', alpha=0.9)
    ax.axvline(d[col].mean(), color='#c0392b', linestyle='--', linewidth=1.4,
               label=f'平均 {d[col].mean():.3f}')
    ax.set_title(titles[col], fontsize=11)
    ax.set_xlabel('割合', fontsize=9)
    ax.set_ylabel('都道府県数', fontsize=9)
    ax.legend(fontsize=9)
    ax.grid(axis='y', alpha=0.3)
axes.ravel()[-1].axis('off')
axes.ravel()[-1].text(0.05, 0.5,
    '【実再現】\nSSDSE-B-2026.csv の人口列から\n表2の計算式で算出\n（2018年度・47都道府県）',
    fontsize=12, va='center', color='#333')
fig1.suptitle('図1：人口構成比の分布（実再現・SSDSE-B 2018年）', fontsize=14, fontweight='bold')
fig1.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2021_H5_2_fig1.png'))
plt.close(fig1)

x = d['女性割合'].values
y = d['65歳以上割合'].values
slope, intercept, r_val, p_val, _ = stats.linregress(x, y)

fig2, ax2 = plt.subplots(figsize=(9, 7))
ax2.scatter(x, y, s=60, color='#4e79c5', alpha=0.85, zorder=3)
for _, row in d.iterrows():
    pref = row['都道府県'].replace('県', '').replace('府', '').replace('都', '').replace('道', '')
    ax2.annotate(pref, (row['女性割合'], row['65歳以上割合']),
                 fontsize=7.5, xytext=(2, 2), textcoords='offset points', color='#555')
xr = np.linspace(x.min(), x.max(), 100)
ax2.plot(xr, intercept + slope * xr, color='#c0392b', linewidth=1.8,
         label=f'回帰線（r={r_val:.3f}, p<0.001）')
ax2.set_xlabel('総人口に占める女性の割合', fontsize=12)
ax2.set_ylabel('総人口に占める65歳以上人口の割合', fontsize=12)
ax2.set_title('図2：女性の割合と高齢者割合の関係（実再現・2018年、N=47）\n'
              '原論文の報告値 r=0.643 を再現 → 偏相関分析が必要な理由',
              fontsize=12, fontweight='bold')
ax2.legend(fontsize=11, loc='upper left')
ax2.grid(alpha=0.25)
fig2.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2021_H5_2_fig2.png'))

POP_VIEWS  = ['男性の割合', '外国人人口割合', '15歳未満人口割合', '15〜64歳人口割合', '65歳以上人口割合']
WASTE_KIND = ['ごみ合計', '生活系ごみ', '事業系ごみ']
REPORTED_R = np.array([
    [-0.218082407,  0.138296204, -0.397686082],   # 男性の割合
    [-0.356373954, -0.028430004, -0.383354428],   # 外国人人口割合
    [-0.311385645, -0.393769041,  0.051500657],   # 15歳未満人口割合
    [-0.384795539, -0.131424138, -0.308414087],   # 15〜64歳人口割合
    [ 0.443669497,  0.250275895,  0.252165197],   # 65歳以上人口割合
])

# 5-2 偏相関係数：単純相関 → 交絡変数を除いた偏相関（女性の割合ベースで報告）
PARTIAL = [
    # (ラベル, 単純相関, 除いた変数, 偏相関)
    ('生活系ごみ × 女性の割合',  -0.148, '65歳以上人口割合', -0.416),
    ('生活系ごみ × 65歳以上割合', 0.250, '女性の人口割合',    0.456),
    ('事業系ごみ × 女性の割合',   0.413, '65歳以上人口割合',  0.339),
    ('事業系ごみ × 65歳以上割合', 0.252, '女性の人口割合',   -0.019),
]

fig4, ax4 = plt.subplots(figsize=(10, 6))
ypos = np.arange(len(labels))
h = 0.36
ax4.barh(ypos + h/2, simple_r,  height=h, color='#b0bec5', label='単純相関係数')
ax4.barh(ypos - h/2, partial_r, height=h, color='#e05c5c', label='偏相関係数')
for i, (s, p) in enumerate(zip(simple_r, partial_r)):
    ax4.text(s + (0.01 if s >= 0 else -0.01), i + h/2, f'{s:.3f}',
             va='center', ha='left' if s >= 0 else 'right', fontsize=10, color='#455a64')
    ax4.text(p + (0.01 if p >= 0 else -0.01), i - h/2, f'{p:.3f}',
             va='center', ha='left' if p >= 0 else 'right', fontsize=10, color='#c0392b')
ax4.axvline(0, color='black', linewidth=0.8)
for thr in (0.3, -0.3):
    ax4.axvline(thr, color='#888', linestyle=':', linewidth=1)
ax4.set_yticks(ypos)
ax4.set_yticklabels(labels, fontsize=11)
ax4.set_xlabel('相関係数（原論文の報告値）', fontsize=12)
ax4.set_title('図4：単純相関 → 偏相関で見え方が変わる【原論文の報告値】\n'
              '点線は |r|=0.3（相関の有無の目安）', fontsize=12, fontweight='bold')
ax4.legend(fontsize=11, loc='lower right')
ax4.grid(axis='x', alpha=0.25)
fig4.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2021_H5_2_fig4.png'))