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

(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. • The matrix Z should be formed in blocks corresponding to the data for each subject with 0's in all columns except for that which corresponds to the random effects for the ith subject and 0's in all rows except for that which corresponds to the$n_i$observations for the ith subject. If you have N subjects with$n_i$observations per subject and$q$random effects, it should have$Nn_i$rows and$qN$columns. – hard2fathom Oct 23 '19 at 14:12 • Regarding the edit to your question, the Wikipedia article is poorly written. When they say "random effect," they really mean "random intercept". Their model is not a general representation of a linear model with mixed effects; i.e., slopes. That's because the slopes would have to multiply not only the factor matrix$Z$but the model matrix$X$as well; see my answer below. – Josh Oct 23 '19 at 15:03 • Also, if I understand your question correctly, I think including the calls to lmer() is only adding confusion. You are assuming that the parameters have already been estimated and are simply trying to represent the prediction operator using linear algebra; the estimation of the parameters doesn't seem relevant here. – Josh Oct 23 '19 at 15:07 • Apologies for the multiple comments -- in the first line of the question you state that you want to simulate a mixed model, but my understanding of your question (based on your comment to my answer below) is that you seek a linear algebra representation of the prediction operator for a mixed effects model. These are not the same thing and you might get better answers by clarifying exactly what you are trying to achieve/learn here. – Josh Oct 23 '19 at 15:38 • @Josh thanks for your comment, I've edited my post based on your suggestions. – mat Oct 23 '19 at 18:22 ## 2 Answers 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)

• Thank you for the code, but it would really help me to pin down the last equation from my post (y <- X%*%b + Z%*%u + e) using linear matrix algebra (matrix multiplication and addition only). – mat Oct 21 '19 at 9:28
• @mat have a look at my edited answer. – Dimitris Rizopoulos Oct 23 '19 at 19:04

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.

• Thanks for your reply, it is much closer to what I was expecting. I will still wait a bit to see if someone else comes with a different approach. – mat Oct 23 '19 at 14:24