都道府県別政治参加（投票率）の決定要因分析
SSDSE-B で迫る社会参加・市民活動の地域格差

import os
import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import statsmodels.api as sm
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'
os.makedirs(FIG_DIR, exist_ok=True)

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)

# 2022年度断面データ
df = df_b[df_b['年度'] == 2022].copy()
df = df.reset_index(drop=True)

df['高齢化率']     = df['65歳以上人口'] / df['総人口'] * 100
df['少子化率']     = df['15歳未満人口'] / df['総人口'] * 100
df['転入超過率']   = (df['転入者数（日本人移動者）'] - df['転出者数（日本人移動者）']) / df['総人口'] * 1000
df['婚姻率']       = df['婚姻件数'] / df['総人口'] * 1000
df['旅券発行率']   = df['一般旅券発行件数'] / df['総人口'] * 1000   # 社会参加代理変数（目的変数）
df['教養娯楽費率'] = df['教養娯楽費（二人以上の世帯）'] / df['消費支出（二人以上の世帯）'] * 100
df['高校進学率']   = df['高等学校卒業者のうち進学者数'] / df['高等学校卒業者数'] * 100
df['リサイクル率'] = df['ごみのリサイクル率']
df['保育所等充足率'] = df['保育所等数'] / df['総人口'] * 10000

region_map = {
    '北海道': '北海道・東北', '青森県': '北海道・東北', '岩手県': '北海道・東北',
    '宮城県': '北海道・東北', '秋田県': '北海道・東北', '山形県': '北海道・東北',
    '福島県': '北海道・東北', '茨城県': '関東', '栃木県': '関東', '群馬県': '関東',
    '埼玉県': '関東', '千葉県': '関東', '東京都': '関東', '神奈川県': '関東',
    '新潟県': '中部', '富山県': '中部', '石川県': '中部', '福井県': '中部',
    '山梨県': '中部', '長野県': '中部', '岐阜県': '中部', '静岡県': '中部', '愛知県': '中部',
    '三重県': '近畿', '滋賀県': '近畿', '京都府': '近畿', '大阪府': '近畿',
    '兵庫県': '近畿', '奈良県': '近畿', '和歌山県': '近畿',
    '鳥取県': '中国・四国', '島根県': '中国・四国', '岡山県': '中国・四国',
    '広島県': '中国・四国', '山口県': '中国・四国', '徳島県': '中国・四国',
    '香川県': '中国・四国', '愛媛県': '中国・四国', '高知県': '中国・四国',
    '福岡県': '九州・沖縄', '佐賀県': '九州・沖縄', '長崎県': '九州・沖縄',
    '熊本県': '九州・沖縄', '大分県': '九州・沖縄', '宮崎県': '九州・沖縄',
    '鹿児島県': '九州・沖縄', '沖縄県': '九州・沖縄'
}
region_colors = {
    '北海道・東北': '#4e9af1',
    '関東':         '#e05c5c',
    '中部':         '#f0a500',
    '近畿':         '#5cb85c',
    '中国・四国':   '#9b59b6',
    '九州・沖縄':   '#f39c12'
}

df['地域区分'] = df['都道府県'].map(region_map)
df['地域色']   = df['地域区分'].map(region_colors)

corr_cols = ['旅券発行率', '高齢化率', '転入超過率', '婚姻率', '教養娯楽費率', '高校進学率', '少子化率', 'リサイクル率']
corr_mat  = df[corr_cols].corr()

fig, ax = plt.subplots(figsize=(9, 7))
im = ax.imshow(corr_mat.values, cmap='RdBu_r', vmin=-1, vmax=1)
plt.colorbar(im, ax=ax, shrink=0.8)
ax.set_xticks(range(len(corr_cols)))
ax.set_yticks(range(len(corr_cols)))
ax.set_xticklabels(corr_cols, rotation=40, ha='right', fontsize=10)
ax.set_yticklabels(corr_cols, fontsize=10)
for i in range(len(corr_cols)):
    for j in range(len(corr_cols)):
        val = corr_mat.values[i, j]
        color = 'white' if abs(val) > 0.6 else 'black'
        ax.text(j, i, f'{val:.2f}', ha='center', va='center', fontsize=9, color=color)
ax.set_title('図2：社会参加・地域特性変数の相関行列（Pearson相関係数, 2022年）',
             fontsize=12, fontweight='bold', pad=12)
fig.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2020_U4_fig2.png'), dpi=150)
plt.close(fig)
print("fig2 saved")

beta_vals = betas.values
beta_pvals = pvals.values

bar_colors = []
for b, p in zip(beta_vals, beta_pvals):
    if p < 0.05:
        bar_colors.append('#c0392b' if b > 0 else '#2980b9')
    else:
        bar_colors.append('#aaaaaa')

