リサイクル活動に対する地域・政策要因の研究

print("\n" + "=" * 70)
print("■ Step3. Ward法 階層的クラスタリング（47都道府県）")
print("=" * 70)

Z = linkage(X_scaled, method='ward')

# k=4 でクラスター割り当て
N_CLUSTERS = 4
labels_ward = fcluster(Z, t=N_CLUSTERS, criterion='maxclust')
df_ana['cluster'] = labels_ward

print(f"\n【クラスター別の都道府県（k={N_CLUSTERS}）】")
for cl in sorted(df_ana['cluster'].unique()):
    prefs_cl = df_ana[df_ana['cluster'] == cl]['都道府県'].tolist()
    print(f"  クラスター{cl} (n={len(prefs_cl)}): {', '.join(prefs_cl)}")

print(f"\n【クラスター別 リサイクル率の平均】")
cl_recycling = df_ana.groupby('cluster')[TARGET].agg(['mean', 'std', 'min', 'max'])
for cl, row in cl_recycling.iterrows():
    print(f"  クラスター{cl}: 平均={row['mean']:.2f}%  SD={row['std']:.2f}%  "
          f"[{row['min']:.1f}%~{row['max']:.1f}%]")

ax2a = axes2[0]

# クラスター色でサンプルを散布
for cl in sorted(df_ana['cluster'].unique()):
    mask = df_ana['cluster'] == cl
    ax2a.scatter(X_pca[mask, 0], X_pca[mask, 1],
                 color=C_CLUSTER[cl], s=60, alpha=0.85,
                 edgecolors='white', linewidth=0.5, zorder=3,
                 label=f'クラスター{cl}')

# 都道府県名ラベル
for i, pref in enumerate(prefs):
    ax2a.annotate(pref, (X_pca[i, 0], X_pca[i, 1]),
                  fontsize=5.5, ha='center', va='bottom',
                  xytext=(0, 3), textcoords='offset points',
                  color='#333333', zorder=4)

# ローディングベクトル（矢印）
scale = 2.5
for j, feat in enumerate(FEAT_SHORT):
    ax2a.annotate('', xy=(loadings[0, j] * scale, loadings[1, j] * scale),
                  xytext=(0, 0),
                  arrowprops=dict(arrowstyle='->', color='#C62828',
                                  lw=1.5, shrinkA=0, shrinkB=0))
    ax2a.text(loadings[0, j] * scale * 1.12, loadings[1, j] * scale * 1.12,
              feat, fontsize=8, color='#C62828', fontweight='bold',
              ha='center', va='center')

ax2a.axhline(0, color='gray', linewidth=0.6, linestyle='--', alpha=0.5)
ax2a.axvline(0, color='gray', linewidth=0.6, linestyle='--', alpha=0.5)
ax2a.set_xlabel(f'PC1（寄与率: {explained_var[0]*100:.1f}%）', fontsize=10)
ax2a.set_ylabel(f'PC2（寄与率: {explained_var[1]*100:.1f}%）', fontsize=10)
ax2a.set_title(f'バイプロット（PC1 vs PC2）\n累積寄与率: {cumulative_var[1]*100:.1f}%',
               fontsize=10, fontweight='bold')
ax2a.legend(fontsize=8, loc='upper right')
ax2a.grid(True, alpha=0.25)

print("\n" + "=" * 70)
print("✓ 全図の生成完了（4枚）")
print("  fig1_corr.png       : 相関ヒートマップ（リサイクル率と説明変数）")
print("  fig2_pca.png        : PCA バイプロット（PC1-PC2・都道府県ラベル）")
print("  fig3_dendrogram.png : Ward法デンドログラム（47都道府県）")
print("  fig4_cluster.png    : クラスター別プロファイル比較（棒グラフ）")
print("=" * 70)

print("\n【最終サマリー】")
print(f"\n  ■ リサイクル率との相関（上位）")
for var, r in recycling_corr.items():
    if abs(r) >= 0.3:
        print(f"    {var:<14}: r = {r:+.4f}")

print(f"\n  ■ PCA 結果")
print(f"    PC1 寄与率: {explained_var[0]*100:.1f}%（主にごみ量・消費支出系）")
print(f"    PC2 寄与率: {explained_var[1]*100:.1f}%")
print(f"    PC1+PC2 累積: {cumulative_var[1]*100:.1f}%")

