You can have a look at the following piece of code on how to simulate data from a linear mixed model:
n <- 100 # number of subjects
K <- 8 # number of measurements per subject
t_max <- 15 # maximum follow-up time
# we construct a data frame with the design:
DF <- data.frame(id = rep(seq_len(n), each = K),
time = rep(seq(0, t_max, length.out = K), n),
sex = rep(gl(2, n/2, labels = c("male", "female")), each = K))
X <- model.matrix(~ sex * time, data = DF)
Z <- model.matrix(~ time, data = DF)
betas <- c(-2.13, 0.5, 1, -0.5) # fixed effects coefficients
sigma <- 1.5 # standard deviation error terms
D11 <- 2 # variance of random intercepts
D22 <- 1 # variance of random slopes
D12 <- 0.8 # covariance random intercepts random slopes
D <- matrix(c(D11, D12, D12, D22), 2, 2)
# we simulate random effects
b <- MASS::mvrnorm(n, rep(0, ncol(Z)), D)
# linear predictor
eta_y <- drop(X %*% betas + rowSums(Z * b[DF$id, ]))
# we simulate normal longitudinal data
DF$y <- rnorm(n * K, mean = eta_y, sd = sigma)
library("lme4")
lmer(y ~ sex * time + (time | id), data = DF)
EDIT: To simulate directly with matrix algebra, the Z matrix needs to become block diagonal. In R you could calculate the linear predictor in this manner using, for example, this syntax:
library("Matrix")
Z2 <- as.matrix(bdiag(lapply(split(Z, DF$id), matrix, ncol = 2)))
b_vec <- c(t(b)) # vector of random effects of all subjects
eta_y2 <- drop(X %*% betas + Z2 %*% b_vec)
all.equal(eta_y, eta_y2)