Different ICC between gmler and brm

Question

I am trying to estimate the same multilevel model with a categorical (dummy) dependent variable using glmer() from lme4 and brm() from brms, but I receive very different ICC: $0.112$ for glmer() and $0.00$ for brm(). For constancy with the other models I am estimating for the same project I have a preference towards brm(), but the $0.00$ ICC cannot justify the use of multilevel modeling. Could you explain the difference in ICC and suggest if it is meaningful to still use brm()?

Reproducible example:

loading data

data <- read.csv("https://www.dropbox.com/s/jmid86jkgzy9tlh/sample.csv?dl=1")

glmer

model.glmer <- glmer(SocialNetworkD ~ 1 + (1|MSOA) + (1|Sector), data=data, family="binomial", control=glmerControl(optimizer="bobyqa"), na.action = "na.exclude") summary(model.glmer) performance::icc(model.glmer)

brm

model.brm <- brm(SocialNetworkD~1+(1|MSOA) + (1|Sector), data=data, family="bernoulli"(link = "logit"), seed = 123, iter = 2000) summary(model.brm) performance::icc(model.brm, ppd = T)

Answer 1

I think the problem is not with glmer or with brms. Instead it is with performance::icc. Having used this before, it is a handy function but seems to break down once you have more complicated structures as yours. You will need to calculate your ICCs the old fashioned way, and with a cross-classified model, you have a few different ICCs.

We can just use the glmer output for this, although you can do the same with brms, just make sure to square the standard deviation estimates of the random effects:

> summary(model.glmer) 
Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
 Family: binomial  ( logit )
Formula: SocialNetworkD ~ 1 + (1 | MSOA) + (1 | Sector)
   Data: dataICC
Control: glmerControl(optimizer = "bobyqa")

     AIC      BIC   logLik deviance df.resid 
  5200.3   5219.2  -2597.2   5194.3     4043 

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.1506 -1.0012  0.5633  0.7938  1.6456 

Random effects:
 Groups Name        Variance Std.Dev.
 MSOA   (Intercept) 0.06483  0.2546  
 Sector (Intercept) 0.34970  0.5914  
Number of obs: 4046, groups:  MSOA, 2339; Sector, 13

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept)   0.3567     0.1685   2.117   0.0343 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

You can get three ICCs from this output:

tot <- .06483+.3497+(pi^2/3) # total variance
3.704398
(msoa <- .06483/tot) #same msoa, different sector
0.01750082
(sector <- .3497/tot) #same sector, different msoa
0.0944013
(sector_msoa <- (.06483+.3497)/tot) #same sector and same msoa
0.1119021

Based on this, performance::ICC is giving you something like sector_msoa, which is the similarity of observations from the same sector and msoa. Although I'm not positive about that given I'm not totally sure what assumption performance::ICC makes about the level 1 residual in a binomial/Bernoulli GLMM. I am went with the latent formuation for the variance at level 1, $\pi^2/3$ because otherwise calculating an ICC is difficult. For more information, see StasK's answer to this question.

Note that the glmer variances are somewhat similar to the estimates you get from brms if you convert them to a variance with brms giving a larger estimate for sector. Obviously they are different in that brms is reporting on a distribution of variances, but if you look at the mean estimate, it's in the ballpark:

Group-Level Effects: 
~MSOA (Number of levels: 2339) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.28      0.11     0.04     0.48 1.01      424      552

~Sector (Number of levels: 13) 
              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept)     0.69      0.17     0.45     1.09 1.00     1296     1819

BTW, I saw a similar pattern of estimates of the variance when using MCMC methods from rstanarm and R2MLwiN. The MCMC procedures tend to give greater variation to sector than maximum likelihood.