#!/usr/bin/env python3 # coding: utf-8 # ============================================================ # 【データの準備】実行前に以下のデータファイルを用意してください # # 必要ファイル: # ・SSDSE-B-2026.csv # → data/raw/SSDSE-B-2026.csv に配置 # ・SSDSE-E-2026.csv # → data/raw/SSDSE-E-2026.csv に配置 # # ダウンロード先: # https://www.nstac.go.jp/use/literacy/ssdse/ # （SSDSE-B（社会・人口統計体系 都道府県データ） の CSV をダウンロード） # （SSDSE-E（社会・人口統計体系 都道府県の指標2） の CSV をダウンロード） # # フォルダ配置（プロジェクトルートからの相対パス）: # code/ ← このスクリプトの場所 # data/raw/ ← CSV ファイルをここに配置 # html/figures/ ← 図の出力先（自動生成） # # 実行方法（ファイルを一切編集せず実行可能）: # python3 code/2025_U1_daijin.py # ============================================================ """ ========================================================================== 2025年度 統計データ分析コンペティション 総務大臣賞（大学生・一般の部） 「何が女性議員比率を左右するのか―東京特別区議会を対象にしたパネルデータ分析―」 著者：永渕真緒（津田塾大学総合政策学部総合政策学科） 【このスクリプトの目的】 論文の手法をトレースすることで、パネルデータ分析の基礎を学ぶ。 - Pooled OLS（プーリング回帰） - 固定効果モデル (FE: Fixed Effects) - 変量効果モデル (RE: Random Effects) - Hausman検定によるモデル選択 【データ】 data/2025_U1/2025_U1_panel.csv（2025_U1_data_prep.py で生成） SSDSE-B/E から構築した 47都道府県×2014〜2022年の実データパネル（423行） ※東京特別区データは非公開のため、同手法を都道府県パネルで再現 出典: 情報・システム研究機構「教育用標準データセット SSDSE-B・E」 目的変数: 転入率（千人あたり）＝「人口移動を左右するもの」 【必要ライブラリ】 pip install pandas numpy matplotlib scipy statsmodels linearmodels ========================================================================== """ 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}') # ========================================================================== # 1. データ読み込み # ========================================================================== 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() # ========================================================================== # 2. 記述統計の表示 # ========================================================================== 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() # ========================================================================== # 3. データの可視化 # ========================================================================== print("=" * 60) print("【Step 3】データの可視化") print("=" * 60) # ---- 図1：転入率の時系列推移（都道府県別） ---- 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') # ---- 図2：固定効果の可視化（OLS/FE/RE 係数比較の前置き） ---- # パネル分析準備 df_panel = df.copy() df_panel['pref_id'] = pd.Categorical(df_panel[ENTITY_COL]).codes + 1 df_panel = df_panel.set_index([ENTITY_COL, TIME_COL]) x_vars = [ '高齢化率', '年少人口比率', '合計特殊出生率', '婚姻率', '保育所等数', '年平均気温', '教育費', '人口密度', ] def standardize(df_in, cols): """変数を標準化（平均0, 標準偏差1）する関数 標準化により単位の異なる変数を比較可能にする。 係数 = 1標準偏差の変化に対する転入率の変化量（標準偏差） """ df_out = df_in.copy() for col in cols: mu = df_out[col].mean() sig = df_out[col].std() df_out[col] = (df_out[col] - mu) / sig return df_out all_vars = [Y_VAR] + x_vars df_sub = df_panel[all_vars].dropna() df_std = standardize(df_sub, x_vars) y = df_std[Y_VAR] X = df_std[x_vars] # Pooled OLS X_ols = sm.add_constant(X) ols_res = sm.OLS(y, X_ols).fit(cov_type='HC3') # Fixed Effects fe_res = PanelOLS(y, X, entity_effects=True, time_effects=False, drop_absorbed=True).fit(cov_type='robust') # Random Effects re_res = RandomEffects(y, X).fit(cov_type='robust') print() print("=" * 60) print("【Step 4】3つのパネルデータモデルの推定") print("=" * 60) print(f" OLS: N={int(ols_res.nobs)}, R²={ols_res.rsquared:.3f}") print(f" FE: N={int(fe_res.nobs)}, R²={fe_res.rsquared:.3f}") print(f" RE: N={int(re_res.nobs)}, R²={re_res.rsquared:.3f}") # ---- 図2: 3モデルの係数比較プロット ---- 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') # ========================================================================== # Hausman 検定 # ========================================================================== 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() # ---- 図3：Hausman 検定の概念図 + 検定結果 ---- 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') # ---- 図4：採用モデルの係数プロット（有意性マーク付き） ---- 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') # ========================================================================== # 結果サマリー # ========================================================================== print() print("=" * 60) print("【Step 6】主要結果のサマリー") print("=" * 60) print() print(f" ─ {adopted_name}（採用モデル）の主要結果 ─") print(f" {'変数名':<25} {'係数':>8} {'p値':>8} {'有意'}") print(" " + "─" * 55) for var in plot_vars: b = adopted_res.params[var] p = adopted_res.pvalues[var] sig = '***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else '.' if p < 0.1 else '' lbl = var_labels.get(var, var)[:25] print(f" {lbl:<25} {b:>8.4f} {p:>8.4f} {sig}") print() print(f" Hausman検定: χ²={H_stat:.2f}, df={df_deg}, p={p_val:.4f}") print(f" → {'FEモデルが統計的に適切' if adopt_fe else 'REモデルが統計的に適切（FEも参考として掲載）'}") print() # ========================================================================== # 学習ポイント # ========================================================================== print("=" * 60) print("【学習ポイント】この論文から学べるデータサイエンスの技術") print("=" * 60) print(""" 1. パネルデータの構造 ・同じ個体（都道府県）の繰り返し観察 → 時系列と横断面を同時活用 ・set_index(['都道府県', '年度']) で MultiIndex を作成 2. 変数の標準化 ・単位が異なる変数（率・温度・件数など）を比較可能にする ・(x - mean) / std で標準化 → 係数 = 1SD変化あたりの効果 3. 3つのモデルの使い分け ・Pooled OLS：シンプルだが不観測の地域差を無視 ・FE（固定効果）：地域固有の不観測要因を統制 ・RE（変量効果）：FEとOLSの中間。個体効果が無相関の場合に有効 4. Hausman 検定 ・FE と RE どちらが適切か統計的に検証 ・χ² 検定: p < 0.05 なら FE を採用 5. ロバスト標準誤差 ・誤差項の不均一分散・自己相関に対処（HC3 / cluster-robust） 6. 実データの活用 ・SSDSE-B：47都道府県×12年の多変数パネル（全て実測値） ・SSDSE-E：都道府県別の面積・県民所得など断面データ ・論文手法を実データで再現 → 教育的再現性の確保 """) print("=" * 60) print("すべての図を保存しました。HTML ファイルで確認してください。") print(" fig1_trend.png : 転入率の時系列推移（47都道府県）") print(" fig2_fe.png : OLS/FE/RE 係数比較プロット") print(" fig3_hausman.png : Hausman検定の概念図と検定結果") print(" fig4_coef.png : 採用モデルの係数（有意性付き）") print("=" * 60)