地方創生目標指標に関する変化要因ネットワークの推定と
それに基づく地域間連携策の提案

import os
import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import networkx as nx
from sklearn.covariance import GraphicalLasso

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

FIG_DIR = 'html/figures'
DATA_A = 'data/raw/SSDSE-A-2025.csv'
os.makedirs(FIG_DIR, exist_ok=True)

# ===== ステップ1: 図1 原論文 図表3 の相関係数（報告値）の可視化 =====
# 原論文 図表3: 2010年総人口と各目標指標の増加率（2010→2015年）の相関係数。
# 以下の数値は原論文 図表3 に記載の報告値の転記であり、再計算ではない。
paper_corr = {
    '就業者増加率': 0.127, '女性就業者増加率': 0.149,
    '完全失業者増加率': 0.097, '正規雇用者増加率': 0.059,
    '転出増加率': 0.036, '転入増加率': 0.068,
    '転出超過増加率': -0.047, '婚姻件数増加率': -0.0046,
    '出生数増加率': 0.027, '保育所等増加率': 0.012,
}
fig, ax = plt.subplots(figsize=(9, 4.8))
names = list(paper_corr.keys())
vals = [paper_corr[k] for k in names]
colors = ['#1565C0' if v >= 0 else '#C62828' for v in vals]
ax.barh(range(len(names))[::-1], vals, color=colors, height=0.6)
ax.set_yticks(range(len(names))[::-1])
ax.set_yticklabels(names, fontsize=10)
ax.axvline(0, color='#555', lw=0.8)
ax.axvline(0.2, color='#999', lw=0.8, ls='--')
ax.axvline(-0.2, color='#999', lw=0.8, ls='--')
ax.set_xlim(-0.25, 0.25)
ax.set_xlabel('2010年総人口との相関係数（原論文 図表3 の報告値）')
ax.set_title('分析1: 目標指標の増加率と人口規模の相関\n'
             '（原論文 図表3 の報告値の可視化。再計算ではない）', fontsize=11)
for i, v in zip(range(len(names))[::-1], vals):
    ax.text(v + (0.008 if v >= 0 else -0.008), i, f'{v:.3f}',
            va='center', ha='left' if v >= 0 else 'right', fontsize=9)
fig.tight_layout()
fig.savefig(f'{FIG_DIR}/2019_U1_fig1.png', bbox_inches='tight')
plt.close(fig)
print('【図1】原論文 図表3 の報告値（2010年総人口と増加率の相関係数）')
print('  最大: 女性就業者増加率 r=0.149 / 最小: 転出超過増加率 r=-0.047')
print('  → いずれも |r|<0.2 で、人口規模は増加率にほぼ影響しない（原論文の結論）')

# ===== ステップ2: SSDSE-A-2025 の読み込み（手法デモ用の現行データ） =====
raw = pd.read_csv(DATA_A, encoding='cp932', header=None, skiprows=3)
raw.columns = pd.read_csv(DATA_A, encoding='cp932', header=None, nrows=3).iloc[2]
print(f'\n【SSDSE-A-2025】{len(raw)}市区町村を読み込み'
      f'（原論文の SSDSE-2019A ベース1722市町村とは別の年次・版）')

# ===== ステップ3: 図2 前処理デモ — ベキ分布を人口あたりに変換 =====
# 原論文 図表2: 就業者数の「変化数」がベキ分布 → 2010年人口で割って正規分布に近づけた。
# ここでは現行データの「水準値」で同じ変換の効果を体験する（変化数のデモではない）。
emp = raw['就業者数'].astype(float)
emp_pc = emp / raw['総人口'].astype(float)
fig, axes = plt.subplots(1, 2, figsize=(9.6, 3.8))
axes[0].hist(emp / 1e4, bins=120, color='#26A69A')
axes[0].set_title('(a) 就業者数（万人）: 裾の重い分布', fontsize=10.5)
axes[0].set_xlabel('就業者数（万人）')
axes[0].set_ylabel('市区町村数')
axes[1].hist(emp_pc, bins=60, color='#26A69A')
axes[1].set_title('(b) 就業者数 ÷ 総人口: 対称な分布に近づく', fontsize=10.5)
axes[1].set_xlabel('就業者数 ÷ 総人口')
fig.suptitle('前処理デモ: 原論文 図表2 と同じ「人口で割る」変換（SSDSE-A-2025 の水準値。'
             '原論文は2010→2015年の変化数に適用）', fontsize=10)
