SDGs指標の地域比較
都道府県別の持続可能性スコア分析

df_raw = pd.read_csv(DATA_B, encoding='shift-jis', header=0)
# 行0がラベル行（日本語名）、行1以降が実データ
df_all = df_raw.iloc[1:].copy()
df_all.columns = df_raw.columns

# 必要な数値列に変換
NUM_COLS = ['A1101', 'A4103', 'E4601', 'E4602',
            'F3102', 'F3103', 'H5610', 'H5614',
            'J2503', 'L3221']
for c in NUM_COLS:
    df_all[c] = pd.to_numeric(df_all[c], errors='coerce')

df_all['年度']    = df_all['SSDSE-B-2026'].astype(int)
df_all['都道府県'] = df_all['Prefecture'].astype(str)

print("=" * 65)
print("■ SSDSE-B-2026 読み込み完了")
print(f"  年度範囲: {df_all['年度'].min()}〜{df_all['年度'].max()}")
print(f"  都道府県数（各年）: {df_all.groupby('年度')['都道府県'].count().iloc[0]}件")
print("=" * 65)

df = df_all[df_all['年度'] == 2022].copy().reset_index(drop=True)

df['消費支出_log']    = np.log(df['L3221'])                         # 経済: 生活水準（対数変換）
df['求人倍率']        = df['F3103'] / df['F3102']                   # 経済: 雇用充実度
df['ごみリサイクル率'] = df['H5614']                                  # 環境: リサイクル率(%)
df['ごみ排出量_inv']  = 1.0 / df['H5610']                           # 環境: 排出少＝良い → 逆数
df['大学進学率']      = df['E4602'] / df['E4601'] * 100             # 社会: 教育到達度(%)
df['保育所密度']      = df['J2503'] / df['A1101'] * 10000           # 社会: 人口1万対保育所数
df['合計特殊出生率']  = df['A4103']                                   # 人口持続性

INDICATOR_COLS = [
    '消費支出_log', '求人倍率', 'ごみリサイクル率', 'ごみ排出量_inv',
    '大学進学率', '保育所密度', '合計特殊出生率'
]
INDICATOR_LABELS = [
    '消費支出\n(log)', '求人倍率', 'ごみ\nリサイクル率', 'ごみ排出量\n(逆数)',
    '大学進学率', '保育所密度', '合計特殊\n出生率'
]

# 欠損値がある行を除外
df_sdgs = df[['都道府県'] + INDICATOR_COLS].dropna().reset_index(drop=True)
prefectures = df_sdgs['都道府県'].values
X_raw = df_sdgs[INDICATOR_COLS].values

print(f"\n■ 2022年 SDGs指標データ: {len(df_sdgs)}都道府県")
print(f"  7指標の基本統計:")
for col in INDICATOR_COLS:
    print(f"    {col}: 平均={df_sdgs[col].mean():.3f}, SD={df_sdgs[col].std():.3f}")

REGION_MAP = {
    '北海道': '北海道・東北',
    '青森県': '北海道・東北', '岩手県': '北海道・東北', '宮城県': '北海道・東北',
    '秋田県': '北海道・東北', '山形県': '北海道・東北', '福島県': '北海道・東北',
    '茨城県': '関東', '栃木県': '関東', '群馬県': '関東', '埼玉県': '関東',
    '千葉県': '関東', '東京都': '関東', '神奈川県': '関東',
    '新潟県': '中部', '富山県': '中部', '石川県': '中部', '福井県': '中部',
    '山梨県': '中部', '長野県': '中部', '岐阜県': '中部', '静岡県': '中部', '愛知県': '中部',
    '三重県': '近畿', '滋賀県': '近畿', '京都府': '近畿', '大阪府': '近畿',
    '兵庫県': '近畿', '奈良県': '近畿', '和歌山県': '近畿',
    '鳥取県': '中国・四国', '島根県': '中国・四国', '岡山県': '中国・四国',
    '広島県': '中国・四国', '山口県': '中国・四国',
    '徳島県': '中国・四国', '香川県': '中国・四国', '愛媛県': '中国・四国', '高知県': '中国・四国',
    '福岡県': '九州・沖縄', '佐賀県': '九州・沖縄', '長崎県': '九州・沖縄',
    '熊本県': '九州・沖縄', '大分県': '九州・沖縄', '宮崎県': '九州・沖縄',
    '鹿児島県': '九州・沖縄', '沖縄県': '九州・沖縄',
}
REGION_ORDER  = ['北海道・東北', '関東', '中部', '近畿', '中国・四国', '九州・沖縄']
REGION_COLORS = {
    '北海道・東北': '#1565C0',
    '関東':         '#E65100',
    '中部':         '#2E7D32',
    '近畿':         '#6A1B9A',
    '中国・四国':   '#795548',
    '九州・沖縄':   '#00838F',
}
regions = np.array([REGION_MAP.get(p, '不明') for p in prefectures])
df_sdgs['地方'] = regions

