1

我正在使用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

4

1 回答 1

4

我相信 thar lm_robust 对象是 lm 的扩展，因此它使用 emmeans 对 lm 的支持。反过来，这意味着估计是通过 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), ")")
    )
}

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

> 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()在未来的更新中添加这个比例转换选项。

于 2020-05-29T22:18:45.070 回答