fig.tight_layout(rect=[0, 0, 1, 0.90])
fig.savefig(f'{FIG_DIR}/2019_U1_fig2.png', bbox_inches='tight')
plt.close(fig)
print('\n【図2】前処理デモ（就業者数）')
print(f'  変換前の歪度: {emp.skew():.2f} → 人口あたり変換後: {emp_pc.skew():.2f}')

# ===== ステップ4: 変数の構築（人口あたり変換と符号反転） =====
# 原論文にならい、規模の変数は総人口で割り、比率（％）の変数はそのまま使う。
goal_vars = ['就業者数', '就業者数（女）', '完全失業者数',
             '転入者数（日本人移動者）', '転出者数（日本人移動者）',
             '婚姻件数', '出生数', '保育所等数（基本票）']  # 目標指標（8個）
other_cnt = ['15～64歳人口', '65歳以上人口', '死亡数', '事業所数（民営）',
             '従業者数（民営）', '歳入決算総額（市町村財政）',
             '歳出決算総額（市町村財政）', '地方税（市町村財政）',
             '民生費（市町村財政）', '教育費（市町村財政）',
             '災害復旧費（市町村財政）']                     # 規模の変数
ratio_vars = ['経常収支比率（市町村財政）', '実質公債費比率（市町村財政）']  # 比率

X = pd.DataFrame(index=raw.index)
pop = raw['総人口'].astype(float)
for c in goal_vars + other_cnt:
    X[c] = raw[c].astype(float) / pop
for c in ratio_vars:
    X[c] = raw[c].astype(float)
X = X.dropna()
X = X[pop.loc[X.index] > 0]

# 多重共線性の確認（原論文は歳入・歳出・投資的経費が r>=0.9 → 歳出等を除外）
r_in_out = X['歳入決算総額（市町村財政）'].corr(X['歳出決算総額（市町村財政）'])
print(f'\n【変数の構築】n={len(X)}市区町村, {X.shape[1]}変数')
print(f'  歳入決算総額 と 歳出決算総額 の相関: r={r_in_out:.3f}'
      ' → 原論文と同様ほぼ同一の変数とみなし歳出を除外')
X = X.drop(columns=['歳出決算総額（市町村財政）'])

# 符号の反転: 地方創生にとって「減ることが望ましい」変数に -1 を掛けて
# 解釈を揃える（原論文 2.3節・3節の処理。原論文は変化数に適用）
flip = ['完全失業者数', '転出者数（日本人移動者）', '死亡数',
        '経常収支比率（市町村財政）', '実質公債費比率（市町村財政）']
X[flip] = -X[flip]
Z = (X - X.mean()) / X.std()  # 標準化（平均0・分散1）
short = {c: c.replace('（市町村財政）', '').replace('（日本人移動者）', '')
             .replace('（基本票）', '').replace('（民営）', '') for c in X.columns}
print(f'  符号反転: {len(flip)}変数（完全失業者数・転出者数・死亡数・経常収支比率・実質公債費比率）')

# ===== ステップ5: glasso のハイパーパラメータを AIC で探索 =====
# 原論文は λ を 0〜0.5 の 0.001 刻みで探索し AIC 最小を選択。
# 本デモは 0.005〜0.5 の 0.005 刻み（計算時間短縮のため）。
n, p = Z.shape
results = []
for alpha in np.arange(0.005, 0.5001, 0.005):
    try:
        gl = GraphicalLasso(alpha=alpha, max_iter=500).fit(Z.values)
    except Exception:
        continue
    k = (np.abs(gl.precision_[np.triu_indices(p, 1)]) > 1e-8).sum() + p
    aic = -2 * n * gl.score(Z.values) + 2 * k
    results.append((alpha, aic, k - p, gl))
