# Rcode: Longitudinal Analysis

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 &lt;- x1[subscripts] y1 &lt;- 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)