text: Q:
  
  Which ICC estimator to use under the fixed effects model?
  
  In the random effects model $$y_{ij} = \mu + \alpha_j + x_{ij} + \epsilon_{ij}$$ the intra-class correlation coefficient is given by $$ICC = \frac{\sigma_{\alpha}}{\sigma_{\alpha}+\sigma_{\epsilon}}$$ and can be estimated by plugging in the variance estimates of the random effects model.
  When the number of clusters / groups is small, let's say below 20, we typically do not want to use a random effects model and model $\alpha_j$ as fixed effects instead.
  What is an accepted estimator for the ICC under the fixed efects model? I am searching the literature but cannot seem to find it.
  
  A:
  
  When the number of clusters / groups is small, let's say below 20, we typically do not want to use a random effects model and model $\alpha_j$ as fixed effects instead.
  
  I think most people would disagree that 20 is too few clusters to use a random effects model. While there is no hard and fast rule, 6 seems to be a reasonable consensus.
  
  What is an accepted estimator for the ICC under the fixed efects model? I am searching the literature but cannot seem to find it.
  
  On way to do this is to estimate a model without the grouping factor (call it m0) and then another model with the grouping factor as a fixed effect (call it m1). Then compute the difference in the residual variance and divide it by the residual variance of m0. The idea behind this is that the fixed effects for the grouping factor absorb a certain amount of variance in the response. When they are normally distributed (as is the assumption in linear mixed models) this should be the same amount of variance as that estimated by the random intercepts.
  Here is how it can be done in R:
  We simulate clustered data with an expected ICC of 0.8 (variance of random intercepts of 4, and residual variance of 1):
  > set.seed(2)
  > dt <- expand.grid(hospID = 1:10, patientID = 1:20)
  > dt$Y <- 1
  > X <- model.matrix(~ 1, data = dt)
  > myFormula <- "Y ~ 1 + (1 | hospID)"
  > foo <- lFormula(eval(myFormula), dt)
  > Z <- t(as.matrix(foo$reTrms$Zt))  # design matrix for random effects
  > betas <- 10    # fixed effects (intercept only in this case)
  > b <- rnorm(10, 0, 2)  # random effects (standard deviation of 2, variance of 4)
  
  > dt$Y <- X %*% betas + Z %*% b + rnorm(nrow(dt))
  
  Now we fit the linear mixed model:
  > library(lme4)
  > (lm0 <- lmer(eval(myFormula), dt)) %>% summary()
  Random effects:
   Groups   Name        Variance Std.Dev.
   hospID   (Intercept) 4.011    2.003   
   Residual             1.188    1.090   
  Number of obs: 200, groups:  hospID, 10
  
  And we see the estimated variance components are as expected.
  Now we fit the models m0 and m1 as described above and compute the ICC from the mixed model and also from the linear models:
  > m0 <- lm(Y ~ 1, dt)
  > m1 <- lm(Y ~ 1 + as.factor(hospID), dt)
  
  > dt.vc <- as.data.frame(VarCorr(lm0))  # extract the variance components
  > (ICC.lmm <- dt.vc[1, 4] / (dt.vc[1, 4] + dt.vc[2, 4]))
  [1] 0.7715357
  > (ICC.lm <- (var(residuals(m0)) - var(residuals(m1))) / var(residuals(m0)) )
  [1] 0.7645219
  
  and these seem to agree well. You can change the seed, change the simulated variances add other fixed effect etc as you see fit.
  
  