(See edit at the bottom for the bounty)
I am trying to learn how to simulate LMM data with matrix linear algebra. So far I've managed to simulate a simple model with a random intercept:
library(data.table) library(lmerTest) # Parameters Ngroups <- 3 NperGroup <- 5 N <- Ngroups * NperGroup groups <- factor(rep(1:Ngroups, each = NperGroup)) b0 <- 2 b1 <- 3 x <- rnorm(N) e <- rnorm(N, sd = .1) # Random intercept u0 <- rnorm(Ngroups, sd = .7) y <- b0 + u0[groups] + b1*x + e # Random intercept [matrix algebra] X <- cbind(intercept = 1, x) b <- rbind(b0, b1) Z <- diag(Ngroups)[rep(1:Ngroups, each = NperGroup), ] y <- X%*%b + Z%*%u0 + e
I created an other model with a random intercept and slope to the model as follow:
# Random intercept and slope u0 <- rnorm(Ngroups, sd = .7) u1 <- rnorm(Ngroups, sd = .4) DT$y <- b0 + u0[groups] + (b1 + u1[groups])*x + e
However, I cannot find the way to generate the same data using a linear matrix algebra approach, this is what I have so far:
u <- cbind(u0, u1) y <- X%*%b + Z%*%u + e
What would the formula be? How can I also incorporate the var-cov between random factors?
To clarify, I'm looking for a neat linear algebra representation of the prediction operator for a mixed effects model with a random intercept and a random slope. I am seeking something similar to the equation from Wikipedia (see below), although, as pointed out by @Josh, it doesn't take into account random slopes.