how to simulate data for linear mixed model I am trying to replicate the paper on linear mixed model by Hossain et al. (2018). But I am unable to replicate their simulation study in which they have generated data by using linear mixed model. I have generated data for fixed effects but it becomes very difficult for me to generate it for random effects. Please suggest some data structure for random effects. I have attached the link of the paper mentioned. [https://link.springer.com/article/10.1007/s00184-018-0656-1]
 A: A mixed model for a random intercept has the form
$$ y_{i, j} = \beta x_{i} + \beta_{0, j} $$
Here, subject $j$ has intercept $\beta_{0, j}$.  The random effects are considered to come from a normal distribution, so $\beta_{0, j}$ can be written as
$$ \beta_{0, j} = \beta_0 + z_{0, j}\sigma_0$$
where $z_{0, j}$ is a standard normal random variable and $\beta_0$ is the population mean intercept.  We can simulate this in R using
library(tidyverse)

set.seed(0)
n_subjects = 10
sigma_gamma = 0.25
n_obs = 10 # Number of times each subject is observed
n = n_obs * n_subjects # Total number of observations

subjects = rep(1:n_subjects, n_obs)
sigma_0 = 0.25
z_0 = rnorm(n_subjects, 0, sigma_0) #Each subject gets their own random effect

x = rnorm(n)
b = 2
b0 = 1
sigma_0 = 0.25
y = b*x + (b0 + z_0[subjects]*sigma_0)

We can extend the model to include random effects for the effect of $x$.  In a similar way, let $\beta_j$ be the effect of $x$ for subejct $j$.  We can rewrite $\beta_j$ to be $\beta + z_{1, j}\sigma_1$ where $z_{1, j}$ is a standard normal random variable

sigma_0 = 0.25
sigma_1 = 1.0
z_0 = rnorm(n_subjects, 0, sigma_0) #Each subject get's their own random effect
z_1 = rnorm(n_subjects, 0, sigma_1) #Each subject get's their own random effect

x = rnorm(n)
b = 2
b0 = 1
sigma_0 = 0.25
y = (b + z_1[subjects]*sigma_1)*x + (b0 + z_0[subjects]*sigma_0)

Lastly, we can simulate correlated random effects by drawing the $z$ jointly from an anisotropic gaussian

sigma_0 = 0.25
sigma_1 = 1.0
Sigma = matrix(c(sigma_0^2, 0.3*sigma_0*sigma_1, 0.3*sigma_0*sigma_1, sigma_1^2), nrow = 2)
Z = MASS::mvrnorm(n_subjects, c(0,0), Sigma)
z0 = Z[,1]
z2 = Z[,2]

x = rnorm(n)
b = 2
b0 = 1
sigma_0 = 0.25
y = (b + z_1[subjects]*sigma_1)*x + (b0 + z_0[subjects]*sigma_0)
```

