How to simulate R data for a random effects model set-up?

Question

Suppose we have $m$ schools chosen randomly from among thousands in a large country. Suppose also that $n$ students of the same age are chosen at each school. Let $Y_{ij}$ be the score of the $j$th student at the $i$th school.

One model we can come up with is:

$$ Y_{ij} = \mu + U_i + W_{ij} $$

where the $\mu$ variable is the average score for the entire population. From this model, $U_i$ is a random effect with $W_{ij}$ being the error term. Suppose that I want $U_i$ to be normally distributed. That is, each school has a different normal distribution parameter: $U_i$.

I am wondering how a scenario like this can be simulated in R, for say $m = 5$ schools, $n = 100$ students?

Answer 1

You could use something like this:

# Simulate random effect model.

m <- 5 #number of schools
mu <- 12.56 # global score average
sdU <- 5.67 # standard deviation for the random effect
sdW <- 3.25 # standard deviation for the noise
sU <- sda^2 
sW <- sdf^2

# values for the random effect
U <- rnorm(m, mean=0, sd = sdU)

# Number of observations in each school
n <- 100

N <- m*n

data <- data.frame(school = rep(1:m, each=n))
data$U <- U[data$school] # random factor value for each school
data$val <- rnorm(m*n, mu + data$U, sdW)

library(lme4)
library(lmerTest)

# Estimate parameters set earlier
m1 <- lmer(val ~ 1 + (1|school), data = data)
summary(m1)

Check also this https://aosmith.rbind.io/2018/04/23/simulate-simulate-part-2/ blog post for extra fixed effects.