在R中复现Stata的margins：获取treat两水平的AME及差值与95%CI

你想要复现Stata中regress bwt age i.smoke后执行margins smoke得到各分组平均预测值的功能，但当前R代码仅输出treat单一水平的平均边际效应（AME），以下是你的问题的解决方案：

问题1：显示treat两个水平的AME及95%置信区间 #

首先明确：Stata的margins smoke输出的是各分组的平均预测值，而非边际效应；你当前用margins(mod, variables = "treat")得到的是Treat相对于参考水平Control的边际效应（即从Control转为Treat时score的平均变化）。

方案1：获取各水平的平均预测值（对应Stata margins smoke功能） #

使用at参数指定treat的所有水平，直接计算每个分组的平均预测值及置信区间：

# 计算treat每个水平的平均预测值
margins_pred <- margins(mod, at = list(treat = c("Treat", "Control")))
# 查看结果（包含95%置信区间）
summary(margins_pred)

方案2：获取各水平的平均边际效应（AME） #

对于二分类因子，参考水平（默认是Control）的AME为0（无基准变化），若要显式展示两个水平的结果，推荐用emmeans包实现：

library(emmeans)
# 提取treat的边际效应及置信区间
emmeans(mod, ~ treat)

问题2：计算两个水平的均值差及95%置信区间 #

方法1：直接用emmeans包（最简便） #

emmeans可以直接计算分组均值差及置信区间：

# 计算Treat与Control的均值差及95%置信区间
emmeans(mod, ~ treat) %>% contrast(method = "pairwise")

方法2：基于margins的预测值计算 #

如果坚持用margins包，可从之前的预测结果中手动计算差值：

# 提取预测结果
preds <- summary(margins_pred)
# 计算均值差、标准误及95%置信区间
diff_result <- data.frame(
  estimate = preds$estimate[preds$at == "treat=Treat"] - preds$estimate[preds$at == "treat=Control"],
  se = sqrt(preds$std.error[preds$at == "treat=Treat"]^2 + preds$std.error[preds$at == "treat=Control"]^2),
  ci_low = preds$estimate[preds$at == "treat=Treat"] - 1.96*preds$std.error[preds$at == "treat=Treat"],
  ci_high = preds$estimate[preds$at == "treat=Treat"] + 1.96*preds$std.error[preds$at == "treat=Treat"]
)
print(diff_result)

另外补充：你最初用summary(margins(mod, variables = "treat"))得到的结果，本质就是Treat相对于Control的均值差（AME），其输出中的置信区间就是该差值的95%置信区间。

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