如何用R的marginaleffects包复现Stata的margins at命令解释交互效应

问题场景 #

我需要在R里复现Stata中基于平均边际效应解释交互效应的操作，以mtcars数据集为例：

• 模型目标：分析马力（hp）对每加仑英里数（mpg）的影响如何随气缸数（cyl）变化
• 先在R中构建交互模型：

mod1 <- lm(mpg ~ hp:as.factor(cyl), data=mtcars)

• 对应Stata中计算不同cyl取值下hp平均边际效应的代码：

margins, dydx(hp) at(cyl == 4) at(cyl == 6) at(cyl == 8)

尝试实现与疑问 #

我用marginaleffects包的avg_slopes函数搭配by参数做了尝试：

library(marginaleffects)

avg_slopes(mod1, variables=c("hp"), by=c("cyl"))

得到输出：

Term    Contrast cyl Estimate Std. Error      z Pr(>|z|)    S   2.5 %   97.5 %
   hp mean(dY/dX)   4  0.00835     0.0347  0.241   0.8097  0.3 -0.0596  0.07633
   hp mean(dY/dX)   6 -0.04492     0.0248 -1.810   0.0704  3.8 -0.0936  0.00373
   hp mean(dY/dX)   8 -0.04736     0.0135 -3.496   <0.001 11.0 -0.0739 -0.02081

但不确定这里的by参数是否和Stata的at参数作用一致，也没太看明白by参数的文档说明。

解答 #

你的用法完全正确，avg_slopes的by参数和Stata的at参数在这个场景下功能等价：

• Stata的at(cyl == x)是固定cyl为某一取值，在对应子样本中计算hp对mpg的平均边际效应
• R中avg_slopes(..., by="cyl")会自动按cyl的不同取值对数据分组，在每个分组内计算hp的平均边际效应，本质逻辑和Stata的操作完全一致

从输出结果也能验证：三个cyl取值对应的边际效应估计值，就是各自子样本中hp每变动1单位时mpg的平均变化量，和Stataat参数的输出逻辑完全匹配。

额外补充：如果你的模型包含主效应（比如mpg ~ hp + as.factor(cyl) + hp:as.factor(cyl)），by参数的作用依然是按cyl分组计算边际效应，结果仍能和Stata的at参数输出对应。

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