print(f"\n  ■ Ward法クラスタリング（k={N_CLUSTERS}）")
for cl in sorted(df_ana['cluster'].unique()):
    n_cl  = (df_ana['cluster'] == cl).sum()
    m_rec = df_ana[df_ana['cluster'] == cl][TARGET].mean()
    prefs_cl = df_ana[df_ana['cluster'] == cl]['都道府県'].tolist()
    print(f"    クラスター{cl}: n={n_cl}, リサイクル率平均={m_rec:.1f}%")
    print(f"      → {', '.join(prefs_cl)}")

print("図1: 相関ヒートマップを作成中...")

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

# 相関行列（リサイクル率を含む全変数）
all_vars_short = [TARGET] + FEAT_SHORT
corr_mat = df_ana[[TARGET] + FEAT_SHORT].corr()

im = ax1.imshow(corr_mat.values, cmap='RdBu_r', vmin=-1, vmax=1, aspect='auto')
plt.colorbar(im, ax=ax1, fraction=0.046, pad=0.04, label='相関係数')

tick_labels = ['リサイクル率', 'ごみ総排出量', '1人1日排出量',
               '消費支出', '食料費', '光熱・水道費', '教養娯楽費', '年平均気温', '高齢化率']
ax1.set_xticks(range(len(tick_labels)))
ax1.set_yticks(range(len(tick_labels)))
ax1.set_xticklabels(tick_labels, rotation=40, ha='right', fontsize=9)
ax1.set_yticklabels(tick_labels, fontsize=9)

# 数値を各セルに表示
for i in range(len(tick_labels)):
    for j in range(len(tick_labels)):
        val = corr_mat.values[i, j]
        color = 'white' if abs(val) > 0.6 else 'black'
        ax1.text(j, i, f'{val:.2f}', ha='center', va='center',
                 fontsize=7.5, color=color, fontweight='bold')

# リサイクル率の行・列をハイライト
for pos in range(len(tick_labels)):
    ax1.add_patch(plt.Rectangle((pos - 0.5, -0.5), 1, 1,
                                 fill=False, edgecolor='#FFA000', linewidth=0, alpha=0))
ax1.axhline(y=0.5, color='#FFA000', linewidth=2)
ax1.axhline(y=-0.5, color='#FFA000', linewidth=2)
ax1.axvline(x=0.5, color='#FFA000', linewidth=2)
ax1.axvline(x=-0.5, color='#FFA000', linewidth=2)

ax1.set_title('相関係数行列ヒートマップ\n（2022年度 47都道府県・SSDSE-B）',
              fontsize=12, fontweight='bold', pad=14)
fig1.text(0.5, 0.01,
          '出典: 統計数理研究所 SSDSE-B-2026（実データ）',
          ha='center', fontsize=8, color='#666')

plt.tight_layout(rect=[0, 0.03, 1, 1])
save_fig('fig1_corr')

ax2b = axes2[1]
n_pcs = min(8, len(explained_var))
pc_labels = [f'PC{i+1}' for i in range(n_pcs)]
bars = ax2b.bar(pc_labels, explained_var[:n_pcs] * 100,
                color='#1565C0', alpha=0.8, edgecolor='white', label='寄与率')
ax2b_twin = ax2b.twinx()
ax2b_twin.plot(pc_labels, cumulative_var[:n_pcs] * 100,
               'o-', color='#E65100', linewidth=2.2, markersize=7,
               markerfacecolor='white', markeredgewidth=2, label='累積寄与率')
ax2b_twin.axhline(80, color='#E65100', linestyle='--', linewidth=1, alpha=0.5)
ax2b_twin.set_ylabel('累積寄与率（%）', fontsize=10, color='#E65100')
ax2b_twin.set_ylim(0, 105)
ax2b_twin.tick_params(axis='y', labelcolor='#E65100')

for i, (bar, ev) in enumerate(zip(bars, explained_var[:n_pcs])):
    ax2b.text(bar.get_x() + bar.get_width() / 2, bar.get_height() + 0.3,
              f'{ev*100:.1f}%', ha='center', va='bottom', fontsize=8, color='#1565C0')
for i, cv in enumerate(cumulative_var[:n_pcs]):
    ax2b_twin.annotate(f'{cv*100:.1f}%', (i, cv * 100),
                       textcoords='offset points', xytext=(6, 4), fontsize=7.5, color='#E65100')

