Rで統計的因果推論:自然実験-差分の差分法

# 必要なパッケージを読み込み
library(tidyverse)
library(knitr)

# 結果を再現可能にするための乱数シード設定
seed <- 20250825
set.seed(seed)

# シミュレーションのパラメータ設定
n <- 500 # 各市のサンプルサイズ
city_a_effect <- 0.5 # A市が元々持つ読書時間の優位性（時間）
time_trend <- -0.2 # 時間経過による自然な変化（時間）
treatment_effect <- 1.5 # 図書館建設の真の因果効果（時間）

# データ生成
sim_data <- tibble(
  # 子供のID
  child_id = 1:(n * 2 * 2),
  # 市（A市が介入群、B市が対照群）
  city = rep(c("A市（介入群）", "B市（対照群）"), each = n * 2),
  # 時期（介入前：2021年、介入後：2023年）
  time = rep(rep(c("介入前 (2021年)", "介入後 (2023年)"), each = n), 2)
) |>
  # 介入ダミー変数と時期ダミー変数を追加
  mutate(
    # 介入群（A市）なら1、対照群（B市）なら0
    treated = if_else(city == "A市（介入群）", 1, 0),
    # 介入後なら1、介入前なら0
    post = if_else(time == "介入後 (2023年)", 1, 0)
  ) |>
  # 結果変数（週平均読書時間）を生成
  mutate(
    reading_time = 5 + # 全体のベースライン読書時間
      city_a_effect * treated + # A市の元々の特性
      time_trend * post + # 時間経過によるトレンド
      treatment_effect * treated * post + # 介入(treated=1)かつ介入後(post=1)の因果効果
      rnorm(n * 2 * 2, mean = 0, sd = 1.5), # 個人のばらつき
    .keep = "all"
  ) |>
  # 時期を因子型に変換（プロットの順序のため）
  mutate(time = factor(time, levels = c("介入前 (2021年)", "介入後 (2023年)")))

# 生成したデータの最初の数行を表示
cat("##### 生成されたデータの一部 #####\n")
kable(head(sim_data))

# グループごとの平均読書時間を計算
summary_stats <- sim_data |>
  group_by(city, time) |>
  summarise(
    mean_reading_time = mean(reading_time),
    .groups = "drop"
  )

cat("##### グループごとの平均読書時間 #####\n")
kable(summary_stats, digits = 2)

# DID推定値の計算
# 1. 各市の介入前後の平均値を取得
y_a_pre <- summary_stats$mean_reading_time[summary_stats$city == "A市（介入群）" & summary_stats$time == "介入前 (2021年)"]
y_a_post <- summary_stats$mean_reading_time[summary_stats$city == "A市（介入群）" & summary_stats$time == "介入後 (2023年)"]
y_b_pre <- summary_stats$mean_reading_time[summary_stats$city == "B市（対照群）" & summary_stats$time == "介入前 (2021年)"]
y_b_post <- summary_stats$mean_reading_time[summary_stats$city == "B市（対照群）" & summary_stats$time == "介入後 (2023年)"]

# 2. 差分1（時間変化）の計算
diff_a <- y_a_post - y_a_pre
diff_b <- y_b_post - y_b_pre

cat("\n##### 差分1：各市の時間変化 #####\n")
cat("A市（介入群）の変化量:", round(diff_a, 3), "時間\n")
cat("B市（対照群）の変化量:", round(diff_b, 3), "時間\n")

# 3. 差分2（差分の差）の計算
did_estimate <- diff_a - diff_b

cat("\n##### 差分2：差分の差（DID推定値） #####\n")
cat(
  "DID推定値 (A市の変化量 - B市の変化量):",
  round(did_estimate, 3),
  "時間\n"
)
cat(
  "この値が、図書館建設による読書時間への因果効果の推定値です。\n"
)