# How to calculate "by hand" between group variance and the ICC for a three-level hierarchical model?

## Introduction:

This question is related to this and this questions that were answered before, but I would like a more detailed answer on how between-groups variance are calculated in a (null) three-level multilevel model. These details would include how those variance estimations could be approximately reproduced with simple algebra.

In the second question linked above, the author claims to have reached a similar answer to my question for a two-level model. However, he does not explain what he did in a generalizable formula.

### Example:

Following the example given in the first link above, let's assume that our outcome variable is "test scores" and we have:

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^2_j$$ is the between class variance, $$\sigma^2_k$$ 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.

### Data sample and code:

Now I am providing some sample data (150 rows and 3 columns of my dataset). I first calculate between-group variances and the ICC with the lme4 package and then I try (and failed) to reproduce these calculations with simple algebra.

My question is: how to reproduce the calculation of the between groups variances and the ICC with simple algebra (or "by hand")

df <- structure(list(school = structure(c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 21, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 43, 57, 57, 57, 57, 57, 57, 57, 57, 57, 57, 57, 57, 57, 57, 57, 57, 57, 57, 57, 57, 57, 57)), class = structure(c(16, 16, 16, 16, 16, 16, 16, 16, 32, 32, 32, 32, 32, 32, 32, 32, 36, 36, 36, 36, 36, 36, 36, 36, 624, 624, 624, 624, 624, 624, 624, 624, 96403, 96403, 96403, 96403, 96403, 96403, 96403, 96403, 112680, 112680, 112680, 112680, 112680, 112680, 112680, 112680, 112680, 864, 864, 864, 864, 864, 864, 864, 864, 864, 864, 60128, 60128, 60128, 60128, 60128, 60128, 60128, 60128, 90214, 90214, 90214, 90214, 90214, 90214, 90214, 90214, 90214, 1886, 1886, 1886, 1886, 1886, 1886, 1886, 1886, 1925, 1925, 1925, 1925, 1925, 1925, 1925,
1925, 64746, 64746, 64746, 64746, 64746, 64746, 64746, 64746, 2534, 2534, 2534, 2534, 2534, 2534, 2534, 2534, 2534, 2534, 2538, 2538, 2538, 2538, 2538, 2538, 2538, 2538, 2538, 54520, 54520, 54520, 54520, 54520, 54520, 54520, 54520, 54520, 19392, 19392, 19392, 19392, 19392, 19392, 19392, 19392, 22007, 22007, 22007, 22007, 22007, 22007, 22007, 22007, 99370, 99370, 99370, 99370, 99370, 99370)), score = c(11.844183, 11.86921, 11.46746, 11.6552581, 9.887791, 9.630146, 11.58856, 11.6230743, 9.519587, 10.00369, 10.05935, 10.4259816, 10.03577, 8.313779, 10.55169, 9.19653872, 11.15101, 10.54217, 10.87741, 11.2772524, 10.66766, 10.58322, 9.378927, 11.4241531, 10.25175, 10.78538, 10.12276, 9.82222516, 10.56482, 9.968682, 9.967876, 10.5387201, 10.33726, 8.983298, 10.94701, 10.6756158, 10.22041, 10.47773, 9.496219, 10.3640281, 11.03855, 10.98184, 11.13349, 10.9017617, 11.12337, 10.57085, 10.86873, 10.46775, 11.183895, 10.95111, 9.671874, 10.2109093, 11.20039, 10.7974, 11.24101, 10.4301442, 11.29554, 9.684299, 11.22843, 10.0653824, 10.20983, 10.23446, 9.689934, 9.90960004, 8.839349, 9.71245, 9.228954, 10.8699612, 10.57655, 11.19091, 10.72089, 10.3339444, 9.841154, 10.86726, 10.65764, 11.782758, 11.71146, 11.35271, 11.03756, 11.5441474, 11.73609, 11.6451, 11.43047, 11.3758962, 11.15836, 11.40763, 11.15246, 11.1065933, 12.02961, 11.35606, 11.76986, 11.5538704, 11.02344, 11.09335, 11.35377, 11.0149371, 11.36095, 11.10366, 11.10921, 11.7059107, 12.00548, 10.8639, 10.97149, 11.0972835, 11.23018, 10.64617, 11.29394, 9.88477143, 11.72005, 11.0903, 11.69961, 11.3036283, 11.46173, 10.96258, 11.25487, 11.5900851, 10.98824, 11.48482, 11.30747, 11.2701052, 10.81601, 11.00135, 11.21832, 11.0314793, 10.73993, 11.19909, 10.82092, 10.1296468, 11.01964, 11.21567, 10.38432, 10.7111348, 10.88975, 9.463143, 10.44034, 11.1747141, 10.77896, 11.50726, 10.75442, 11.1716335, 11.67676, 10.62803, 11.04481, 10.2000019, 7.867141, 11.01964, 9.605369, 10.2336048, 9.865208, 0.04218)), row.names = c(NA, -150L), groups = structure(list(school = structure(c(1, 7, 10, 21, 43, 57)),
.rows = structure(list(1:24, 25:49, 50:76, 77:100, 101:128, 129:150), ptype = integer(0), class = c("vctrs_list_of", "vctrs_vctr", "list"))), row.names = c(NA, 6L), class = c("tbl_df", "tbl", "data.frame"), .drop = TRUE), class = c("grouped_df", "tbl_df", "tbl", "data.frame"))



Calculating between-group variances and ICC with lme4:

library(lme4)

m <- lmer(score ~ 1 + (1 | school/class), data = df)
summary(m)


Which gives me as output:

