""" 教育用再現コード: 2022年 統計データ分析コンペティション 優秀賞（大学生・一般の部） ================================================================================= 論文タイトル：社会保障政策と犯罪の関係 ―都道府県パネルデータによる実証分析― 手法：パネルデータ分析（固定効果・変量効果）、Hausman検定 【分析概要】 犯罪統計（刑法犯認知件数）はSSDSE-Bに含まれないため、 「求職圧力指数（月間有効求職者数/総人口）」を経済的犯罪リスクの代理変数として使用。 社会保障指標（保健医療費・保育所等数・一般病院数）がこの経済的脆弱性指標に どう関連するかをPanelOLS（固定効果モデル）で実証分析する。 - 目的変数: 求職圧力指数 = 月間有効求職者数（一般）/ 総人口 - 説明変数: 保健医療費（二人以上の世帯）、保育所等数（人口1万対）、 一般病院数（人口1万対）、65歳以上人口割合 - モデル: PanelOLS（entity_effects=True）、clustered標準誤差 【使用データ（実データ）】 data/raw/SSDSE-B-2026.csv — 47都道府県 パネル 2012-2023 出典: 政府統計の総合窓口（e-Stat）、SSDSE（教育用標準データセット） 【出力図】 html/figures/2022_U2_fig1_corr.png — 相関ヒートマップ（主要変数） html/figures/2022_U2_fig2_ts.png — 時系列: 求職圧力指数と保健医療費の推移 html/figures/2022_U2_fig3_fe_coef.png — 固定効果係数プロット html/figures/2022_U2_fig4_scatter.png — 散布図: 求職圧力指数 vs 社会保障指標 ================================================================================= """ # ============================================================ # 【データの準備】実行前に以下のデータファイルを用意してください # # 必要ファイル: # ・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/2022_U2_yushu.py # ============================================================ import os import warnings 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 scipy import stats warnings.filterwarnings('ignore') # ────────────────────────────────────────────────────────────── # 共通設定 # ────────────────────────────────────────────────────────────── plt.rcParams['font.family'] = 'Hiragino Sans' plt.rcParams['axes.unicode_minus'] = False plt.rcParams['figure.dpi'] = 150 FIGURE_DIR = os.path.normpath('html/figures') DATA_DIR = os.path.normpath('data/raw') os.makedirs(FIGURE_DIR, exist_ok=True) COLORS = { 'primary': '#1565C0', 'secondary': '#E65100', 'success': '#2E7D32', 'danger': '#C62828', 'purple': '#6A1B9A', 'teal': '#00695C', 'gray': '#616161', } # ────────────────────────────────────────────────────────────── # データ読み込み（SSDSE-B-2026） # ────────────────────────────────────────────────────────────── print("=" * 70) print("■ データ読み込み: SSDSE-B-2026.csv（47都道府県、2012-2023年）") print("=" * 70) DATA_B = os.path.join(DATA_DIR, 'SSDSE-B-2026.csv') 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() # 必要カラムを数値変換 num_cols = [ '総人口', '65歳以上人口', '月間有効求職者数（一般）', '月間有効求人数（一般）', '保健医療費（二人以上の世帯）', '消費支出（二人以上の世帯）', '保育所等数', '保育所等定員数', '一般病院数', ] for col in num_cols: df_b[col] = pd.to_numeric(df_b[col], errors='coerce') # ────────────────────────────────────────────────────────────── # 変数作成 # ────────────────────────────────────────────────────────────── # 求職圧力指数（月間有効求職者数 / 総人口） df_b['求職圧力指数'] = df_b['月間有効求職者数（一般）'] / df_b['総人口'] # 保育所等数（人口1万人対） df_b['保育所等数_万対'] = df_b['保育所等数'] / df_b['総人口'] * 10000 # 一般病院数（人口1万人対） df_b['一般病院数_万対'] = df_b['一般病院数'] / df_b['総人口'] * 10000 # 65歳以上人口割合 df_b['高齢化率'] = df_b['65歳以上人口'] / df_b['総人口'] # 保健医療費（万円単位に正規化、スケール調整） df_b['保健医療費_千円'] = df_b['保健医療費（二人以上の世帯）'] / 1000 # 分析対象変数リスト VARS = ['求職圧力指数', '保健医療費_千円', '保育所等数_万対', '一般病院数_万対', '高齢化率'] # 欠損値除去 df_panel = df_b[['年度', '地域コード', '都道府県'] + VARS].dropna().copy() df_panel = df_panel.sort_values(['地域コード', '年度']).reset_index(drop=True) print(f"分析対象: {df_panel['地域コード'].nunique()}都道府県 × {df_panel['年度'].nunique()}年") print(f"年度: {sorted(df_panel['年度'].unique())}") print(f"総観測数: {len(df_panel)}") print() print(df_panel[VARS].describe().round(4)) # ────────────────────────────────────────────────────────────── # PanelOLS 固定効果モデル # ────────────────────────────────────────────────────────────── print("\n" + "=" * 70) print("■ PanelOLS 固定効果モデル（entity_effects=True）") print("=" * 70) from linearmodels import PanelOLS # マルチインデックス設定（地域コード × 年度） df_fe = df_panel.set_index(['地域コード', '年度']) # 目的変数・説明変数 y_fe = df_fe['求職圧力指数'] X_vars = ['保健医療費_千円', '保育所等数_万対', '一般病院数_万対', '高齢化率'] X_fe = sm.add_constant(df_fe[X_vars]) # 固定効果モデル推定 model_fe = PanelOLS(y_fe, X_fe, entity_effects=True) result_fe = model_fe.fit(cov_type='clustered', cluster_entity=True) print("\n【固定効果モデル 推定結果】") print(result_fe.summary) # 係数・標準誤差・p値の抽出 fe_params = result_fe.params fe_stderr = result_fe.std_errors fe_pvalues = result_fe.pvalues fe_tstat = result_fe.tstats print("\n【係数サマリー（定数項を除く）】") for var in X_vars: sig = ('***' if fe_pvalues[var] < 0.01 else '**' if fe_pvalues[var] < 0.05 else '*' if fe_pvalues[var] < 0.10 else 'n.s.') print(f" {var:22s}: β={fe_params[var]:+.6f} SE={fe_stderr[var]:.6f}" f" t={fe_tstat[var]:+.3f} p={fe_pvalues[var]:.4f} {sig}") # ────────────────────────────────────────────────────────────── # 変量効果モデル & Hausman 検定（手動実装） # ────────────────────────────────────────────────────────────── print("\n" + "=" * 70) print("■ 変量効果モデル & Hausman 検定") print("=" * 70) from linearmodels import RandomEffects model_re = RandomEffects(y_fe, X_fe) result_re = model_re.fit(cov_type='robust') print("\n【変量効果モデル 推定結果（係数のみ）】") for var in X_vars: sig = ('***' if result_re.pvalues[var] < 0.01 else '**' if result_re.pvalues[var] < 0.05 else '*' if result_re.pvalues[var] < 0.10 else 'n.s.') print(f" {var:22s}: β={result_re.params[var]:+.6f} " f"p={result_re.pvalues[var]:.4f} {sig}") # Hausman 検定（手動） # H₀: 変量効果モデルが一致推定量（ランダム効果と固定効果の差がゼロ） # H₁: 固定効果モデルが適切 b_fe = fe_params[X_vars].values b_re = result_re.params[X_vars].values # 分散共分散行列の差 cov_fe = result_fe.cov.loc[X_vars, X_vars].values cov_re = result_re.cov.loc[X_vars, X_vars].values cov_diff = cov_fe - cov_re # 固有値が負になりうるため、擬似逆行列を使用 try: cov_diff_inv = np.linalg.pinv(cov_diff) diff = b_fe - b_re hausman_stat = float(diff @ cov_diff_inv @ diff) hausman_df = len(X_vars) hausman_pval = 1 - stats.chi2.cdf(hausman_stat, df=hausman_df) except np.linalg.LinAlgError: hausman_stat = float('nan') hausman_df = len(X_vars) hausman_pval = float('nan') print(f"\n【Hausman 検定】") print(f" H₀: 変量効果モデルが適切（FE = RE）") print(f" H₁: 固定効果モデルが適切（FE ≠ RE）") print(f" χ²統計量 = {hausman_stat:.4f} df = {hausman_df}") print(f" p値 = {hausman_pval:.4f}") if hausman_pval < 0.05: print(" → 固定効果モデルが支持される（p < 0.05）") else: print(" → 変量効果モデルが支持される（p ≥ 0.05）") # ────────────────────────────────────────────────────────────── # 全国年度平均（時系列プロット用） # ────────────────────────────────────────────────────────────── annual = df_panel.groupby('年度')[['求職圧力指数', '保健医療費_千円', '保育所等数_万対']].mean() annual = annual.sort_index() print("\n" + "=" * 70) print("■ 全国年度平均（主要変数）") print("=" * 70) print(annual.round(5)) # ================================================================ # 図1：相関ヒートマップ（主要変数） # ================================================================ print("\n図1: 相関ヒートマップを作成中...") CORR_VARS = ['求職圧力指数', '保健医療費_千円', '保育所等数_万対', '一般病院数_万対', '高齢化率'] CORR_LABELS = [ '求職圧力\n指数', '保健医療費\n（千円）', '保育所等数\n（万対）', '一般病院数\n（万対）', '高齢化率', ] corr_mat = df_panel[CORR_VARS].corr() n_vars = len(CORR_VARS) fig1, ax1 = plt.subplots(figsize=(8, 6.5)) im = ax1.imshow(corr_mat.values, cmap='RdBu_r', vmin=-1, vmax=1, aspect='auto') plt.colorbar(im, ax=ax1, shrink=0.8, label='相関係数 r') ax1.set_xticks(range(n_vars)) ax1.set_yticks(range(n_vars)) ax1.set_xticklabels(CORR_LABELS, fontsize=10) ax1.set_yticklabels(CORR_LABELS, fontsize=10) ax1.tick_params(axis='x', labelrotation=0) # 相関係数の数値をセルに表示 for i in range(n_vars): for j in range(n_vars): r = corr_mat.values[i, j] # 有意性検定 if i != j: x_arr = df_panel[CORR_VARS[j]].dropna().values y_arr = df_panel[CORR_VARS[i]].dropna().values mask = ~(np.isnan(x_arr) | np.isnan(y_arr)) _, pval = stats.pearsonr(x_arr[mask], y_arr[mask]) sig_mark = '***' if pval < 0.001 else '**' if pval < 0.01 else '*' if pval < 0.05 else '' else: sig_mark = '' text_color = 'white' if abs(r) > 0.6 else 'black' ax1.text(j, i, f'{r:.2f}{sig_mark}', ha='center', va='center', fontsize=9.5, color=text_color, fontweight='bold') ax1.set_title('主要変数の相関ヒートマップ\n（47都道府県 × 2012-2023年、プールドデータ）\n' '※ *p<0.05、**p<0.01、***p<0.001', fontsize=11, fontweight='bold', pad=14) ax1.set_xlabel('データ出典: SSDSE-B-2026（e-Stat）', fontsize=9, color='gray') # 枠線 for spine in ax1.spines.values(): spine.set_visible(False) ax1.set_xticks(np.arange(-0.5, n_vars, 1), minor=True) ax1.set_yticks(np.arange(-0.5, n_vars, 1), minor=True) ax1.grid(which='minor', color='white', linewidth=2) ax1.tick_params(which='minor', bottom=False, left=False) plt.tight_layout() fig1.savefig(os.path.join(FIGURE_DIR, '2022_U2_fig1_corr.png'), bbox_inches='tight', dpi=150) plt.close(fig1) print(" → 2022_U2_fig1_corr.png 保存完了") # ================================================================ # 図2：時系列推移（求職圧力指数 × 保健医療費） # ================================================================ print("図2: 時系列推移プロットを作成中...") years_ts = annual.index.values fig2, ax2a = plt.subplots(figsize=(10, 5)) ax2b = ax2a.twinx() l1 = ax2a.plot(years_ts, annual['求職圧力指数'] * 1000, color=COLORS['primary'], linewidth=2.5, marker='o', markersize=7, markerfacecolor='white', markeredgewidth=2, label='求職圧力指数（×1000）') l2 = ax2b.plot(years_ts, annual['保健医療費_千円'], color=COLORS['secondary'], linewidth=2.5, marker='s', markersize=7, markerfacecolor='white', markeredgewidth=2, linestyle='--', label='保健医療費（千円/世帯）') ax2a.set_xlabel('年度', fontsize=12) ax2a.set_ylabel('求職圧力指数（×1000）', fontsize=12, color=COLORS['primary']) ax2b.set_ylabel('保健医療費（千円/二人以上世帯）', fontsize=12, color=COLORS['secondary']) ax2a.tick_params(axis='y', labelcolor=COLORS['primary']) ax2b.tick_params(axis='y', labelcolor=COLORS['secondary']) ax2a.set_title('全国平均 求職圧力指数と保健医療費の推移（2012-2023年）\n' 'データ出典: SSDSE-B-2026（e-Stat）', fontsize=12, fontweight='bold') ax2a.set_xlim(2011.5, 2023.5) ax2a.set_xticks(years_ts) ax2a.grid(True, alpha=0.3) # 注釈: COVID-19 if 2020 in annual.index: ax2a.axvspan(2019.5, 2021.5, alpha=0.08, color='red', label='COVID-19期') ax2a.annotate('COVID-19\n（2020-21年）', xy=(2020, annual.loc[2020, '求職圧力指数'] * 1000), xytext=(2020.8, annual.loc[2020, '求職圧力指数'] * 1000 * 1.05), arrowprops=dict(arrowstyle='->', color='gray', lw=1.2), fontsize=8.5, color='gray') lines = l1 + l2 labels = [l.get_label() for l in lines] ax2a.legend(lines, labels, loc='upper right', fontsize=10) plt.tight_layout() fig2.savefig(os.path.join(FIGURE_DIR, '2022_U2_fig2_ts.png'), bbox_inches='tight', dpi=150) plt.close(fig2) print(" → 2022_U2_fig2_ts.png 保存完了") # ================================================================ # 図3：固定効果係数プロット # ================================================================ print("図3: 固定効果係数プロットを作成中...") # 表示用変数名マッピング var_label_map = { '保健医療費_千円': '保健医療費\n（千円/世帯）', '保育所等数_万対': '保育所等数\n（人口1万対）', '一般病院数_万対': '一般病院数\n（人口1万対）', '高齢化率': '高齢化率', 'const': '定数項', } # CI（95%） ci_lo = fe_params - 1.96 * fe_stderr ci_hi = fe_params + 1.96 * fe_stderr # 定数項を除いて表示 plot_vars = X_vars # ['保健医療費_千円', '保育所等数_万対', '一般病院数_万対', '高齢化率'] plot_labels = [var_label_map[v] for v in plot_vars] plot_coefs = [fe_params[v] for v in plot_vars] plot_cilo = [ci_lo[v] for v in plot_vars] plot_cihi = [ci_hi[v] for v in plot_vars] plot_pvals = [fe_pvalues[v] for v in plot_vars] # p値に応じた色 bar_colors = [] for p in plot_pvals: if p < 0.01: bar_colors.append(COLORS['danger']) elif p < 0.05: bar_colors.append(COLORS['secondary']) elif p < 0.10: bar_colors.append(COLORS['success']) else: bar_colors.append(COLORS['gray']) fig3, ax3 = plt.subplots(figsize=(9, 5.5)) y_pos = np.arange(len(plot_vars)) ax3.barh(y_pos, plot_coefs, color=bar_colors, alpha=0.85, edgecolor='white', linewidth=0.8, height=0.5) # 95% CI エラーバー for i, (lo, hi, coef) in enumerate(zip(plot_cilo, plot_cihi, plot_coefs)): ax3.plot([lo, hi], [i, i], color='black', linewidth=2.0, solid_capstyle='round') ax3.plot(lo, i, '|', color='black', markersize=10, markeredgewidth=2) ax3.plot(hi, i, '|', color='black', markersize=10, markeredgewidth=2) ax3.axvline(0, color='black', linewidth=1.0, linestyle='--', alpha=0.7) ax3.set_yticks(y_pos) ax3.set_yticklabels(plot_labels, fontsize=11) ax3.set_xlabel('推定係数（β）', fontsize=12) ax3.set_title('固定効果モデル（PanelOLS）の推定係数\n' '目的変数: 求職圧力指数、エラーバー: 95%信頼区間', fontsize=12, fontweight='bold') ax3.grid(axis='x', alpha=0.3) ax3.invert_yaxis() # 係数値・有意水準のラベル for i, (coef, p) in enumerate(zip(plot_coefs, plot_pvals)): sig = ('***' if p < 0.01 else '**' if p < 0.05 else '*' if p < 0.10 else 'n.s.') x_text = max(plot_cihi[i], abs(coef)) * 1.05 + 0.00001 ax3.text(plot_cihi[i] + abs(max(plot_cihi) - min(plot_cilo)) * 0.03, i, f'{coef:+.5f} {sig}', va='center', fontsize=9, fontweight='bold') # 凡例 legend_patches = [ mpatches.Patch(color=COLORS['danger'], label='p < 0.01 ***'), mpatches.Patch(color=COLORS['secondary'], label='p < 0.05 **'), mpatches.Patch(color=COLORS['success'], label='p < 0.10 *'), mpatches.Patch(color=COLORS['gray'], label='n.s.（p ≥ 0.10）'), ] ax3.legend(handles=legend_patches, loc='lower right', fontsize=9) # R² など r2_text = (f"R²（within）= {result_fe.rsquared:.3f}\n" f"観測数 N = {int(result_fe.nobs)}\n" f"クラスタード標準誤差（都道府県）") ax3.text(0.02, 0.97, r2_text, transform=ax3.transAxes, fontsize=9, va='top', color='gray', bbox=dict(boxstyle='round', facecolor='#F8F9FA', alpha=0.8)) plt.tight_layout() fig3.savefig(os.path.join(FIGURE_DIR, '2022_U2_fig3_fe_coef.png'), bbox_inches='tight', dpi=150) plt.close(fig3) print(" → 2022_U2_fig3_fe_coef.png 保存完了") # ================================================================ # 図4：散布図（求職圧力指数 vs 社会保障指標） # ================================================================ print("図4: 散布図を作成中...") # 3つの代表年度 years_sc = [2014, 2018, 2022] df_sc = df_panel[df_panel['年度'].isin(years_sc)].copy() scatter_pairs = [ ('保健医療費_千円', '保健医療費（千円/世帯）'), ('保育所等数_万対', '保育所等数（人口1万対）'), ] fig4, axes4 = plt.subplots(2, 3, figsize=(15, 9)) fig4.suptitle('求職圧力指数 vs 社会保障指標（散布図・複数年比較）\n' 'データ出典: SSDSE-B-2026（e-Stat）', fontsize=13, fontweight='bold', y=1.01) scatter_colors_yr = { 2014: COLORS['primary'], 2018: COLORS['secondary'], 2022: COLORS['success'], } for row_i, (xvar, xlabel) in enumerate(scatter_pairs): for col_i, yr in enumerate(years_sc): ax = axes4[row_i, col_i] df_yr = df_sc[df_sc['年度'] == yr].dropna(subset=['求職圧力指数', xvar]) x = df_yr[xvar].values y = df_yr['求職圧力指数'].values * 1000 # ×1000 で見やすく col = scatter_colors_yr[yr] ax.scatter(x, y, color=col, alpha=0.7, s=55, edgecolors='white', linewidth=0.7, zorder=3) # 回帰直線 if len(x) > 3: X_fit = sm.add_constant(x) fit_yr = sm.OLS(y, X_fit).fit() x_line = np.linspace(x.min(), x.max(), 100) y_line = fit_yr.params[0] + fit_yr.params[1] * x_line ax.plot(x_line, y_line, color=col, linewidth=2.0, linestyle='-', alpha=0.9, zorder=2) r2 = fit_yr.rsquared pv = fit_yr.pvalues[1] sig = '***' if pv < 0.01 else '**' if pv < 0.05 else '*' if pv < 0.10 else 'n.s.' stat_text = f'β={fit_yr.params[1]:.4f}{sig}\nR²={r2:.3f}' else: stat_text = 'データ不足' ax.set_xlabel(xlabel, fontsize=9.5) ax.set_ylabel('求職圧力指数（×1000）', fontsize=9.5) ax.set_title(f'{yr}年度\n{stat_text}', fontsize=10.5, fontweight='bold') ax.grid(True, alpha=0.3) # 上位・下位都道府県にラベル df_yr_sorted = df_yr.assign(y_plot=y).sort_values('y_plot') for _, row in df_yr_sorted.tail(3).iterrows(): ax.annotate(row['都道府県'], (row[xvar], row['求職圧力指数'] * 1000), fontsize=7, alpha=0.8, xytext=(3, 3), textcoords='offset points') for _, row in df_yr_sorted.head(2).iterrows(): ax.annotate(row['都道府県'], (row[xvar], row['求職圧力指数'] * 1000), fontsize=7, alpha=0.8, xytext=(3, -9), textcoords='offset points') plt.tight_layout() fig4.savefig(os.path.join(FIGURE_DIR, '2022_U2_fig4_scatter.png'), bbox_inches='tight', dpi=150) plt.close(fig4) print(" → 2022_U2_fig4_scatter.png 保存完了") # ────────────────────────────────────────────────────────────── # 最終サマリー # ────────────────────────────────────────────────────────────── print("\n" + "=" * 70) print("✓ 全図の生成完了") print("=" * 70) print("\n【主要結果サマリー】") print(f" 固定効果モデル R²（within）= {result_fe.rsquared:.4f}") print(f" 観測数: {int(result_fe.nobs)}") print() for var in X_vars: sig = ('***' if fe_pvalues[var] < 0.01 else '**' if fe_pvalues[var] < 0.05 else '*' if fe_pvalues[var] < 0.10 else 'n.s.') print(f" [{sig}] {var}: β={fe_params[var]:+.6f} p={fe_pvalues[var]:.4f}") print() print(f" Hausman検定: χ²={hausman_stat:.4f} df={hausman_df} p={hausman_pval:.4f}") if hausman_pval < 0.05: print(" → 固定効果モデル採用を支持") else: print(" → 変量効果モデルも許容される")