emmeans 在计算转换结果变量的置信区间时是否使用 lm_robust 集群稳健标准误差？

Question

我正在使用该包检查两个连续预测变量之间的相互作用emmeans。我使用包lm_robust()中的estimatr内容来执行线性回归并获得集群稳健的标准误差。结果变量以 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包创建者网站上的一个简短示例表明lm_robust对象可以与 一起使用，但我在“emmeans 支持的模型”小插图页面或包文档中看emmeans不到将其列为受支持的模型。lm_robust

Answer 1

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

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

附加信息

首先，我发现lm_robust不继承自lm; 相反，estimatr包包含了它自己对emmeans的支持。没有给出太多细节，但是estimatr的开发人员必须相信所提供的内容必须是适当的。

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

解决方法是创建emmeans()需要反转转换的对象及其相关对象；stats::make.link()这是or返回形式的函数列表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), ")")
    )
}

要使用此功能，您需要手动指定转换，因为它不会自动检测。warpbreaks以下是使用R 中已有数据的示例：

> 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()我会考虑在未来的更新中添加此比例转换选项。