margins包计算含交互项logistic回归边际效应问题求解

问题背景 #

构建三个解释变量均为分类变量的Logistic回归模型，部分变量为二分类，其余为3个及以上水平的多分类变量：

• 模型1：基准模型，仅纳入所有变量主效应
• 模型2：纳入多组两向交互项
• 模型3：在模型2基础上额外纳入三向交互项

模型构建代码：

model1 <- glm(y ~ x + z + a + b + c + d, data=df,
              family=binomial(link="logit"), na.action=na.omit)

model2 <- glm(y ~ x*d + z*d + a*d + b*d + c*d, data=df,
              family=binomial(link="logit"), na.action=na.omit)

model3 <- glm(y ~ x*d + z*d + a*d + b*d + c*d + x*z*d, data=df, 
              family=binomial(link="logit"), na.action=na.omit)

使用margins包计算边际效应时，基准模型输出正常，但含交互项的模型2、模型3未输出交互项对应的边际效应，仅显示和基准模型一致的主变量系数。计算代码如下：

library(margins)

mod0 <- margins(model1)
mod1 <- margins(model2)
mod2 <- margins(model3)

library(modelsummary)
modelsummary(list(mod0, mod1, mod3))

输出结果截图：

尝试使用at参数调整时存在两个问题：一是部分系数正负方向和原始交互项回归系数符号不一致，二是结果缺少p值与置信区间；尝试dydx参数时无法明确指定所有交互项，尤其是三向交互项的计算设置。

核心原因与正确实现 #

为什么没有单独的交互项边际效应？ #

margins默认计算的是平均边际效应（AME），衡量的是某一个解释变量变动1单位时，因变量预测概率的绝对变动幅度。交互项本身不存在独立的“边际效应”——交互项的统计含义是某变量的边际效应会随其他变量的取值变化而产生异质性，因此不会像回归系数表那样单独输出一行交互项的效应值，你需要计算的是「调节变量取不同值时，核心变量的边际效应」，这才是交互项对应的实际可解释结果。

具体操作步骤 #

1. 两向交互项的边际效应计算
   不需要单独指定交互项，通过at参数设定调节变量的所有分类取值水平，即可输出核心变量在调节变量不同分组下的边际效应，结果自动包含标准误、p值和置信区间。以模型2中x*d、z*d等所有和d交互的项为例，代码如下：

   # 先确保分类变量为因子型，若为数值存储请提前转换：df$d <- as.factor(df$d)
   me_model2 <- margins(model2, at = list(d = levels(df$d)))
   # 查看完整结果，含p值、置信区间
   summary(me_model2)
   

2. 三向交互项的边际效应计算
   三向交互的统计含义是「两个变量的交互效应大小，会随第三个变量的取值变化而变动」，只需要在at参数中同时指定两个调节变量的所有分类水平组合即可。以模型3中x*z*d三向交互为例，代码如下：

   me_model3 <- margins(model3, at = list(z = levels(df$z), d = levels(df$d)))
   summary(me_model3)
   

常见问题说明 #

• 边际效应符号和原始回归系数符号不一致是正常现象：原始glm输出的是对数发生比（log-odds）尺度的系数，margins输出的是预测概率尺度的边际效应，二者量纲完全不同，符号、显著性出现差异不属于计算错误。
• 结果缺少p值、置信区间的解决方法：直接调用margins()返回的对象仅存储点估计值，对边际效应对象调用summary()方法才会输出标准误、检验统计量、p值和置信区间。modelsummary可以自动识别这些统计量，直接传入边际效应对象即可生成规范结果表：

  modelsummary(list(
    "基准模型" = margins(model1),
    "两向交互模型" = margins(model2, at = list(d = levels(df$d))),
    "三向交互模型" = margins(model3, at = list(z = levels(df$z), d = levels(df$d)))
  ), statistic = "p.value", conf.int = TRUE)
  

注意：如果分类变量以数值型格式存储，必须先通过as.factor()转换为因子类型，否则at参数无法正确识别所有分类水平，会导致计算结果偏差。

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