CODEWITHSUNDEEP
rlme4standard-errorrandom-effects
Matthias Herrmann
Matthias Herrmann

Reputation: 81

Cluster-Robust Standard Errors for Lmer and Glmer in Stargazer (lme4 package)

I have an experimental data set in which subjects were assigned to a specific treatment. Each treatment consisted of 5 groups. I want to estimate a model that that includes random effects on subject level and then cluster the standard errors by the assigned group.

Does anyone know how to get stargazer to display clustered SEs on group level for i) lmer and ii) glmer models?

A similar question was asked some time ago for plm models

individual random effects model with standard errors clustered on a different variable in R (R-project)

Cluster-robust errors for a plm with clustering at different level as fixed effects

Upvotes: 2

Views: 3004

Answers (1)

Vincent
Vincent

Reputation: 17853

There are two main problems.

First, I do not think that stargazer supports this. And if the feature is not currently supported, chances are good that it will not be supported in the near term. That package has not received a major update in several years.

Second, the two main packages to compute robust-cluster standard errors are sandwich and clubSandwich. sandwich does not support lme4 models. clubSandwich supports lmer models but not glmer.

This means that you can get “half” of what you want by if you are willing to consider a more “modern” alternative to stargazer, such as the modelsummary package. (Disclaimer: I am the author.)

(Note that in the example below, I am using version 1.0.1 of the package. A small bug was introduced in 1.0.0 which slowed down mixed-effects models. It is fixed in the latest version.)

install.packages("modelsummary", type = "source")

In this model, I print standard errors clustered by the cyl variable for the linear model:

library(modelsummary)
library(lme4)

mod1 <- lmer(mpg ~ hp + drat + (1 | cyl), data = mtcars)
mod2 <- glmer(am ~ hp + drat + (1 | cyl), data = mtcars, family = binomial)
#> boundary (singular) fit: see help('isSingular')

varcov1 <- vcovCR(mod1, cluster = mtcars$cyl, type = "CR0")
varcov2 <- vcov(mod2)

# converting the variance-covariance matrices to "standard" matrices
varcov1 <- as.matrix(varcov1)
varcov2 <- as.matrix(varcov2)

modelsummary(
    list(mod1, mod2),
    vcov = list(varcov1, varcov2))
Model 1 Model 2
(Intercept) 12.790 -29.076
(5.104) (12.418)
hp -0.045 0.011
(0.012) (0.009)
drat 3.851 7.310
(1.305) (3.047)
SD (Intercept cyl) 1.756 0.000
SD (Observations) 3.016 1.000
Num.Obs. 32 32
R2 Marg. 0.616 0.801
R2 Cond. 0.713
AIC 175.5 28.1
BIC 182.8 34.0
ICC 0.3
RMSE 2.81 0.33

Upvotes: 3

Related Questions