How to simulate data to demonstrate mixed effects with R (lme4)?

Question

As a counterpart to this post, I worked on simulating data with continuous variables, lending themselves to correlated intercepts and slopes.

Although there are great posts on this topic on the site, and outside the site, I had difficulties in coming across a beginning-to-end example with simulated data that paralleled a simple, real-life scenario.

So the question is how to simulate these data, and "test" it with lmer. Nothing new to many, but possibly useful to many others searching to understand mixed models.

Answer 1

If you prefer a blog article format, Hierarchical linear models and lmer is an article I wrote that features a simulation with random slopes and intercepts. Here's the simulation code I used:

rm(list = ls())
set.seed(2345)

N <- 30
unit.df <- data.frame(unit = c(1:N), a = rnorm(N))

head(unit.df, 3)
unit.df <-  within(unit.df, {
  E.alpha.given.a <-  1 - 0.15 * a
  E.beta.given.a <-  3 + 0.3 * a
})
head(unit.df, 3)

library(mvtnorm)
q = 0.2
r = 0.9
s = 0.5
cov.matrix <- matrix(c(q^2, r * q * s, r * q * s, s^2), nrow = 2,
                     byrow = TRUE)
random.effects <- rmvnorm(N, mean = c(0, 0), sigma = cov.matrix)
unit.df$alpha <- unit.df$E.alpha.given.a + random.effects[, 1]
unit.df$beta <- unit.df$E.beta.given.a + random.effects[, 2]
head(unit.df, 3)

J <- 30
M = J * N  #Total number of observations
x.grid = seq(-4, 4, by = 8/J)[0:30]

within.unit.df <-  data.frame(unit = sort(rep(c(1:N), J)), j = rep(c(1:J),
                              N), x =rep(x.grid, N))
flat.df = merge(unit.df, within.unit.df)

flat.df <-  within(flat.df, y <-  alpha + x * beta + 0.75 * rnorm(n = M))
simple.df <-  flat.df[, c("unit", "a", "x", "y")]
head(simple.df, 3)

library(lme4)
my.lmer <-  lmer(y ~ x + (1 + x | unit), data = simple.df)
cat("AIC =", AIC(my.lmer))
my.lmer <-  lmer(y ~ x + a + x * a + (1 + x | unit), data = simple.df)
summary(my.lmer)

Answer 2

The data is completely fictional and the code that I used to generate it can be found here.

The idea is that we would take measurements on glucose concentrations on a group of 30 athletes at the completion of 15 races in relation to the concentration of the made-up amino acid A (AAA) in these athletes blood.

The model is: lmer(glucose ~ AAA + (1 + AAA | athletes)

There is a fixed effect slope (glucose ~ amino acid A concentration); however, the slopes also vary between different athletes with a mean = 0 and sd = 0.5, while the intercepts for the different athletes are spread a random effects around 0 with sd = 0.2. Further there is a correlation between intercepts and slopes of 0.8 within the same athlete.

These random effects are added to a chosen intercept = 1 for fixed effects, and slope = 2.

The values of the concentration of glucose were calculated as alpha + AAA * beta + 0.75 * rnorm(observations), meaning the intercept for every athlete (i.e. 1 + random effects changes in the intercept) $+$ the concentration of amino acid, AAA $*$ the slope for every athlete (i.e. 2 + random effect changes in slopes for each athlete) $+$ noise ($\epsilon$), which we set up to have a sd = 0.75.

So the data look like:

      athletes races      AAA   glucose
    1        1     1  51.79364 104.26708
    2        1     2  49.94477 101.72392
    3        1     3  45.29675  92.49860
    4        1     4  49.42087 100.53029
    5        1     5  45.92516  92.54637
    6        1     6  51.21132 103.97573
    ...

Unrealistic levels of glucose, but still...

The summary returns:

Random effects:
 Groups   Name        Variance Std.Dev. Corr
 athletes (Intercept) 0.006045 0.07775      
          AAA         0.204471 0.45218  1.00
 Residual             0.545651 0.73868      
Number of obs: 450, groups:  athletes, 30

Fixed effects:
             Estimate Std. Error        df t value Pr(>|t|)    
(Intercept)   1.31146    0.35845 401.90000   3.659 0.000287 ***
AAA           1.93785    0.08286  29.00000  23.386  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

The random effects correlation is 1 instead of 0.8. The sd = 2 for the random variation in intercepts is interpreted as 0.07775. The standard deviation of 0.5 for random changes in slopes among athletes is calculated as 0.45218. The noise set up with a standard deviation 0.75 was returned as 0.73868.

The intercept of fixed effects was supposed to be 1, and we got 1.31146. For the slope it was supposed to be 2, and the estimate was 1.93785.

Fairly close!