print("\n" + "=" * 65)
print("■ 図1: PCAスコアプロット（PC1 vs PC2）")
print("=" * 65)

fig1, ax1 = plt.subplots(figsize=(12, 8))

for region in REGION_ORDER:
    mask = regions == region
    if mask.sum() == 0:
        continue
    ax1.scatter(scores[mask, 0], scores[mask, 1],
                color=REGION_COLORS[region], s=80, alpha=0.85,
                edgecolors='white', linewidth=0.7,
                label=f'{region} (n={mask.sum()})', zorder=3)

# 都道府県ラベル（全都道府県に短縮名を表示）
for i, pref in enumerate(prefectures):
    short = pref.replace('県', '').replace('都', '').replace('道', '').replace('府', '')
    ax1.annotate(short, (scores[i, 0], scores[i, 1]),
                 xytext=(3, 3), textcoords='offset points',
                 fontsize=6.5, color='#333', alpha=0.9)

ax1.axhline(0, color='gray', linestyle='--', linewidth=0.8, alpha=0.5)
ax1.axvline(0, color='gray', linestyle='--', linewidth=0.8, alpha=0.5)

ax1.set_xlabel(f'PC1（寄与率 {ev_ratio[0]*100:.1f}%）── 総合持続可能性スコア',
               fontsize=12)
ax1.set_ylabel(f'PC2（寄与率 {ev_ratio[1]*100:.1f}%）── 経済 vs 環境・人口',
               fontsize=12)
ax1.set_title(
    'SDGs指標の主成分スコアプロット（2022年, 都道府県別）\n'
    '〜PC1: 総合持続可能性スコア。右ほど持続可能性が高い傾向〜',
    fontsize=13, fontweight='bold'
)
ax1.legend(loc='lower left', fontsize=9, framealpha=0.88,
           bbox_to_anchor=(0.0, 0.0))
ax1.grid(True, alpha=0.25)
plt.tight_layout()
fig1.savefig(os.path.join(FIG_DIR, '2022_H5_15_fig1_pca_score.png'),
             bbox_inches='tight', dpi=150)
plt.close(fig1)
print("図1保存: 2022_H5_15_fig1_pca_score.png")

# PC1スコア上位・下位を表示
df_pca = df_sdgs.copy()
df_pca['PC1'] = scores[:, 0]
df_pca['PC2'] = scores[:, 1]
top5  = df_pca.nlargest(5, 'PC1')[['都道府県', '地方', 'PC1', 'PC2']]
bot5  = df_pca.nsmallest(5, 'PC1')[['都道府県', '地方', 'PC1', 'PC2']]
print(f"\n  PC1スコア TOP5（持続可能性が高い都道府県）:")
print(top5.to_string(index=False))
print(f"\n  PC1スコア BOTTOM5（持続可能性が低い都道府県）:")
print(bot5.to_string(index=False))

print("\n" + "=" * 65)
print("■ 図2: 主成分負荷量ヒートマップ")
print("=" * 65)

N_PC_SHOW = 4   # PC1〜PC4を表示
load_show = loadings[:, :N_PC_SHOW]

fig2, ax2 = plt.subplots(figsize=(9, 6))

im = ax2.imshow(load_show, cmap='RdBu_r', vmin=-1, vmax=1, aspect='auto')
plt.colorbar(im, ax=ax2, label='因子負荷量', shrink=0.85)

