使用R的marginaleffects包绘制动态处理效应的问题

核心思路 #

要让avg_comparisons()或plot_comparisons()的结果与模型中treated*periods交互项完全匹配，关键是确保工具正确计算每个时间点内处理组与控制组的平均差异，同时与模型的固定效应设定对齐。

步骤1：构建正确的双向固定效应DiD模型 #

首先确保你的模型是标准的双向固定效应设定（个体+时间固定效应），且periods为因子型（保证交互项对应每个时间点的处理效应）。这里以fixest包为例（DiD分析常用，效率更高）：

# 加载依赖包
library(marginaleffects)
library(fixest)
library(tibble)

# 模拟符合要求的DiD数据（2组、7个时间点）
set.seed(123)
dat <- tibble(
  id = rep(1:40, each = 7),          # 40个个体，每个7个时间点
  periods = factor(rep(1:7, 40)),    # 时间点转为因子型，确保每个点对应单独交互项
  treated = factor(rep(c(0, 1), each = 20*7)), # 处理组/对照组标记
  y = rnorm(280, 
            mean = 1 + 0.5*as.numeric(treated) + 
              0.3*as.numeric(periods) + 
              c(0, 0, 0, 0.8, 1.2, 1.5, 1.8)*as.numeric(treated), # 前3期无效应，后4期递增
            sd = 1)
)

# 拟合双向固定效应DiD模型
# 公式说明：主效应包含treated*periods交互项，固定效应为个体(id)+时间(periods)
model <- feols(y ~ treated*periods | id + periods, data = dat)
summary(model)

步骤2：用avg_comparisons()匹配交互项结果 #

通过指定variables = "treated"和by = "periods"，让工具按时间点分组计算处理组与控制组的平均差异，这会直接对应模型中treated*periods的交互项系数：

# 计算每个时间点的动态处理效应
dyn_effects <- avg_comparisons(
  model,
  variables = list(treated = c("0", "1")), # 指定比较：控制组(0) vs 处理组(1)
  by = "periods",                          # 按时间点分组计算
  re.form = ~0                             # 适配固定效应模型，排除随机效应（fixest可省略，但保留更稳妥）
)

# 查看结果：estimate列即为每个时间点的处理效应，与模型交互项系数完全匹配
print(dyn_effects)

步骤3：用plot_comparisons()绘制动态效应图 #

直接基于上述逻辑绘制可视化结果：

# 绘制动态处理效应折线图
plot_comparisons(
  model,
  variables = "treated",
  by = "periods",
  re.form = ~0
) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "darkred") +
  labs(x = "时间点", y = "处理效应（处理组 - 控制组）", title = "双向固定效应DiD动态处理效应") +
  theme_minimal()

常见问题排查（如果结果仍不匹配） #

1. 检查periods类型：必须是因子型，若为数值型，模型交互项会变成线性趋势，而非每个时间点的单独效应，导致结果无法匹配。
2. 核对参考组：确保treated的参考组（模型中默认是0）与avg_comparisons()中指定的比较组一致。
3. 模型固定效应重复：不要同时在主效应和固定效应中加入periods（比如y ~ treated*periods + factor(periods) | id），这会引发多重共线性，导致交互项系数失真。
4. 包版本问题：确保marginaleffects和fixest是最新版本，旧版本可能存在固定效应处理的bug。

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