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")
Loading the dataset: 150 rows, 3 variables: "school", "class", "score":
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?
