emmeans 在计算转换后的结果变量的置信区间时是否使用 lm_robust 聚类稳健标准误差?

Does emmeans use lm_robust cluster-robust standard errors when calculating confidence intervals for transformed outcome variables?

我正在使用 emmeans 包检查两个连续预测变量之间的相互作用。我正在使用 estimatr 包中的 lm_robust() 来执行线性回归并获得聚类稳健标准误差。结果变量居中并缩放为 SD 单位方差。例如:

fit <- lm_robust(scale(Y) ~ X1 * X2 + X3 + X4, data = mydata, cluster = school, se_type = 'CR2')

然后我可以使用类似于以下的代码在 X2 的三个级别执行成对对比或可视化线条:

emmip(fit, X2 ~ X1, CIs = TRUE, at = list(X2 = c(mean(X2) - sd(X2),
                                                 mean(X2),
                                                 mean(X2) + sd(X2))))

我不希望将结果变量反向转换为其原始比例。

我的问题是 emmeans 是否使用聚类稳健标准误差来计算其报告的置信区间或 p 值,这种行为是否取决于结果变量是否在其原始尺度上或转变? estimatr 包创建者 website suggests that lm_robust objects can be used with emmeans, but I can't see lm_robust listed as a supported model on the "Models supported by emmeans" vignette page 或包文档的简短示例。

我相信 lm_robust 对象是 lm 的扩展,因此它使用 lm 的 emmeans 支持。反过来,这意味着估计值是通过 coef(model) 获得的,而它们的 SE 是使用 vcov(model) 得出的。因此,如果 vcov() returns 您需要稳健的方差,emmeans 将使用它们。

对于大多数转换,它将按照转换小插图中的描述进行工作。特别是,指定 type = "response" 会导致对估计值和置信限进行反向转换,保留 P 值,并通过 delta 方法计算 SE(但 不会 在 CI 和测试中使用)。

附加信息

首先,我发现lm_robust并没有继承自lm;相反,estimatr 包包含了它自己对 emmeans 的支持。没有给出很多细节,但是 estimatr 的开发者必须相信所提供的必须是合适的。

scale() 转换不是内置的,因为它很复杂。只是说我们使用了 "scale" 并不像说它是 "log" 那么简单,比如说,因为要处理 scale() 结果,我们需要知道用什么来居中和划分结果.

解决方法是创建emmeans()及其亲属需要反转转换的对象;这是 stats::make.link()emmeans::make.tran() 返回的函数列表。这是一个可以达到这个目的的函数:

make.scaletran = function(y, ...) {
    sy = scale(y, ...)
    if(is.null(m <- attr(sy, "scaled:center")))
       m = 0
    if(is.null(s <- attr(sy, "scaled:scale")))
        s = 1
    list(
        linkfun = function(mu) (mu - m) / s,
        linkinv = function(eta) s * eta + m,
        mu.eta = function(eta) s,
        valideta = function(eta) TRUE,
        name = paste0("scale(", signif(m, 3), ", ", signif(s, 3), ")")
    )
}

要使用它,您需要手动指定转换,因为它不会被自动检测到。这是一个使用 R 中已有的 warpbreaks 数据的示例:

> warp.lmr = lm_robust(scale(breaks) ~ tension, cluster = wool, 
+     se_type = 'CR2', data = warpbreaks)

> tran = make.scaletran(warpbreaks$breaks)

> emmeans(warp.lmr, "tension", tran = tran)
 tension emmean    SE df lower.CL upper.CL
 L        0.624 0.619 51   -0.618   1.8666
 M       -0.133 0.181 51   -0.497   0.2301
 H       -0.491 0.219 51   -0.930  -0.0517

Results are given on the scale(28.1, 13.2) (not the response) scale. 
Confidence level used: 0.95 

> emmeans(warp.lmr, "tension", tran = tran, type = "response")
 tension response   SE df lower.CL upper.CL
 L           36.4 8.17 51     20.0     52.8
 M           26.4 2.39 51     21.6     31.2
 H           21.7 2.89 51     15.9     27.5

Confidence level used: 0.95 
Intervals are back-transformed from the scale(28.1, 13.2) scale 

OP 中 emmip() 调用的代码不正确,因为它使用 emmeans() 的规范,而不是 emmip()

我会考虑在以后的更新中 emmeans::make.tran() 添加这个比例转换选项。