何が女性議員比率を左右するのか

import numpy as np
import pandas as pd
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
from matplotlib.gridspec import GridSpec
import statsmodels.api as sm
from linearmodels.panel import PooledOLS, PanelOLS, RandomEffects
from scipy import stats
import warnings
import os

warnings.filterwarnings('ignore')

# 日本語フォント設定（macOS）
plt.rcParams['font.family'] = 'Hiragino Sans'
plt.rcParams['axes.unicode_minus'] = False
plt.rcParams['font.size'] = 10

# 出力先ディレクトリ
OUT_DIR  = 'html/figures'
DATA_DIR = 'data/2025_U1'
os.makedirs(OUT_DIR, exist_ok=True)

def save_fig(name):
    path = os.path.join(OUT_DIR, f'2025_U1_{name}.png')
    plt.savefig(path, dpi=150, bbox_inches='tight', facecolor='white')
    plt.close()
    print(f'  → 保存: {path}')

print("=" * 60)
print("【Step 1】パネルデータの読み込み")
print("  出典: SSDSE-B-2026.csv / SSDSE-E-2026.csv（実データ）")
print("=" * 60)

_csv_path = os.path.join(DATA_DIR, '2025_U1_panel.csv')
if not os.path.exists(_csv_path):
    raise FileNotFoundError(
        f"CSVが見つかりません: {_csv_path}\n"
        "先に code/2025_U1_data_prep.py を実行してください。"
    )
df = pd.read_csv(_csv_path)

# 変数名の整理
ENTITY_COL = '都道府県'
TIME_COL   = '年度'
Y_VAR      = '転入率'

PREFS = df[ENTITY_COL].unique().tolist()
YEARS = sorted(df[TIME_COL].unique().tolist())
N_PREFS = len(PREFS)
N_YEARS = len(YEARS)

print(f"  データ形状: {df.shape}")
print(f"  都道府県数: {N_PREFS}, 年数: {N_YEARS} ({min(YEARS)}〜{max(YEARS)})")
print()

print("=" * 60)
print("【Step 2】記述統計")
print("  （SSDSE-B/E 実データに基づく 47都道府県×9年パネル）")
print("=" * 60)

vars_to_describe = [
    '転入率', '高齢化率', '年少人口比率', '合計特殊出生率',
    '婚姻率', '保育所等数', '年平均気温', '教育費', '人口密度',
]
desc = df[vars_to_describe].describe().T[['mean', 'std', 'min', 'max', 'count']]
desc.columns = ['平均', '標準偏差', '最小値', '最大値', 'N']
print(desc.round(3))
print()

print()
print("=" * 60)
print("【Step 5】Hausman 検定：FE か RE か？")
print()
print("  【Hausman 検定の考え方】")
print("  ・帰無仮説 H0：個別効果と説明変数は「無相関」（RE が適切）")
print("  ・対立仮説 H1：個別効果と説明変数は「有相関」（FE が適切）")
print("  ・p値 < 0.05 → H0 棄却 → FE を採用")
print("=" * 60)

def hausman_test(fe_result, re_result):
    """Hausman 検定の実装
    H = (b_FE - b_RE)' × [Var(b_FE) - Var(b_RE)]^{-1} × (b_FE - b_RE)
    H は漸近的にカイ二乗分布に従う（自由度 = 係数の数）
    固有値分解により数値安定性を確保する。
    """
    fe_params = fe_result.params
    re_params = re_result.params
    common = [v for v in fe_params.index if v in re_params.index]

    b_fe = fe_params[common].values
    b_re = re_params[common].values

    V_fe  = fe_result.cov.loc[common, common].values
    V_re  = re_result.cov.loc[common, common].values
    V_diff = V_fe - V_re
    diff   = b_fe - b_re

    eigvals, eigvecs = np.linalg.eigh(V_diff)
    pos_mask = eigvals > 1e-10
    if pos_mask.sum() == 0:
        H_stat = float(diff @ np.linalg.pinv(V_fe + V_re) @ diff)
        df_deg  = len(common)
    else:
        V_inv = (eigvecs[:, pos_mask]
                 @ np.diag(1.0 / eigvals[pos_mask])
                 @ eigvecs[:, pos_mask].T)
        H_stat = float(diff @ V_inv @ diff)
        df_deg  = int(pos_mask.sum())

    H_stat = max(H_stat, 0.0)
    p_val  = 1 - stats.chi2.cdf(H_stat, df=df_deg)
    return H_stat, df_deg, p_val