best_alpha, best_aic, best_edges, best_gl = min(results, key=lambda t: t[1])
print(f'\n【glasso デモ】λ を {len(results)}通り探索（AIC 基準）')
print(f'  AIC 最小: λ={best_alpha:.3f}, AIC={best_aic:.1f}, エッジ数={best_edges}')
print('  ※原論文は 0.001 刻みで探索し、変化数データに対して AIC 最小の λ を採用')

# ===== ステップ6: 図3 推定されたネットワークの可視化（手法デモ） =====
# 偏相関 r_jk = -ω_jk / sqrt(ω_jj・ω_kk) をエッジの重みとする
prec = best_gl.precision_
d = np.sqrt(np.diag(prec))
pcorr = -prec / np.outer(d, d)
np.fill_diagonal(pcorr, 0.0)
cols = list(Z.columns)
G = nx.Graph()
G.add_nodes_from(range(p))
for i in range(p):
    for j in range(i + 1, p):
        if abs(pcorr[i, j]) > 1e-4:
            G.add_edge(i, j, w=pcorr[i, j], dist=1.0 / abs(pcorr[i, j]))

# 正・負それぞれの係数ネットワークで媒介中心性を計算し、差をとる（原論文 2.3節）
Gp = nx.Graph()
Gn = nx.Graph()
Gp.add_nodes_from(range(p))
Gn.add_nodes_from(range(p))
for u, v, d_ in G.edges(data=True):
    (Gp if d_['w'] > 0 else Gn).add_edge(u, v, dist=d_['dist'])
bc_p = nx.betweenness_centrality(Gp, weight='dist')
bc_n = nx.betweenness_centrality(Gn, weight='dist')
score = {i: bc_p.get(i, 0.0) - bc_n.get(i, 0.0) for i in range(p)}

pos = nx.spring_layout(G, seed=42, k=1.4)
fig, ax = plt.subplots(figsize=(9.6, 7.6))
sizes = [300 + 2600 * max(score[i], 0) for i in range(p)]
node_col = ['#FFC107' if cols[i] in goal_vars else '#90CAF9' for i in range(p)]
top3 = sorted(score, key=score.get, reverse=True)[:3]
for i in top3:
    node_col[i] = '#F57F17' if cols[i] in goal_vars else '#1565C0'
epos = [(u, v) for u, v, d_ in G.edges(data=True) if d_['w'] > 0]
eneg = [(u, v) for u, v, d_ in G.edges(data=True) if d_['w'] < 0]
nx.draw_networkx_edges(G, pos, edgelist=epos, edge_color='#78909C',
                       width=[3 * abs(G[u][v]['w']) + 0.3 for u, v in epos], ax=ax)
nx.draw_networkx_edges(G, pos, edgelist=eneg, edge_color='#E57373', style='dashed',
                       width=[3 * abs(G[u][v]['w']) + 0.3 for u, v in eneg], ax=ax)
nx.draw_networkx_nodes(G, pos, node_size=sizes, node_color=node_col,
                       edgecolors='#455A64', ax=ax)
nx.draw_networkx_labels(G, pos, labels={i: short[cols[i]] for i in range(p)},
                        font_size=8, font_family='Hiragino Sans', ax=ax)
ax.set_title('glasso により推定した変数間ネットワーク（手法デモ: SSDSE-A-2025 の水準値, '
             f'λ={best_alpha:.3f}, AIC最小）\n'
             '黄=目標指標系 / 青=その他 / 濃色=媒介中心性スコア上位3 / '
             '実線=正・破線=負の偏相関（原論文 図表4(a) の再計算ではない）',
             fontsize=10)
ax.axis('off')
fig.tight_layout()
fig.savefig(f'{FIG_DIR}/2019_U1_fig3.png', bbox_inches='tight')
plt.close(fig)
print(f'\n【図3】ネットワーク描画: ノード{p}個, 正エッジ{len(epos)}本, 負エッジ{len(eneg)}本')

