地方創生における三つの「鍵」

# ── ライブラリの読み込みと、SSDSE-B から都道府県指標を作る ──────────────────
import os
import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from scipy import stats
from sklearn.preprocessing import StandardScaler
from sklearn.cluster import KMeans
from sklearn.decomposition import PCA

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)

# 先頭行が列コード、2行目が日本語名。skiprows=[1] で日本語名の行を落とし、
# 年度列 'SSDSE-B-2026' で 2023 年の47都道府県だけを取り出す。
df = pd.read_csv(DATA_B, encoding='cp932', skiprows=[1])
d = df[df['SSDSE-B-2026'] == 2023].copy()

# 原論文の指標定義に対応する形で、SSDSE-B に存在する列だけで指標を作る。
d['自然増減率'] = (d['A4101'] - d['A4200']) / d['A1102'] * 1000     # (出生-死亡)/日本人人口×1000
d['社会増減率'] = (d['A5101'] - d['A5102']) / d['A1102'] * 1000     # (転入-転出)/日本人人口×1000
d['人口増減率'] = d['自然増減率'] + d['社会増減率']
d['大学等進学率'] = d['E4602'] / d['E4601'] * 100                   # 高卒者のうち進学者/高卒者×100
d['高校進学率'] = d['E3702'] / d['E3701'] * 100                     # 中卒者のうち進学者/中卒者×100
d['高齢化率'] = d['A1303'] / (d['A1301'] + d['A1302'] + d['A1303']) * 100
d['合計特殊出生率'] = d['A4103']
d['pref'] = d['Prefecture'].str.replace('[都府県]$', '', regex=True)  # 北海道はそのまま残す

print("=== [1] SSDSE-B から作った都道府県指標（2023年・実データ） ===")
print(f"対象: {len(d)} 都道府県（原論文は市区町村1737・2015年）")
print(d[['pref', '人口増減率', '大学等進学率', '高校進学率', '高齢化率', '合計特殊出生率']]
      .describe().round(2).to_string())

# ── 図1【実再現】：大学等進学率の都道府県ランキング ────────────────────────
rank = d.sort_values('大学等進学率')
mean_adv = rank['大学等進学率'].mean()
colors = ['#e8792b' if v >= mean_adv else '#3f7cc0' for v in rank['大学等進学率']]

fig, ax = plt.subplots(figsize=(9, 11))
ax.barh(rank['pref'], rank['大学等進学率'], color=colors)
ax.axvline(mean_adv, color='#555', ls='--', lw=1)
ax.text(mean_adv + 0.3, 0.5, f'全国平均 {mean_adv:.1f}%', color='#555', fontsize=11)
ax.set_xlabel('大学等進学率（％）＝ 高卒者のうち進学者 / 高卒者 ×100')
ax.set_title('図1【実再現】大学等進学率の都道府県ランキング（SSDSE-B・2023年）\n'
             '※原論文の3つの「鍵」のうち SSDSE-B に残る唯一の変数。市区町村・3年平均の原論文値とは一致しない',
             fontsize=11)
ax.margins(y=0.01)
plt.tight_layout()
plt.savefig(f'{FIG_DIR}/2018_U1_fig1.png', bbox_inches='tight')
plt.close()

print("=== [2] 大学等進学率ランキング（実再現・2023年） ===")
print("上位5:", [f"{p} {v:.1f}%" for p, v in
                 zip(rank['pref'].tail(5)[::-1], rank['大学等進学率'].tail(5)[::-1])])
print("下位5:", [f"{p} {v:.1f}%" for p, v in
                 zip(rank['pref'].head(5), rank['大学等進学率'].head(5))])
print(f"全国平均 {mean_adv:.1f}%")

# ── 図2【報告値の可視化】：原論文 表6 の相関行列を転記して描く ──────────────
# 下の数値は原論文 表6（市区町村1737・2015年・ピアソン積率相関）の報告値そのもの。
# ここで新たに計算はしていない（再計算ではない）。
report_corr = pd.DataFrame(
    {
        '自然増減率': [0.088, 0.278, 0.087, 0.046, -0.679, -0.164, -0.003, 0.149],
        '社会増減率': [-0.023, 0.265, 0.041, 0.026, -0.861, -0.048, 0.062, 0.068],
        '人口増減率': [0.083, 0.288, 0.059, 0.034, -0.864, -0.119, 0.045, 0.089],
    },
    index=['完全失業率', '労働力人口率', '人口10万対の医師数', '人口10万対の病院病床数',
           '老齢人口千対の老人ホーム定員数', '乳幼児人口千対の幼稚園･保育所利用者数',
           '三年間の高校進学率', '三年間の大学進学率'],
)

fig, ax = plt.subplots(figsize=(7.5, 7))
im = ax.imshow(report_corr.values, cmap='RdBu_r', vmin=-1, vmax=1, aspect='auto')
ax.set_xticks(range(report_corr.shape[1]))
ax.set_xticklabels(report_corr.columns)
ax.set_yticks(range(report_corr.shape[0]))
ax.set_yticklabels(report_corr.index, fontsize=10)
for i in range(report_corr.shape[0]):
    for j in range(report_corr.shape[1]):
        v = report_corr.values[i, j]
        ax.text(j, i, f'{v:+.3f}', ha='center', va='center',
                color='white' if abs(v) > 0.5 else '#222', fontsize=11)
