robust standard errors for mixed-effects models in lme4 package

Question

I am using the lme4 package for linear mixed effect modeling

The mixed-effect model is:

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

The above command returns an S4 object called fm01.

This object includes fixed effects and their OLS standard errors (below)

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"

However, I need to provide robust standard errors.

How can I do this with an S4 object such as returned by lme4?

Horne

something like this? – Transcript

Excatly, but for a mixed-effects regression. Unfortunately, vcovHC(model, type="HC0") does not work on those model outputs. You can obtain vcov(model) but you cannot obtain vcovHC(model). – Horne

you (or someone) would need to look at vignette("sandwich-OOP",package="sandwich") and figure out how to write estfun.merMod and bread.merMod functions, starting from sandwich:::estfun.lm and sandwich:::bread.lm and adapting as necessary. – Hearst

merDeriv package and clubSandwich package would do some help in extracting some components for sandwich robust standard errors and hypothesis test. – Psychographer

Answer 1

It looks like robust SEs for lmerMod objects are available via the merDeriv and clubSandwich packages:

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

using merDeriv

(From the replication materials of the merDeriv JSS paper, thanks to @AchimZeileis for the tip)

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

clubSandwich

(all possible types: I don't know enough to know which is 'best' in any given case)

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)

combine results

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's results match type="CR0" (merDeriv provides robust Wald estimates for all components, including the random effect parameters; it's up to you to decide if Wald estimates for RE parameters are reliable enough)

Hearst

In the replication materials of the merDeriv paper in JSS (doi.org/10.18637/jss.v087.c01), the authors show how to combine their package with sandwich. You have to indicate the specific flavor of bread() and estfun() explicitly: sandwich(m, bread = bread(m, full = TRUE), meat = meat(m, level = 2)) – Soll

Answer 2

I think this is what you're looking for: https://cran.r-project.org/web/packages/robustlmm/vignettes/rlmer.pdf

It's the robustlmm package, which has the rlmer function.

"The structure of the objects and the methods are implemented to be as similar as possible to the ones of lme4 with robustness specific extensions where needed."

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

Surtax

This is a different meaning of "robust" than the OP means ... – Hearst