Code for ‘Data Analysis And Graphics Using R’, 3rd edn, CUP, 2010
John Maindonald and John Braun
Oct 3, 2014
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: Rcppif(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()