ax2.set_xticks(range(N_PC_SHOW))
ax2.set_xticklabels(
    [f'PC{i+1}\n({ev_ratio[i]*100:.1f}%)' for i in range(N_PC_SHOW)],
    fontsize=11
)
ax2.set_yticks(range(len(INDICATOR_LABELS)))
ax2.set_yticklabels(INDICATOR_LABELS, fontsize=11)

# セル内に数値を表示
for i in range(len(INDICATOR_COLS)):
    for j in range(N_PC_SHOW):
        val = load_show[i, j]
        color = 'white' if abs(val) > 0.55 else '#222'
        ax2.text(j, i, f'{val:.2f}',
                 ha='center', va='center', fontsize=10,
                 fontweight='bold', color=color)

ax2.set_title(
    'SDGs指標の主成分負荷量ヒートマップ（PC1〜PC4）\n'
    '〜赤：正の負荷（その指標がPCと同方向）、青：負の負荷〜',
    fontsize=13, fontweight='bold'
)
ax2.set_xlabel('主成分（説明した分散の割合）', fontsize=11)
ax2.set_ylabel('SDGs指標', fontsize=11)
plt.tight_layout()
fig2.savefig(os.path.join(FIG_DIR, '2022_H5_15_fig2_loadings.png'),
             bbox_inches='tight', dpi=150)
plt.close(fig2)
print("図2保存: 2022_H5_15_fig2_loadings.png")

print(f"\n  因子負荷量（PC1〜PC4）:")
df_load = pd.DataFrame(load_show, index=INDICATOR_COLS,
                        columns=[f'PC{i+1}' for i in range(N_PC_SHOW)])
print(df_load.round(3).to_string())

print("\n" + "=" * 65)
print("■ 図3: Ward法クラスタリング デンドログラム")
print("=" * 65)

# Ward法でクラスタリング（PCA後のスコア全7成分を使用）
Z = linkage(X_std, method='ward', metric='euclidean')

# 都道府県の短縮名リスト
pref_short = [p.replace('県', '').replace('都', '').replace('道', '').replace('府', '')
              for p in prefectures]

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

dn = dendrogram(Z,
                labels=pref_short,
                ax=ax3,
                color_threshold=0.5 * Z[-5:, 2].mean(),
                leaf_font_size=8.5,
                leaf_rotation=75)

ax3.set_title(
    'Ward法クラスタリング デンドログラム（SDGs7指標, 47都道府県）\n'
    '〜類似した持続可能性プロファイルを持つ都道府県がグループ化される〜',
    fontsize=13, fontweight='bold'
)
ax3.set_xlabel('都道府県', fontsize=11)
ax3.set_ylabel('Ward距離（ユークリッド距離）', fontsize=11)
ax3.yaxis.grid(True, alpha=0.3)
plt.tight_layout()
fig3.savefig(os.path.join(FIG_DIR, '2022_H5_15_fig3_dendrogram.png'),
             bbox_inches='tight', dpi=150)
plt.close(fig3)
print("図3保存: 2022_H5_15_fig3_dendrogram.png")

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

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

import os
FIG_DIR = 'html/figures'
DATA_B  = 'data/raw/SSDSE-B-2026.csv'
os.makedirs(FIG_DIR, exist_ok=True)

scaler = StandardScaler()
X_std = scaler.fit_transform(X_raw)

pca = PCA(n_components=7)
pca.fit(X_std)
scores = pca.transform(X_std)          # 主成分スコア (47 × 7)
loadings = pca.components_.T           # 因子負荷行列 (7変数 × 7成分)
ev_ratio = pca.explained_variance_ratio_

print(f"\n■ PCA 寄与率:")
for i, ev in enumerate(ev_ratio[:5]):
    cum = ev_ratio[:i+1].sum()
    print(f"  PC{i+1}: {ev*100:.1f}%  (累積 {cum*100:.1f}%)")

print("\n" + "=" * 65)
print("■ 図4: 地域別レーダーチャート（6地域 × 4指標）")
print("=" * 65)

