`g6.1` <- function(device="") { if(device!="")hardcopy(width=2.25, height=2.25, device=device) oldpar<-par(mar=c(4.1,4.1,2.1,1.1), mgp=c(2.5,0.75,0), pty="s") on.exit(par(oldpar)) xy<-par()$usr[c(1,4)] attach(allbacks) xlim <- range(volume) xlim <- xlim+c(-.075,.075)*diff(xlim) ## Plot of weight vs volume: data frame allbacks (DAAG) plot(weight ~ volume, data=allbacks, pch=c(16,1)[unclass(cover)], lwd=1.25, xlim=xlim) ## unclass(cover) gives the integer codes that identify levels ## As text() does not accept the parameter data, use with() ## to specify the data frame. with(allbacks,text(weight ~ volume, labels=paste(1:15), cex=0.75, offset=0.35, pos=c(2,4)[unclass(cover)])) par(xpd=TRUE) legend(x=660, y=1700, pch=c(16,1), legend=c("hardback ","softback"), horiz=T, bty="n", xjust=0.5, x.intersp=0.75) detach(allbacks) if(device!="")dev.off() } `g6.10` <- function(device="") { if(device!="")hardcopy(width=5.5, height=3.5, device=device) oldpar <- par(mfrow=c(2,4), mar=c(4.1,4.1,3.1,0.6), mgp=c(2.0,.5,0),oma=c(0,0,1,2), pty="s") on.exit(par(oldpar)) hills.lm <- lm(log(time) ~ log(climb) + log(dist), data = hills) plot.new() par(mfg=c(1,1)) mtext(side=3, line=2.5, "A: Include observation 18 (Knock Hill)", adj=0) par(new=T) plot(hills.lm, sub="", caption=c("Resids vs Fitted", "Normal Q-Q", "Scale-Location", "", "Resids vs Leverage")) lmobj <- lm(log(time) ~ log(climb) + log(dist), data = hills[-18,]) # plot(hills.lm, caption="A: Resids vs Fitted", which=1, sub="") plot.new() mtext(side=3, line=2.5, "B: Omit observation 18 (Knock Hill)", adj=0) par(new=T) plot(lmobj, sub="", caption=c("Resids vs Fitted", "Normal Q-Q", "Scale-Location", "", "Resids vs Leverage"), cex.main=1, cook.levels=c(.25,.5)) if(device!="")dev.off() invisible() } `g6.11` <- function(lmobj = NULL, device="") { if(device!="")hardcopy(width=4, height=2, device=device) oldpar <- par(mfrow=c(1,2), mar=c(4.1,4.1,2.1,1.1), mgp=c(2.25,.75,0),oma=c(0,0,0,1), pty="s") on.exit(par(oldpar)) hills.loglm <- lm(log(time) ~ log(dist) + log(climb), data = hills[-18, ]) lqsobj <- lqs(log(time) ~ log(climb) + log(dist), data = hills[-18,]) hat <- fitted(lqsobj) res <- resid(lqsobj) plot(resid(hills.loglm), res, xlab="Residuals, lm model", ylab="Residuals, lqs model") abline(0,1, col=2) mtext(side=3, line=0.5, "A", adj=0) plot(hat, res, xlab="Fitted values, lqs model", ylab="Residuals, lqs model") mtext(side=3, line=0.5, "B", adj=0) if(device!="")dev.off() invisible() } `g6.12` <- function(y = hills$time[-18], pred = T, device="") { if(device!="")hardcopy(width=2, height=2, device=device) here<-rep(T,35) here[18]<-F # All elements except the 18th are set to T hills.logm <- lm(log(time) ~ log(dist) + log(climb), data = hills, subset=here) hills.ci<-exp(predict(hills.logm,interval="confidence")) low <- hills.ci[,"lwr"] hi <- hills.ci[,"upr"] x<-hills.ci[,"fit"] ylim <- range(c(y, low, hi)) titl <- paste("Observed versus predicted values, for the", "\nhill race data. The line y=x is shown.", "\nThe horizontal limits of the bounding lines", "\nare 95% pointwise confidence limits for", "\npredicted times.",sep="") oldpar <- par(mar = c(4.1,4.1,1.1,1.1), mgp = c(2.5, 0.75, 0) ) on.exit(par(oldpar)) r <- cor(x, y) plot(x, y, xlab = "Predicted time", ylab = "Observed time", type="n", ylim = ylim) abline(0, 1) ord <- order(x) lines(low[ord], x[ord], col = "grey60") lines(hi[ord], x[ord], col = "grey60") points(x, y, pch=16, cex=0.5) topleft <- par()$usr[c(1, 4)] chw <- par()$cxy[1] chh <- par()$cxy[2] cat(titl,"\n") if(device!="")dev.off() invisible(hills.logm) } `g6.13` <- function(device=""){ if(device!="")hardcopy(width=3.25, height=3.25, device=device, pointsize=9) logoddbooks <- log(oddbooks) pairs(log(oddbooks), labels=c("thick\n\nlog(mm)", "breadth\n\nlog(cm)", "height\n\nlog(cm)", "weight\n\nlog(g)")) figtxt<-paste("Scatterplot matrix for the oddbooks data.") print(figtxt) if(device!="")dev.off() invisible() } `g6.14` <- function(device=""){ if(device!="")hardcopy(width=5, height=1.7, device=device) logoddbooks <- log(oddbooks) oldpar <- par(mfrow=c(1,3), mar=c(4.1,4.1,1.6,2.1), mgp=c(2.5,.75,0), oma=c(0,0,0,0), pty="s") on.exit(par(oldpar)) attach(logoddbooks) u1 <- lm(weight~I(thick+breadth+height)) plot(weight~I(thick+breadth+height), xlab="log(thick)+log(breadth)+log(height)", ylab="log(weight)") abline(u1) mtext(side=3,line=0.25, "A", adj=0) u2 <- lm(weight~thick+I(breadth+height)) plot(weight~fitted(u2), xlab="Predicted log(weight)", ylab="log(weight)") abline(0,1) mtext(side=3,line=0.25, "B", adj=0) u3 <- lm(weight~thick+breadth+height) plot(weight~fitted(u2), xlab="Predicted log(weight)", ylab="log(weight)") abline(0,1) mtext(side=3,line=0.25, "C", adj=0) if(device!="")dev.off() } `g6.15` <- function(device=""){ if(device!="")hardcopy(width=2.75, height=2.75, device=device, pointsize=9) attach(oddbooks) area <- breadth*height volume <- area*thick density <- weight/volume pairs(cbind(log(density), log(area), log(thick)), labels=c("log(density)", "log(area)", "log(thick)")) detach(oddbooks) figtxt<-paste("Data derived from the oddbooks data.") print(figtxt) if(device!="")dev.off() invisible() } `g6.16` <- function(device=""){ if(device!="")hardcopy(width=4.5, height=4.5, device=device, pointsize=8) if(!exists("carprice")){ require(MASS) data(Cars93) US <- Cars93$Origin=="USA" carprice <- Cars93[US, c("Type", "Min.Price","Price","Max.Price")] carprice$Range.Price <- carprice$Max.Price - carprice$Min.Price carprice$Range.Price <- round(carprice$Range.Price, 1) carprice$RoughRange <- carprice$Range.Price+rnorm(48, sd=0.01) carprice$RoughRange <- round(carprice$RoughRange, 2) avMPG <- apply(Cars93[US, c("MPG.city","MPG.highway")], 1, mean) carprice$gpm100 <- round(100/avMPG, 1) cat("Data set has been reconstructed. Calculations involving", "\nRoughRange will give different results from those obtained", "\n using the data set carprice in the DAAG package.","\n") } pairs(carprice[,2:7], oma=rep(2,4)) if(device!="")dev.off() } `g6.17` <- function(mu=20, n=100,a=15, b=2.5, sigma=12.5, device=""){ if(device!="")hardcopy(width=3.75, height=2.75, device=device, pointsize=8) if(F){ errorsINx <- function(mu = 20, n = 100, a = 15, b = 2.5, sigma = 12.5, timesSigma=(1:5)/2.5){ mat <- matrix(0, nrow=n, ncol=length(timesSigma)+2) x0 <- mu*exp(rnorm(n,1,0.15)) y <- a + b*x0+rnorm(n,0,sigma) mat[, length(timesSigma)+2] <- y mat[,1] <- x0 sx <- sd(x0) k <- 1 for(i in timesSigma){ k <- k+1 x1 <- x0+rnorm(n, 0, sx*i) mat[, k] <- x1 } mat } } oldpar <- par(mar=c(3.6,3.1,1.6,0.6), mgp=c(2.5,0.75,0), oma=c(1,1,0.6,1), mfrow=c(2,3), pty="s") mu <- 20; n <- 100; a <- 15; b <- 2.5; sigma <- 12.5; timesSigma<-(1:5)/2.5 mat <- errorsINx(mu = 20, n = 100, a = 15, b = 2.5, sigma = 5, timesSDx=(1:5)/2.5) beta <- numeric(dim(mat)[2]-1) sx <- sd(mat[,2]) y <- mat[, 1] for(j in 1:length(beta)){ xj <- mat[,j+1] plot(y ~ xj, xlab="", ylab="", col="gray30") if(j==1) mtext(side=3, line=0.5, "No error in x") else{ xm <- timesSigma[j-1] mtext(side=3, line=0.5, substitute(tau == xm*s[z], list(xm=xm))) } if(j>=4)mtext(side=1, line=2, "x") if(j%%3 == 1)mtext(side=2, line=2, "y") errors.lm <- lm(y ~ xj) abline(errors.lm) beta[j] <- coef(errors.lm)[2] bigsigma <- summary(errors.lm)$sigma print(bigsigma/sigma) abline(a, b, lty=2) } print(round(beta, 3)) if(device!="") dev.off() } `g6.2` <- function(device="", stats=F){ if(device!="")hardcopy(width=5.5, height=1.4, device=device) oldpar <- par(mfrow=c(1,4), mar=c(4.1,4.1,2.1,0.6), mgp=c(2.0,.5,0),oma=c(0,0,0,1), pty="s") on.exit(par(oldpar)) allbacks.lm<-lm(weight~-1+volume+area,data=allbacks) plot(allbacks.lm,caption= c("A: Resids vs Fitted", "B: Normal Q-Q", "C: Scale-Location", "D: Cook's distance"), which=1:4) if(stats) print(summary(allbacks.lm)) if(device!="")dev.off() invisible(allbacks.lm) } `g6.3` <- function(device="", smooth=T){ if(device!="")hardcopy(width=3, height=3, device=device, pointsize=9) if(smooth)here.panel <- panel.smooth else here.panel <- points pairs(litters,cex.labels=1.2, labels=c("lsize\n\n(litter size)","bodywt\n\n(Body Weight)", "brainwt\n\n(Brain Weight)"), oma=rep(2,4), lower.panel=here.panel, upper.panel=here.panel) figtxt<-paste("Scatterplot matrix for the litters data set.", sep="") print(figtxt) if(device!="")dev.off() invisible() } `g6.4` <- function(device="") { if(device!="")hardcopy(width=4, height=4, device=device, pointsize=8) titl <- "Brain Weight as a Function of Litter Size and Body Weight" oldpar <- par(mfrow = c(2, 2), mar = c(5.1,4.1,2.1,2.1), mgp = c(2, 0.5, 0), xpd=FALSE) on.exit(par(oldpar)) plot(litters[, 1], litters[, 3], xlab = "Litter Size", ylab = "Brain Weight (g)", pch = 16) u1 <- lm(brainwt ~ lsize, data = litters) cat("\n") print(summary(u1)$coef) cat("\n") se <- summary(u1)$coef[2, 2] mtext(side = 3, line = 0.25, text = paste("brainwt =", round(u1$ coef[1], 2), "+", round(u1$coef[2], 5), "[SE =", round(se, 5), "]", "lsize"), cex = 0.65) plot(litters[, 2], litters[, 3], xlab = "Body Weight (g)", ylab = "Brain Weight (g)", pch = 16) u2 <- lm(brainwt ~ bodywt, data = litters) print(summary(u2)$coef) cat("\n") mtext(side = 3, line = 0.25, text = paste("brainwt =", round(u2$ coef[1], 2), "+", round(u2$coef[2], 4), "bodywt"), cex = 0.65) plot(litters[, 1], litters[, 2], xlab = "Litter Size", ylab = "Body Weight (g)", pch = 16, mkh = 0.04) r3 <- cor(litters[, 1], litters[, 2]) mtext(side = 3, line = 0.25, text = paste("Correlation =", round(r3, 2) ), cex = 0.65) u <- lm(brainwt ~ lsize + bodywt, data = litters) print(summary(u)$coef) hat <- fitted(u) plot(hat, litters[, 3], xlab = "Fitted Weight (g)", ylab = "Brain Weight (g)", pch = 1, lwd=2) se <- summary(u)$coef[2, 2] se1 <- summary(u)$coef[3, 2] mtext(side = 3, line = 0.5, text = paste("brainwt =", round(u$ coef[1], 2), "+", round(u$coef[2], 4), "[SE =", round(se, 4), "] ", "lsize \n+", round(u$coef[3], 4), "[SE =", round(se1, 4), "] ", "bodywt",sep=""), cex = 0.65) abline(0, 1, lty = 2) cat(titl, "\n") if(device!="")dev.off() invisible(u) } `g6.5` <- function(device="", col=1){ if(device!="")hardcopy(width=5.5, height=1.4, device=device) oldpar <- par(mfrow=c(1,4), mar=c(4.1,4.1,2.1,0.6), mgp=c(2.0,.5,0),oma=c(0,0,0,1), pty="s") on.exit(par(oldpar)) litters.lm<-lm(brainwt~bodywt+lsize,data=litters) plot(litters.lm,caption= c("A: Resids vs Fitted", "B: Normal Q-Q", "C: Scale-Location", "D: Cook's distance"), col=col) figtxt<-paste("Diagnostic plots for the regression of brain", " weight \non body weight and litter size",".", sep="") cat(figtxt, "\n") if(device!="")dev.off() invisible() } `g6.6` <- function(df = litters, device="") { if(device!="")hardcopy(width=4, height=2.15, device=device, pointsize=8) titl <- paste("termplots from the Regression of", "\nof Brain Weight on Litter Size & Body Weight.", sep = "") oldpar <- par(mfrow = c(1, 2), mar = c(4.6, 4.1, 2.1, 2.1), mgp = c(2.5, 0.75, 0), pty="s") par(mex = 1, cex = 1) on.exit(par(oldpar)) u <- lm(brainwt ~ lsize+bodywt, data=df) termplot(u, partial.resid=TRUE, smooth=panel.smooth, col.res="gray30") cat(titl,"\n") if(device!="")dev.off() invisible() } `g6.7` <- function(device="") { if(device!="")hardcopy(width=2.25, height=2.25, device=device) oldpar<-par(mar=c(4.1,4.1,2.1,1.1), mgp=c(2.5,0.75,0), pty="s") on.exit(par(oldpar)) plot(mice12.lm, which=5, sub.caption="") if(device!="")dev.off() } `g6.8` <- function(device="") { if(device!="")hardcopy(width=4.25, height=1.5, device=device) oldpar<-par(mar=c(5.1,3.6,2.1,1.1), mgp=c(2.25,0.75,0), las=1) on.exit(par(oldpar)) allbacks.lm0 <- lm(weight ~ -1+volume+area, data=allbacks) z <- dfbetas(allbacks.lm0) plot(range(z), c(0.25,1), xlab="Standardized change in coefficient", ylab="", axes=F, type="n") chh <- par()$cxy[2] par(xpd=T) axis(1) abline(h=par()$usr[3]) m <- dim(z)[1] points(z[,1], rep(0.5,m), pch="|") points(z[,2], rep(1,m), pch="|") idrows <- (1:m)[apply(z,1,function(x)any(abs(x)>2))] if(length(idrows)>0)text(z[idrows,],rep(c(.5,1)-0.25*chh,rep(length(idrows),2)), rep(paste(idrows), 2), pos=1, cex=0.8) par(family="mono") text(rep(par()$usr[1],2), c(.5,1), c("volume", "area"), pos=2) title(sub="dfbetas(allbacks.lm0)") par(family="sans") if(device!="")dev.off() } `g6.9` <- function(device=""){ if(device!="")hardcopy(width=5.2, height=3.1, device=device, pointsize=c(6,4), trellis=TRUE) oldpar <- par(fig=c(0,0.5,0,1), mar=c(0.5,0.5,2.6,0.5)) on.exit(par(oldpar)) plot.new() mtext(side=3, line=1.5, "A: Untransformed scales:", adj=0, cex=0.9) par(fig=c(0.5,1,0,1), mar=c(0.5,0.5,2.6,0.5),new=T) mtext(side=3, line=1.5, "B: Logarithmic scales", adj=0, cex=0.9) hills.splom <- splom(~hills, cex.labels=1.2, varnames=c("dist\n(miles)","climb\n(feet)", "time\n(hours)"), newpage=FALSE, xlab="") print(hills.splom, position=c(0, 0, 0.5, 1), newpage=FALSE) lhills.splom <- splom(~log(hills), cex.labels=1.2, varnames=c("dist\n(log miles)", "climb\n(log feet)", "time\n(log hours)"), xlab="") print(lhills.splom, position=c(0.5, 0, 1, 1), newpage=FALSE) if(device!="")dev.off() invisible() } `gdump` <- function(fnam=NULL, prefix="~/r-book/ed2/figures/figs", xtras=c("hardcopy","renum.fun","renum.files"), splitchar="/ch"){ if(is.null(fnam)){ path <- getwd() pathtag <- strsplit(path, "/ch", fixed=TRUE)[[1]] fnam <- paste(prefix, pathtag[length(pathtag)], ".R", sep="") } else fnam <- paste(prefix, fnam, sep="/") objnames <- c(objects(pattern="^g", envir=sys.frame(0)), xtras) cat("\nDump to file:", fnam, "\n") print(objnames) dump(objnames, fnam) } `gfocus.demo` <- function(device=""){ library(lattice) library(grid) if(device!="")hardcopy(device=device, width=4, height=2.25, trellis=TRUE, color=TRUE) trellis.par.set(layout.heights=list(key.top=0.25, axis.top=0.5)) hp1.plt1 <- xyplot(o2 ~ wattsPerKg, groups=id, data=humanpower1, panel=function(x,y,subscripts,groups,...){ u <- lm(y~groups*x); hat <- fitted(u) panel.superpose(x,y,subscripts,groups) panel.superpose(x,hat,subscripts,groups, type="l") }, ## key=simpleKey(text=rep("",5), lines=TRUE, columns=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)))) print(hp1.plt1, position=c(0,0,.535,1)) u <- lme(o2 ~ wattsPerKg, random=~wattsPerKg|id, data=humanpower1) hp1.plt2 <- xyplot(o2 ~ wattsPerKg, groups=id, data=humanpower1, panel=function(x,y,subscripts,groups,...){ u <- lm(y~groups*x); hat <- fitted(u) panel.superpose(x,hat,subscripts,groups, type="l") }, xlab="Watts per kilogram", ylab="", legend=list(top=list(fun=textGrob, args=list(label="B", x=0)))) hat <- fitted(u) print(hp1.plt2, position=c(.465, 0,1,1), newpage=FALSE) trellis.focus("panel", row=1, column=1) arglist <- trellis.panelArgs() panel.superpose(x=arglist$x,y=hat,subscripts=arglist$subscripts, groups=arglist$groups, , type="l", lty=2) trellis.unfocus() if(device!="")dev.off() } `gsave` <- function(fnam=NULL, prefix="~/r-book/ed2/figures/figs", splitchar="/ch", xtras=c("renum.fun","renum.files","hardcopy")){ if(is.null(fnam)){ path <- getwd() pathtag <- strsplit(path, "/ch", fixed=TRUE)[[1]] fnam <- paste(prefix, pathtag[length(pathtag)], ".RData", sep="") } else fnam <- paste(prefix, fnam, sep="/") objnames <- c(objects(pattern="^g", envir=sys.frame(0)), xtras) cat("\nDump to file:", fnam, "\n") print(objnames) save(list=objnames, file=fnam) } `gtest` <- function(){ trellis.device(postscript, file="test.eps") trellis.par.set(layout.heights=list(key.top=0.5)) zz.d <- dotplot(variety ~ yield, data = barley, legend=list(top=list(fun=grid.text, args=list(label="ABC", x=0)))) pushViewport(viewport(layout=grid.layout(2, 1))) pushViewport(viewport(layout.pos.row=1)) print(zz.d,newpage=FALSE) upViewport() pushViewport(viewport(layout.pos.row=2)) print(zz.d, newpage=FALSE) popViewport(2) dev.off() } `hardcopy` <- function(width=3.75, height=3.75, color=FALSE, trellis=FALSE, device=c("","pdf","ps"), path=getwd(), file=NULL, format=c("nn-nn", "name"), split="\\.", pointsize=c(8,4), fonts=NULL, horiz=FALSE, ...){ if(!trellis)pointsize <- pointsize[1] funtxt <- sys.call(1) nam <- strsplit(as.character(funtxt), "(", fixed=TRUE)[[1]][1] suffix <- switch(device, ps=".eps", pdf=".pdf") if(is.character(path) & nchar(path)>1 & substring(path, nchar(path))!="/") path <- paste(path, "/", sep="") if(is.null(file)) if(format[1]=="nn-nn"){ if(!is.null(split))dotsplit <- strsplit(nam, split)[[1]] else dotsplit <- nam if(length(dotsplit)==1)dotsplit <- c("", dotsplit) nn2 <- paste(if(nchar(dotsplit[2])==1)"0" else "", dotsplit[2], sep="") if(nchar(dotsplit[1])>0){ numstart <- which(unlist(strsplit(dotsplit[1], "")) %in% paste(0:9))[1] nn1 <- substring(dotsplit[1], numstart) nn1 <- paste(if(nchar(nn1) == 1) "0" else "", nn1, "-", sep="") } else nn1 <- "" file <- paste(nn1, nn2, sep="") } else file <- nam if(nchar(file)>4 & substring(file, nchar(file)-nchar(suffix)+1)==suffix) suffix <- "" file <- paste(path, file, suffix, sep="") print(paste("Output will be directed to file:", file)) dev.out <- device[1] dev.fun <- switch(dev.out, pdf=pdf, ps=postscript) if(trellis){ library(lattice) if(device=="ps") trellis.device(file=file, device=dev.fun, color = color, horiz=horiz, fonts=fonts, width=width, height=height, ...) else trellis.device(file=file, device=dev.fun, fonts=fonts, color = color, width=width, height=height, ...) trellis.par.set(list(fontsize=list(text=pointsize[1], points=pointsize[2]))) } else if (dev.out!=""){ print(c(width, height)) if(device=="ps") dev.fun(file=file, paper="special", horiz=horiz, fonts=fonts, width=width, height=height, pointsize=pointsize[1], ...) else dev.fun(file=file, paper="special", fonts=fonts, width=width, height=height, pointsize=pointsize[1], ...) } if(trellis)trellis.par.set(list(fontsize=list(text=pointsize[1], points=pointsize[2]))) } `renum.fun` <- function(from.prefix="g", to.prefix="g",from=20:7, to=21:8, doit=F){ path <- getwd() pathtag <- strsplit(path, "/ch", fixed=TRUE)[[1]] endbit <- pathtag[length(pathtag)] from.prefix <- paste(from.prefix, endbit, sep="") to.prefix <- paste(to.prefix, endbit, sep="") for(i in 1:length(to)) {txt<-paste(to.prefix,".",to[i]," <- ", from.prefix,".",from[i],sep="") if(doit)eval(parse(text=txt),envir=sys.frame(0)) print(txt) } } `renum.files` <- function(from.prefix="~/r-book/ed2/Art/", to.prefix="~/r-book/ed2/Art/", from=20:7, to=21:8, doit=F){ path <- getwd() pathtag <- strsplit(path, "/ch", fixed=TRUE)[[1]] endbit <- pathtag[length(pathtag)] if(nchar(endbit)==2)chap <- paste(endbit) else chap <- paste("0",endbit,sep="") from.prefix <- paste(from.prefix, chap, "-", sep="") to.prefix <- paste(to.prefix, chap, "-", sep="") for(i in 1:length(from)){ if (from[i]<=9) ltext <- paste("0",from[i],sep="") else ltext <- paste(from[i]) if (to[i]<=9) rtext <- paste("0",to[i],sep="") else rtext <- paste(to[i]) txt<-paste("mv ", from.prefix, ltext, ".eps", " ", to.prefix, rtext, ".eps", sep="") backup<-paste("cp ", from.prefix, ltext, ".eps", " ", "archive", sep="") if(doit)system(backup) if(doit)system(txt) print(backup) print(txt) } }