H_stat, df_deg, p_val = hausman_test(fe_res, re_res)
adopt_fe = p_val < 0.05

print(f"\n  Hausman検定: χ²={H_stat:.2f}, 自由度={df_deg}, p={p_val:.4f}")
if adopt_fe:
    print("  → p < 0.05: H0 棄却 → 固定効果モデル(FE)を採用")
else:
    print("  → p ≥ 0.05: H0 採択 → 変量効果モデル(RE)を採用（FEより効率的）")
print()

print("  図1：転入率の時系列推移（都道府県別）")

fig, ax = plt.subplots(figsize=(12, 6))
colors = plt.cm.tab20(np.linspace(0, 1, N_PREFS))
pref_avg_by_year = df.groupby(TIME_COL)[Y_VAR].mean()

# 上位・下位都道府県を目立たせる（東京・大阪など高転入率）
top_prefs = df[df[TIME_COL] == max(YEARS)].nlargest(3, Y_VAR)[ENTITY_COL].tolist()
bottom_prefs = df[df[TIME_COL] == max(YEARS)].nsmallest(3, Y_VAR)[ENTITY_COL].tolist()

for i, pref in enumerate(PREFS):
    pref_data = df[df[ENTITY_COL] == pref].sort_values(TIME_COL)
    if pref in top_prefs:
        ax.plot(pref_data[TIME_COL], pref_data[Y_VAR],
                color='#E53935', alpha=0.9, linewidth=1.4, zorder=4)
        last_row = pref_data.iloc[-1]
        ax.text(last_row[TIME_COL] + 0.1, last_row[Y_VAR], pref, fontsize=7, color='#E53935')
    elif pref in bottom_prefs:
        ax.plot(pref_data[TIME_COL], pref_data[Y_VAR],
                color='#1565C0', alpha=0.9, linewidth=1.4, zorder=4)
        last_row = pref_data.iloc[-1]
        ax.text(last_row[TIME_COL] + 0.1, last_row[Y_VAR], pref, fontsize=7, color='#1565C0')
    else:
        ax.plot(pref_data[TIME_COL], pref_data[Y_VAR],
                color=colors[i], alpha=0.25, linewidth=0.8)

ax.plot(pref_avg_by_year.index, pref_avg_by_year.values,
        color='black', linewidth=2.5, linestyle='--', label='全国平均', zorder=5)

ax.set_xlabel('年度', fontsize=12)
ax.set_ylabel('転入率（千人あたり）', fontsize=12)
ax.set_title('図1：都道府県別 転入率の時系列推移（2014〜2022年）\n'
             '出典: SSDSE-B（実データ）  赤=高転入率、青=低転入率、破線=全国平均', fontsize=11)
ax.legend(fontsize=9, loc='upper right')
ax.set_xlim(min(YEARS) - 0.3, max(YEARS) + 1.5)
ax.grid(True, alpha=0.3)
save_fig('fig1_trend')

print()
print("  図2：OLS / FE / RE の係数比較（95% CI付き）")

var_labels = {
    '高齢化率':       '高齢化率（%）',
    '年少人口比率':   '年少人口比率（%）',
    '合計特殊出生率': '合計特殊出生率',
    '婚姻率':         '婚姻率（千人）',
    '保育所等数':     '保育所等数',
    '年平均気温':     '年平均気温（℃）',
    '教育費':         '教育費（世帯）',
    '人口密度':       '人口密度（人/km²）',
}

fig, axes = plt.subplots(1, 3, figsize=(15, 6), sharey=True)
models_info = [
    ('Pooled OLS',   ols_res, 'steelblue'),
    ('固定効果(FE)', fe_res,  'darkorange'),
    ('変量効果(RE)', re_res,  'forestgreen'),
]