Random effects:
Groups       Name        Variance Std.Dev.
class:school (Intercept) 0.43589  0.6602
school       (Intercept) 0.08417  0.2901
Residual                 0.89498  0.9460
Number of obs: 150, groups:  class:school, 18; school, 6

Fixed effects:
Estimate Std. Error t value
(Intercept)  10.6461     0.2104    50.6


So, the output shows the between-classes variance (0.43589), the between-schools variances (0.08417) and the residual variance (0.89498).

To get the total variance I just need to sum up all the variances given in the output. To compute the ICC I just apply the formulas that were given above (in the example section of my question). So:

class_var_m <- 0.43589
school_var_m <- 0.08417
residual_var_m <- 0.89498
total_var_m <- class_var_m + school_var_m + residual_var_m
ICC_class_m <- (class_var_m + school_var_m)/total_var_m


The total variance of scores is 1.41504 and the model ICC for classes (or ICC_class_m) is 0.3675 .

Calculating between-group variances with algebra:

Now, I want to reach the same results by my own calculations, to try to understand better the ICC.

For the total variance, I calculate the variance of the score for the whole sample:

total_var_h <- var(df$score) # "h" stands for "hand calculations"  Which yields: 1.31584, a different but relatively close value to the model estimate for the total variance (1.41504). I wonder if this difference could be due to some sort of weighting performed by the model to account for different number of cases in each group. Then I try to compute a between group variance for classes. First, I create a dataframe with classes' means: df_class <- df %>% group_by(class) %>% summarise(class_score = mean(score))  Second, I compute classes' variance around their mean: class_var_h <- var(df_class$class_score)


This yields a quite different result from what I get from the model: 0.6662 (against 0.43589). Why is this so different?

I repeated the same computations for school level:

df_school <- df %>%
group_by(school) %>%
summarise(school_score = mean(score))

school_var_h <- var(df_school$school_score)  Then I got 0.2322 (against 0.08417 from the model). Could someone clarify what I am doing wrong? Last comment: When I do the same calculations above with my original dataframe (28k observations), the results are somewhat different. The total variances are almost identical in both calculations. However, the between-groups variances in my "calculations by hand" are almost exactly the twice the size of the between groups variance that I get from the model. Therefore, if I divide by two the between-groups variances I computed, I get roughly the same result as the lme4 model. Any clues on why so? ## 1 Answer I wonder if this difference could be due to some sort of weighting performed by the model to account for different number of cases in each group. The difference between the results of lmer and your manual computations is due to that you are manually computing fixed effects-estimates of the group means (and their sample variance) a.k.a. BLUEs (best linear unbiased estimator), while lmer computes random-effects estimates, a.k.a. BLUPs (best linear unbiased predictor). BLUEs have higher variance than BLUPs, as also indicated by your results. My question is: how to reproduce the calculation of the between groups variances and the ICC with simple algebra (or "by hand") You cannot do that by 'simple' algebra, that's why we use lmer: in order to obtain REML estimates of the variances. In a mixed-effects model, the variances of the random components are not computed by taking the sample variance over the group means. The variance of the random effects is estimated, first. Given those variances, predictions of group means (a.k.a. BLUPs) can then be computed from the data. The variances determine the amount of shrinkage of the predicted group means (BLUPs) towards the overall mean, with higher variance yielding less shrinkage: grand_mean <- mean(df$score)

## Plot differences between grand mean, and BLUEs and BLUPs for classes:
n_per_class <- table(df$$class) # students per class BLUE <- tapply(df$$score, df$$class, mean) BLUP <- coef(m)$$class[,1]
par(mfrow = c(1, 2))
plot(BLUE - grand_mean, BLUP - grand_mean, main = "Classes", type = "n")
text(BLUE - grand_mean, BLUP - grand_mean, label = n_per_class, cex = .7)
abline(0, 1)

## Plot differences between grand mean, and BLUEs and BLUPs for schools:
n_per_school <- colSums(table(df$$class, df$$school) > 0) # classes per school
BLUE <- tapply(df$$score, df$$school, mean)
BLUP <- coef(m)\$school[,1]
plot(BLUE - grand_mean, BLUP - grand_mean, main = "Schools", type = "n")
text(BLUE - grand_mean, BLUP - grand_mean, label = n_per_school, cex = .7)
abline(0, 1)


Note that I used the within-group sample size as a plotting symbol. The straight line indicates where the BLUP and BLUE estimates would be equal. What we see from the plots:

• BLUPs are shrunken towards the grand mean, compared to BLUEs (i.e., values on the $$y$$ axis are closer to 0 than the values on the $$x$$ axis).

• The further the BLUE group mean deviates from the grand mean, the more the BLUP will be shrunken towards the grand mean (i.e., shrunken towards a value of 0 on the $$y$$ axis in the plot).

• The smaller the variance of the random effect, the more the BLUPs will be shrunken towards the grand mean: There is stronger shrinkage for the means of schools than for the means of classes, because the variance of that component is smaller.

• The group sizes also factor into the BLUP estimates. If the group sizes are equal (for school), all group means are shrunken in a linear fashion. If group sizes differ (for class), the amount of shrinkage increases as a function of 1) deviation from the group mean, and 2) group size.

• @Socrates There is no closed-form solution for the variances. Function lmer iterates between estimating the fixed- and the random-effect parameters (variances) in order to optimize the log-likelihood for the full mixed-effects model. You can track this optiomization, for example, if you type: m <- lmer(score ~ 1 + (1 | school/class), data = df, verbose = 1) Mar 31, 2021 at 10:29
• You'll find an example of how to 'manually' reproduce in R the iteratively-estimated results of lmer here: sia.webpopix.org/EMlme.html Mar 31, 2021 at 13:07