如何在R中绘制双重差分（DID）安慰剂检验图？

核心思路 #

安慰剂检验的逻辑是：多次随机重新分配「处理组」标签，每次拟合DID模型并提取交互项（treatment*after）的系数，最终把这些安慰剂系数的分布可视化，再将真实模型的交互项系数标注在图中，以此判断真实效应是否显著异于随机分配产生的效应。

具体步骤&代码 #

1. 加载依赖包

我们需要用tidyverse做数据处理和绘图，broom提取模型结果：

library(tidyverse)
library(broom)

2. 定义安慰剂模拟函数

写一个函数，每次随机打乱原数据的treatment变量，拟合DID模型并提取交互项系数：

placebo_sim <- function(data) {
  # 随机重分配处理组
  data_sim <- data %>% 
    mutate(treatment = sample(treatment, n(), replace = FALSE))
  
  # 拟合DID模型
  model <- lm(score ~ treatment + after + treatment*after, data = data_sim)
  
  # 提取交互项系数
  tidy(model) %>% 
    filter(term == "treatment:after") %>% 
    pull(estimate)
}

3. 运行多次模拟

这里我们跑1000次模拟（次数可按需调整），得到所有安慰剂系数：

set.seed(123) # 设置随机种子保证结果可重复
placebo_coefs <- replicate(1000, placebo_sim(data))

4. 提取真实模型的交互项系数

从你已经拟合好的did模型中提取真实的交互项系数：

real_coef <- tidy(did) %>% 
  filter(term == "treatment:after") %>% 
  pull(estimate)

5. 绘制安慰剂检验图

用密度图+直方图结合的方式展示，同时用竖线标出真实系数：

tibble(placebo_coef = placebo_coefs) %>% 
  ggplot(aes(x = placebo_coef)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30, fill = "lightblue", alpha = 0.7) +
  geom_density(color = "darkblue", linewidth = 1) +
  geom_vline(xintercept = real_coef, color = "red", linewidth = 1.5, linetype = "dashed") +
  labs(
    x = "安慰剂交互项系数",
    y = "密度",
    title = "DID安慰剂检验结果",
    subtitle = paste0("真实交互项系数: ", round(real_coef, 3))
  ) +
  theme_minimal()

结果解读 #

如果红色虚线（真实系数）明显落在安慰剂系数分布的尾部（比如95%置信区间之外），说明你的真实处理效应不太可能是随机产生的，DID估计结果更可靠。

内容的提问来源于stack exchange，提问作者Ajern