# ===== ステップ7: 図4 媒介中心性スコア（正−負）の可視化（手法デモ） =====
order = sorted(range(p), key=lambda i: score[i])
fig, ax = plt.subplots(figsize=(8.6, 6.4))
vals4 = [score[i] for i in order]
cols4 = ['#1565C0' if v >= 0 else '#C62828' for v in vals4]
ax.barh(range(p), vals4, color=cols4, height=0.62)
ax.set_yticks(range(p))
ax.set_yticklabels([short[cols[i]] for i in order], fontsize=9)
ax.axvline(0, color='#555', lw=0.8)
ax.set_xlabel('媒介中心性スコア（正の係数ネットワーク − 負の係数ネットワーク）')
ax.set_title('手法デモ: 各変数の媒介中心性スコア（SSDSE-A-2025 の水準値）\n'
             '原論文 図表4(b) と同じ手順（距離=1/|係数|）。変化数分析の再計算ではない',
             fontsize=10.5)
fig.tight_layout()
fig.savefig(f'{FIG_DIR}/2019_U1_fig4.png', bbox_inches='tight')
plt.close(fig)
print('\n【図4】媒介中心性スコア（デモ）上位5変数:')
for i in sorted(score, key=score.get, reverse=True)[:5]:
    print(f'  {short[cols[i]]}: {score[i]:+.4f}')

print('\n【原論文 4.2節の報告】媒介中心性が最も高かった変数（変化数, 順に）:')
print('  市町村の一般財源, 15〜64歳人口, 課税対象所得, 地方税, 歳入決算額')
print('  負の値が大きい変数: 災害復旧費, 実質公債費比率, 経常収支比率')
print('  → 経済規模の変化に関する変数群 =「稼ぐ力」が地方創生の出発点（原論文の結論）')

print('\n=== 完了 ===')
for i in range(1, 5):
    print(f'図{i}: html/figures/2019_U1_fig{i}.png')

# ===== ステップ7: 図4 媒介中心性スコア（正−負）の可視化（手法デモ） =====
order = sorted(range(p), key=lambda i: score[i])
fig, ax = plt.subplots(figsize=(8.6, 6.4))
vals4 = [score[i] for i in order]
cols4 = ['#1565C0' if v >= 0 else '#C62828' for v in vals4]
ax.barh(range(p), vals4, color=cols4, height=0.62)
ax.set_yticks(range(p))
ax.set_yticklabels([short[cols[i]] for i in order], fontsize=9)
ax.axvline(0, color='#555', lw=0.8)
ax.set_xlabel('媒介中心性スコア（正の係数ネットワーク − 負の係数ネットワーク）')
ax.set_title('手法デモ: 各変数の媒介中心性スコア（SSDSE-A-2025 の水準値）\n'
             '原論文 図表4(b) と同じ手順（距離=1/|係数|）。変化数分析の再計算ではない',
             fontsize=10.5)
fig.tight_layout()
fig.savefig(f'{FIG_DIR}/2019_U1_fig4.png', bbox_inches='tight')
plt.close(fig)
print('\n【図4】媒介中心性スコア（デモ）上位5変数:')
for i in sorted(score, key=score.get, reverse=True)[:5]:
    print(f'  {short[cols[i]]}: {score[i]:+.4f}')

print('\n【原論文 4.2節の報告】媒介中心性が最も高かった変数（変化数, 順に）:')
print('  市町村の一般財源, 15〜64歳人口, 課税対象所得, 地方税, 歳入決算額')
print('  負の値が大きい変数: 災害復旧費, 実質公債費比率, 経常収支比率')
print('  → 経済規模の変化に関する変数群 =「稼ぐ力」が地方創生の出発点（原論文の結論）')

print('\n=== 完了 ===')
for i in range(1, 5):
    print(f'図{i}: html/figures/2019_U1_fig{i}.png')