R – How to Simulate Data for a Random Intercept Effect Model

Question

I want to simulate data of the following random intercept effect model in R
$$Y_{ij}=\alpha+\beta x_{ij}+u_{0,i}+\epsilon_{ij}$$
$$u_{0,i} \sim N(0,\tau^2)$$
$$\epsilon_{ij} \sim N(0,\sigma^2)$$

Here the intercept $\alpha=10$ and slope $\beta=5$, and I know $\tau=10,\sigma=1$.

I wonder how to simulate the data of this model in R? I found some tutorials online but they were not that helpful.

Edit: I forgot to mention that there are 20 individuals and 25 observations for each individual!

Answer 1

Fairly straightforward.

One way to express a random intercept model is in matrix notation like

$$ y = X\beta + Zu$$

Here, $X$ is a design matrix, $\beta$ are regression coefficeints, $Z$ is an indicator matrix for group membership, and $u$ are the random effects.

Using R, this made really easy for your problem

a = 10
b = 5
beta = c(a, b)
j = rep(1:10, 10)
x = rnorm(length(j))
tau = 10
sigma = 1

X = model.matrix(~x)
Z = model.matrix(~factor(j)-1)
u = rnorm(ncol(Z), 0, tau)
eps = rnorm(length(j), 0, sigma)

y = X %*% beta + Z %*% u  + eps