社会保障政策と犯罪の関係

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))

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 保存完了")

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 保存完了")

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 保存完了")

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',
}

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}")

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))

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("  → 変量効果モデルも許容される")