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!

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.
