#' --- #' title: "Code for 'Data Analysis And Graphics Using R', 3rd edn, CUP, 2010" #' author: "John Maindonald and John Braun" #' date: "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) 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()