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

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