g5.1 <- function(y = roller$depression, x = roller$weight, device="") { if(device!="")hardcopy(width=2.25, height=2.25, device=device) oldpar <- par(mar = c(4.1, 4.1, 1.1, 1.1), mgp=c(2.5,0.75,0)) on.exit(par(oldpar)) titl <- paste("Lawn depression (mm) versus roller weight (t).", sep = "") roller.obj <- lm(y ~ x) yhat <- predict(roller.obj) ymax <- max(c(y, yhat)) plot(x, y, xlab = "Roller weight (t)", ylab = "Depression in lawn (mm)", pch = 4, xlim=c(0, max(x)), ylim=c(0, ymax)) abline(roller.obj) b <- summary(roller.obj)$coef options(digits=3) print(anova(roller.obj)) cat("\n\nCoefficients\n\n") print(b) topleft <- par()$usr[c(1, 4)] chw <- par()$cxy[1] chh <- par()$cxy[2] legend(topleft[1], topleft[2]+0.25*chh,pch=c(1,4), legend=c("Fitted values", "Data values"), adj=0,cex=0.8, x.intersp=0.8, y.intersp=0.8, bty="n") r <- cor(x, y) bottomright <- par()$usr[c(2,3)] text(bottomright[1] - chw/2, bottomright[2]+0.5*chh, paste("r =", format(round(r, 3))), adj = 1,cex=0.8) here <- y > yhat z <- as.vector(rbind(y[here], yhat[here], rep(NA, sum(here)))) zx <- as.vector(rbind(x[here], x[here], x[here])) lines(zx, z) here <- y < yhat z <- as.vector(rbind(y[here], yhat[here], rep(NA, sum(here)))) zx <- as.vector(rbind(x[here], x[here], x[here])) lines(zx, z, lty = 2) n <- length(y) ns <- min((1:n)[y - yhat >= 0.75*max(y - yhat)]) ypos <- 0.5 * (y[ns] + yhat[ns]) chw <- par()$cxy[1] text(x[ns] - 0.25*chw, ypos, "+ve residual", adj = 1,cex=0.8) points(x, yhat, pch = 1) ns <- (1:n)[y - yhat == min(y - yhat)][1] ypos <- 0.5 * (y[ns] + yhat[ns]) text(x[ns] + 0.4*chw, ypos, "-ve residual", adj = 0,cex=0.8) titl <- paste("Lawn depression for various weights of roller,", "with fitted line.") cat("\n", titl, "\n") if(device!="")dev.off() invisible() } g5.10 <- function(sd = 2, npoints=5, nrep=4, slope=0.8, confint = F, stats = F, device="", type = "") { if(device!="")hardcopy(width=4.5, height=2.6, device=device) figtxt <- paste("Test for Linear Trend, versus Multiple Comparisons. The", "\nsix panels are six simulated results from a straight", "\nline model with slope ", slope,", sd=", sd, ", and ", nrep, " reps per level.", sep="") oldpar <- par(mfrow = c(2, 3), mar = par()$mar - c(1, 0, 1, 0), mgp = c(2, 0.5, 0), oma=c(0,0,.5,0)) on.exit(par(oldpar)) x <- rep(1:npoints, rep(nrep, npoints)) for(i in 1:6) { y <- 100 + slope * x + rnorm(20, 0, sd) stripchart(y~x, xlab = "Treatment level", ylab = "Response", vertical=TRUE, cex=1.25, xlim=c(.5,5.5)) u <- lm(y ~ factor(x)) z <- summary.aov(u) p.aov <- z[[1]][1,"Pr(>F)"] u <- lm(y ~ x) z1 <- summary(u) p.slope <- z1$coef[2, 4] if(stats) print(z1) if(i==1)mtext(side=3, line=2, "p-values are -", adj=0, cex=0.8) mtext(side = 3, line = .25, paste("Linear trend:", format(p.slope, digits = 2), "\naov: ", format(p.aov, digits = 2), sep = ""), adj = 1,cex=.8) } cat(figtxt, "\n") if(device!="")dev.off() invisible() } g5.11 <- function(device=""){ if(device!="")hardcopy(width=2.25, height=2.25, device=device) simulate.linear(sd = 2, npoints=5, nrep=4, nsets=200, seed=21, col="gray30") if(device!="")dev.off() invisible() } g5.12 <- function(df = houseprices, m = 3, x = "area", y = "sale.price", seed=29, dots=F, device="") { if(device!="")hardcopy(width=2.25, height=2.25, device=device) if(!is.null(seed))set.seed(seed) oldpar<-par(mar=c(4.6, 4.6, 2.1, 1.1), mgp=c(2.5, 0.75, 0)) on.exit(par(oldpar)) coltypes <- c(2, 3, 6) ltypes <- 1:3 ptypes <- 2:4 options(digits=3) n <- dim(df)[1] rand <- sample(n)%%m+1 xv <- df[, x] yv <- df[, y] plot(xv, yv, xlab = "Floor area", ylab = "Sale price", type = "n") xval <- pretty(xv, n = 20) houseprices.lm <- lm(yv ~ xv, data=df) print(anova(houseprices.lm)) cat("\n") sumss<-0 sumdf<-0 par(lwd=2) for(i in sort(unique(rand))) { cat("\nfold", i, "\n") n.in <- (1:n)[rand != i] n.out <- (1:n)[rand == i] cat("Observations in test set:", n.out, "\n") ab <- lm(yv ~ xv, subset = n.in)$coef z <- xv[n.out] points(xv[n.out], yv[n.out], col=coltypes[i], pch = ptypes[i], cex = 1.1) if(dots) points(xv[n.out], yv[n.out], col = coltypes[i], pch = 16) pred <- ab[1] + ab[2] * z resid <- yv[n.out] - pred xy <- data.frame(rbind(z,pred, yv[n.out], resid ), row.names=c("Floor area","Predicted price", "Observed price","Residual")) yval <- ab[1] + ab[2] * xval lines(xval, yval, lwd = 2, col = coltypes[i], lty=ltypes[i]) num <- length(n.out) print(xy,collab=rep("",num)) ss <- sum(resid^2) sumss<-sumss+ss sumdf<-sumdf+num ms <- ss/num cat("\nSum of squares =", round(ss, 2), " Mean square =", round(ms, 2), " n =", num, "\n") } print(c("Overall ms"=sumss/sumdf)) topleft<-par()$usr[c(1,4)] + c(-1, 1.75)*par()$cxy par(lwd=1, xpd=T) legend(topleft[1],topleft[2],legend=paste("Fold",1:m," "),pch=ptypes, lty=ltypes,col=coltypes, cex=0.75, horiz=T, bty="n", xjust=0) par(col = 1, xpd=F) if(device!="")dev.off() } g5.13 <- function(device=""){ if(device!="")hardcopy(width=4, height=2, device=device) oldpar <- par(mar = c(4.1,4.1,1.6,2.1), mgp=c(2.5,0.75,0)) on.exit(par(oldpar)) par(mfrow=c(1,2)) houseprices2.fn<-function (houseprices,index) { house.resample<-houseprices[index,] house.lm<-lm(sale.price~area,data=house.resample) houseprices$sale.price-predict(house.lm,houseprices) # resampled prediction errors } require(boot) set.seed(1111) n<-length(houseprices$area) R<-200 houseprices.lm<-lm(sale.price~area,data=houseprices) houseprices2.boot<-boot(houseprices,R=R,statistic=houseprices2.fn) house.fac<-factor(rep(1:n,rep(R,n))) plot(house.fac,as.vector(houseprices2.boot$t), ylab="", xlab="House") mtext(side=2, line=3, "Prediction Errors", cex=0.9) mtext(side = 3, line = 0.25, "A", adj = 0) boot.se <- apply(houseprices2.boot$t,2,sd) model.se <- predict.lm(houseprices.lm,se.fit=T)$se.fit plot(boot.se/model.se, ylab="", xlab="House",pch=16) mtext(side=2, line=2.5, "Ratio of SEs\nBootstrap to Model-Based", cex=0.9) mtext(side = 3, line = 0.25, "B", adj = 0) abline(1,0) if(device!="")dev.off() invisible() } g5.14 <- function(tcol = 2, col="gray40", device="") { if(device!="")hardcopy(width=4.4, height=3.2, device=device) titl <- paste("The above panels show (as solid curves) some alternative", "\npatterns of response. The dotted line shows the power", "\nfamily transformation that makes the response linear.", "\nThe formula for this power family transformation appears", "\nabove the graph." ) oldpar <- par(mfrow = c(2, 3), mar = par()$mar - c( 1, 1, 1.0, 1), mgp = c(1.5, 0.5, 0), oma=c(0,1,0,1)) on.exit(par(oldpar)) powerplot(expr="sqrt(x)", xlab="", col=col) powerplot(expr="x^0.25", xlab="", ylab="", col=col) powerplot(expr="log(x)", xlab="", ylab="", col=col) powerplot(expr="x^2", col=col) powerplot(expr="x^4", ylab="", col=col) powerplot(expr="exp(x)", ylab="", col=col) if(device!="")dev.off() par(mfrow=c(1,1)) invisible() } g5.15 <- function(yvar = "heart", xvar="weight", species = "Cape fur seal", dset = cfseal, stats = F, device="") { if(device!="")hardcopy(width=2.25, height=2.25, device=device) library(DAAG); data(cfseal) pretext <- c(heart = "Heart weight", brain = "Brain weight", liver = "Liver weight", weight="Body weight")[yvar] xtext <- c(heart = "Heart weight", brain = "Brain weight", liver = "Liver weight", weight="Body weight")[xvar] oldpar <- par(mar = c(4.1,4.1,2.1,1.1), mgp=c(2.5,0.75,0)) on.exit(par(oldpar)) x <- log(dset[,xvar]) y <- log(dset[, yvar]) ylim <- range(y) if (yvar == "heart") ylim<- log(c(82.5,1100)) xlim <- range(x) if(xvar == "weight") xlim <- log(c(17,180)) ylab <- switch(yvar, heart = "Heart weight (g, log scale)", brain = "Brain weight (g, log scale)", liver = "Liver weight (g, log scale)", weight = "Body weight (kg, log scale)") xlab <- switch(xvar, heart = "Heart weight (g, log scale)", brain = "Brain weight (g, log scale)", liver = "Liver weight (g, log scale)", weight = "Body weight (kg, log scale)") xtik <- pretty(exp(x),10) ytik <- pretty(exp(y), 10) xtik <- xtik[xtik >= exp(xlim[1]) & xtik <= exp(xlim[2])] ytik <- ytik[ytik >= exp(ylim[1]) & ytik <= exp(ylim[2])] plot(x, y, xlab = xlab, ylab = ylab, axes = F, xlim = xlim, ylim = ylim, pch = 16, cex=0.85) axis(1, at = log(xtik), labels = paste(xtik)) axis(2, at = log(ytik), labels = paste(ytik)) box() form1 <- formula(y ~ x) u <- lm(form1, data = dset) abline(u$coef[1], u$coef[2]) usum <- summary(u)$coef options(digits=3) print(usum) if(yvar == "brain") { u2 <- lm(form1, data = dset, subset = brain < log(1.7)) redrange <- range(x[y < max(y)]) yval <- u2$coef[1] + u2$coef[2] * redrange lines(redrange, yval, lty = 2) if(stats) print(summary(u2)) } cwh <- par()$cxy eqn <- paste("log y =", round(usum[1, 1], 2), " [SE ", round(usum[1, 2], 3), "] + ", round(usum[2, 1], 3), " [", round(usum[2, 2], 2), "] log x") mtext(side=3, line=0.25, eqn, adj = 0.5, cex = 0.75) figtxt <- paste(pretext, " versus ", xtext, " for ", length(y), " members \n of the species ", species, ".", collapse = "") if(yvar == "logbrain") figtxt <- paste(figtxt, " The dotted line shows the", "\neffect of omitting the data point", "for the largest animal.") cat(figtxt, "\n") if(device!="")dev.off() if(stats) summary(u) } g5.16 <- function(dset1 = pair65, dset2 = leafshape, device="", col="gray30") { if(device!="")hardcopy(width=4.2, height=2.1, device=device) oldpar <- par(mfrow = c(1, 2), mar = c(4.1,4.1,2.1,2.1), mgp=c(2.5,.75,0), pty="s") library(DAAG); data(leafshape) on.exit(par(oldpar)) x <- dset1$ambient y <- dset1$heated ylab <- "Stretch (heated band)" plot(x, y, xlab = "Stretch (band held at ambient)", ylab = ylab, pch = 16) topleft <- par()$usr[c(1, 4)] + c(0.5, -0.5) * par()$cxy text(topleft[1], topleft[2], paste("r =", round(cor(x, y), 2)), adj = 0) mtext(side = 3, line = 0.23, "A", adj = 0) u1 <- lm(heated ~ ambient, data = dset1) abline(u1$coef[1], u1$coef[2]) u2 <- lm(ambient ~ heated, data = dset1) abline( - u2$coef[1]/u2$coef[2], 1/u2$coef[2], lty = 2) x <- dset2$logpet y <- dset2$loglen plot(x, y, xlab = "Petiole length (mm)", ylab = "Leaf length (mm)", axes = F, type="n") points(x, y, cex=0.65, col=col) topleft <- par()$usr[c(1, 4)] + c(0.5, -0.5) * par()$cxy text(topleft[1], topleft[2], paste("r =", round(cor(x, y), 3)), adj = 0) mtext(side = 3, line = 0.23, "B", adj = 0) xlab <- pretty(exp(x)) xlab <- c(0.5, 1, xlab) xlab <- xlab[xlab > 0] ylab <- pretty(exp(y)) ylab <- ylab[ylab > 0] axis(1, at = log(xlab), labels = paste(xlab)) axis(2, at = log(ylab), labels = paste(ylab)) box() u1 <- lm(y ~ x) u2 <- lm(x ~ y) abline(u1$coef[1], u1$coef[2]) abline( - u2$coef[1]/u2$coef[2], 1/u2$coef[2], lty = 2) figtxt <- paste("Plot (A) is for a dataset that compared the stretch", "\nof an elastic band that was heated with that for an", "\nelastic band that was held at ambient temperature,", "\nwhile (B) is for a leaf dataset. In each plot we", "\nshow the regression line for y on x (solid line),", "\nand the regression line for x on y (dotted line).", "\nIn (A) the lines are quite similar, while in (B)", "\nwhere the correlation is smaller, the lines are", "\nquite different.") cat(figtxt, "\n") if(device!="")dev.off() } g5.2 <- function(y = roller$depression, x = roller$weight, curve = c("reg"), fignum=5.2, device="") { if(device!="")psfile(fignum=fignum, width=4, height=2.2, device=device) titl <- paste("Diagnostic plots for the regression of Figure 5.1:", "\nPanel (A) plots residuals against fitted values.", "\nPanel (B) is a normal probability plot of residuals.", "\nIf residuals follow a normal distribution", "the points \nshould fall close to a line.") oldpar <- par(mfrow = c(1, 2), mgp = c(2.25, 0.5, 0), mar = par()$mar - c(1, 0.5, 0.5, 0)) on.exit(par(oldpar)) u <- lm(depression ~ weight, data = roller) res <- resid(u) hat <- fitted(u) plot(hat, res, xlab = "Fitted values", ylab = "Residuals", pch = 1) points(hat, res, pch = 1, cex = 0.8, lwd = 2) mtext(side = 3, line = 0.25, "A", adj = 0) abline(h = 0, lty = 2) qqnorm(res, xlab = "Quantiles of Normal Distribution", ylab = "Ordered data values", pch = 16, cex.main=1.15) mtext(side = 3, line = 0.25, "B", adj = 0) cat(titl, "\n") if(device!="")dev.off() invisible(u) } g5.3 <- function(fignum=5.3, device=""){ if(device!="")psfile(fignum=fignum, width=4.8, height=2.75, device=device, pointsize=9) if(!exists("roller.lm"))roller.lm <- lm(depression ~ weight, data=roller) qreference(residuals(roller.lm), nrep=8, nrows=2, xlab="") if(device!="")dev.off() } g5.4 <- function(y = ironslag$chemical, x = ironslag$magnetic, device="") { if(device!="")hardcopy(width=3.75, height=3.25, device=device, pointsize=9) oldpar <- par(mfrow = c(2, 2), mar = c(4.8, 4.1, 2.6,2.1), mgp = c(2.5, 0.75, 0), cex = 0.85, mex = 0.85) on.exit(par(oldpar)) mtext3 <- function(item="A", txt=leg[1], xleft=par()$usr[1]+1.25*par()$cxy[1]){ mtext(side = 3, line = 0.25, item, adj = 0) mtext(side = 3, line = 0.25, txt, cex = 0.7, at=xleft, adj = 0) } titl <- paste(reg = "Chemical vs Magnetic Test of Iron Content in Slag.", "\nThe fitted curves used the Splus loess smooth.", "\nIn (d) the downward slope suggests a lower variance", "\nfor larger fitted values." ) leg <- c("Fitted line and loess curve", "Residuals from line, with smooth", "Observed vs predicted, for fitted line", "Is variance constant about line?") plot(x, y, xlab = "Magnetic", ylab = "Chemical", type="n") u <- lm(y ~ x) abline(u$coef[1], u$coef[2]) yhat <- predict(u) lines(x, yhat) mtext3(xleft=par()$usr[1]+1.25*par()$cxy[2]) print(panel.smooth(x, y, span = 0.95, lty = 2, lwd = 1.5, pch=0)) res <- residuals(u) plot(x, res, xlab = "Magnetic", ylab = "Residual", type = "n") mtext3(item="B", txt=leg[2], xleft=par()$usr[1]+1.25*par()$cxy[1]) print(panel.smooth(x, res, span = 0.95)) points(x, res, pch = 1, cex = 0.9, lwd = 1) hat <- fitted(u) plot(hat, y, xlab = "Predicted chemical", ylab = "Chemical", type = "n") print(panel.smooth(hat, y, span = 0.95)) mtext3(item="C", txt=leg[3], xleft=par()$usr[1]+1.25*par()$cxy[1]) yabs <- sqrt(abs(res)) plot(hat, yabs, xlab = "Predicted chemical", ylab = "", type = "n") mtext(side = 2, line = 2.5, expression(sqrt(abs(residual)))) mtext3(item="D", txt=leg[4], xleft=par()$usr[1]+1.25*par()$cxy[1]) print(panel.smooth(hat, yabs, span = 0.95)) cat(titl, "\n") if(device!="")dev.off() invisible() } g5.5 <- function(y = ironslag$chemical, x = ironslag$magnetic, device="") { if(device!="")hardcopy(width=4, height=2, device=device) oldpar <- par(mfrow = c(1, 2), mar = c(4.6, 4.1,2.1,1.6), mgp = c( 2.5, 0.5, 0)) on.exit(par(oldpar)) leg <- c("A. Residuals from fitted loess curve.", "B. Is variance about curve constant?") u <- loess(y ~ x) resval <- residuals(u) yhat <- predict(u) yabs <- sqrt(abs(resval)) plot(x, resval, xlab = "Magnetic", ylab = "Residual", pch = 1) points(x, resval, cex = 0.8, type="n") abline(h=0,lty=2) mtext(side = 3, line = 0.25, leg[1], adj = 0, cex=0.7) plot(yhat, yabs, xlab = "Predicted chemical", ylab = expression(sqrt(abs(residual))), type = "n") panel.smooth(yhat, yabs, span = 1.1, cex = 1.1) mtext(side = 3, line = 0.25, leg[2], adj = 0, cex=0.7) if(device!="")dev.off() invisible() } g5.6 <- function(y = softbacks$weight, x = softbacks$volume, curve = c("reg"), show.fits = T, device="") { if(device!="")hardcopy(width=2.25, height=2.25, device=device) titl <- switch(curve[1], reg = paste("Weight versus volume for softcover books,", "\nwith fitted line."), lo = paste( "Weight versus volume for softcover books,", "\nwith fitted line and S-PLUS loess smooth curve.")) oldpar <- par(mar = c(4.1,4.1,1.1,1.1), mgp = c(2.5, 0.5, 0)) on.exit(par(oldpar)) u <- lm(weight ~ volume, data = softbacks) cat("\nCoefficients\n\n") options(digits=3) print(summary(u)$coef) cat("\n\n") print(anova(u)) yhat <- predict(u) r <- cor(x, y) xlim <- range(x) xlim[2] <- xlim[2]+diff(xlim)*0.08 plot(x, y, xlab = "Volume (cc)", xlim=xlim, ylab = "Weight (g)", pch = 1, ylim = range(c(y, yhat))) if(match("reg", curve, nomatch = 0)) { abline(u$coef[1], u$coef[2], lty = 1) z <- summary(u)$coef if(show.fits) { points(x, yhat, pch = 1, cex = 0.75) res <- resid(u) for(i in 1:length(res)) { resi <- res[i] izzy <- as.numeric(resi > 0) xi <- x[i] yhati <- yhat[i] yi <- y[i] lines(rep(xi, 2), c(yhati, yi), col=2-izzy, lty=2-izzy) eps <- par()$cxy[1] * 0.2 if(i == 6) { adji <- 1 eps <- - eps } else adji <- 0 text(xi + eps, yhati + resi/2, paste(round(resi, 1)), adj = adji, cex = 0.65) } } bottomright <- par()$usr[c(2, 3)] chw <- par()$cxy[1] chh <- par()$cxy[2] btxt <- c(paste("a =", format(round(z[1, 1], 3)), " SE =", format(round(z[1, 2], 3))), paste("b =", format(round(z[2, 1], 3)), " SE =", format(round(z[2, 2], 3)))) legend(bottomright[1], bottomright[2], legend=btxt, xjust=1, yjust=0, cex=0.8, bty="n") } cat("\n",titl,"\n") if(device!="")dev.off() invisible(u) } g5.7 <- function(device=""){ 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") require(DAAG); data(softbacks) softbacks.lm<-lm(weight~volume,data=softbacks) on.exit(par(oldpar)) plot(softbacks.lm, which=1:4, caption= c("A: Residuals vs Fitted", "B: Normal Q-Q", "C: Scale-Location", "D: Cook's distance")) figtxt <- "Diagnostic plots for previous figure." cat(figtxt, "\n") if(device!="")dev.off() invisible(softbacks.lm) } g5.8 <- function(y = roller$depression, x = roller$weight, pred = T, device="") { if(device!="")psfile(width=2, height=2, device=device) titl <- paste("Lawn Depression, for Various Weights of Roller,", "\nwith fitted line and showing 95% pointwise", "\nconfidence bounds for points on the fitted line.") oldpar <- par(mar = c(4.1, 4.1, 1.1, 1.1), mgp = c(2.5, 0.5, 0) ) on.exit(par(oldpar)) xlim <- c(0, max(x)) ylim <- c(0, max(y)) plot(x, y, xlab = "Weight of Roller (tonnes)", xlim=xlim, ylim=ylim, ylab = "Depression in Lawn (mm)", pch = 16) u <- lm(depression ~ weight, data = roller) abline(u$coef[1], u$coef[2], lty = 1) z <- summary(u)$coef y12 <- par()$usr[3:4] if(pred) { assign("xy", data.frame(weight = pretty(roller$weight, 20))) hat <- predict(u, newdata = xy, interval="confidence") ci<-data.frame(lower=hat[,"lwr"],upper=hat[,"upr"]) here <- ci$lower >= y12[1] lines(xy$weight[here], ci$lower[here], lty = 2,lwd=2,col="grey") here <- ci$upper < y12[2] lines(xy$weight[here], ci$upper[here], lty = 2,lwd=2,col="grey") } cat(titl, "\n") if(device!="")dev.off() invisible(u) } g5.9 <- function(color = F, device="") { if(device!="")hardcopy(width=3.4, height=2, device=device) figtxt<-paste("Two rubber band experiments, with different ranges", "\nof x-values. The dashed curves are pointwise 95%", "\nconfidence bounds for the fitted line.") panelci<-function(data,...) { nrows<-list(...)$nrows ncols<-list(...)$ncols if(ncols==1)axis(2) if(ncols==1)axis(1) else axis(3) x<-data$stretch; y<-data$distance u <- lm(y ~ x) upred <- predict(u, interval="confidence") ci <- data.frame(fit=upred[,"fit"],lower=upred[,"lwr"], upper=upred[,"upr"]) ord<-order(x) lines(x[ord], ci$fit[ord], lty=1, lwd=2) lines(lowess(x[ord], ci$upper[ord]), lty=2, lwd=2, col="grey") lines(lowess(x[ord], ci$lower[ord]), lty=2, lwd=2, col="grey") } xy<-rbind(elastic2,elastic1) nam <- c("Range of stretch 30-60 mm","Range of stretch 42-54 mm") trial<-rep(nam, c(dim(elastic2)[1],dim(elastic1)[1])) xlim<-range(elastic2$stretch) ylim<-range(elastic2$distance) xy<-split(xy,trial) xy<-lapply(1:length(xy),function(i){c(as.list(xy[[i]]), list(xlim=xlim, ylim=ylim))}) # xy<-lapply(xy,function(z){z$xlim<-xlim; z$ylim<-ylim; z}) names(xy) <- nam panelplot(xy,panel=panelci,totrows=1,totcols=2, par.strip.text=list(cex=.8), oma=c(4,4,2.5,2)) cat(figtxt, "\n") mtext(side = 2, line = 2.75, "Distance moved (cm)") mtext(side=1,line=3.5,"Amount of stretch (mm)") if(device!="")dev.off() } gdump <- function(fnam=NULL, prefix="~/r-book/ed2/figures/figs", 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)), "hardcopy") cat("\nDump to file:", fnam, "\n") print(objnames) dump(objnames, fnam) } gfile <- function(width=3.75, height=3.75, color=F, trellis=F, device=c("","pdf","ps"), path="", pointsize=c(8,5), horiz=F){ ## 1 x 1: 2.25" x 2.25" ## 2 x 2: 2.75" x 2.75" ## 3 x 3: 3.75" x 3.75" or 3.25" x 3.25" for simple scatterplots ## 1 x 2: 4" x 2.25" ## 2 x 3: 4" x 2.8" ## 3 x 4: 4.5" x 3.25 if(!trellis)pointsize <- pointsize[1] funtxt <- sys.call(1) fnam <- strsplit(as.character(funtxt), "(", fixed=T)[[1]][1] dotsplit <- strsplit(fnam, "\\.")[[1]] dotsplit[1] <- substring(dotsplit[1], 2) prefix1 <- paste(if(nchar(dotsplit[1])==1)"0" else "", dotsplit[1], sep="") prefix2 <- paste(if(nchar(dotsplit[2])==1)"0" else "", dotsplit[2], sep="") if(device=="")stop("No device has been specified") suffix <- switch(device, ps=".eps", pdf=".pdf") fnam <- paste("~/r-book/second/Art/",prefix1,"-",prefix2, suffix, sep="") print(fnam) dev.out <- device[1] dev.fun <- switch(dev.out, pdf=pdf, ps=postscript) if(trellis){ library(lattice) trellis.device(file=fnam, device=dev.fun, color = color, width=width, height=height, horiz=horiz) trellis.par.set(fontsize=list(text=pointsize[1], points=pointsize[2])) } else if (dev.out!=""){ print(c(width, height)) dev.fun(file=fnam, paper="special", enc="MacRoman", horiz=horiz, width=width, height=height, pointsize=pointsize[1]) } } 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]))) }