ax2b.set_xlabel('主成分', fontsize=10)
ax2b.set_ylabel('寄与率（%）', fontsize=10, color='#1565C0')
ax2b.set_title('主成分の寄与率・累積寄与率', fontsize=10, fontweight='bold')
ax2b.grid(axis='y', alpha=0.3)
ax2b.tick_params(axis='y', labelcolor='#1565C0')

# 因子負荷量の凡例表として追加
load_text = '【因子負荷量 PC1 / PC2】\n'
for j, feat in enumerate(FEAT_SHORT):
    load_text += f'{feat}: {loadings[0,j]:+.2f} / {loadings[1,j]:+.2f}\n'
ax2b.text(0.01, 0.01, load_text.strip(), transform=ax2b.transAxes,
          fontsize=7, va='bottom', ha='left', fontfamily='monospace',
          bbox=dict(boxstyle='round', facecolor='#EFF3FF', alpha=0.9, edgecolor='#1565C0'))

plt.tight_layout()
save_fig('fig2_pca')

print("図3: Ward法デンドログラムを作成中...")

fig3, ax3 = plt.subplots(figsize=(16, 7))

# クラスター色対応のデンドログラム
cl_colors_ward = {}
for i, pref in enumerate(prefs):
    cl = df_ana.loc[i, 'cluster']
    cl_colors_ward[pref] = C_CLUSTER[cl]

# デンドログラムを描画
ddata = dendrogram(
    Z,
    labels=prefs,
    ax=ax3,
    leaf_rotation=90,
    leaf_font_size=8,
    color_threshold=Z[-N_CLUSTERS + 1, 2],   # k=4 のカット高さ
    above_threshold_color='#9E9E9E',
)

# カット線
cut_height = (Z[-N_CLUSTERS, 2] + Z[-N_CLUSTERS + 1, 2]) / 2
ax3.axhline(y=cut_height, color='red', linestyle='--', linewidth=2,
            label=f'カット高さ（k={N_CLUSTERS}クラスター）')

ax3.set_title(f'Ward法 階層的クラスタリング デンドログラム\n'
              f'47都道府県 ごみ・消費・人口・気温変数（2022年度・SSDSE-B）',
              fontsize=12, fontweight='bold')
ax3.set_xlabel('都道府県', fontsize=10)
ax3.set_ylabel('Ward距離（クラスター内分散の増分）', fontsize=10)
ax3.legend(fontsize=9, loc='upper right')
ax3.grid(axis='y', alpha=0.25)

# クラスター凡例
patch_list = [mpatches.Patch(color=C_CLUSTER[cl], label=f'クラスター{cl}')
              for cl in sorted(df_ana['cluster'].unique())]
ax3.legend(handles=patch_list + [
    mpatches.Patch(color='red', label=f'カット高さ（k={N_CLUSTERS}）')
], fontsize=9, loc='upper left')

fig3.text(0.5, 0.01,
          '出典: 統計数理研究所 SSDSE-B-2026（実データ）| Ward法: scipy.cluster.hierarchy',
          ha='center', fontsize=8, color='#666')

plt.tight_layout(rect=[0, 0.03, 1, 1])
save_fig('fig3_dendrogram')

import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from scipy.cluster.hierarchy import dendrogram, linkage, fcluster
import warnings
warnings.filterwarnings('ignore')

plt.rcParams['font.family'] = 'Hiragino Sans'
plt.rcParams['axes.unicode_minus'] = False
plt.rcParams['figure.dpi'] = 150

import os
FIGURE_DIR = 'html/figures'
os.makedirs(FIGURE_DIR, exist_ok=True)

def save_fig(name):
    path = os.path.join(FIGURE_DIR, f'2022_U3_{name}.png')
    plt.savefig(path, bbox_inches='tight', dpi=150)
    plt.close()
    print(f"  → {path} 保存完了")

# クラスターカラー（k=4）
C_CLUSTER = {1: '#E53935', 2: '#1E88E5', 3: '#43A047', 4: '#FB8C00'}

DATA_B = 'data/raw/SSDSE-B-2026.csv'
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_2022 = df_b[df_b['年度'] == 2022].copy()

