How compute the Intra-Class Correlation for a Negative Binomial Mixed Model in lme4

Question

I have a negative binomial mixed model with counties nested within states. I've read that the formula for the ICC for a negative binomial model is:

found here: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3916583/

When I run VarCorr(m10), on my random intercept model, I get:
Groups Name Std.Dev. state (Intercept) 0.22709

Could someone explain how to extract the other information from my glmer() model so that I can write a function to compute the ICC? e.g., fixef(m10) gets me the intercept term.

It seems that in your model you did not estimate several of the terms that are referenced in the ICC equation. The equation refers to several different variance components but you only have a single random intercept term. — Jake Westfall, Sep 01 '15 at 21:48

timothy.s.lau · Answer 1 · 2015-09-01T21:02:10.240

2

This is my functional code (using glmer.nb()):

ICC.NB <- function(model){
  sigma_a2 <- as.numeric(exp(data.frame(VarCorr(model))["vcov"]))
  beta <- as.numeric(fixef(model)["(Intercept)"])
  r <- getME(model,"glmer.nb.theta")
  icc <- (exp(sigma_a2) - 1) / ((exp(sigma_a2) - 1) + ((exp(sigma_a2) - 1) / r) + (exp(-beta - (sigma_a2 / 2))))
  icc
}
ICC.NB(model = m12)

edited Sep 01 '15 at 21:02

answered Sep 01 '15 at 20:11

timothy.s.lau

1,113
3
14
28

1

I think there are some ways to make this less fragile and more readable, e.g. use ["vcov"] in place of [4] in the sigma extraction; use "(Intercept)" instead of [1] to get the intercept; use getME(model,"glmer.nb.theta")) instead of the big ugly gsub call ... – Ben Bolker Sep 01 '15 at 20:55
@BenBolker , I edited the post and implemented your suggested changes. Thank you. – timothy.s.lau Sep 01 '15 at 21:03
1

don't forget to test them ... :-) – Ben Bolker Sep 01 '15 at 21:09
1

PS I think the as.numeric()s are unnecessary? – Ben Bolker Sep 01 '15 at 21:10
Good idea. Now, I'm trying to add a "var" parameter to make it work with a three level model by setting what the sigma should be. – timothy.s.lau Sep 01 '15 at 21:12
a) I wondered if the "exp" in "sigma_a2" is required? b) I have been trying getME(model,"glmer.nb.theta")), but no success even after updating lme4? – Hassan Nov 26 '15 at 22:58
@Hassan, you should give more details/post an issue. – Ben Bolker Dec 20 '15 at 16:53

How compute the Intra-Class Correlation for a Negative Binomial Mixed Model in lme4

1 Answers1

Linked