for ax, (name, res, color) in zip(axes, models_info):
    coefs, ci_lo, ci_hi, labels = [], [], [], []

    for var in x_vars:
        if var not in res.params.index:
            continue
        b = res.params[var]
        se = res.std_errors[var] if hasattr(res, 'std_errors') else res.bse[var]
        coefs.append(b)
        ci_lo.append(b - 1.96 * se)
        ci_hi.append(b + 1.96 * se)
        labels.append(var_labels.get(var, var))

    y_pos = range(len(coefs))
    ax.barh(y_pos, coefs, color=color, alpha=0.7, height=0.6)
    ax.errorbar(coefs, y_pos, xerr=[
        np.array(coefs) - np.array(ci_lo),
        np.array(ci_hi) - np.array(coefs),
    ], fmt='none', color='black', capsize=3, linewidth=1)
    ax.axvline(0, color='black', linewidth=0.8, linestyle='--')
    ax.set_yticks(y_pos)
    ax.set_yticklabels(labels, fontsize=9)
    ax.set_title(f'{name}\n(N={int(res.nobs)}, R²={res.rsquared:.3f})', fontsize=10)
    ax.set_xlabel('標準化偏回帰係数（95% CI）', fontsize=9)
    ax.grid(True, axis='x', alpha=0.3)

plt.suptitle('図2：3モデルの回帰係数比較（標準化係数＋95%信頼区間）\n'
             '被説明変数: 転入率（千人あたり）｜出典: SSDSE-B 実データ', fontsize=11, y=1.02)
plt.tight_layout()
save_fig('fig2_fe')

print("  図3：Hausman検定の概念図と検定結果")

fig, axes = plt.subplots(1, 2, figsize=(13, 6))

# (a) FE vs RE の係数比較（実データ）
ax = axes[0]
common_vars = [v for v in x_vars if v in fe_res.params.index]
fe_coefs = [fe_res.params[v] for v in common_vars]
re_coefs = [re_res.params[v] for v in common_vars]
labs     = [var_labels.get(v, v) for v in common_vars]

y_pos = np.arange(len(common_vars))
ax.scatter(fe_coefs, y_pos + 0.15, color='darkorange', s=60, zorder=4,
           label='FE 係数', marker='o')
ax.scatter(re_coefs, y_pos - 0.15, color='forestgreen', s=60, zorder=4,
           label='RE 係数', marker='s')
for i, (fe_c, re_c) in enumerate(zip(fe_coefs, re_coefs)):
    ax.plot([fe_c, re_c], [y_pos[i] + 0.15, y_pos[i] - 0.15],
            color='gray', linewidth=0.8, linestyle=':')

ax.axvline(0, color='black', linewidth=0.8, linestyle='--', alpha=0.5)
ax.set_yticks(y_pos)
ax.set_yticklabels(labs, fontsize=9)
ax.set_xlabel('標準化偏回帰係数', fontsize=10)
ax.set_title('FE と RE の係数比較\n（乖離が大きい→個体効果と変数に相関→FE採用）', fontsize=10)
ax.legend(fontsize=9)
ax.grid(True, axis='x', alpha=0.3)
ax.invert_yaxis()

# (b) Hausman 検定結果サマリー
ax2 = axes[1]
ax2.set_xlim(0, 10); ax2.set_ylim(0, 10); ax2.axis('off')

# カイ二乗分布の参考図
x_chi = np.linspace(0, max(H_stat * 1.5, 20), 300)
y_chi = stats.chi2.pdf(x_chi, df=df_deg)
crit  = stats.chi2.ppf(0.95, df=df_deg)

ax_ins = fig.add_axes([0.56, 0.55, 0.38, 0.32])
ax_ins.plot(x_chi, y_chi, 'b-', linewidth=1.5)
ax_ins.fill_between(x_chi, y_chi, where=(x_chi > crit), color='red', alpha=0.3, label='棄却域(5%)')
ax_ins.axvline(H_stat, color='darkred', linewidth=2, linestyle='--',
               label=f'H 統計量={H_stat:.2f}')
ax_ins.set_xlabel(f'χ²(df={df_deg})', fontsize=8)
ax_ins.set_title('Hausman 統計量の位置', fontsize=8)
ax_ins.legend(fontsize=7)
ax_ins.grid(True, alpha=0.3)