print("=" * 70)
print("■ SSDSE-B-2026 実データ（2022年度、N=47都道府県）")
print("  出典: 統計数理研究所 SSDSE（社会・人口統計体系）")
print("=" * 70)
print(f"  行数: {len(df_2022)}, 列数: {df_2022.shape[1]}")

# 分析変数の定義
TARGET = 'ごみのリサイクル率'
FEATURES = {
    'ごみ総排出量（総量）':         'ごみ総排出量',
    '1人1日当たりの排出量':          '1人1日排出量',
    '消費支出（二人以上の世帯）':    '消費支出',
    '食料費（二人以上の世帯）':      '食料費',
    '光熱・水道費（二人以上の世帯）': '光熱・水道費',
    '教養娯楽費（二人以上の世帯）':  '教養娯楽費',
    '年平均気温':                     '年平均気温',
}

# 高齢化率を計算（65歳以上人口 / 総人口）
df_2022 = df_2022.copy()
df_2022['高齢化率'] = df_2022['65歳以上人口'] / df_2022['総人口'] * 100
FEATURES['高齢化率'] = '高齢化率'

# 分析用データフレームを作成
feat_cols = list(FEATURES.keys())
all_cols  = [TARGET] + feat_cols
df_ana = df_2022[['都道府県'] + all_cols].copy().reset_index(drop=True)
df_ana.columns = ['都道府県', TARGET] + [FEATURES[c] for c in feat_cols]

FEAT_SHORT = list(FEATURES.values())  # 短縮名リスト

print(f"\n【分析変数一覧（N=47都道府県、2022年度）】")
print(f"  目的変数: {TARGET}")
print(f"  説明変数: {', '.join(FEAT_SHORT)}")
print(f"\n【基本統計量】")
print(df_ana[FEAT_SHORT + [TARGET]].describe().round(2))

# 標準化
scaler = StandardScaler()
X_raw    = df_ana[FEAT_SHORT].values
X_scaled = scaler.fit_transform(X_raw)

prefs = df_ana['都道府県'].tolist()

print("\n" + "=" * 70)
print("■ Step1. 相関分析（リサイクル率 vs 説明変数）")
print("=" * 70)

corr_df = df_ana[[TARGET] + FEAT_SHORT].corr()
print(f"\n【リサイクル率との相関係数】")
recycling_corr = corr_df[TARGET].drop(TARGET).sort_values(ascending=False)
for var, r in recycling_corr.items():
    direction = '+' if r > 0 else '-'
    strength  = '強' if abs(r) >= 0.5 else '中' if abs(r) >= 0.3 else '弱'
    print(f"  {var:<14}: r = {r:+.4f}  ({direction} {strength})")

print("\n" + "=" * 70)
print("■ Step2. 主成分分析（PCA）")
print("=" * 70)

pca = PCA()
pca.fit(X_scaled)
X_pca = pca.transform(X_scaled)

explained_var  = pca.explained_variance_ratio_
cumulative_var = np.cumsum(explained_var)
loadings       = pca.components_  # shape: (n_components, n_features)

print(f"\n【主成分の寄与率】")
print(f"  {'PC':<6} {'固有値':>8} {'寄与率':>9} {'累積寄与率':>11}")
print("  " + "-" * 38)
eigenvalues = pca.explained_variance_
for i in range(min(6, len(explained_var))):
    print(f"  PC{i+1:<4} {eigenvalues[i]:>8.4f} {explained_var[i]*100:>8.2f}% {cumulative_var[i]*100:>10.2f}%")

print(f"\n【PC1・PC2 因子負荷量（loadings）】")
print(f"  {'変数':<14} {'PC1 負荷量':>11} {'PC2 負荷量':>11}")
print("  " + "-" * 38)
for j, feat in enumerate(FEAT_SHORT):
    print(f"  {feat:<14} {loadings[0, j]:>+11.4f} {loadings[1, j]:>+11.4f}")

print(f"\n  PC1 寄与率: {explained_var[0]*100:.1f}%  PC2 寄与率: {explained_var[1]*100:.1f}%")
print(f"  PC1+PC2 累積: {cumulative_var[1]*100:.1f}%")

print("\n" + "=" * 70)
print("■ Step4. クラスター別プロファイル比較（標準化平均）")
print("=" * 70)

