How to simulate random effects models? It is quite simple to simulate linear models:
set.seed(42)    
years <- rnorm(100, 12, 8)
work_hours <- rnorm(100, 8, 2)
income <- 2*years + 0.5*work_hours + 2000 + rnorm(100, 0, 10)
plot(work_hours, income2)
lmmodel <- lm(income ~ years + work_hours)
summary(lmmodel)

Or logistic models:
set.seed(42)
x1 <-  rnorm(100) 
x2 <-  rnorm(100)
z <- 1 + 2*x1 + 3*x2   
pr <- 1/(1+exp(-z)) 

y = rbinom(100,1,pr) 

df <- data.frame(y=y,x1=x1,x2=x2)
logitmodel <- glm( y~x1+x2,data=df,family="binomial")
summary(logitmodel)

So, how does one simulate random effects models?
I mean, there are lots of "flavors" with this class of models. Looking at Faraway's [book][1] there are:


*

*Blocks as Random Effects

*Split Plots

*Nested Effects

*Crossed Effects

*Multilevel Models

*Repeated Measures

*Longitudinal/Panel Data

*Mixed effect models for nonnormal Responses


How would I simulate them so I can toy with them?
[1]: Extending the Linear Model with R - John Faraway
 A: Just write down an (algebraic) formula for the model, and simulate from that description. I will give a very simple example, a model with multiple observations of the same subjects, with an exchangeable covariance structure. Such a structure can be represented with a random intercept for each subject. Also an subject-level covariate:
$$ y_{ij}=\mu + \alpha x_i + \epsilon_i + \epsilon_{ij} $$ for $i=1,2,\dotsc,n$ and $j=1,\dotsc,k$ within each subject. So this is a balanced model. The same principle is used for unbalanced situations, but that gives more programming. Then we must specify values for fixed parameters  and distributions for random effects $\epsilon_i, \epsilon_{ij}$.   Some simple R code is:
N <- 20 # Number subjects
k <- 4  # Number obs within subject
set.seed(7*11*13) # My public seed

id <- as.factor(1:N)
x <-  runif(N, 1, 5)
idran <- rnorm(N, 0, 1)
obsran <- rnorm(N*k, 0, 2)
mu <- 10.
alpha <- 1.

X <- rep(x, each=k)
Y <- mu + alpha*X + rep(idran, each=k) + obsran

A plot of this simulated data is:

For more complex situations it would help with some preprogrammed package, there is a package simstudy on CRAN which can help. See also Model Matrices for Mixed Effects Models and https://stackoverflow.com/questions/30896540/extract-raw-model-matrix-of-random-effects-from-lmer-objects-lme4-r,   https://stackoverflow.com/questions/55199251/how-to-create-a-simulation-of-a-small-data-set-in-r.