# テキスト注記
result_text = (
    f"【Hausman 検定結果】\n\n"
    f"  H 統計量  = {H_stat:.3f}\n"
    f"  自由度    = {df_deg}\n"
    f"  p 値      = {p_val:.4f}\n\n"
    f"  {'→ FE モデルを採用 **' if adopt_fe else '→ RE モデルを採用'}\n\n"
    f"【判断フロー】\n"
    f"  p < 0.05 → 個体効果と説明変数が相関\n"
    f"              → RE はバイアスあり\n"
    f"              → FE を使うべき\n\n"
    f"  p ≥ 0.05 → 相関なしの仮定OK\n"
    f"              → RE が一致かつ有効（効率的）"
)
ax2.text(5, 5, result_text, ha='center', va='center', fontsize=10,
         bbox=dict(boxstyle='round', facecolor='lightyellow', alpha=0.9, edgecolor='#F9A825'),
         family='monospace')
ax2.set_title('Hausman検定によるモデル選択', fontsize=10, y=0.98)

plt.suptitle('図3：Hausman検定 ― FE と RE どちらを使うべきか？\n'
             '被説明変数: 転入率（千人あたり）｜出典: SSDSE-B 実データ', fontsize=11, y=1.02)
plt.tight_layout()
save_fig('fig3_hausman')

print("  図4：採用モデルの係数プロット（p値・有意性付き）")

adopted_res   = fe_res if adopt_fe else re_res
adopted_name  = '固定効果モデル (FE)' if adopt_fe else '変量効果モデル (RE)'
adopted_color = 'darkorange' if adopt_fe else 'forestgreen'

fig, ax = plt.subplots(figsize=(10, 6))
plot_vars = [v for v in x_vars if v in adopted_res.params.index]
coefs = [adopted_res.params[v]       for v in plot_vars]
pvals = [adopted_res.pvalues[v]      for v in plot_vars]
se_s  = [adopted_res.std_errors[v]   for v in plot_vars]
ci_lo = [b - 1.96 * s for b, s in zip(coefs, se_s)]
ci_hi = [b + 1.96 * s for b, s in zip(coefs, se_s)]
labs  = [var_labels.get(v, v)        for v in plot_vars]

y_pos = np.arange(len(plot_vars))
bar_colors = [adopted_color if p < 0.05 else '#BDBDBD' for p in pvals]
ax.barh(y_pos, coefs, color=bar_colors, alpha=0.8, height=0.6)
ax.errorbar(coefs, y_pos,
            xerr=[np.array(coefs) - np.array(ci_lo),
                  np.array(ci_hi) - np.array(coefs)],
            fmt='none', color='black', capsize=4, linewidth=1.2)
ax.axvline(0, color='black', linewidth=1.0, linestyle='--')
ax.set_yticks(y_pos)
ax.set_yticklabels(labs, fontsize=10)
ax.set_xlabel('標準化偏回帰係数（95% CI）\n※ 有彩色: p<0.05 有意、グレー: 非有意', fontsize=10)
ax.set_title(f'図4：{adopted_name} の回帰係数（転入率モデル）\n'
             f'N={int(adopted_res.nobs)}, R²={adopted_res.rsquared:.3f}'
             f'  ｜出典: SSDSE-B 実データ（47都道府県×9年）', fontsize=11)
ax.grid(True, axis='x', alpha=0.3)

# 有意性マーク
for i, (b, se, p) in enumerate(zip(coefs, se_s, pvals)):
    sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else ''
    if sig:
        x_text = b + 1.96 * se + 0.02
        ax.text(x_text, i, sig, va='center', fontsize=10, color='#C62828', fontweight='bold')

# 凡例
sig_patch   = mpatches.Patch(color=adopted_color, alpha=0.8, label='p<0.05 有意')
nosig_patch = mpatches.Patch(color='#BDBDBD',     alpha=0.8, label='p≥0.05 非有意')
ax.legend(handles=[sig_patch, nosig_patch], fontsize=9, loc='lower right')
plt.tight_layout()
save_fig('fig4_coef')