3

I'm using lme4 in R, and I have a model set up that uses a three-level hierarchy for a negative binomial regression.

There is previously a question (How compute the Intra-Class Correlation for a Negative Binomial Mixed Model in lme4) that addresses how to calculate ICC for a 2-level negative binomial model, but it's unclear to me how to include the third level.

Community
  • 1
alxlvt
  • 145
  • 1
  • 6

1 Answers1

6

I don't know if you still need the answer for this, but I'll try anyway.

The ICC for a two level negative binomial model (Tseloni and Pease, 2003) can be easily calculated by:

$$ \rho = \frac{\sigma_{j}^2}{\sigma_{j}^2 + \alpha} $$

where $\sigma_{j}^2$ is the variance of between-group differences (level 2), and $\alpha$ is the variance at level 1, though the parameter reported by lme4 as the overdispersion parameter ($\theta$) is $\alpha^{-1}$.

So, in a three level model, there are two intra-class correlations that can be calculated: individuals within level-2 groups, and level-2 groups within level-3 groups.

Using a hypothetical example, if we had pupils nested in classes ($j$) nested in schools ($k$), the formulas for the ICC are:

$$ \rho_{class} = \frac{\sigma_{j}^2 + \sigma_{k}^2}{\sigma_{j}^2 +\sigma_k^2 + \alpha} $$

$$ \rho_{school} = \frac{\sigma_{k}^2}{\sigma_{j}^2 +\sigma_k^2 + \alpha} $$

where $\sigma_j^2$ is the between class variance, $\sigma_k^2$ is the between school variance, and $\alpha$ is the between pupil variance. $\rho_{class}$ would be the correlation between two pupils in the same class, and $\rho_{school}$ the correlation between two classes in the same school.

Now, to calculate these ICCs using lme4 we need to access the specific estimates stored in the glmerMod object generated by glmer.nb.

So assuming you have a three level model named m

class(m)   ## "glmerMod"

### store the intercepts variance, which
### frustratingly, is also called theta in lme4 models
### and it's stored as the sd, hence the need to square it
var_k <- as.numeric(getME(m, "theta")[2]^2) # level 3 variance
var_j <- as.numeric(getME(m, "theta")[1]^2) # level 2 variance

### store the alpha value (which lme4 stores as theta = 1/alpha)
alpha <- 1/getME(m, "glmer.nb.theta")

### ICC for level 2
ICC_l2 <- (var_k + var_j)/(var_k + var_j + alpha) 

### ICC for level 3
ICC_l3 <- var_k/(var_k + var_j + alpha)

References

Tseloni, A., & Pease, K. (2003). Repeat personal victimization. ‘Boosts’ or ‘flags’? British Journal of Criminology, 43(1), 196-212. doi:10.1093/bjc/43.1.19

prestevez
  • 191
  • 2
  • 6