R: linear algebra representation of the prediction operator for a mixed effects model

Question

(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?

---

Edit
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.

Answer 1

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)

Answer 2

Replace your last line with:

y <- (X %*% b %x% t(c(1,1,1)) + (X %*% t(u)) + e) * Z

And you'll see that

> all.equal(c(y[1:5, 1], y[6:10, 2], y[11:15, 3]), DT$y)
[1] TRUE

Alternatively, you can just look at y itself, and you'll see that the nonzero elements correspond to DT$y; i.e.

> all.equal(rowSums(y), DT$y)
[1] TRUE

Remember, you can't just multiply the random effect coefficients u by only the factor matrix Z because the columns of u must have different units, and all the elements of a factor matrix have either the same units or are unitless, depending on context. The columns of u are analogous to the rows of b and therefore must multiply the model matrix X. The only reason you were able to get the desired result without doing this in the intercept-only example is that the random intercept would only be multiplying the intercept column of X anyway, i.e., 1.

I don't know of a matrix operation that performs the column-wise selection of elements; if such an operation exists perhaps someone can edit this answer to include it.

Edit: In a comment above, you state that you want a formula giving the predictions using only matrix multiplication (%*%) and addition. In this answer I use the tensor and Hadamard products (%x% and *, respectively) as well, but these are legitimate matrix operations so I see no reason to avoid their use.

Edit 2: As for the covariance of the random effects, that has to be considered only when estimating the coefficients. Once they're estimated, generating predictions is a straightforward matter of multiplying coefficients with predictors and summing.