fig.colorbar(im, ax=ax, shrink=0.8, label='ピアソン相関係数（原論文 表6 の報告値）')
ax.set_title('図2【報告値の可視化】人口増減と8指標の相関行列\n'
             '※原論文 表6（市区町村1737・2015年）の報告値を転記（再計算ではない）', fontsize=11)
plt.tight_layout()
plt.savefig(f'{FIG_DIR}/2018_U1_fig2.png', bbox_inches='tight')
plt.close()

print("=== [3] 原論文 表6 の相関行列（報告値・転記） ===")
print(report_corr.to_string())
print("※ 病院病床数・高校進学率は人口増減率と5%水準で有意でなかった（原論文）。")

# ── 図3【実再現】：人口増減率 × 大学等進学率（原論文の主張を都道府県で検証）──
x = d['大学等進学率'].values
y = d['人口増減率'].values
r, p = stats.pearsonr(x, y)
slope, intercept = np.polyfit(x, y, 1)
xs = np.linspace(x.min(), x.max(), 100)

fig, ax = plt.subplots(figsize=(9, 6.5))
ax.scatter(x, y, s=45, color='#3f7cc0', edgecolor='white', zorder=3)
ax.plot(xs, slope * xs + intercept, color='#e8792b', lw=2,
        label=f'回帰直線 y={slope:.3f}x{intercept:+.1f}')
for _, row in d.iterrows():
    if row['pref'] in ['東京', '沖縄', '秋田', '京都', '北海道', '愛知']:
        ax.annotate(row['pref'], (row['大学等進学率'], row['人口増減率']),
                    fontsize=10, xytext=(4, 4), textcoords='offset points')
ax.axhline(0, color='#999', lw=0.8)
ax.set_xlabel('大学等進学率（％）')
ax.set_ylabel('人口増減率（‰＝日本人人口千対）')
ax.set_title(f'図3【実再現】人口増減率 × 大学等進学率（都道府県・2023年）\n'
             f'実データで計算した相関 r={r:+.3f}（p={p:.4g}）'
             f'／原論文 表6 の報告値は r=+0.089（市区町村・2015年）', fontsize=11)
ax.legend()
plt.tight_layout()
plt.savefig(f'{FIG_DIR}/2018_U1_fig3.png', bbox_inches='tight')
plt.close()

print("=== [4] 人口増減率 × 大学等進学率（実再現・2023年） ===")
print(f"相関係数 r = {r:+.4f}  (p = {p:.4g})   回帰の傾き = {slope:+.4f}")
print("原論文 表6 の報告値（市区町村・2015年）は r=+0.089（1%有意）。")
print("→ 単位も年次も違うが、大学進学率が高い地域ほど人口が増えやすい向き（正の関連）は一致。")

# ── 図4【実再現】：都道府県の類型化（k-means＋主成分分析）────────────────────
# 原論文の類型化は市区町村を3指標の「人口増加地域の平均」で上回るか否かで9分類する
# ルールベース。ここでは SSDSE-B に残る指標だけを標準化し、k-means でデータ駆動に
# 4類型へまとめ、主成分分析(PCA)で2次元に落として散布図にする（原論文の地図の代替）。
feats = ['人口増減率', '大学等進学率', '高齢化率', '合計特殊出生率']
Z = StandardScaler().fit_transform(d[feats])
km = KMeans(n_clusters=4, random_state=42, n_init=10).fit(Z)
d['cluster'] = km.labels_

# クラスタを「人口増減率の平均」が高い順に並べ替えてラベルを安定させる
order = d.groupby('cluster')['人口増減率'].mean().sort_values(ascending=False).index
relabel = {c: i for i, c in enumerate(order)}
d['type'] = d['cluster'].map(relabel)
names = {0: '成長・都市型', 1: '中間型', 2: '高齢・地方型', 3: '過疎・低出生型'}

pca = PCA(n_components=2).fit(Z)
XY = pca.transform(Z)
palette = ['#2e7d32', '#1565c0', '#e8792b', '#8e24aa']

fig, ax = plt.subplots(figsize=(9.5, 7.5))
for t in range(4):
    m = d['type'].values == t
    ax.scatter(XY[m, 0], XY[m, 1], s=70, color=palette[t],
               edgecolor='white', label=f'類型{t+1}: {names[t]}（{m.sum()}県）', zorder=3)
labels = d['pref'].tolist()
for i in range(len(d)):
    ax.annotate(labels[i], (XY[i, 0], XY[i, 1]), fontsize=8,
                xytext=(3, 3), textcoords='offset points', color='#333')
ax.set_xlabel(f'第1主成分（寄与率 {pca.explained_variance_ratio_[0]*100:.0f}%）')
ax.set_ylabel(f'第2主成分（寄与率 {pca.explained_variance_ratio_[1]*100:.0f}%）')
ax.set_title('図4【実再現】都道府県の類型化（k-means＋PCA・SSDSE-B・2023年）\n'
             '※原論文の市区町村・ルールベース9分類とは手法も単位も異なる現代的再現', fontsize=11)
ax.legend(loc='best')
plt.tight_layout()
plt.savefig(f'{FIG_DIR}/2018_U1_fig4.png', bbox_inches='tight')
plt.close()

print("=== [5] 都道府県の類型化（実再現・k-means4類型） ===")
prof = d.groupby('type')[feats].mean().round(2)
prof.index = [f'類型{t+1}:{names[t]}' for t in prof.index]
print(prof.to_string())
for t in range(4):
    members = d[d['type'] == t]['pref'].tolist()
    print(f'類型{t+1} {names[t]}（{len(members)}県）: ' + '、'.join(members))