How to simulate data of a random intercept effect model in R? 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!
 A: 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

A: 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.