df_ana_z = df_ana.copy()
for j, feat in enumerate(FEAT_SHORT):
    df_ana_z[feat + '_z'] = X_scaled[:, j]

profile = df_ana_z.groupby('cluster')[[f + '_z' for f in FEAT_SHORT]].mean()
print(f"\n【クラスター別 標準化変数平均値（高いほど全国平均より上）】")
print(f"  {'変数':<14}", end='')
for cl in sorted(df_ana['cluster'].unique()):
    print(f" {'CL' + str(cl):>8}", end='')
print()
print("  " + "-" * (14 + 9 * N_CLUSTERS))
for feat in FEAT_SHORT:
    print(f"  {feat:<14}", end='')
    for cl in sorted(df_ana['cluster'].unique()):
        v = profile.loc[cl, feat + '_z']
        print(f" {v:>+8.3f}", end='')
    print()

print("\n\n" + "=" * 70)
print("■ 図の生成（4枚）")
print("=" * 70)

print("図2: PCA バイプロットを作成中...")

fig2, axes2 = plt.subplots(1, 2, figsize=(16, 6.5))
fig2.suptitle('Step2. 主成分分析（PCA）バイプロット\n'
              '2022年度 47都道府県 ごみ・消費・人口・気温変数（SSDSE-B）',
              fontsize=12, fontweight='bold')

print("図4: クラスター別プロファイル棒グラフを作成中...")

fig4, axes4 = plt.subplots(2, 4, figsize=(16, 9))
fig4.suptitle(f'Step4. Ward法クラスタリング（k={N_CLUSTERS}）クラスター別プロファイル比較\n'
              '2022年度 47都道府県（SSDSE-B-2026 実データ）',
              fontsize=12, fontweight='bold')

axes4_flat = axes4.flatten()

# リサイクル率をプロファイル変数に追加
profile_vars = [TARGET] + FEAT_SHORT
n_vars_p = len(profile_vars)

cl_order = sorted(df_ana['cluster'].unique())
x_pos = np.arange(len(cl_order))

for ax_idx, var in enumerate(profile_vars):
    if ax_idx >= len(axes4_flat):
        break
    ax = axes4_flat[ax_idx]

    means = [df_ana[df_ana['cluster'] == cl][var].mean() for cl in cl_order]
    stds  = [df_ana[df_ana['cluster'] == cl][var].std()  for cl in cl_order]

    bars = ax.bar(x_pos, means, yerr=stds, capsize=4,
                  color=[C_CLUSTER[cl] for cl in cl_order],
                  alpha=0.85, edgecolor='white',
                  error_kw={'elinewidth': 1.2, 'capthick': 1.2})

    ax.set_xticks(x_pos)
    ax.set_xticklabels([f'CL{cl}' for cl in cl_order], fontsize=9)
    ax.set_title(var, fontsize=9, fontweight='bold')
    ax.grid(axis='y', alpha=0.3)

    # 全国平均との差を示す水平線
    overall_mean = df_ana[var].mean()
    ax.axhline(overall_mean, color='black', linestyle='--', linewidth=1.2,
               alpha=0.7, label=f'全国平均 ({overall_mean:.1f})')
    ax.legend(fontsize=6.5, loc='upper right')

    # 値ラベル
    for bar, m, s in zip(bars, means, stds):
        ax.text(bar.get_x() + bar.get_width() / 2,
                m + s + (max(means) - min(means)) * 0.03,
                f'{m:.1f}', ha='center', va='bottom', fontsize=7, fontweight='bold')

# 余った軸を非表示
for ax_idx in range(n_vars_p, len(axes4_flat)):
    axes4_flat[ax_idx].set_visible(False)

# クラスター凡例（図全体）
patch_list4 = [mpatches.Patch(color=C_CLUSTER[cl], label=f'クラスター{cl}')
               for cl in cl_order]
fig4.legend(handles=patch_list4, fontsize=10,
            loc='lower right', bbox_to_anchor=(0.98, 0.02),
            title='クラスター', title_fontsize=10)

fig4.text(0.5, 0.005,
          '出典: 統計数理研究所 SSDSE-B-2026 | Ward法（scipy）| エラーバー: ±1 SD',
          ha='center', fontsize=8, color='#666')

plt.tight_layout(rect=[0, 0.03, 1, 0.96])
save_fig('fig4_cluster')