Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

robust standard errors for mixed-effects models in lme4 package of R

Tags:

r

mixed-models

I am using the lme4 package for linear mixed effect modeling

the mixed-effect model is below:

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?

like image 588
Adrienne Avatar asked Oct 16 '14 19:10

Adrienne


2 Answers

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)

like image 168
Ben Bolker Avatar answered Nov 10 '22 11:11

Ben Bolker


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))
like image 25
R me matey Avatar answered Nov 10 '22 13:11

R me matey