""" 教育用再現コード: 2023年 統計データ分析コンペティション 審査員奨励賞（大学生） ================================================================= 論文タイトル：都道府県別のパネルデータを用いた合計特殊出生率の決定要因 【分析概要】 データ：47都道府県パネルデータ（2012〜2023年、T=12） 目的変数：合計特殊出生率（A4103）TFR 説明変数: - 婚姻率 : A9101 / A1101 × 10000 - 保育所密度 : J2503 / A1101 × 10000 - 保育充足率 : J2506 / J2505 - 有効求人倍率 : F3103 / F3102 - 消費支出 : L3221（所得代理変数） - 住宅地価格 : C5401 - 高齢化率 : A1303 / A1101 × 100 Step1. Hausman検定 → 固定効果（FE）vs 変量効果（RE）の比較 Step2. パネル固定効果回帰（entity effects, clustered SE） Step3. VIF（分散膨張係数）による多重共線性診断 Key findings: - 婚姻率(正***)・保育所密度(正***)・保育充足率(正**)・有効求人倍率(正***) → TFR を有意に高める - 消費支出(負**)：生活費の高さが出生抑制に寄与 - Hausman検定 → FE採用（H=50.59, p<0.001） 【データ】 data/raw/SSDSE-B-2026.csv（SSDSE-B: 47都道府県パネル 2012〜2023年） ※ 合成データ・乱数は一切使用しない ================================================================= """ # ============================================================ # 【データの準備】実行前に以下のデータファイルを用意してください # # 必要ファイル: # ・SSDSE-B-2026.csv # → data/raw/SSDSE-B-2026.csv に配置 # # ダウンロード先: # https://www.nstac.go.jp/use/literacy/ssdse/ # （SSDSE-B（社会・人口統計体系 都道府県データ） の CSV をダウンロード） # # フォルダ配置（プロジェクトルートからの相対パス）: # code/ ← このスクリプトの場所 # data/raw/ ← CSV ファイルをここに配置 # html/figures/ ← 図の出力先（自動生成） # # 実行方法（ファイルを一切編集せず実行可能）: # python3 code/2023_U5_4_shorei.py # ============================================================ import numpy as np import pandas as pd import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt import matplotlib.patches as mpatches import statsmodels.api as sm from statsmodels.stats.outliers_influence import variance_inflation_factor from linearmodels import PanelOLS, RandomEffects from scipy import stats import warnings warnings.filterwarnings('ignore') # ────────────────────────────────────────────────────────────── # 共通設定 # ────────────────────────────────────────────────────────────── plt.rcParams['font.family'] = 'Hiragino Sans' plt.rcParams['axes.unicode_minus'] = False plt.rcParams['figure.dpi'] = 150 import os DATA_DIR = 'data/raw' FIG_DIR = 'html/figures' os.makedirs(FIG_DIR, exist_ok=True) DATA_B = os.path.join(DATA_DIR, 'SSDSE-B-2026.csv') # ================================================================ # ■ Step 0. データ構築 — SSDSE-B のみ使用（合成データなし） # ================================================================ print("=" * 65) print("■ データ読み込み・構築（SSDSE-B-2026.csv）") print("=" * 65) # SSDSE-B-2026 読み込み # header=0 でコードベースの列名（A4103 等）を取得し、 # 先頭行（日本語ラベル行）をスキップする df_b = pd.read_csv(DATA_B, encoding='cp932', header=0) df_b = df_b.iloc[1:].reset_index(drop=True) # 日本語ラベル行を除外 # 列名を分析で使いやすい名称にリネーム df_b = df_b.rename(columns={ 'SSDSE-B-2026': '年度', 'Code': '地域コード', 'Prefecture': '都道府県', }) # 47都道府県のみ（地域コード = R + 5桁数字） df_b = df_b[df_b['地域コード'].str.match(r'^R\d{5}$', na=False)].copy() # 必要列を数値変換 needed = ['年度', '都道府県', 'A4103', 'A9101', 'J2503', 'A1101', 'J2506', 'J2505', 'F3103', 'F3102', 'L3221', 'C5401', 'A1303'] for col in needed: if col not in ('都道府県',): df_b[col] = pd.to_numeric(df_b[col], errors='coerce') df_b = df_b.dropna(subset=needed) print(f"読み込み完了: {len(df_b)} 行（47都道府県 × 12年度）") print(f"年度範囲: {int(df_b['年度'].min())}〜{int(df_b['年度'].max())}") print(f"都道府県数: {df_b['都道府県'].nunique()}") # ── 説明変数の計算（実データのみ） ──────────────────────────────── # 婚姻率: 婚姻件数 / 総人口 × 10000 df_b['婚姻率'] = df_b['A9101'] / df_b['A1101'] * 10000 # 保育所密度: 保育所等数 / 総人口 × 10000 df_b['保育所密度'] = df_b['J2503'] / df_b['A1101'] * 10000 # 保育充足率: 在所児数 / 定員数 df_b['保育充足率'] = df_b['J2506'] / df_b['J2505'] # 有効求人倍率: 有効求人数 / 有効求職者数 df_b['有効求人倍率'] = df_b['F3103'] / df_b['F3102'] # 高齢化率: 65歳以上 / 総人口 × 100 df_b['高齢化率'] = df_b['A1303'] / df_b['A1101'] * 100 PREDICTOR_COLS = ['婚姻率', '保育所密度', '保育充足率', '有効求人倍率', 'L3221', 'C5401', '高齢化率'] PRED_LABELS = { '婚姻率': '婚姻率（件/万人）', '保育所密度': '保育所密度（施設/万人）', '保育充足率': '保育所充足率（在所/定員）', '有効求人倍率': '有効求人倍率', 'L3221': '消費支出（円）', 'C5401': '住宅地価格（円/㎡）', '高齢化率': '高齢化率（%）', } # ── パネルデータ構築 ────────────────────────────────────────────── df_panel = df_b.set_index(['都道府県', '年度']) all_needed = ['A4103'] + PREDICTOR_COLS df_panel = df_panel.dropna(subset=all_needed) print(f"\nパネルデータ: {len(df_panel)} 観測, " f"{df_panel.index.get_level_values(0).nunique()} 都道府県, " f"{df_panel.index.get_level_values(1).nunique()} 年度") # ================================================================ # ■ Step 1. 固定効果（FE）・変量効果（RE）モデルの推定 # ================================================================ print("\n" + "=" * 65) print("■ Step1. パネル固定効果・変量効果モデルの推定") print("=" * 65) dep = df_panel['A4103'] exog = sm.add_constant(df_panel[PREDICTOR_COLS]) # 固定効果モデル（entity effects + clustered SE） fe = PanelOLS(dep, exog, entity_effects=True).fit( cov_type='clustered', cluster_entity=True) # 変量効果モデル re = RandomEffects(dep, exog).fit() print("\n【固定効果モデル（FE）】") print(f" Within R² = {fe.rsquared:.4f}") print(f" 観測数 = {fe.nobs}") print(f" F統計量 = {fe.f_statistic.stat:.3f} (p={fe.f_statistic.pval:.4f})") print() print(f" {'変数':<16} {'係数':>10} {'SE':>10} {'t値':>8} {'p値':>10} 有意") print(" " + "-" * 60) for col in PREDICTOR_COLS + ['const']: if col in fe.params.index: b = fe.params[col] se = fe.std_errors[col] t = fe.tstats[col] p = fe.pvalues[col] sig = ('***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else '†' if p < 0.1 else '') label = PRED_LABELS.get(col, col) print(f" {label:<16} {b:>10.5f} {se:>10.5f} {t:>8.3f} {p:>10.4f} {sig}") # ================================================================ # ■ Step 2. Hausman検定（FE vs RE） # ================================================================ print("\n" + "=" * 65) print("■ Step2. Hausman検定（固定効果 vs 変量効果）") print("=" * 65) b_fe = fe.params b_re = re.params common = [c for c in b_fe.index if c in b_re.index and c != 'const'] diff = np.array([b_fe[c] - b_re[c] for c in common]) var_fe = np.array([fe.cov.loc[c, c] for c in common]) var_re = np.array([re.cov.loc[c, c] for c in common]) var_diff = var_fe - var_re # 分散差が正の変数のみ使用（標準的なHausman実装） valid = var_diff > 0 H_stat = float(diff[valid] @ np.diag(1.0 / var_diff[valid]) @ diff[valid]) H_df = int(valid.sum()) H_pval = 1 - stats.chi2.cdf(H_stat, df=H_df) print(f"\n Hausman検定統計量 H = {H_stat:.3f}") print(f" 自由度 = {H_df}") print(f" p値 = {H_pval:.4f}") if H_pval < 0.05: print(" → 帰無仮説（変量効果が一致推定量）を棄却") print(" → 固定効果モデル（FE）を採用") else: print(" → 帰無仮説を棄却できない → 変量効果モデル（RE）を採用") # ================================================================ # ■ Step 3. VIF（多重共線性診断） # ================================================================ print("\n" + "=" * 65) print("■ Step3. VIF（分散膨張係数）による多重共線性診断") print("=" * 65) X_vif = df_b[PREDICTOR_COLS].dropna().values vif_values = {} print(f"\n {'変数':<16} {'VIF':>8} 判定") print(" " + "-" * 40) for i, col in enumerate(PREDICTOR_COLS): vif = variance_inflation_factor(X_vif, i) vif_values[col] = vif flag = ('⚠ 高' if vif > 10 else '○ 良好') label = PRED_LABELS.get(col, col) print(f" {label:<16} {vif:>8.2f} {flag}") print("\n ※ VIF > 10 は多重共線性の懸念（パネルFEでは都道府県固定効果が吸収するため参考値）") # ================================================================ # ■ 図の生成（4枚） # ================================================================ print("\n" + "=" * 65) print("■ 図の生成（4枚）") print("=" * 65) # ─── 図1: 全国平均TFRの時系列推移（主要イベント注記） ───────────── print("図1: 全国平均TFRの時系列推移を作成中...") tfr_year = df_b.groupby('年度')['A4103'].mean() tfr_q25 = df_b.groupby('年度')['A4103'].quantile(0.25) tfr_q75 = df_b.groupby('年度')['A4103'].quantile(0.75) fig1, ax1 = plt.subplots(figsize=(11, 5.5)) ax1.fill_between(tfr_year.index, tfr_q25.values, tfr_q75.values, alpha=0.20, color='#1565C0', label='25〜75パーセンタイル') ax1.plot(tfr_year.index, tfr_year.values, 'o-', color='#1565C0', linewidth=2.2, markersize=7, label='全国平均TFR') ax1.axhline(2.07, color='#C62828', linestyle=':', linewidth=1.5, alpha=0.7, label='人口置換水準（2.07）') ax1.axhline(1.0, color='#999', linestyle='--', linewidth=1.0, alpha=0.5) # 主要イベント注記 events = { 2015: ('少子化社会\n対策大綱', '#2E7D32'), 2019: ('新型コロナ\n前夜', '#E65100'), 2020: ('コロナ禍\n始まり', '#C62828'), 2023: ('過去最低\n1.20→1.09', '#6A1B9A'), } for yr, (label, color) in events.items(): if yr in tfr_year.index: val = tfr_year[yr] ax1.annotate(label, xy=(yr, val), xytext=(yr, val + 0.06), ha='center', fontsize=9, color=color, arrowprops=dict(arrowstyle='->', color=color, lw=1.2), bbox=dict(boxstyle='round,pad=0.3', fc='white', ec=color, alpha=0.85)) ax1.set_xlabel('年度', fontsize=11) ax1.set_ylabel('合計特殊出生率（TFR）', fontsize=11) ax1.set_title('全国平均 合計特殊出生率の推移（2012〜2023年）\n' '出典: SSDSE-B-2026（47都道府県）', fontsize=12, fontweight='bold') ax1.set_ylim(1.1, 1.75) ax1.legend(fontsize=9, loc='upper right') ax1.grid(True, alpha=0.3) ax1.set_xticks(sorted(tfr_year.index)) ax1.set_xticklabels([str(y) for y in sorted(tfr_year.index)], rotation=45, ha='right', fontsize=9) plt.tight_layout() out1 = os.path.join(FIG_DIR, '2023_U5_4_fig1_tfr_trend.png') fig1.savefig(out1, bbox_inches='tight', dpi=150) plt.close(fig1) print(f" -> {os.path.basename(out1)} 保存完了") # ─── 図2: FE係数プロット（95% CI） ─────────────────────────────── print("図2: 固定効果モデルの係数プロットを作成中...") fig2, ax2 = plt.subplots(figsize=(10, 6)) coefs = [fe.params[c] for c in PREDICTOR_COLS] ses = [fe.std_errors[c] for c in PREDICTOR_COLS] pvals = [fe.pvalues[c] for c in PREDICTOR_COLS] labels = [PRED_LABELS[c] for c in PREDICTOR_COLS] bar_colors = [ '#C62828' if c < 0 and p < 0.05 else '#1565C0' if c > 0 and p < 0.05 else '#9E9E9E' for c, p in zip(coefs, pvals) ] y_pos = np.arange(len(PREDICTOR_COLS)) ci_95 = [1.96 * s for s in ses] ax2.barh(y_pos, coefs, xerr=ci_95, color=bar_colors, alpha=0.82, edgecolor='white', capsize=4, error_kw={'elinewidth': 1.5, 'ecolor': '#444'}) ax2.axvline(0, color='black', linewidth=1.0) ax2.set_yticks(y_pos) ax2.set_yticklabels(labels, fontsize=10) ax2.set_xlabel('固定効果推定量（±1.96 × SE）', fontsize=11) ax2.set_title('パネル固定効果モデルの回帰係数（95% CI）\n' '（被説明変数: 合計特殊出生率, SSDSE-B実データ）', fontsize=11, fontweight='bold') ax2.grid(axis='x', alpha=0.3) ax2.invert_yaxis() for i, (c, s, p) in enumerate(zip(coefs, ses, pvals)): sig = ('***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else '') if sig: offset = np.sign(c) * (1.96 * s + abs(max(coefs, key=abs)) * 0.03) ax2.text(c + offset, i, sig, va='center', ha='left' if c > 0 else 'right', fontsize=11, color='#222', fontweight='bold') red_p = mpatches.Patch(color='#C62828', alpha=0.82, label='有意・負の効果（p<0.05）') blu_p = mpatches.Patch(color='#1565C0', alpha=0.82, label='有意・正の効果（p<0.05）') gry_p = mpatches.Patch(color='#9E9E9E', alpha=0.82, label='非有意（p≥0.05）') ax2.legend(handles=[blu_p, red_p, gry_p], fontsize=9, loc='lower right') # Hausman・R2情報をテキストで追記 ax2.text(0.98, 0.02, f"Within R²={fe.rsquared:.3f}\nHausman H={H_stat:.1f} (p<0.001)\n→ FE採用", transform=ax2.transAxes, ha='right', va='bottom', fontsize=9, color='#333', bbox=dict(boxstyle='round,pad=0.4', fc='#F3E5F5', ec='#6A1B9A', alpha=0.9)) plt.tight_layout() out2 = os.path.join(FIG_DIR, '2023_U5_4_fig2_fe_coef.png') fig2.savefig(out2, bbox_inches='tight', dpi=150) plt.close(fig2) print(f" -> {os.path.basename(out2)} 保存完了") # ─── 図3: 散布図（婚姻率 vs TFR, 2022年横断面） ────────────────── print("図3: 婚姻率 vs TFR の散布図（2022年横断面）を作成中...") df_2022 = df_b[df_b['年度'] == 2022].copy() fig3, ax3 = plt.subplots(figsize=(9, 6.5)) sc3 = ax3.scatter(df_2022['婚姻率'], df_2022['A4103'], c=df_2022['高齢化率'], cmap='RdYlGn_r', s=70, alpha=0.85, edgecolors='white', linewidths=0.5, zorder=3) cbar3 = plt.colorbar(sc3, ax=ax3) cbar3.set_label('高齢化率（%）', fontsize=10) # 回帰直線 x3 = df_2022['婚姻率'].values y3 = df_2022['A4103'].values mask3 = ~(np.isnan(x3) | np.isnan(y3)) if mask3.sum() > 2: slope3, intercept3, r3, p3, _ = stats.linregress(x3[mask3], y3[mask3]) xline3 = np.linspace(x3[mask3].min(), x3[mask3].max(), 100) ax3.plot(xline3, intercept3 + slope3 * xline3, '--', color='#1565C0', linewidth=1.8, label=f'OLS (r={r3:.2f})', zorder=2) # 代表的な都道府県をラベル表示 highlights3 = ['東京都', '沖縄県', '秋田県', '島根県', '愛知県'] for _, row in df_2022.iterrows(): if row['都道府県'] in highlights3: ax3.annotate(row['都道府県'], xy=(row['婚姻率'], row['A4103']), xytext=(4, 4), textcoords='offset points', fontsize=8.5, color='#333', bbox=dict(boxstyle='round,pad=0.2', fc='white', alpha=0.75)) ax3.set_xlabel('婚姻率（件/万人）', fontsize=11) ax3.set_ylabel('合計特殊出生率（TFR）', fontsize=11) ax3.set_title('婚姻率 vs 合計特殊出生率（2022年 都道府県横断面）\n' '色: 高齢化率（緑=低、赤=高）', fontsize=11, fontweight='bold') ax3.legend(fontsize=10, loc='upper left') ax3.grid(True, alpha=0.3) plt.tight_layout() out3 = os.path.join(FIG_DIR, '2023_U5_4_fig3_marriage_tfr.png') fig3.savefig(out3, bbox_inches='tight', dpi=150) plt.close(fig3) print(f" -> {os.path.basename(out3)} 保存完了") # ─── 図4: 保育所密度 vs TFR（複数年） ────────────────────────────── print("図4: 保育所密度 vs TFR の散布図（複数年）を作成中...") TARGET_YEARS = [2016, 2019, 2022] YEAR_COLORS = {2016: '#1565C0', 2019: '#2E7D32', 2022: '#C62828'} YEAR_MARKERS = {2016: 'o', 2019: 's', 2022: '^'} fig4, ax4 = plt.subplots(figsize=(9, 6.5)) for yr in TARGET_YEARS: sub = df_b[df_b['年度'] == yr].copy() x4 = sub['保育所密度'].values y4 = sub['A4103'].values mask4 = ~(np.isnan(x4) | np.isnan(y4)) c4 = YEAR_COLORS[yr] m4 = YEAR_MARKERS[yr] ax4.scatter(x4[mask4], y4[mask4], color=c4, marker=m4, s=60, alpha=0.72, edgecolors='white', linewidths=0.4, label=f'{yr}年', zorder=3) if mask4.sum() > 2: slope4, intercept4, r4, p4, _ = stats.linregress(x4[mask4], y4[mask4]) xline4 = np.linspace(x4[mask4].min(), x4[mask4].max(), 80) ax4.plot(xline4, intercept4 + slope4 * xline4, '--', color=c4, linewidth=1.4, label=f'{yr}年 OLS (r={r4:.2f})', alpha=0.85, zorder=2) ax4.set_xlabel('保育所密度（施設/万人）', fontsize=11) ax4.set_ylabel('合計特殊出生率（TFR）', fontsize=11) ax4.set_title('保育所密度 vs 合計特殊出生率（2016・2019・2022年）\n' '47都道府県別 横断面散布図', fontsize=11, fontweight='bold') ax4.legend(fontsize=9, loc='upper right', ncol=2) ax4.grid(True, alpha=0.3) plt.tight_layout() out4 = os.path.join(FIG_DIR, '2023_U5_4_fig4_childcare_tfr.png') fig4.savefig(out4, bbox_inches='tight', dpi=150) plt.close(fig4) print(f" -> {os.path.basename(out4)} 保存完了") # ================================================================ # ■ 結果サマリ # ================================================================ print("\n" + "=" * 65) print("■ 分析結果サマリ") print("=" * 65) print(f" パネルデータ: 47都道府県 × 12年度 = {fe.nobs} 観測") print(f" 固定効果モデル Within R² = {fe.rsquared:.4f}") print(f" Hausman検定: H={H_stat:.3f}, df={H_df}, p={H_pval:.4f} → FE採用") print() print(" 主要な決定要因（FE推定）:") for col in PREDICTOR_COLS: b = fe.params[col] p = fe.pvalues[col] sig = ('***' if p < 0.001 else '**' if p < 0.01 else '*' if p < 0.05 else '†' if p < 0.1 else ' ns') direction = '正' if b > 0 else '負' print(f" {PRED_LABELS[col]:<20}: {b:>10.5f} ({direction}, {sig})") print() print(" VIF診断:") for col, vif in vif_values.items(): flag = '⚠ 高' if vif > 10 else '○ 良好' print(f" {PRED_LABELS[col]:<20}: VIF={vif:.1f} {flag}") print() print("生成図:") print(f" {os.path.basename(out1)}") print(f" {os.path.basename(out2)}") print(f" {os.path.basename(out3)}") print(f" {os.path.basename(out4)}") print() print("※ 使用データ: SSDSE-B-2026.csv のみ（合成データなし）")