1

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

Matthias Herrmann
  • 71
  • 7

1 Answers1

1

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
Vincent
  • 15,809
  • 7
  • 37
  • 39
  • Hello Vincent, I implemented your solution for my dataset, however I get the following error when running vcovCR: "Non-nested random effects detected. clubSandwich methods are not available for such models." The problem is that my random effects are for the individual subject while the clustering of SEs are on group level. Any idea how to solve this issue? – Matthias Herrmann Jul 14 '22 at 10:35
  • No, sorry. Maybe you'll find something in the `clubSandwich` documentation: https://cran.r-project.org/web/packages/clubSandwich/clubSandwich.pdf – Vincent Jul 14 '22 at 10:50
  • Alright, thanks anyway! – Matthias Herrmann Jul 14 '22 at 10:58