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")


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)


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)

  • $\begingroup$ Are MSOAs perfectly nested within Sectors or are these crossed. That is, are there some MSOAs that belong to more than one Sector (or vice-versa)? $\endgroup$ – Erik Ruzek Feb 27 at 15:08
  • $\begingroup$ Yes, each row represents an individual firm located within an MSOA and working on one primary Sector. For example there are 70 firms located within E02000001 and these firms focus on 10 different Sectors. $\endgroup$ – et_ Feb 27 at 15:24

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
(msoa <- .06483/tot) #same msoa, different sector
(sector <- .3497/tot) #same sector, different msoa
(sector_msoa <- (.06483+.3497)/tot) #same sector and same msoa

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.

| cite | improve this answer | |
  • $\begingroup$ Thanks, so for the brms the ICC will be calculated as follows: (total <- .28^2 + .69^2 + (pi^2/3)) (msoa <- .28^2/total) (sector <- .69^2/total) (sector <- (.28^2 + .69^2)/total) $\endgroup$ – et_ Feb 27 at 16:32
  • $\begingroup$ That's the way! $\endgroup$ – Erik Ruzek Feb 27 at 16:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.