# レーダーチャートに使う4指標（解釈しやすいものを選択）
RADAR_COLS   = ['求人倍率', 'ごみリサイクル率', '大学進学率', '合計特殊出生率']
RADAR_LABELS = ['求人倍率\n（雇用）', 'ごみリサイクル率\n（環境）',
                '大学進学率\n（教育）', '合計特殊出生率\n（人口）']

# 地域別平均を計算（正規化用に全国平均・標準偏差を使用）
df_radar = df_sdgs[['地方'] + RADAR_COLS].copy()
region_means = df_radar.groupby('地方')[RADAR_COLS].mean()

# 全国平均・SD で0〜1正規化（各列を min-max スケール）
col_min = df_sdgs[RADAR_COLS].min()
col_max = df_sdgs[RADAR_COLS].max()
region_norm = (region_means - col_min) / (col_max - col_min)

print("\n  地域別 正規化スコア（0〜1スケール）:")
print(region_norm.round(3).to_string())

# レーダーチャートの角度計算
N_VARS = len(RADAR_COLS)
angles = np.linspace(0, 2 * np.pi, N_VARS, endpoint=False).tolist()
angles += angles[:1]  # 閉じるために先頭を末尾に追加

fig4, ax4 = plt.subplots(figsize=(10, 8),
                          subplot_kw=dict(projection='polar'))

# 全体の平均をグレーの参照線として描画
national_mean = np.ones(N_VARS) * 0.5  # 正規化後の中央値（参考）
vals_nat = national_mean.tolist() + national_mean[:1].tolist()
ax4.plot(angles, vals_nat, 'k--', linewidth=1.0, alpha=0.4, label='_nolegend_')
ax4.fill(angles, vals_nat, color='gray', alpha=0.08)

# 各地域をプロット
for region in REGION_ORDER:
    if region not in region_norm.index:
        continue
    vals = region_norm.loc[region, RADAR_COLS].values.tolist()
    vals += vals[:1]  # 閉じる
    ax4.plot(angles, vals,
             color=REGION_COLORS[region], linewidth=2.2,
             label=region, alpha=0.9)
    ax4.fill(angles, vals,
             color=REGION_COLORS[region], alpha=0.12)

# 軸ラベル設定
ax4.set_xticks(angles[:-1])
ax4.set_xticklabels(RADAR_LABELS, fontsize=10.5)
ax4.set_ylim(0, 1)
ax4.set_yticks([0.25, 0.5, 0.75])
ax4.set_yticklabels(['25%', '50%', '75%'], fontsize=8, color='#666')
ax4.set_title(
    '地域別 SDGs4指標 レーダーチャート（2022年）\n'
    '〜各指標を全国min-max正規化。グレー点線は中央値（0.5）〜',
    fontsize=13, fontweight='bold', pad=20
)
ax4.legend(loc='lower left', bbox_to_anchor=(-0.18, -0.08),
           fontsize=9.5, framealpha=0.88)
ax4.grid(True, alpha=0.3)
plt.tight_layout()
fig4.savefig(os.path.join(FIG_DIR, '2022_H5_15_fig4_radar.png'),
             bbox_inches='tight', dpi=150)
plt.close(fig4)
print("図4保存: 2022_H5_15_fig4_radar.png")

print("\n" + "=" * 65)
print("■ 全図生成完了（4枚）")
print(f"  fig1_pca_score.png   : PCAスコアプロット（PC1 vs PC2, 地域色分け）")
print(f"  fig2_loadings.png    : 変数の主成分負荷量ヒートマップ")
print(f"  fig3_dendrogram.png  : Ward法クラスタリング デンドログラム")
print(f"  fig4_radar.png       : 地域別レーダーチャート（6地域 × 4指標）")
print("=" * 65)

print(f"\n■ PCA寄与率サマリ:")
for i in range(4):
    print(f"  PC{i+1}: {ev_ratio[i]*100:.1f}% (累積: {ev_ratio[:i+1].sum()*100:.1f}%)")

print(f"\n■ PC1スコア上位5都道府県（持続可能性が高い）:")
print(top5[['都道府県', '地方', 'PC1']].to_string(index=False))
print(f"\n■ PC1スコア下位5都道府県（持続可能性が低い）:")
print(bot5[['都道府県', '地方', 'PC1']].to_string(index=False))