Longitudinal data: repeated measurements of a response (and, perhaps, some covariates) over time on several experimental (or observational) units. Usually the experimental unit is “Subject”

Goal: learn how the response varies over time within a subject, as well as the variation in the time trends between subjects

sleepstudy: A simple longitudinal Data

library(lme4)
library(lattice)

?sleepstudy
data(sleepstudy)
str(sleepstudy)
sleepstudy[1:5,]
attach(sleepstudy)

## 10 days, 18 subjects
table(Days, Subject)
length(unique(Subject))

## order the subjects by their intercepts
xyplot(Reaction ~ Days | Subject, sleepstudy, aspect = "xy",
             layout = c(6,3), type = c("g", "p", "r"),
             index.cond = function(x,y) coef(lm(y ~ x))[1],
             xlab = "Days of sleep deprivation",
             ylab = "Average reaction time (ms)")

Model 0: A simple one-way random-effects model

fm0 = lmer(Reaction ~ 1 + (1 | Subject), sleepstudy)
summary(fm0)

Model 1: A mixed model with random intercepts

fm1 = lmer(Reaction ~ 1 + Days + (1 | Subject), sleepstudy)
summary(fm1)

Model 2: A mixed model with Correlated random intercepts and slopes

fm2 = lmer(Reaction ~ Days + (Days | Subject),   sleepstudy)

## short-hand notation for
## fm2 = lmer(Reaction ~ 1 + Days + (1 + Days | Subject), 
         sleepstudy)

summary(fm2)

## The default is RMLE. To get MLE, we need to set REML = FALSE
fm2.mle = lmer(Reaction ~ Days + (Days | Subject), sleepstudy, REML=FALSE)

## Compare the output from fm2 and fm2.mle
summary(fm2.mle)

Model 3: A mixed model with (independent) random intercepts and slopes

fm3 = lmer(Reaction ~ 1 + Days + (1 | Subject) + 
          (0 + Days | Subject), sleepstudy)

## Note that the expression is 
## NOT "1 + Days + (1 | Subject) + (Days | Subject)"

summary(fm3)

Comparing Model 2 (fm2) with Model 3 (fm3)

anova(fm3, fm2)
fm3.mle = lmer(Reaction ~ Days + (1 | Subject) + (Days | Subject),
         sleepstudy, REML=FALSE)
fm2.mle
fm3.mle

Predict random intercepts and slopes (for Model 3)

ranef(fm3)
dotplot(ranef(fm3), condVar=TRUE)  ## not working; no PIs
dotplot(ranef(fm2), condVar=TRUE)

mycoef=ranef(fm3)[["Subject"]]
mycoef[1:5,]

lm(Reaction ~ Days, sleepstudy, subset=Subject==308)$coef
lm(Reaction ~ Days, sleepstudy)$coef

Visualize how the subject-specific coefficients are shrunk toward the fixed effects estimates.

df = coef(lmList(Reaction ~ Days | Subject, sleepstudy))
fclow = subset(df, `(Intercept)`  251)
cc1 = as.data.frame(coef(fm3)$Subject)
names(cc1) = c("A", "B")
df = cbind(df, cc1)
ff = fixef(fm3)
with(df,
     print(xyplot(`(Intercept)` ~ Days, aspect = 1,
           x1 = B, y1 = A,
           panel = function(x, y, x1, y1, subscripts, ...) {
           panel.grid(h = -1, v = -1)
           x1 <- x1[subscripts]
           y1 <- y1[subscripts]
           larrows(x, y, x1, y1, type = "closed";, length = 0.1,
                angle = 15, ...)
           lpoints(x, y,
           pch = trellis.par.get("superpose.symbol")$pch[2],
           col = trellis.par.get("superpose.symbol")$col[2])
           lpoints(x1, y1,
           pch = trellis.par.get("superpose.symbol")$pch[1],
           col = trellis.par.get("superpose.symbol")$col[1])
           lpoints(ff[2], ff[1], 
           pch = trellis.par.get("superpose.symbol")$pch[3],
           col = trellis.par.get("superpose.symbol")$col[3])
           ltext(fclow[,2], fclow[,1], row.names(fclow),
                 adj = c(0.5, 1.7))
           ltext(fchigh[,2], fchigh[,1], row.names(fchigh),
                 adj = c(0.5, -0.6))
                },
           key = list(space = "top", columns = 3,
                  text = list(c("Mixed model", 
               "Within-group", "Population")),
           points = list(
             col = trellis.par.get("superpose.symbol")$col[1:3],
             pch = trellis.par.get("superpose.symbol")$pch[1:3]))
     )))

Lattice plot for visualizing the shrinkage effect

xyplot(Reaction ~ Days | Subject, sleepstudy, aspect = "xy",
     layout = c(6,3), type = c("g", "p", "r"),
     coef.list = df[,3:4],
     panel = function(..., coef.list) {
     panel.xyplot(...)
     panel.abline(as.numeric(coef.list[packet.number(),]),
     col.line = trellis.par.get("superpose.line")$col[2],lty = 2)
     panel.abline(fixef(fm3),
            col.line = trellis.par.get("superpose.line")$col[4],lty=2) },
     index.cond = function(x,y) coef(lm(y ~ x))[1],
             xlab = "Days of sleep deprivation",
              ylab = "Average reaction time (ms)",
              key = list(space = "top", columns = 3,
              text = list(c("Within-subject", 
                "Mixed model", "Population")),
             lines = list(
                   col = trellis.par.get("superpose.line")$col[c(1:2,4)],
                   lty = c(1,2,2))))

Compare Model 3 (fm3) and Model 1 (fm1)

anova(fm1,fm3)