Rで差分の差分法

library(tidyverse)

# パラメータ設定
n_per_group <- 100 # 各グループの個体数 (サンプルサイズ)
time_points <- 2 # 時間点 (0: 介入前, 1: 介入後)
baseline_mean <- 10 # 介入前の平均アウトカム（対照群）
time_trend_effect <- 2 # 時間経過による共通の効果（両群に影響）
treatment_effect <- 6.5 # 介入（処置）の効果（処置群の介入後にのみ影響）
sd_noise <- 2 # 個体差や測定誤差などのノイズの標準偏差

# サンプルデータを格納するデータフレームを作成
df <- expand.grid(
  id = 1:(n_per_group * 2), # 個体ID
  time = 0:(time_points - 1) # 時間 (0: 前, 1: 後)
) %>%
  as_tibble() %>%
  arrange(id, time)

# 処置群か対照群かを割り当てる (idが101以上を処置群とする)
df <- df %>%
  mutate(
    treatment_group = ifelse(id > n_per_group, 1, 0), # 0: 対照群, 1: 処置群
    # time_post = time # timeが0/1なので time_post = time と同じ。
  )

head(df)

seed <- 20250420
set.seed(seed)

# アウトカム変数の生成
# 介入前のベースラインのアウトカムを生成（個体差を含む）
baseline_outcome <- rnorm(n_per_group * 2, mean = baseline_mean, sd = sd_noise)
df <- df %>%
  mutate(
    # 各個体のベースラインアウトカムを対応させる
    baseline_y = baseline_outcome[id],
    # 時間経過効果と処置効果、ノイズを加味してアウトカムを計算
    outcome = baseline_y + # ベースライン
      time * time_trend_effect + # 時間トレンド効果 (両群共通)
      treatment_group * time * treatment_effect + # 処置効果 (処置群のtime=1のみ)
      rnorm(n(), mean = 0, sd = sd_noise) # 各時点での追加ノイズ
  ) %>%
  # 不要な中間変数を削除
  select(-baseline_y)

head(df)
str(df)

(summary_df <- df %>%
  group_by(treatment_group, time) %>%
  summarise(mean_outcome = mean(outcome), .groups = "drop")) # .groups='drop'でgroup化を解除しています

result_did <- function() {
  mean_control_pre <- summary_df %>%
    filter(treatment_group == 0, time == 0) %>%
    pull(mean_outcome)
  mean_control_post <- summary_df %>%
    filter(treatment_group == 0, time == 1) %>%
    pull(mean_outcome)
  diff_control <- mean_control_post - mean_control_pre
  mean_treatment_pre <- summary_df %>%
    filter(treatment_group == 1, time == 0) %>%
    pull(mean_outcome)
  mean_treatment_post <- summary_df %>%
    filter(treatment_group == 1, time == 1) %>%
    pull(mean_outcome)
  diff_pre <- mean_treatment_pre - mean_control_pre
  diff_treatment <- mean_treatment_post - mean_treatment_pre
  did_estimate <- diff_treatment - diff_control
  list("対照群の介入前平均" = mean_control_pre, "対照群の介入後平均" = mean_control_post, "対照群の変化量" = diff_control, "処置群の介入前平均" = mean_treatment_pre, "処置群の介入後平均" = mean_treatment_post, "処置群の変化量" = diff_treatment, "介入前における、処置群と対照群のアウトカム平均の差" = diff_pre, "DID推定値 (処置群の変化量 - 対照群の変化量)" = did_estimate, "データ生成時の処置効果(再掲)" = treatment_effect)
}
result_did()

# 回帰モデルの実行
# time も treatment_group も 0/1 のダミー変数として投入
# 相互作用項 treatment_group:time の係数が DID 推定値となる
did_model <- lm(outcome ~ treatment_group + time + treatment_group:time, data = df)
# did_model <- lm(outcome ~ treatment_group * time, data = df)
summary(did_model)
did_estimate_regression <- coef(did_model)["treatment_group:time"]
list("回帰モデルから得られたDID推定値 (treatment_group:time の係数)" = did_estimate_regression)