R中lme4包的混合效应模型的鲁棒标准误差

Question

我正在使用lme4包进行线性混合效应建模。

下面是混合效应模型：

fm01 <- lmer(sublat <- goal + (1|userid))

上述命令返回一个名为fm01的S4对象

此对象包括固定效应及其OLS标准误差 (如下所示)

Fixed effects:

            Estimate Std. Error t value
(Intercept)   31.644      3.320   9.530
goaltypeF1    -4.075      3.243  -1.257
goaltypeF2    -9.187      5.609  -1.638
goaltypeF3   -13.935      9.455  -1.474
goaltypeF4   -20.219      8.196  -2.467
goaltypeF5   -12.134      8.797  -1.379"

然而，我需要提供强健的标准误差。

对于像lme4返回的S4对象这样的对象，我该如何做到这一点？

- Adrienne

---

类似于这个的编程相关内容？ - Paulo E. Cardoso

确切地说，这适用于混合效应回归。不幸的是，vcovHC(model, type="HC0")对这些模型输出无效。您可以获得vcov(model)，但无法获得vcovHC(model)。 - Adrienne

你（或他人）需要查看 vignette("sandwich-OOP",package="sandwich") 并从 sandwich:::estfun.lm 和 sandwich:::bread.lm 开始适应必要的调整，以编写 estfun.merMod 和 bread.merMod 函数。 - Ben Bolker

merDeriv包和clubSandwich包可以帮助提取三明治健壮标准误差和假设检验的一些组件。 - J.D.

Answer 1

看起来，merDeriv和clubSandwich包提供了lmerMod对象的强大SEs:

library(lme4)
library(clubSandwich)
m <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy)

使用merDeriv

(来自merDeriv JSS论文的复制材料，感谢@AchimZeileis提供提示)

library(merDeriv)
sand <- sandwich(m, bread = bread(m, full = TRUE),
         mean = meat(m, level = 2))

clubSandwich

（所有可能的类型：我不知道哪种在任何情况下都是“最佳”的）

cstypes <- paste0("CR", c("0", "1", "1p", "1S", "2", "3"))
rob_se_fun <- function(type) sqrt(diag(vcovCR(m, type = type)))
rob_se <- sapply(cstypes, rob_se_fun)

合并结果

std_se <- sqrt(diag(vcov(m)))
cbind(std = std_se, rob_se,
      merDeriv = sqrt(diag(sand)[1:2]))\

                 std      CR0      CR1     CR1p     CR1S      CR2      CR3
(Intercept) 6.824597 6.632277 6.824557 7.034592 6.843700 6.824557 7.022411
Days        1.545790 1.502237 1.545789 1.593363 1.550125 1.545789 1.590604
            merDeriv
(Intercept) 6.632277
Days        1.502237

merDeriv 的结果与 type="CR0" 匹配（merDeriv 提供了所有组件的稳健的 Wald 估计，包括随机效应参数；由您决定是否随机效应参数的 Wald 估计足够可靠）

- Ben Bolker

---

在JSS的merDeriv论文的复制材料中（https://doi.org/10.18637/jss.v087.c01），作者展示了如何将他们的软件包与`sandwich`相结合。您必须明确指出`bread()`和`estfun()`的具体类型：`sandwich(m, bread = bread(m, full = TRUE), meat = meat(m, level = 2))`。 - Achim Zeileis

Answer 2

我认为这就是你所需要的: https://cran.r-project.org/web/packages/robustlmm/vignettes/rlmer.pdf 它是一个名为robustlmm的包，其中包含rlmer函数。

“对象的结构和方法的实现尽可能与 lme4 的相似，并在需要时具有特定的健壮性扩展。”

fm01_rob <- rlmer(sublat <- goal + (1|userid))

- R me matey

---

这里的“健壮”一词与OP所指的含义不同... - Ben Bolker