fig, ax = plt.subplots(figsize=(8, 5))
y_pos = range(len(X_cols))
bars = ax.barh(y_pos, beta_vals, color=bar_colors, height=0.55, edgecolor='white')
ax.axvline(0, color='black', lw=0.8)
ax.set_yticks(list(y_pos))
ax.set_yticklabels(X_cols, fontsize=11)
ax.set_xlabel('標準化偏回帰係数 (β)', fontsize=12)
ax.set_title(f'図3：重回帰分析 標準化偏回帰係数（R²={r2:.3f}, Adj.R²={r2_adj:.3f}）',
             fontsize=12, fontweight='bold')

for i, (b, p) in enumerate(zip(beta_vals, beta_pvals)):
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else ''
    offset = 0.02 if b >= 0 else -0.02
    ha = 'left' if b >= 0 else 'right'
    ax.text(b + offset, i, f'{b:+.3f}{sig}', va='center', ha=ha, fontsize=10)

from matplotlib.patches import Patch
legend_elements = [
    Patch(facecolor='#c0392b', label='有意（正, p<0.05）'),
    Patch(facecolor='#2980b9', label='有意（負, p<0.05）'),
    Patch(facecolor='#aaaaaa', label='非有意 (ns)'),
]
ax.legend(handles=legend_elements, loc='lower right', fontsize=9)
ax.grid(axis='x', alpha=0.3)
fig.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2020_U4_fig3.png'), dpi=150)
plt.close(fig)
print("fig3 saved")

X_cols = ['高齢化率', '転入超過率', '婚姻率', '教養娯楽費率', '高校進学率']
y_col  = '旅券発行率'

y = df[y_col]
X = df[X_cols]
X_std = (X - X.mean()) / X.std()
X_std_const = sm.add_constant(X_std)

model = sm.OLS(y, X_std_const).fit()

betas  = model.params[1:]
pvals  = model.pvalues[1:]
r2     = model.rsquared
r2_adj = model.rsquared_adj
f_stat = model.fvalue
f_pval = model.f_pvalue

print("=" * 60)
print("OLS 重回帰分析結果（目的変数：旅券発行率）")
print("=" * 60)
print(f"R² = {r2:.3f},  Adj.R² = {r2_adj:.3f}")
print(f"F 統計量 = {f_stat:.2f},  p = {f_pval:.2e}")
print()
print("標準化偏回帰係数:")
for col, b, p in zip(X_cols, betas, pvals):
    sig = "***" if p < 0.001 else "**" if p < 0.01 else "*" if p < 0.05 else "ns"
    print(f"  {col:12s}  beta={b:+.3f}  p={p:.3f}  {sig}")
print()

# 相関分析
print("相関分析（旅券発行率 vs 各説明変数）:")
for col in X_cols:
    r, p = stats.pearsonr(df[y_col], df[col])
    sig = "***" if p < 0.001 else "**" if p < 0.01 else "*" if p < 0.05 else "ns"
    print(f"  {col:12s}  r={r:+.3f}  p={p:.3f}  {sig}")

print()
print("記述統計（旅券発行率）:")
print(f"  全国平均: {df['旅券発行率'].mean():.2f}  (SD={df['旅券発行率'].std():.2f})")
print(f"  最大: {df.loc[df['旅券発行率'].idxmax(), '都道府県']}  {df['旅券発行率'].max():.2f}")
print(f"  最小: {df.loc[df['旅券発行率'].idxmin(), '都道府県']}  {df['旅券発行率'].min():.2f}")

df_sorted = df.sort_values('旅券発行率', ascending=True).copy()

fig, ax = plt.subplots(figsize=(10, 12))
colors = [region_colors[r] for r in df_sorted['地域区分']]
bars = ax.barh(df_sorted['都道府県'], df_sorted['旅券発行率'],
               color=colors, height=0.75, edgecolor='white')
ax.axvline(df['旅券発行率'].mean(), color='#333', lw=1.5, ls='--',
           label=f'全国平均 {df["旅券発行率"].mean():.2f}')
ax.set_xlabel('旅券発行率（人口千人あたり）', fontsize=12)
ax.set_title('図4：都道府県別 旅券発行率ランキング（2022年）', fontsize=13, fontweight='bold')

from matplotlib.patches import Patch
legend_elements = [Patch(facecolor=c, label=r) for r, c in region_colors.items()]
legend_elements.append(plt.Line2D([0], [0], color='#333', ls='--', lw=1.5, label=f'全国平均 {df["旅券発行率"].mean():.2f}'))
ax.legend(handles=legend_elements, loc='lower right', fontsize=9, framealpha=0.9)
ax.grid(axis='x', alpha=0.3)
fig.tight_layout()
fig.savefig(os.path.join(FIG_DIR, '2020_U4_fig4.png'), dpi=150)
plt.close(fig)
print("fig4 saved")

print()
print("全図の生成が完了しました。")
print(f"保存先: {os.path.abspath(FIG_DIR)}")