Why estimated population variance differs from estimated $\sigma^2 + \tau^2$ in this random effects ANOVA?

Question

A random effects ANOVA model is typically written as

$Y_{ij} = \gamma_{00} + u_{0j} + \epsilon_{ij}$

. and the total variance of the outcome variable is decomposed into

$var(Y_{ij}) = \tau^2 + \sigma^2$

where $\tau^2$ represents the variance of $u_{0j}$ and $\sigma^2$ the variance of $\epsilon_{ij}$.

I have some code which creates estimates of $\sigma^2$ and $\tau^2$ from artificial data.

Then I run an empty multilevel model and can see that the estimates of $\sigma^2$ and $\tau^2$ exactly match the results lme4 gives me.

However, when I then compute estimated population variance as

estimated_population_variance <- var(my_data$y) * (N-1)/N

it doesn't equal $\sigma^2 + \tau^2$. Why not? Aren't $\sigma^2 + \tau^2$ just a decomposition of the total variance?

I figure it must have something to do with how I'm estimating population variance but I'm not sure what I've done wrong.

Code:

library(lme4)

set.seed(123) 
groups <- rep(1:10, each = 10)
random_numbers <- sample(1:99, length(groups), replace = TRUE)
my_list <- list(group = groups, y = random_numbers)
my_data <- as.data.frame(my_list)

N <- 100
num_groups <- 10
average_group_size <- 10

table_of_groups <- my_data %>%
  group_by(group) %>%
  summarise(mean = mean(y), variance = var(y), n = n())

s2within <- mean(table_of_groups$variance)
estimated_sigma2 <- s2within

s2between <- sum((table_of_groups$mean-mean(my_data$y))^2) * (1/(num_groups-1))
estimated_tau2 <- s2between - (s2within/average_group_size)

empty_model <- lmer("y ~ 1 + (1 | group)", data = my_data) 
summary(empty_model)

estimated_population_variance <- var(my_data$y) * (N-1)/N

estimated_tau2 + estimated_sigma2

Answer 1

The "common" estimate of the variance of $Y$ (obtained with var(Y) in R), namely the total sum of squares divided by $M-1$ ($M$ being the total number of cases in all groups together) can be proven to be equal to $\hat{\sigma}^2 + \frac{n*(N-1)}{M-1} \hat{\tau}^2$, where $n$ is the number of cases in each group (assumed equal for simplicity), 10 in your example), $N$ is the number of groups (also 10 in your example); $\hat{\sigma}^2$ and $\hat{\tau}^2$ are the estimates of the within and between variances, as provided by the random effects one-way Anova (and shown under Random Effects by your summary(empty_model)).

The fraction $\frac{n*(N-1)}{M-1}$ equals 90/99 in your example. In general, this fraction is always less then 1, meaning that the "common" variance formula leads to a lower variance estimate than the estimate provided by the random-effects Anova model, being $\hat{\sigma}^2 + \hat{\tau}^2$. Only when all group means would be equal, and hence $\hat{\tau}^2$=0, the common estimate and the random effect estimate would be equal, both being equal to $\hat{\sigma}^2$.

Because the estimates $\hat{\sigma}^2$ and $\hat{\tau}^2$ are unbiased, so is there sum, meaning that $\hat{\sigma}^2 + \hat{\tau}^2$ is un unbiased estimate of the variance of $Y$, i.e. if REML estimation is used (which is done in R's "lmer" by default). This entails that the "common" variance estimate is downwardly biased, in case the true data are generated by a random-effects Anova model $Y_{ij}=\gamma_{00}+u_{0j}+e_{ij}$. Notice that the "common" variance estimate is in line with the simple empty model $ Y=b_0+e_i$.

EDIT

I would like to add a more intuitive explanation of the downward bias of the common variance estimate. Due to the similarity of $Y$ values within the groups, the common variance estimate is lower than in case such similarity would not exist. This can made plausible by realising that the common sample variance can also be expressed in terms of squared differences between all individual $Y$ values:

$var(Y) = \frac{\sum\limits_{i=1}^N (Y_i-\bar{Y})^2} {N-1} = \frac{\sum\limits_{i=1}^N \sum\limits_{j=1}^N (Y_i-Y_j)^2}{2N(N-1)}$

The closer the two values $Y_i$ and $Y_j$ are within each of the possible pairs, the lower the estimated variance is. If the groups are homogeneous ($Y$ values within the same group are "close") the common variance estimate is lower: if one would have taken $Y$ values completely at random instead of from predefined homogeneous groups, the common variance would have been higher and equal to $\hat{\sigma}^2 + \hat{\tau}^2$.