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!
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
j houses the IDs for subjects. There are 10 subjects observed 10 times (if you want to use 20 subejcts observed 25 times then do j=rep(1:20, 25)). $Z$ is then an indicator matrix telling me which subject is being observed. The columns correspond to subejcts, the rows to observations. So if subject $j$ is the $i^{th}$ observation then $Z_{i, j} = 1$.
$\endgroup$
Commented
Apr 15, 2022 at 15:03
As Xi'an indicates, start with the upper levels. Simulate the random intercepts for each group, $u_{0,i}$. The rnorm(∙) function will work.
Next, you want to generate the $x_{i,j}$ values and the residual error terms at the lowest level. I recommend you make sure these are uncorrelated (as this is an assumption of the model).
The easiest way to do this is to generate two random variables $x$ and $\varepsilon^*$. Run a regression predicting $\varepsilon^*$ from $x$, and then retain the residuals from this regression. These can be your $\varepsilon$ terms (and they will random and uncorrelated with $x$).
This is most easily done using rnorm(∙) as above, and then something like
var.eps <- lm(var.eps ~ x)$resid
Now you can use these values to generate your dependent variables.
