I'm a bit confused about your two different implementations of computing the ICC for neg. bin. models shown here in your GitHub code and at Cross Validated.

The results differ, so I wonder which ICC is the correct one?

Taking the example code from ?lme4::glmer.nb:

fit <- glmer.nb(y ~ f1*f2 + (1|g), data=dd, verbose=TRUE)
# code from cross validated
ICC.NB(fit)
[1] 0.3172471
# code from this repository
ICC.NB(fit, "g")
[1] 8.08294e-12

Could you please give advice on this?

I am afraid the formula used for the ICC of the poisson model is wrong. I have used the “exact” formula suggested by Stryhn et al. in the past which gave reasonable results.

That’s my implementation of the Stryhn “exact” formula:

StryhnICC.POI <- function(model){
  require(lme4)
  sigma2 <- as.numeric(VarCorr(model)[[1]])
  lp <- as.numeric(fixef(model)%*%colMeans(model.matrix(model)))
  icc <- (exp(2*lp+2*sigma2)-exp(2*lp+sigma2)) / (exp(2*lp+2*sigma2)-exp(2*lp+sigma2)+exp(lp+sigma2/2))
  names(icc)<-"ICC"
  return(icc)
}

Let’s have a look at the sample data:

ml <- read.dta("https://stats.idre.ucla.edu/stat/data/hsbdemo.dta")
model <- glmer(formula = awards ~ 1 + (1 | cid), family = poisson, data = ml)
plot(awards~cid, data=ml, pch=20, col="#00000020", cex=1.5)
points(exp(predict(model))~model@frame$cid, pch=3, cex=2)

The dots are observed counts (the darker the more frequent), crosses are the expected average counts per group as predicted by the poisson model.
It’s evident from that figure that most of the variation is attributable to group membership (represented by cid).
However, the two formulae give vastly different results:

ICC.POI(model, "cid")
[1] 0.5902718
StryhnICC.POI(model)
ICC
0.8679164

Let’s consider the extreme scenarios.

1. No variability is due to the group membership.

set.seed(44)
ml$rancid <- ml$cid[sample(nrow(ml))] # permutating group membership, to break group effects
model <- glmer(formula = awards ~ 1 + (1 | rancid), family = poisson, data = ml)
plot(awards~rancid, data=ml, pch=20, col="#00000020", cex=1.5)
points(exp(predict(model))~model@frame$rancid, pch=3, cex=2)

ICC.POI(model, "rancid")
[1] 0.00251346
StryhnICC.POI(model)
ICC
0.004196454

In this case both calculations are close to zero.

2. The other extreme. All variability is explained by group membership.

library(dplyr)
ml <- data.frame(ml %>% group_by(cid) %>% mutate(meanawards=round(mean(awards)))) # fix counts to group means
model <- glmer(formula = meanawards ~ 1 + (1 | cid), family = poisson, data = ml)
plot(meanawards~cid, data=ml, pch=20, col="#00000020", cex=1.5)
points(exp(predict(model))~model@frame$cid, pch=3, cex=2)

ICC.POI(model, "cid")
[1] 0.7308025
StryhnICC.POI(model)
ICC
0.976839

In this case the Stryhn “exact” formula is close to 1. While the other formula severely underestimates the ICC.

Notes on the implementation of the Stryhn “exact” formula

This implementation of the Stryhn “exact” formula also works with models including an offset for modeling event rates.
It can also be used for calculation of the conditional ICC.
For example, if we include honors as an explanatory variable:

model <- glmer(formula = awards ~ honors + 1 + (1 | cid), family = poisson, data = ml)
StryhnICC.POI(model)
ICC
0.4955388

The ICC is reduced to about 0.5. So after accounting for honors 50% of the remaining variance can be explained by group membership.
Please note, in case explanatory variables are included the conditional ICC reported is for mean covariate levels (in this case for a population with 26.5% honors enrolled).