Chapter 10: Multi-level Models, and Repeated Measures
Packages required: “DAAG”, “lattice”, “MEMSS”, “lme4”, “grid”
The script that follows is designed to be executed as it stands. It sets up
functions that create Ch 10 graphs, then runs those functions. If executed by
clicking on the RStudio “Compile Notebook …” command, it will be processed
through R Markdown, generating a document that includes code, output and
graphs generated by executing the functions. If some needed packages are
missing, warning messages will appear. See further:
http://rmarkdown.rstudio.com/r_notebook_format.html
g10.1 <-
function(device=""){
if(device!="")hardcopy(width=3.25, height=2.5,
device=device, trellis=TRUE,
color=FALSE, pointsize=c(9,6))
oldpar <- par(mar = par()$mar - c(1, 0, 1, 0))
on.exit(par(oldpar))
if(!require(lattice))return("Package 'lattice' must be installed.")
Site <- with(ant111b, reorder(site, harvwt, FUN=mean))
print(stripplot(Site ~ harvwt, data=ant111b, scales=list(tck=0.5),
xlab="Harvest weight of corn"))
if(device!="")dev.off()
}
g10.2 <-
function(device=""){
if(device!="")hardcopy(width=5.5, height=1.8, pointsize=c(8,4),
device=device, trellis=TRUE, color=FALSE)
if(!require(DAAG))return("Package 'DAAG' must be installed.")
if(!require(lme4))return("Package 'lme4' must be installed.")
if(!require(lattice))return("Package 'lattice' must be installed.")
trellis.par.set(layout.heights=list(key.top=0.5, axis.top=1.15,
bottom.padding=0.15, main.key.padding=2.5))
ant111b.lmer <- lmer(harvwt ~ 1 + (1|site), data=ant111b)
prof.lmer <- profile(ant111b.lmer)
if(!require(lattice))return("Package 'lattice' must be installed.")
print(xyplot(prof.lmer, conf=c(50, 80, 95, 99)/100, aspect=0.8,
between=list(x=0.35)))
if(device!="")dev.off()
}
g10.3 <-
function(device="")
{
if(device!="")hardcopy(width=3, height=1.85, device=device)
oldpar <- par(mar = c(4.1, 6.1,1.1,1.1), mgp=c(2.5, 0.75,0), las=2)
on.exit(par(oldpar))
require(DAAG)
classmeans <- with(science, aggregate(like,by=list(PrivPub,Class),mean))
names(classmeans) <- c("PrivPub","Class","like")
with(classmeans, {
boxplot(split(like, PrivPub), xlab = "Class average of score",
boxwex = 0.4, horizontal=TRUE)
rug(like[PrivPub == "private"], side = 1)
rug(like[PrivPub == "public"], side = 3)
})
if(device!="")dev.off()
}
g10.4 <-
function(pnum=c(1,3), device="", term = "school:class"){
# leg <- c("Class effect vs #", "W/i class var vs #",
# "qnorm(); site effect", "qnorm(); w/i class resids")
if(!require(DAAG))return("Package 'DAAG' must be installed.")
if(!require(lme4))return("Package 'lme4' must be installed.")
leg <- rep("",4)
mtext3 <- function(item="A", txt=leg[1],
xleft=par()$usr[1]+1.25*par()$cxy[1]){
mtext(side = 3, line = 0.4, item, cex=0.925, adj = 0)
mtext(side = 3, line = 0.4, txt, at=xleft, adj = 0)
}
science$class <- factor(science$class)
science$school <- factor(science$school)
if(device!="")hardcopy(width=3.75, height=3.75,
device=device, pointsize=9)
science <- na.omit(science)
science1.lmer <- lmer(like ~ sex + PrivPub + (1 | school:class),
data = science, na.action=na.exclude)
oldpar <- par(mfrow=c(2,2), mar = c(4.6, 4.1, 2.6, 2.1), mgp=c(2.5,0.75,0))
on.exit(par(oldpar))
ranf <- ranef(obj = science1.lmer, drop=TRUE)[[term]]
flist <- science1.lmer@flist[[term]]
privpub <- science[match(names(ranf), flist), "PrivPub"]
num <- unclass(table(flist))
if(!all(names(num)==names(flist)))stop("!all(names(num)==names(flist))")
plot(sqrt(num), ranf, xaxt="n", pch=pnum[as.numeric(privpub)],
xlab="# in class (square root scale)",
ylab="Estimate of class effect")
lines(lowess(sqrt(num[privpub=="private"]),
ranf[privpub=="private"], f=1.1), lty=2)
lines(lowess(sqrt(num[privpub=="public"]),
ranf[privpub=="public"], f=1.1), lty=3)
numlabs <- pretty(num)
axis(1, at=sqrt(numlabs), labels=paste(numlabs))
mtext3(item="A", txt=leg[1], xleft=par()$usr[1]+1.25*par()$cxy[1])
res <- residuals(science1.lmer)
vars <- tapply(res, INDEX=list(flist),
FUN=var)
vars <- vars*(num-1)/(num-2)
plot(sqrt(num), vars, xaxt="n", pch=pnum[as.numeric(privpub)],
xlab="# in class (square root scale)",
ylab="Within class variance")
lines(lowess(sqrt(num[privpub=="private"]),
as.vector(vars)[privpub=="private"], f=1.1), lty=2)
lines(lowess(sqrt(num[privpub=="public"]),
as.vector(vars)[privpub=="public"], f=1.1), lty=3)
axis(1, at=sqrt(numlabs), labels=paste(numlabs))
mtext3(item="B", txt=leg[2], xleft=par()$usr[1]+1.25*par()$cxy[1])
qqnorm(ranf, ylab="Ordered site effects", main="")
mtext3(item="C", txt=leg[3], xleft=par()$usr[1]+1.25*par()$cxy[1])
qqnorm(res, ylab="Ordered w/i class residuals", main="")
mtext3(item="D", txt=leg[4], xleft=par()$usr[1]+1.25*par()$cxy[1])
par(mfrow=c(1,1))
par(fig=c(0,1,0.8,1), mar=c(0,0,0,0), mgp=c(0,0,0),
oma=c(0,0,0.5,0), new=TRUE)
plot.new()
plot.window(xlim=c(0,1), ylim=c(0,1))
legend(x="top", legend=c("Private ", "Public"), pch=c(1,3),
lwd=c(1,1), lty=2:3,
xjust=0.5, yjust=0, horiz=TRUE, merge=FALSE, bty="n")
if(device!="")dev.off()
}
g10.5 <-
function(device=""){
if(device!="")hardcopy(width=4.25, height=4.25, device=device)
oldpar <- par(mar=c(0.5,1,1,1),mgp=c(0,.5,.5))
on.exit(par(oldpar))
plot(c(0,13),c(0,13),type="n",xlab="",ylab="",
axes=F)
eps <- 0.1
vines<-function(x=1,y=1,subp=0){
lines(c(x,x+1,x+1,x,x),c(y,y,y+1,y+1,y))
points(c(x+.2,x+.8,x+.8,x+.2),c(y+.2,y+.2,y+.8,y+.8),pch=3)
text(x+.5,y+.5,paste(subp))
}
k<-0
for(i in c(1,3,5,7)){k<-k+1; vines(1,i,c(3,1,0,2)[k])}
k<-0
for(i in c(1,3,5,7)){k<-k+1; vines(4,i,c(2,1,0,3)[k])}
k <- 0
for(i in c(1,4,4,1)){k<-k+1
j<-c(9,9,11,11)[k]
vines(i,j,c(3,2,1,0)[k])
}
lines(c(3,3,NA,3,3),c(0,2.85,NA,10.15,13),lty=2)
lines(c(0,0,NA,0,0),c(0,2.85,NA,10.15,13),lty=2)
lines(c(8,8,NA,8,8),c(0,4.5,NA,8.5,13),lty=2)
lines(c(0,1.25,NA,6.75,8),rep(0,5),lty=2)
lines(c(0,1.25,NA,6.75,8),rep(13,5),lty=2)
lines(c(1,5,5,1,1)+c(-eps,eps,eps,-eps,-eps),
c(9,9,12,12,9)+c(-eps,-eps,eps,eps,-eps),lwd=2)
lines(c(1,2,2,1,1)+c(-eps,eps,eps,-eps,-eps),
c(1,1,8,8,1)+c(-eps,-eps,eps,eps,-eps),lwd=2)
lines(c(1,2,2,1,1)+3+c(-eps,eps,eps,-eps,-eps),
c(1,1,8,8,1)+c(-eps,-eps,eps,eps,-eps),lwd=2)
text(0,6.5,"6 meters height artifical shelter belt", srt=90)
text(4,0,"9 meters height shelter belt")
text(4,13,"19 meters height shelter belt")
text(8,6.5,"Willow shelter belt",srt=90)
text(8.5,10.5,"0 Unshaded \n1 Shaded Aug-Dec \n2 Dec-Feb \n3 Feb-May",
adj=0)
text(3,6.5,"16 meters height willow shelter belt", srt=90)
lines(c(8.25,8.75),c(13,13-sqrt(3)*.5))
lines(rep(8.25,2),c(12.6,13))
lines(c(8.25,8.25+sqrt(3)/5),c(13,12.8))
text(8.85,12.65,"N")
if(device!="")dev.off()
}
g10.6 <-
function(device="")
{
if(device!="")hardcopy(width=3.75, height=2.25, device=device)
oldpar<-par(mar=c(4.5,6.1,1.1,1.1), mfrow=c(1,1))
on.exit(par(oldpar))
if(!require(DAAG))return("Package 'DAAG' must be installed.")
kiwimeans <- with(kiwishade, aggregate(yield,by=list(block,shade),mean))
names(kiwimeans) <- c("block","shade","meanyield")
plotmeans.lm <- lm(meanyield~block+shade, data=kiwimeans)
effects <- predict(plotmeans.lm, type="terms")
kiwishade.lm <- lm(yield ~ block*shade, data=kiwishade)
y <- c(effects[,"block"]/sqrt(2) * sqrt(16),
effects[,"shade"]/sqrt(3) * sqrt(12),
residuals(plotmeans.lm)/sqrt(6) * sqrt(4),
residuals(kiwishade.lm)/
sqrt(12))
n <- rep(4:1, c(12, 12, 12, 48))
gps <- rep(c("block effect\n(ms=86.2)", "shade\n(464.8)",
"plot\n(20.9)", "vine\n(12.2)"), c(12, 12, 12, 48))
gps <- factor(gps, levels = rev(c("block effect\n(ms=86.2)",
"shade\n(464.8)", "plot\n(20.9)",
"vine\n(12.2)")))
gps.sd <- sapply(split(y,gps), sd)
gps.av <- sapply(split(y,gps), mean)
plot(range(y), range(n)+c(-0.5, 0.5), xlab="", ylab="", yaxt="n", type="n")
text(y, n+0.15, "|")
un <- 1:4
sapply(un, function(j){lines(gps.av[j]+c(-gps.sd[j], gps.sd[j]),
rep(j-0.15,2), col="gray")
lines(rep(gps.av[j],2)-gps.sd[j],
j-0.15+c(-.06, .06), col="gray")
lines(rep(gps.av[j],2)+gps.sd[j],
j-0.15+c(-.06, .06), col="gray")
})
mtext(side=1,line=3,
text="Variation in Yield (kg)\n(Add to grand mean of yield = 96.53)")
par(las=2)
axis(2, at=1:4, labels=levels(gps))
par(las=0)
if(device!="")dev.off()
}
g10.7 <-
function(device="", path="")
{
if(device!="")hardcopy(width=4, height=4,
device=device, trellis=TRUE, pointsize=c(7, 5))
if(!require(lattice))return("Package 'lattice' must be installed.")
if(!require(grid))return("Package 'grid' must be installed.")
if(!require(lme4))return("Package 'lme4' must be installed.")
trellis.par.set(layout.heights=list(key.top=0.5, axis.top=1.15,
bottom.padding=0.15, main.key.padding=2.5))
if (!exists("kiwishade.lmer"))
kiwishade.lmer <- lmer(yield ~ shade + (1|block/shade), data=kiwishade)
pk1 <- xyplot(residuals(kiwishade.lmer) ~ fitted(kiwishade.lmer)|block,
groups=shade, layout=c(3,1), par.strip.text=list(cex=1.0),
pch=1:4, data=kiwishade, grid=TRUE,
xlab="Fitted values (Treatment + block + plot effects)",
ylab="Residuals",
scales=list(x=list(alternating=FALSE), tck=0.35),
legend=list(top=list(fun=textGrob,
args=list(label="A", x=0,
just="left"))),
key=list(x=0.052, y=1.24, points=list(pch=1:4),
text=list(labels=levels(kiwishade$shade)),columns=4))
print(pk1, position=c(0,.52,1,1))
ploteff <- ranef(kiwishade.lmer, drop=TRUE)[[1]]
nam <- names(sort(ploteff, decreasing=TRUE)[1:4])
trellis.par.set(layout.heights=list(key.top=0.5, axis.top=0.85,
bottom.padding=0.15, main.key.padding=1))
pk2 <- qqmath(ploteff, ylab="Plot effect estimates",
scales=list(tck=0.35),
key=list(x=0, y=0.98, corner=c(0,1),
text=list(labels=nam, cex=0.75)),
legend=list(top=list(fun=textGrob,
args=list(label="B", x=0,
just="left"))), aspect=1,
xlab="Normal quantiles")
print(pk2, newpage=FALSE, position=c(0,0,.5,.5))
simvals <- simulate(kiwishade.lmer, nsim=2)
simeff <- apply(simvals, 2, function(y) scale(ranef(refit(kiwishade.lmer, y),
drop=TRUE)[[1]]))
simeff <- data.frame(v1=simeff[,1], v2=simeff[,2])
pk3 <- qqmath(~ v1+v2, data=simeff,
ylab="Simulated plot effects\n(2 sets, standardized)",
Qs=list(tck=0.35), aspect=1,
legend=list(top=list(fun=textGrob,
args=list(label="C", x=0,
just="left"))),
xlab="Normal quantiles")
print(pk3, newpage=FALSE, position=c(0.48,0,1,.5))
if(device!="")dev.off()
invisible()
}
g10.8 <-
function(device=""){
if(!require(lattice))return("Package 'lattice' must be installed.")
# if(!require(nlme))return("Package 'nlme' must be installed.")
if(!require(lme4))return("Package 'lme4' must be installed.")
if(device!="")hardcopy(device=device, width=3.6, height=3.4,
color=FALSE, pointsize=9)
oldpar <- par(mgp=c(2.0, .5,0), mar=c(4.1,4.1,1.6,1.1),
pty="s", cex=1, cex.lab=1)
on.exit(par(oldpar))
par(fig=c(0, 0.925, 0, 1))
plot(o2 ~ wattsPerKg, data=humanpower1,
pch=(1:5)[unlist(id)],
col=palette()[unlist(id)],
xlab="Watts per kilogram",
ylab=expression("Oxygen intake ("*ml.min^{-1}*.kg^{-1}*")"))
hp.lmList <- lmList(o2 ~ wattsPerKg|id, data=humanpower1)
sapply(1:length(hp.lmList), function(i)abline(hp.lmList[[i]], col=i))
mtext(side=3, line=0.25, "A", adj=-0.165)
par(fig=c(0.575, 1, 0.025, 0.525), new=TRUE, cex=0.8, mgp=c(1.8,0.5,0))
coefs <- data.frame(t(sapply(hp.lmList, coef)))
names(coefs) <- c("Intercept", "Slope")
par(pty="s", tck=-0.025, xpd=TRUE, cex.axis=0.8, cex.lab=0.75,
mgp=c(1.25,0.25,0))
plot(Slope ~ Intercept, data=coefs, pch=1:5, type="n")
xy <- par()$usr
chw <- par()$cxy[1]
chh <- par()$cxy[2]
rect(xy[1]-3.25*chw,xy[3]-2.75*chh,xy[2]+chw,xy[4]+chh, col="gray85")
axis(1)
axis(2)
box()
with(coefs, points(Slope ~ Intercept, cex=1, pch=1:5))
par(xpd=F)
abline(lm(Slope ~ Intercept, data=coefs)$coef)
par(xpd=TRUE, cex=1)
text(xy[1]-2.25*chw, xy[4]+0.25*chh, "B")
mtext(side=1,line=1.5,"Intercept", cex=0.8)
mtext(side=2,line=1.5,"Slope", cex=0.8)
par(fig=c(0,1,0,1), cex=1, xpd=FALSE)
if(device!="")dev.off()
}
g10.9 <-
function(device=""){
if(!require(lattice))return("Package 'lattice' must be installed.")
if(!require(lme4))return("Package 'lme4' must be installed.")
if(device!="")hardcopy(device=device, width=5.5, height=2.75,
trellis=TRUE, color=FALSE)
trellis.par.set(layout.heights=list(key.top=0.25, axis.top=0.65))
hp.lmer <- lmer(o2 ~ wattsPerKg + (wattsPerKg | id),
data=humanpower1)
hat2 <- fitted(hp.lmer)
hp1.plt1 <- xyplot(o2 ~ wattsPerKg, groups=id, data=humanpower1,
panel=function(x,y,subscripts,groups,hat2,...){
u <- lm(y ~ groups*x);
hat <- fitted(u)
panel.superpose(x,hat,subscripts,groups, type="l")
panel.superpose(x, y=hat2, subscripts, groups,
type="l", lty=2)
},
scales=list(tck=0.5),
xlab="Watts per kilogram",
ylab=expression("Oxygen intake ("*ml.min^{-1}*.kg^{-1}*")"),
legend=list(top=list(fun=textGrob, args=list(label="A", x=0, gp=gpar(cex=1.15)))),
hat2=hat2, aspect=1)
print(hp1.plt1, position=c(0,0,0.5,1))
res <- resid(hp.lmer)
hp1.plt2 <- xyplot(res ~ wattsPerKg, groups=id,
data=humanpower1, xlab="Watts per kilogram",
ylab="Residuals from random lines",
scales=list(tck=0.5), aspect=1,
legend=list(top=list(fun=textGrob, args=list(label="B", x=0, gp=gpar(cex=1.15)))), type="l")
print(hp1.plt2, position=c(0.5,0,1,1), newpage=FALSE)
if(device!="")dev.off()
}
g10.10 <-
function(device="", log=2){
if(!require(lattice))return("Package 'lattice' must be installed.")
if(device!="")hardcopy(device=device, width=5.25, height=3.25,
trellis=TRUE, color=FALSE, pointsize=c(7,5))
plt <- xyplot(distance ~ age | Subject, groups=Sex,
data=Orthodont, type=c("p","r"),
par.strip.text=list(cex=0.75), scale=list(log=log),
layout=c(11,3))
print(plt)
if(device!="")dev.off()
}
g10.11 <-
function(device=""){
if(device!="")hardcopy(device=device, width=4.5, height=2.25,
color=TRUE, pointsize=8)
oldpar <- par(mar=c(4.1,4.1,1.1,1.1), mgp=c(2.5,0.75,0), pty="s",
mfrow=c(1,2))
on.exit(par(oldpar))
if(!require(lme4))return("Package 'lme4' must be installed.")
if(!require(MEMSS))return("Package 'MEMSS' must be installed.")
Orthodont$age <- Orthodont$age-11
ab <- coef(lmList(distance ~ age|Subject, data=Orthodont))
ab$sex <- substring(rownames(ab),1,1)
names(ab) <- c("a", "b", "sex")
xlim <- range(ab[,1])
xlim <- xlim+diff(xlim)*c(-.015, .185)
ylim <- range(ab[,2])
ylim <- ylim+diff(ylim)*c(-.015, .015)
plot(ab[,1], ab[,2], col=c(F="gray40", M="black")[ab$sex],
pch=c(F=1, M=3)[ab$sex], xlim=xlim, ylim=ylim,
xlab="Intercept", ylab="Slope")
use <- ab$a %in% range(ab$a) | ab$b %in% range(ab$b) |
ab$b==min(ab$b[ab$sex=="M"])
text(ab[use, 1], ab[use, 2], rownames(ab)[use], pos=4, xpd=TRUE)
mtext(side=3, line=1,"A: Distances", adj=0)
##
Orthodont$logdist <- log(Orthodont$distance)
ab <- coef(lmList(logdist ~ age|Subject, data=Orthodont))
ab$sex <- substring(rownames(ab),1,1)
names(ab) <- c("a", "b", "sex")
xlim <- range(ab[,1])
xlim <- xlim+diff(xlim)*c(-.015, .185)
ylim <- range(ab[,2])
ylim <- ylim+diff(ylim)*c(-.015, .015)
plot(ab[,1], ab[,2], col=c(F="gray40", M="black")[ab$sex],
pch=c(F=1, M=3)[ab$sex], xlim=xlim, ylim=ylim,
xlab="Intercept", ylab="Slope")
text(ab[use, 1], ab[use, 2], rownames(ab)[use], pos=4, xpd=TRUE)
mtext(side=3, line=1,"B: Logarithms of distances", adj=0)
if(device!="")dev.off()
}
g10.11Col <-
function(device=""){
if(device!="")hardcopy(device=device, width=4.5, height=2.25,
color=TRUE, pointsize=8)
oldpar <- par(mar=c(4.1,4.1,1.1,1.1), mgp=c(2.5,0.75,0), pty="s",
mfrow=c(1,2))
on.exit(par(oldpar))
if(!require(lme4))return("Package 'lme4' must be installed.")
if(!require(MEMSS))return("Package 'MEMSS' must be installed.")
Orthodont$age <- Orthodont$age-11
ab <- coef(lmList(distance ~ age|Subject, data=Orthodont))
ab$sex <- substring(rownames(ab),1,1)
names(ab) <- c("a", "b", "sex")
xlim <- range(ab[,1])
xlim <- xlim+diff(xlim)*c(-.015, .185)
ylim <- range(ab[,2])
ylim <- ylim+diff(ylim)*c(-.015, .015)
plot(ab[,1], ab[,2], col=c(F="red", M="blue")[ab$sex],
pch=c(F=1, M=3)[ab$sex], xlim=xlim, ylim=ylim,
xlab="Intercept", ylab="Slope")
use <- ab$a %in% range(ab$a) | ab$b %in% range(ab$b) |
ab$b==min(ab$b[ab$sex=="M"])
text(ab[use, 1], ab[use, 2], rownames(ab)[use], pos=4, xpd=TRUE)
mtext(side=3, line=1,"A: Distances", adj=0)
##
Orthodont$logdist <- log(Orthodont$distance)
ab <- coef(lmList(logdist ~ age|Subject, data=Orthodont))
ab$sex <- substring(rownames(ab),1,1)
names(ab) <- c("a", "b", "sex")
xlim <- range(ab[,1])
xlim <- xlim+diff(xlim)*c(-.015, .185)
ylim <- range(ab[,2])
ylim <- ylim+diff(ylim)*c(-.015, .015)
plot(ab[,1], ab[,2], col=c(F="red", M="blue")[ab$sex],
pch=c(F=1, M=3)[ab$sex], xlim=xlim, ylim=ylim,
xlab="Intercept", ylab="Slope")
text(ab[use, 1], ab[use, 2], rownames(ab)[use], pos=4, xpd=TRUE)
mtext(side=3, line=1,"B: Logarithms of distances", adj=0)
if(device!="")dev.off()
}
pkgs <- c("DAAG", "lattice", "MEMSS", "lme4", "grid")
z <- sapply(pkgs, require, character.only=TRUE, warn.conflicts=FALSE, quietly=TRUE)
## Loading required package: Matrix
## Loading required package: Rcpp
if(any(!z)){
notAvail <- paste(names(z)[!z], collapse=", ")
warning(paste("The following packages should be installed:", notAvail))
}
g10.1()
g10.2()
g10.3()
g10.4()
g10.5()
g10.6()
g10.7()
g10.8()
g10.9()
g10.10()
g10.11()