RDX2 X   g5.1 source3 >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)) K titl <- paste("Lawn depression (mm) versus roller weight (t).", sep = "")  roller.obj <- lm(y ~ x)  yhat <- predict(roller.obj)  ymax <- max(c(y, yhat)) 0 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] 4 legend(topleft[1], topleft[2]+0.25*chh,pch=c(1,4), 2 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)] 6 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)))) 3 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)))) 3 zx <- as.vector(rbind(x[here], x[here], x[here]))  lines(zx, z, lty = 2)  n <- length(y) 2 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) A titl <- paste("Lawn depression for various weights of roller,", $ "with fitted line.")  cat("\n", titl, "\n")  if(device!="")dev.off() invisible() }ў§ y $ roller depressionў xџџ weightў device ў { if != џ ў hardcopy width@ height@ џ џўў <- oldpar par mar c@ffffff@ffffff?ё™™™™™š?ё™™™™™šў mgpџ@?шўўў on.exitџџўўџ titl paste .Lawn depression (mm) versus roller weight (t). sep ўўџ roller.obj lm ~џџўўўџ yhat predictџўўџ ymax maxџџџўўў plotџџ xlab Roller weight (t) ylab Depression in lawn (mm) pch@ xlimџ џџўў ylimџџўў ablineџўџ bџ summaryџў coefўў options digits@ў print anovaџўў cat  Coefficients ў-џ(џўџ topleft [џџў usrўџ?№@ўўўџ chw1џџџў cxyў?№ўўџ chh1џџџў4џў@ўў legend1џ0џ?№ў +1џ0џ@ў *?а5џўў$џџ?№@ў6џџ Fitted values Data valuesў adj cex?щ™™™™™š x.intersp?щ™™™™™š y.intersp?щ™™™™™š bty nўџ r corџџўўџ bottomright1џџџў2џўџ@@ўўў text -1џ@џ?№ў /3џ@ўў7џ1џ@џ@ў8џ?р5џўўџ r = format round>џ@ўўў9џ?№:џ?щ™™™™™šўџ here >џџўўџ z as.vector rbind1џџFџў1џџFџў rep € sumFџўўўўўџ zxIџJџ1џџFџў1џџFџў1џџFџўўўў linesMџHџўџFџ <џџўўџHџIџJџ1џџFџў1џџFџўKџ €LџFџўўўўўџMџIџJџ1џџFџў1џџFџў1џџFџўўўўNџMџHџ lty@ўџ n lengthџўўџ ns min1џ ( :?№Qџўў >=Bџџџў8џ?ш џBџџџўўўўўўўџ ypos8џ?рUџ7џ1џџSџў1џџSџўўўўўџ3џ1џџџў4џў?№ўўAџBџ1џџSџў8џ?а3џўўXџ +ve residual9џ?№:џ?щ™™™™™šў pointsџџ$џ?№ўџSџ1џ1џUџVџ?№Qџўў ==BџџџўTџBџџџўўўў?№ўўџXџ8џ?рUџ7џ1џџSџў1џџSџўўўўўAџ7џ1џџSџў8џ?й™™™™™š3џўўXџ -ve residual9џ:џ?щ™™™™™šўџџџ .Lawn depression for various weights of roller, with fitted line.ўў/џ  џ  ў џ џ џ ў dev.offўў invisibleўў g5.10џ" .function(sd = 2, npoints=5, nrep=4, slope=0.8, 8 confint = F, stats = F, device="", type = "") { > if(device!="")hardcopy(width=4.5, height=2.6, device=device) L figtxt <- paste("Test for Linear Trend, versus Multiple Comparisons. The", K "\nsix panels are six simulated results from a straight", A "\nline model with slope ", slope,", sd=", sd, = ", and ", nrep, " reps per level.", sep="") A oldpar <- par(mfrow = c(2, 3), mar = par()$mar - c(1, 0, 1, 0), 4 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) 5 stripchart(y~x, xlab = "Treatment level", ylab = C "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) C if(i==1)mtext(side=3, line=2, "p-values are -", adj=0, cex=0.8) 7 mtext(side = 3, line = .25, paste("Linear trend:", > format(p.slope, digits = 2), "\naov: ", K format(p.aov, digits = 2), sep = ""), adj = 1,cex=.8)  }  cat(figtxt, "\n")  if(device!="")dev.off() invisible() }ў§ sd@ npoints@ nrep@ slope?щ™™™™™š confint F statscџ џ  type ў џ џ џ џ ў џџ@џ@ЬЬЬЬЬЭ џ џўўџ figtxtџ 7Test for Linear Trend, versus Multiple Comparisons. The 5 six panels are six simulated results from a straight  line model with slope aџ , sd=^џ , and `џ  reps per level.џ ўўџџџ mfrowџ@@ўџBџџџўџўџ?№?№ўўџџ@?рў omaџ?рўўўџџџўўџџKџVџ?№_џўKџ`џ_џўўў for iVџ?№@ў џџџ7џ7џ@Y8џaџџўў rnorm@4^џўўў stripchartџџџў"џ Treatment level#џ Response vertical :џ?є%џџ?р@ўўџ uџџџ factorџўўўўџHџ summary.aovnџўўџ p.aov1џ [[Hџ?№ў?№ Pr(>F)ўўџnџџџџџўўўџ z1)џnџўўџ p.slope1џџsџ*џў@@ўў џdџ-џsџўў џZџjџ?№ў mtext side@ line@ p-values are -9џ:џ?щ™™™™™šўўuџvџ@wџ?аџ Linear trend:Dџtџ,џ@ў  aov: Dџqџ,џ@ўџ ў9џ?№:џ?щ™™™™™šўўў/џfџ  ў џ џ џ ў[џўў\џўў g5.11џ function(device=""){ B if(device!="")hardcopy(width=2.25, height=2.25, device=device) P simulate.linear(sd = 2, npoints=5, nrep=4, nsets=200, seed=21, col="gray30")  if(device!="")dev.off()  invisible()  }ў§ џ ў џ џ џ џ ў џџ@џ@ џ џўў simulate.linear^џ@_џ@`џ@ nsets@i seed@5 col gray30ў џ џ џ ў[џўў\џўў g5.12џ9 Hfunction(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] D 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] 1 cat("Observations in test set:", n.out, "\n") ) ab <- lm(yv ~ xv, subset = n.in)$coef  z <- xv[n.out] M 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 3 xy <- data.frame(rbind(z,pred, yv[n.out], resid I ), 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)) 4 topleft<-par()$usr[c(1,4)] + c(-1, 1.75)*par()$cxy  par(lwd=1, xpd=T) G legend(topleft[1],topleft[2],legend=paste("Fold",1:m," "),pch=ptypes, F lty=ltypes,col=coltypes, cex=0.75, horiz=T, bty="n", xjust=0)  par(col = 1, xpd=F)  if(device!="")dev.off() }ў§ df houseprices m@џ areaџ sale.price{џ@= dotscџ џ ў џ џ џ џ ў џџ@џ@ џ џўў џ ! is.null{џўў set.seed{џўўџџџџџ@ffffff@ffffff@ЬЬЬЬЬЭ?ё™™™™™šўџџ@?шўўўџџџўўџ coltypesџ@@@ўўџ ltypesVџ?№@ўўџ ptypesVџ@@ўў+џ,џ@ўџQџ1џ dim~џў?№ўўџ rand7џ %% sampleQџў€џў?№ўўџ xv1џ~џћџўўџ yv1џ~џћџўў!џŒџџ"џ Floor area#џ Sale priceeџ nўџ xval prettyŒџQџ@4ўўџ houseprices.lmџџџŒџў data~џўў-џ.џџўў/џ  ўџ sumssўџ sumdfўџ lwd@ўiџjџ sort unique‰џўў џ/џ  foldjџ  ўџ n.in1џUџVџ?№Qџўў џ‰џjџўўўџ n.out1џUџVџ?№QџўўZџ‰џjџўўў/џ Observations in test set:˜џ  ўџ abџџџџŒџў subset—џў*џўўџHџ1џŒџ˜џўўYџ1џŒџ˜џў1џџ˜џў|џ1џ…џjџў$џ1џ‡џjџў:џ?ё™™™™™šў џџYџ1џŒџ˜џў1џџ˜џў|џ1џ…џjџў$џ@0ўўџ pred7џ1џ™џ?№ў8џ1џ™џ@ўHџўўўџ residBџ1џџ˜џў›џўўџ xy data.frameJџHџ›џ1џџ˜џўœџў row.namesџ Floor area Predicted price Observed price Residualўўўџ yval7џ1џ™џ?№ў8џ1џ™џ@ўŽџўўўNџŽџ џ”џ@|џ1џ…џjџўPџ1џ†џjџўўџ numRџ˜џўў-џџ collabKџ Ёџўўџ ssLџ ^œџ@ўўўџ’џ7џ’џЃџўўџ“џ7џ“џЁџўўџ msCџЃџЁџўў/џ  Sum of squares =EџЃџ@ў  Mean square =EџЅџ@ў  n =Ёџ  ўўў-џџ Overall msCџ’џ“џўўўџ0џ7џ1џџџў2џўџ?№@ўў8џџBџ?№ў?ќўџџў4џўўўўџ”џ?№ xpd Tў6џ1џ0џ?№ў1џ0џ@ў6џџ FoldVџ?№€џў  ў$џ‡џPџ†џ|џ…џ:џ?ш horizЈџ=џ n xjustўџ|џ?№Їџcџў џ џ џ ў[џўўў 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)) 1 houseprices2.fn<-function (houseprices,index)  { + house.resample<-houseprices[index,] 9 house.lm<-lm(sale.price~area,data=house.resample) < houseprices$sale.price-predict(house.lm,houseprices) E # resampled prediction errors  }  require(boot)  set.seed(1111)  n<-length(houseprices$area) R<-200 8 houseprices.lm<-lm(sale.price~area,data=houseprices) F houseprices2.boot<-boot(houseprices,R=R,statistic=houseprices2.fn) ( house.fac<-factor(rep(1:n,rep(R,n))) 2 plot(house.fac,as.vector(houseprices2.boot$t),  ylab="", xlab="House") 7 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 8 plot(boot.se/model.se, ylab="", xlab="House",pch=16) E 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()  }ў§ џ ў џ џ џ џ ў џџ@џ@ џ џўўџџџџџ@ffffff@ffffff?љ™™™™™š@ЬЬЬЬЬЭўџџ@?шўўўџџџўўџgџџ?№@ўўџ houseprices2.fn functionџћ indexћў џџ house.resample1џџЎџћўўџ house.lmџџ sale.price areaў‘џЏџўўBџџџБџўџАџџўўў function (houseprices,index)  { + house.resample<-houseprices[index,] 9 house.lm<-lm(sale.price~area,data=house.resample) < houseprices$sale.price-predict(house.lm,houseprices) E # resampled prediction errors  }ўў require bootў„џ@‘\ўџQџRџџџВџўўўџ R@iўџџџџБџВџў‘џџўўџ houseprices2.bootДџџЕџЕџ statisticЌџўўџ house.facoџKџVџ?№QџўKџЕџQџўўўў!џИџIџџЖџ tўў#џ "џ Houseўuџvџ@wџ@ Prediction Errors:џ?ьЬЬЬЬЬЭўuџvџ@wџ?а A9џўџ boot.se applyџЖџЙџў@^џўўџ model.seџ predict.lmџ se.fitЈџўОџўў!џCџКџМџў#џ "џ House$џ@0ўuџvџ@wџ@ %Ratio of SEs Bootstrap to Model-Based:џ?ьЬЬЬЬЬЭўuџvџ@wџ?а B9џў'џ?№ў џ џ џ ў[џўў\џўў g5.14џ +function(tcol = 2, col="gray40", device="") { > if(device!="")hardcopy(width=4.4, height=3.2, device=device) K titl <- paste("The above panels show (as solid curves) some alternative", J "\npatterns of response. The dotted line shows the power", K "\nfamily transformation that makes the response linear.", M "\nThe formula for this power family transformation appears", $ "\nabove the graph."  ) oldpar <- par(mfrow = c(2, 3), 3 mar = par()$mar - c( 1, 1, 1.0, 1), 5 mgp = c(1.5, 0.5, 0), oma=c(0,1,0,1))  on.exit(par(oldpar)) - powerplot(expr="sqrt(x)", xlab="", col=col) 5 powerplot(expr="x^0.25", xlab="", ylab="", col=col) 5 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() }ў§ tcol@|џ gray40 џ ў џ џ џ џ ў џџ@™™™™™šџ@ ™™™™™š џ џўўџџџ 8The above panels show (as solid curves) some alternative 6 patterns of response. The dotted line shows the power 6 family transformation that makes the response linear. 9 The formula for this power family transformation appears  above the graph.ўўџџџgџџ@@ўџBџџџўџўџ?№?№?№?№ўўџџ?ј?рўhџџ?№?№ўўўџџџўў powerplot expr sqrt(x)"џ |џ|џўСџТџ x^0.25"џ #џ |џ|џўСџТџ log(x)"џ #џ |џ|џўСџТџ x^2|џ|џўСџТџ x^4#џ |џ|џўСџТџ exp(x)#џ |џ|џў џ џ џ ў[џўўџgџџ?№?№ўў\џўў g5.15џD 'function(yvar = "heart", xvar="weight", 4 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", B 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, 7 heart = "Heart weight (g, log scale)", 7 brain = "Brain weight (g, log scale)", 7 liver = "Liver weight (g, log scale)", 8 weight = "Body weight (kg, log scale)")  xlab <- switch(xvar, 7 heart = "Heart weight (g, log scale)", 7 brain = "Brain weight (g, log scale)", 7 liver = "Liver weight (g, log scale)", 8 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])] 8 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 8 eqn <- paste("log y =", round(usum[1, 1], 2), " [SE ", B round(usum[1, 2], 3), "] + ", round(usum[2, 1], 3), 5 " [", round(usum[2, 2], 2), "] log x") 6 mtext(side=3, line=0.25, eqn, adj = 0.5, cex = 0.75) - figtxt <- paste(pretext, " versus ", xtext, D " 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) }ў§ yvar heart xvar weight species Cape fur seal dset cfsealdџcџ џ ў џ џ џ џ ў џџ@џ@ џ џўў library DAAGў‘џШџўџ pretext1џџ heart Heart weight brain Brain weight liver Liver weightџ Body weightўФџўўџ xtext1џџЬџ Heart weightЭџ Brain weightЮџ Liver weightџ Body weightўХџўўџџџџџ@ffffff@ffffff@ЬЬЬЬЬЭ?ё™™™™™šўџџ@?шўўўџџџўўџџ log1џЧџћХџўўўџџаџ1џЧџћФџўўўџ&џ rangeџўў џZџФџ heartўџ&џаџџ@T @‘0ўўўўџ%џбџџўў џZџХџ weightўџ%џаџџ@1@f€ўўўўџ#џ switchФџЬџ Heart weight (g, log scale)Эџ Brain weight (g, log scale)Юџ Liver weight (g, log scale)џ Body weight (kg, log scale)ўўџ"џвџХџЬџ Heart weight (g, log scale)Эџ Brain weight (g, log scale)Юџ Liver weight (g, log scale)џ Body weight (kg, log scale)ўўџ xtikџ expџў@$ўўџ ytikџдџџў@$ўўџгџ1џгџ &Wџгџдџ1џ%џ?№ўўў <=гџдџ1џ%џ@ўўўўўўџеџ1џеџжџWџеџдџ1џ&џ?№ўўўзџеџдџ1џ&џ@ўўўўўў!џџџ"џ"џ#џ#џ axescџ%џ%џ&џ&џ$џ@0:џ?ы333333ў axis?№ atаџгџў labelsџгџўўйџ@кџаџеџўлџџеџўў boxўџ form1 formulaџџџўўўџnџџнџ‘џЧџўў'џ1џџnџ*џў?№ў1џџnџ*џў@ўўџ usumџ)џnџў*џўў+џ,џ@ў-џпџў џZџФџ brainў џџ u2џнџ‘џЧџšџOџЭџаџ?ћ333333ўўўўџ redrangeбџ1џџOџџ џџўўўўўџ џ7џ1џџрџ*џў?№ў8џ1џџрџ*џў@ўсџўўўNџсџ џPџ@ў џdџ-џ)џрџўўўўўџ cwhџџў4џўўџ eqnџ log y =Eџ1џпџ?№?№ў@ў  [SE Eџ1џпџ?№@ў@ў ] + Eџ1џпџ@?№ў@ў  [Eџ1џпџ@@ў@ў ] log xўўuџvџ@wџ?ауџ9џ?р:џ?шўџfџџЫџ  versus Яџ  for Rџџў  members of the species Цџ . collapse ўў џZџФџ logbrainўџfџџfџ  The dotted line shows the " effect of omitting the data point for the largest animal.ўўў/џfџ  ў џ џ џ ў[џўў џdџ)џnџўўў g5.16џ5 Dfunction(dset1 = pair65, dset2 = leafshape, device="", col="gray30") { > if(device!="")hardcopy(width=4.2, height=2.1, device=device) L 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)" 6 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 A 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) 9 abline( - u2$coef[1]/u2$coef[2], 1/u2$coef[2], lty = 2)  x <- dset2$logpet  y <- dset2$loglen F 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 A 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]) 9 abline( - u2$coef[1]/u2$coef[2], 1/u2$coef[2], lty = 2) H figtxt <- paste("Plot (A) is for a dataset that compared the stretch", J "\nof an elastic band that was heated with that for an", I "\nelastic band that was held at ambient temperature,", G "\nwhile (B) is for a leaf dataset. In each plot we", H "\nshow the regression line for y on x (solid line),", H "\nand the regression line for x on y (dotted line).", G "\nIn (A) the lines are quite similar, while in (B)", F "\nwhere the correlation is smaller, the lines are", ' "\nquite different.")  cat(figtxt, "\n")  if(device!="")dev.off() }ў§ dset1 pair65 dset2 leafshape$ џ |џ gray30ў џ џ џ џ ў џџ@ЬЬЬЬЬЭџ@ЬЬЬЬЬЭ џ џўўџџџgџџ?№@ўџџ@ffffff@ffffff@ЬЬЬЬЬЭ@ЬЬЬЬЬЭўџџ@?шў pty sўўЩџЪџў‘џщџўџџџўўџџџцџ ambientўўџџџцџ heatedўўџ#џ Stretch (heated band)ў!џџџ"џ Stretch (band held at ambient)#џ#џ$џ@0ўџ0џ7џ1џџџў usrўџ?№@ўў8џџ?рBџ?рўўџџў cxyўўўўAџ1џ0џ?№ў1џ0џ@ўџ r =Eџ?џџџў@ўў9џўuџvџ@wџ?ЭpЃз =q A9џўџ u1џџ heated ambientў‘џцџўў'џ1џџыџ coefў?№ў1џџыџ coefў@ўўџрџџџэџьџў‘џцџўў'џCџBџ1џџрџ coefў?№ўў1џџрџ coefў@ўўCџ?№1џџрџ coefў@ўўPџ@ўџџџшџ logpetўўџџџшџ loglenўў!џџџ"џ Petiole length (mm)#џ Leaf length (mm)иџcџeџ nўYџџџ:џ?фЬЬЬЬЬЭ|џ|џўџ0џ7џ1џџџў usrўџ?№@ўў8џџ?рBџ?рўўџџў cxyўўўўAџ1џ0џ?№ў1џ0џ@ўџ r =Eџ?џџџў@ўў9џўuџvџ@wџ?ЭpЃз =q B9џўџ"џџдџџўўўџ"џџ?р?№"џўўџ"џ1џ"џGџ"џўўўџ#џџдџџўўўџ#џ1џ#џGџ#џўўўйџ?№кџаџ"џўлџџ"џўўйџ@кџаџ#џўлџџ#џўўмџўџыџџџџџўўўџрџџџџџўўў'џ1џџыџ coefў?№ў1џџыџ coefў@ўў'џCџBџ1џџрџ coefў?№ўў1џџрџ coefў@ўўCџ?№1џџрџ coefў@ўўPџ@ўџfџџ 3Plot (A) is for a dataset that compared the stretch 4 of an elastic band that was heated with that for an 3 elastic band that was held at ambient temperature, 1 while (B) is for a leaf dataset. In each plot we 2 show the regression line for y on x (solid line), 2 and the regression line for x on y (dotted line). 1 In (A) the lines are quite similar, while in (B) 0 where the correlation is smaller, the lines are  quite different.ўў/џfџ  ў џ џ џ ў[џўўў g5.2џ Dfunction(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) G titl <- paste("Diagnostic plots for the regression of Figure 5.1:", H "\nPanel (A) plots residuals against fitted values.", K "\nPanel (B) is a normal probability plot of residuals.", @ "\nIf residuals follow a normal distribution", > "the points \nshould fall close to a line.") 9 oldpar <- par(mfrow = c(1, 2), mgp = c(2.25, 0.5, 0), 6 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) G plot(hat, res, xlab = "Fitted values", ylab = "Residuals", pch = 1) 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) B 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) }ў§џџџџўџџџџў curveџ regў fignum@ЬЬЬЬЬЭ џ ў srcref  srcfileђ encoding native.enc timestampAвГD^€ class POSIXt POSIXctў filename /tmp/johnm.g5.2.R wd /Users/johnm/r/ch5ўўіџ srcfileўіџ srcrefў ''ђџѓџіџ srcrefў  >>ђџѓџіџ srcrefў   66ђџѓџіџ srcrefў ђџѓџіџ srcrefў //ђџѓџіџ srcrefў ђџѓџіџ srcrefў ђџѓџіџ srcrefў GGђџѓџіџ srcrefў 11ђџѓџіџ srcrefў ..ђџѓџіџ srcrefў ђџѓџіџ srcrefў ::ђџѓџіџ srcrefў ..ђџѓџіџ srcrefў ђџѓџіџ srcrefў ђџѓџіџ srcrefў ђџѓџіџ srcrefўђџѓџў џ џ џ џ ў psfile№џ№џџ@џ@™™™™™š џ џўўџџџ 2Diagnostic plots for the regression of Figure 5.1: 1 Panel (A) plots residuals against fitted values. 5 Panel (B) is a normal probability plot of residuals. * If residuals follow a normal distribution (the points should fall close to a line.ўўџџџgџџ?№@ўџџ@?рўџBџџџўџўџ?№?р?рўўўўџџџўўџnџџџџџў‘џџўўџ resœџnџўўџ hat fittednџўў!џћџњџ"џ Fitted values#џ Residuals$џ?№ўYџћџњџ$џ?№:џ?щ™™™™™š”џ@ўuџvџ@wџ?а A9џў'џ hPџ@ў qqnormњџ"џ Quantiles of Normal Distribution#џ Ordered data values$џ@0 cex.main?ђffffffўuџvџ@wџ?а B9џў/џџ  ў џ џ џ ў[џўў\џnџўў g5.3џ function(fignum=5.3, device=""){ C if(device!="")psfile(fignum=fignum, width=4.8, height=2.75, 8 device=device, pointsize=9) Q if(!exists("roller.lm"))roller.lm <- lm(depression ~ weight, data=roller) C qreference(residuals(roller.lm), nrep=8, nrows=2, xlab="")  if(device!="")dev.off()  }ў§№џ@333333 џ ўёџ $$$$ђџђєџ native.encѕџAвГDРіџ POSIXt POSIXctўїџ /tmp/johnm.g5.3.Rјџ /Users/johnm/r/ch5ўўіџ srcfileўіџ srcrefў  8 8ђџџіџ srcrefў  Q Qђџџіџ srcrefў  B Bђџџіџ srcrefў   ђџџіџ srcrefўђџџў џ џ џ џ ўљџ№џ№џџ@333333џ@ џ џ pointsize@"ўў џ‚џ exists roller.lmўўџ roller.lmџџџџў‘џџўўў qreference residualsџў`џ@  nrows@"џ ў џ џ џ ў[џўўў g5.4џ0 Afunction(y = ironslag$chemical, x = ironslag$magnetic, device="") { 1 if(device!="")hardcopy(width=3.75, height=3.25, 4 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) C mtext(side = 3, line = 0.25, txt, cex = 0.7, at=xleft, adj = 0)  }  titl <- paste(reg = E "Chemical vs Magnetic Test of Iron Content in Slag.", D "\nThe fitted curves used the Splus loess smooth.", H "\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", 5 "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]) C 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") D 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)) D mtext3(item="C", txt=leg[3], xleft=par()$usr[1]+1.25*par()$cxy[1])  yabs <- sqrt(abs(res)) E plot(hat, yabs, xlab = "Predicted chemical", ylab = "", type = "n") > mtext(side = 2, line = 2.5, expression(sqrt(abs(residual)))) D 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() }ў§џџ ironslag chemicalўџџ џ magneticў$ џ ў џ џ џ џ ў џџ@џ@  џ џџ@"ўўџџџgџџ@@ўџџ@333333@ffffff@ЬЬЬЬЬЭ@ЬЬЬЬЬЭўџџ@?шў:џ?ы333333 mex?ы333333ўўџџџўўџ mtext3­џ$ item A$ txt1џ leg?№ў$ xleft7џ1џџџў2џў?№ў8џ?є1џџџў4џў?№ўўўў џuџvџ@wџ?а џ9џўuџvџ@wџ?а џ:џ?цffffffкџџ9џўў function(item="A", txt=leg[1], ; xleft=par()$usr[1]+1.25*par()$cxy[1]){ / mtext(side = 3, line = 0.25, item, adj = 0) C mtext(side = 3, line = 0.25, txt, cex = 0.7, at=xleft, adj = 0)  }ўўџџџ reg 2Chemical vs Magnetic Test of Iron Content in Slag. / The fitted curves used the Splus loess smooth. 4 In (d) the downward slope suggests a lower variance  for larger fitted values.ўўџџџ Fitted line and loess curve Residuals from line, with smooth &Observed vs predicted, for fitted line Is variance constant about line?ўў!џџџ"џ Magnetic#џ Chemicaleџ nўџnџџџџџўўў'џ1џџnџ coefў?№ў1џџnџ coefў@ўўџџџnџўўNџџџў џџ7џ1џџџў usrў?№ў8џ?є1џџџў cxyў@ўўўў-џ panel.smoothџџ span?юffffffPџ@”џ?ј$џўўџњџџnџўў!џџњџ"џ Magnetic#џ Residualeџ nў џ џ B џ1џџ@ўџ7џ1џџџў usrў?№ў8џ?є1џџџў cxyў?№ўўўў-џџџњџџ?юffffffўўYџџњџ$џ?№:џ?ьЬЬЬЬЬЭ”џ?№ўџћџќџnџўў!џћџџ"џ Predicted chemical#џ Chemicaleџ nў-џџћџџџ?юffffffўў џ џ C џ1џџ@ўџ7џ1џџџў usrў?№ў8џ?є1џџџў cxyў?№ўўўўџ yabs sqrt absњџўўў!џћџџ"џ Predicted chemical#џ eџ nўuџvџ@wџ@ expressionџџ residualўўўў џ џ D џ1џџ@ўџ7џ1џџџў usrў?№ў8џ?є1џџџў cxyў?№ўўўў-џџћџџџ?юffffffўў/џџ  ў џ џ џ ў[џўў\џўў g5.5џ Afunction(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)) 4 leg <- c("A. Residuals from fitted loess curve.", 2 "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) 8 mtext(side = 3, line = 0.25, leg[1], adj = 0, cex=0.7) 7 plot(yhat, yabs, xlab = "Predicted chemical", ylab = 3 expression(sqrt(abs(residual))), type = "n") 1 panel.smooth(yhat, yabs, span = 1.1, cex = 1.1)  8 mtext(side = 3, line = 0.25, leg[2], adj = 0, cex=0.7)  if(device!="")dev.off() invisible() }ў§џџ џ chemicalўџџ џ magneticў$ џ ў џ џ џ џ ў џџ@џ@ џ џўўџџџgџџ?№@ўџџ@ffffff@ffffff@ЬЬЬЬЬЭ?љ™™™™™šўџџ@?рўўўџџџўўџџџ %A. Residuals from fitted loess curve. $B. Is variance about curve constant?ўўџnџ loessџџџўўўџ resvalџnџўўџџџnџўўџџџџџўўў!џџџ"џ Magnetic#џ Residual$џ?№ўYџџџ:џ?щ™™™™™šeџ nў'џ§џPџ@ўuџvџ@wџ?а1џџ?№ў9џ:џ?цffffffў!џџџ"џ Predicted chemical#џџџџџўўўeџ nўџџџџ?ё™™™™™š:џ?ё™™™™™šўuџvџ@wџ?а1џџ@ў9џ:џ?цffffffў џ џ џ ў[џўў\џўў g5.6џ? Gfunction(y = softbacks$weight, x = softbacks$volume, curve = c("reg"), $ show.fits = T, device="") { 1 if(device!="")hardcopy(width=2.25, height=2.25, % device=device)  titl <- switch(curve[1], I reg = paste("Weight versus volume for softcover books,", * "\nwith fitted line."),  lo = paste( @ "Weight versus volume for softcover books,", H "\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] 5 btxt <- c(paste("a =", format(round(z[1, 1], 3)), 9 " SE =", format(round(z[1, 2], 3))), 5 paste("b =", format(round(z[2, 1], 3)), 9 " 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) }ў§џџ softbacksџўџџџ volumeўяџџ regў show.fitsЈџ џ ў џ џ џ џ ў џџ@џ@ џ џўўџџвџ1џяџ?№ўџџ )Weight versus volume for softcover books,  with fitted line.ў loџ )Weight versus volume for softcover books, 0 with fitted line and S-PLUS loess smooth curve.ўўўџџџџџ@ffffff@ffffff?ё™™™™™š?ё™™™™™šўџџ@?рўўўџџџўўџnџџџџџў‘џџўў/џ  Coefficients ў+џ,џ@ў-џџ)џnџў*џўў/џ  ў-џ.џnџўўџџџnџўўџ>џ?џџџўўџ%џбџџўўџ1џ%џ@ў7џ1џ%џ@ў8џ diff%џў?ДzсGЎ{ўўў!џџџ"џ Volume (cc)%џ%џ#џ Weight (g)$џ?№&џбџџџџўўў џ match regяџ nomatchў џ'џ1џџnџ*џў?№ў1џџnџ*џў@ўPџ?№ўџHџџ)џnџў*џўў џџ џYџџџ$џ?№:џ?шўџњџœџnџўўiџjџVџ?№Rџњџўў џџ resi1џњџjџўўџ izzy as.numericGџ#џўўўџ xi1џџjџўўџ yhati1џџjџўўџ yi1џџjџўўNџKџ&џ@ўџ'џ(џў|џBџ@$џўPџBџ@$џўўџ eps8џ1џџџў4џў?№ў?Щ™™™™™šўў џZџjџ@ў џџ adji?№ўџ)џBџ)џўўўџ*џўўAџ7џ&џ)џў7џ'џCџ#џ@ўўџEџ#џ?№ўў9џ*џ:џ?фЬЬЬЬЬЭўўўўўџ@џ1џџџў2џўџ@@ўўўџ3џ1џџџў4џў?№ўўџ5џ1џџџў4џў@ўўџ btxtџџ a =DџEџ1џHџ?№?№ў@ўў  SE =DџEџ1џHџ?№@ў@ўўўџ b =DџEџ1џHџ@?№ў@ўў  SE =DџEџ1џHџ@@ў@ўўўўў6џ1џ@џ?№ў1џ@џ@ў6џ+џЊџ?№ yjust:џ?щ™™™™™š=џ nўўў/џ  џ  ў џ џ џ ў[џўў\џnџўў g5.7џ function(device=""){ @ if(device!="")hardcopy(width=5.5, height=1.4, device=device) 7 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) 2 softbacks.lm<-lm(weight~volume,data=softbacks)  on.exit(par(oldpar)) * plot(softbacks.lm, which=1:4, caption= % c("A: Residuals vs Fitted", 0 "B: Normal Q-Q", "C: Scale-Location", ! "D: Cook's distance")) 5 figtxt <- "Diagnostic plots for previous figure."  cat(figtxt, "\n")  if(device!="")dev.off()  invisible(softbacks.lm)  }ў§$ џ ў џ џ џ џ ў џџ@џ?іffffff џ џўўџџџgџџ?№@ўџџ@ffffff@ffffff@ЬЬЬЬЬЭ?у333333ўџџ@?рўhџџ?№ўъџ sўўГџЪџў‘џџўџ softbacks.lmџџџџў‘џџўўџџџўў!џ.џ whichVџ?№@ў captionџ A: Residuals vs Fitted B: Normal Q-Q C: Scale-Location D: Cook's distanceўўџfџ %Diagnostic plots for previous figure.ў/џfџ  ў џ џ џ ў[џўў\џ.џўў g5.8џ >function(y = roller$depression, x = roller$weight, pred = T,  device="") { 8 if(device!="")psfile(width=2, height=2, device=device) B titl <- paste("Lawn Depression, for Various Weights of Roller,", ? "\nwith fitted line and showing 95% pointwise", E "\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)) F plot(x, y, xlab = "Weight of Roller (tonnes)", xlim=xlim, ylim=ylim, 2 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") 7 ci<-data.frame(lower=hat[,"lwr"],upper=hat[,"upr"])  here <- ci$lower >= y12[1] D lines(xy$weight[here], ci$lower[here], lty = 2,lwd=2,col="grey")  here <- ci$upper < y12[2] D lines(xy$weight[here], ci$upper[here], lty = 2,lwd=2,col="grey")  }  cat(titl, "\n")  if(device!="")dev.off()  invisible(u) }ў§џџџџўџџџџў›џЈџ џ ў џ џ џ џ ўљџџ@џ@ џ џўўџџџ /Lawn Depression, for Various Weights of Roller, + with fitted line and showing 95% pointwise 1 confidence bounds for points on the fitted line.ўўџџџџџ@ffffff@ffffff?ё™™™™™š?ё™™™™™šўџџ@?рўўўџџџўўџ%џџ џџўўўџ&џџ џџўўў!џџџ"џ Weight of Roller (tonnes)%џ%џ&џ&џ#џ Depression in Lawn (mm)$џ@0ўџnџџџџџў‘џџўў'џ1џџnџ*џў?№ў1џџnџ*џў@ўPџ?№ўџHџџ)џnџў*џўўџ y121џџџў2џўVџ@@ўўў џ›џ џ assign xyžџџџџџџў@4ўўўџћџџnџ newdataџ interval confidenceўўџ cižџ lower1џћџћ lwrў upper1џћџћ uprўўўџFџWџџ6џ7џў1џ2џ?№ўўўNџ1џџџџўFџў1џџ6џ7џўFџўPџ@”џ@|џ greyўџFџOџџ6џ8џў1џ2џ@ўўўNџ1џџџџўFџў1џџ6џ8џўFџўPџ@”џ@|џ greyўўў/џџ  ў џ џ џ ў[џўў\џnџўў g5.9џ' function(color = F, device="") { < if(device!="")hardcopy(width=3.4, height=2, device=device) E figtxt<-paste("Two rubber band experiments, with different ranges", E "\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) 0 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) D lines(lowess(x[ord], ci$upper[ord]), lty=2, lwd=2, col="grey") D lines(lowess(x[ord], ci$lower[ord]), lty=2, lwd=2, col="grey")  }  xy<-rbind(elastic2,elastic1) C nam <- c("Range of stretch 30-60 mm","Range of stretch 42-54 mm") 7 trial<-rep(nam, c(dim(elastic2)[1],dim(elastic1)[1]))  xlim<-range(elastic2$stretch) ylim<-range(elastic2$distance)  xy<-split(xy,trial) I 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 1 panelplot(xy,panel=panelci,totrows=1,totcols=2, : par.strip.text=list(cex=.8), oma=c(4,4,2.5,2))  cat(figtxt, "\n") 5 mtext(side = 2, line = 2.75, "Distance moved (cm)") 1 mtext(side=1,line=3.5,"Amount of stretch (mm)")  if(device!="")dev.off() }ў§ colorcџ џ ў џ џ џ џ ў џџ@ 333333џ@ џ џўўџfџџ 2Two rubber band experiments, with different ranges 1 of x-values. The dashed curves are pointwise 95% ' confidence bounds for the fitted line.ўўџ panelci­џ‘џћ ...ћў џџџџ list<џўџўўџ ncolsџ=џ<џў>џўў џZџ>џ?№ўйџ@ўў џZџ>џ?№ўйџ?№ўйџ@ўўџџџ‘џ stretchўўџџџ‘џ distanceўўџnџџџџџўўўџ upredџnџ5џ confidenceўўџ6џžџ fit1џAџћ fitў7џ1џAџћ lwrў8џ1џAџћ uprўўўџ ord orderџўўNџ1џџCџў1џџ6џBџўCџўPџ?№”џ@ўNџ lowess1џџCџў1џџ6џ8џўCџўўPџ@”џ@|џ greyўNџEџ1џџCџў1џџ6џ7џўCџўўPџ@”џ@|џ greyўў 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) 0 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) D lines(lowess(x[ord], ci$upper[ord]), lty=2, lwd=2, col="grey") D lines(lowess(x[ord], ci$lower[ord]), lty=2, lwd=2, col="grey")  }ўўџџJџ elastic2 elastic1ўўџ namџ Range of stretch 30-60 mm Range of stretch 42-54 mmўўџ trialKџHџџ1џˆџFџў?№ў1џˆџGџў?№ўўўўџ%џбџџFџ?џўўўџ&џбџџFџ@џўўўџџ splitџIџўўџџ lapplyVџ?№Rџџўў­џjџћў џџ as.listrџџjџўў=џ%џ%џ&џ&џўўў /function(i){c(as.list(xy[[i]]), list(xlim=xlim, ? ylim=ylim))}ўўўџ namesџўHџў panelplotџ panel;џ totrows?№ totcols@ par.strip.text=џ:џ?щ™™™™™šўhџџ@@@@ўў/џfџ  ўuџvџ@wџ@ Distance moved (cm)ўuџvџ?№wџ@  Amount of stretch (mm)ў џ џ џ ў[џўўў gdumpџ 7function(fnam=NULL, prefix="~/r-book/ed2/figures/figs",  splitchar="/ch"){  if(is.null(fnam)){  path <- getwd() 7 pathtag <- strsplit(path, "/ch", fixed=TRUE)[[1]] 5 fnam <- paste(prefix, pathtag[length(pathtag)], ! ".R", sep="")  } - else fnam <- paste(prefix, fnam, sep="/") H objnames <- c(objects(pattern="^g", envir=sys.frame(0)), "hardcopy") & cat("\nDump to file:", fnam, "\n")  print(objnames)  dump(objnames, fnam)  }ў§ fnamў prefix ~/r-book/ed2/figures/figs splitchar /chў џ џƒџTџў џџ path getwdўўџ pathtagrџ strsplitWџ /ch fixed ў?№ўўџTџџUџ1џYџRџYџўў .Rџ ўўўџTџџUџTџџ /ўўўџ objnamesџ objects pattern ^g envir sys.frameўў hardcopyўў/џ  Dump to file:Tџ  ў-џ\џў dump\џTџўў gfileџ& 5function(width=3.75, height=3.75, color=F, trellis=F, H 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" D ## 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) G prefix1 <- paste(if(nchar(dotsplit[1])==1)"0" else "", dotsplit[1],  sep="") G prefix2 <- paste(if(nchar(dotsplit[2])==1)"0" else "", dotsplit[2],  sep="") 6 if(device=="")stop("No device has been specified") 3 suffix <- switch(device, ps=".eps", pdf=".pdf") = fnam <- paste("~/r-book/second/Art/",prefix1,"-",prefix2, ! suffix, sep="")  print(fnam)  dev.out <- device[1] 6 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) L 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, C width=width, height=height, pointsize=pointsize[1])  }  }ў§$џ@$џ@:џcџ trelliscџ$ џџ  pdf psўWџ џџ@ @ўЉџcџў џ џ‚џcџўџџ1џџ?№ўўўџ funtxt sys.call?№ўўџTџ1џrџZџ as.characterdџў ([џЈџў?№ў?№ўўџ dotsplitrџZџTџ \.ў?№ўўџ1џgџ?№ў substring1џgџ?№ў@ўўџ prefix1џ џZџ nchar1џgџ?№ўў?№ў 0 ў1џgџ?№ўџ ўўџ prefix2џ џZџjџ1џgџ@ўў?№ў 0 ў1џgџ@ўџ ўў џZџ џ ў stop No device has been specifiedўўџ suffixвџ џ ps .eps pdf .pdfўўџTџџ ~/r-book/second/Art/iџ -kџmџџ ўў-џTџўџ dev.out1џ џ?№ўўџ dev.funвџpџoџoџnџ postscriptўў џcџ џЩџ latticeў trellis.device fileTџ џqџ:џ:џџџџџЉџЉџў trellis.par.set fontsize=џAџ1џџ?№ўYџ1џџ@ўўўў џ џpџ ў џ-џџџџўўqџuџTџ paper special enc MacRomanЉџЉџџџџџџ1џџ?№ўўўўўў gfocus.demoџ$ function(device=""){  library(lattice)  library(grid) Q if(device!="")hardcopy(device=device, width=4, height=2.25, trellis=TRUE, ' color=TRUE) H trellis.par.set(layout.heights=list(key.top=0.25, axis.top=0.5)) H hp1.plt1 <- xyplot(o2 ~ wattsPerKg, groups=id, data=humanpower1, 9 panel=function(x,y,subscripts,groups,...){ ' u <- lm(y~groups*x); # hat <- fitted(u) 9 panel.superpose(x,y,subscripts,groups) E panel.superpose(x,hat,subscripts,groups, type="l")  }, F## key=simpleKey(text=rep("",5), lines=TRUE, columns=5), ) xlab="Watts per kilogram", K ylab=expression("Oxygen intake ("*ml.min^{-1}*.kg^{-1}*")"), N legend=list(top=list(fun=textGrob, args=list(label="A", x=0)))) / print(hp1.plt1, position=c(0,0,.535,1)) J u <- lme(o2 ~ wattsPerKg, random=~wattsPerKg|id, data=humanpower1) H hp1.plt2 <- xyplot(o2 ~ wattsPerKg, groups=id, data=humanpower1, 9 panel=function(x,y,subscripts,groups,...){ ' u <- lm(y~groups*x); # hat <- fitted(u) E panel.superpose(x,hat,subscripts,groups, type="l")  }, 2 xlab="Watts per kilogram", ylab="", N 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() H panel.superpose(x=arglist$x,y=hat,subscripts=arglist$subscripts, A groups=arglist$groups, , type="l", lty=2) ! trellis.unfocus()  if(device!="")dev.off()  }ў§ џ ў џЩџsџўЩџ gridў џ џ џ ў џ џ џџ@џ@cџ :џ ўўvџ layout.heights=џ key.top?а axis.top?рўўџ hp1.plt1 xyplotџ o2 wattsPerKgў groups id‘џ humanpower1Oџ­џџћџћ subscriptsћƒџћ<џћў џџnџџџџ8џƒџџўўўўџћџќџnџўў panel.superposeџџ†џƒџў‡џџћџ†џƒџeџ lўў $function(x,y,subscripts,groups,...){ ' u <- lm(y~groups*x); # hat <- fitted(u) 9 panel.superpose(x,y,subscripts,groups) E panel.superpose(x,hat,subscripts,groups, type="l")  }ў"џ Watts per kilogram#џџ8џ8џ8џ Oxygen intake (Єџ ml.min џBџ?№ўўўўЄџ .kg џBџ?№ўўўў )ўў6џ=џ top=џ fun textGrob args=џ label Aџўўўўў-џџ positionџ?сИQы…?№ўўџnџ lmeџџ‚џў randomџ |‚џ„џўў‘џ…џўўџ hp1.plt2€џџџ‚џўƒџ„џ‘џ…џOџ­џџћџћ†џћƒџћ<џћў џџnџџџџ8џƒџџўўўўџћџќџnџўў‡џџћџ†џƒџeџ lўў $function(x,y,subscripts,groups,...){ ' u <- lm(y~groups*x); # hat <- fitted(u) E panel.superpose(x,hat,subscripts,groups, type="l")  }ў"џ Watts per kilogram#џ 6џ=џŠџ=џ‹џŒџџ=џŽџ Bџўўўўўџћџќџnџўў-џ“џџџ?нТ\(ѕУ?№?№ў newpage ў trellis.focus panel row?№ column?№ўџ arglist trellis.panelArgsўў‡џџџ˜џџўџћџ†џџ˜џ†џўƒџџ˜џƒџўћeџ lPџ@ў trellis.unfocusў џ џ џ ў[џўўў gsaveџ 7function(fnam=NULL, prefix="~/r-book/ed2/figures/figs", K splitchar="/ch", xtras=c("renum.fun","renum.files","hardcopy")){  if(is.null(fnam)){  path <- getwd() 7 pathtag <- strsplit(path, "/ch", fixed=TRUE)[[1]] 5 fnam <- paste(prefix, pathtag[length(pathtag)], % ".RData", sep="")  } - else fnam <- paste(prefix, fnam, sep="/") C objnames <- c(objects(pattern="^g", envir=sys.frame(0)), xtras) & cat("\nDump to file:", fnam, "\n")  print(objnames) " save(list=objnames, file=fnam)  }ў§TџўUџ ~/r-book/ed2/figures/figsVџ /ch xtrasџ renum.fun renum.files hardcopyўў џ џƒџTџў џџWџXџўўџYџrџZџWџ /ch[џ ў?№ўўџTџџUџ1џYџRџYџўў .RDataџ ўўўџTџџUџTџџ /ўўўџ\џџ]џ^џ ^g_џ`џўўœџўў/џ  Dump to file:Tџ  ў-џ\џў save=џ\џuџTџўў gtestџ function(){ / trellis.device(postscript, file="test.eps") 5 trellis.par.set(layout.heights=list(key.top=0.5)) 6 zz.d <- dotplot(variety ~ yield, data = barley, B legend=list(top=list(fun=grid.text, I args=list(label="ABC", x=0)))) 4 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()  }ў§ў џtџrџuџ test.epsўvџ|џ=џ}џ?рўўџ zz.d dotplotџ variety yieldў‘џ barley6џ=џŠџ=џ‹џ grid.textџ=џŽџ ABCџўўўўў pushViewport viewport layout grid.layout@?№ўўўЅџІџ layout.pos.row?№ўў-џŸџ”џ ў upViewportўЅџІџЉџ@ўў-џŸџ”џ ў popViewport@ў[џўў renum.funџ Dfunction(from.prefix="g", to.prefix="g",from=20:7, to=21:8, doit=F){  path <- getwd() 5 pathtag <- strsplit(path, "/ch", fixed=TRUE)[[1]] & endbit <- pathtag[length(pathtag)] 5 from.prefix <- paste(from.prefix, endbit, sep="") 1 to.prefix <- paste(to.prefix, endbit, sep="")  for(i in 1:length(to)) M {txt<-paste(to.prefix,".",to[i]," <- ", from.prefix,".",from[i],sep="") 7 if(doit)eval(parse(text=txt),envir=sys.frame(0))  print(txt)  }  }ў§ from.prefix g to.prefix g fromVџ@4@ў toVџ@5@ ў doitcџў џџWџXџўўџYџrџZџWџ /ch[џ ў?№ўўџ endbit1џYџRџYџўўўџ­џџ­џВџџ ўўџЎџџЎџВџџ ўўiџjџVџ?№RџАџўў џџ џџЎџ .1џАџjџў  <- ­џ .1џЏџjџўџ ўў џБџ eval parseAџ џў_џ`џўўў-џ џўўўў renum.filesџ Hfunction(from.prefix="~/r-book/ed2/Art/", to.prefix="~/r-book/ed2/Art/", ' from=20:7, to=21:8, doit=F){  path <- getwd() 5 pathtag <- strsplit(path, "/ch", fixed=TRUE)[[1]] & endbit <- pathtag[length(pathtag)] 2 if(nchar(endbit)==2)chap <- paste(endbit) else ' chap <- paste("0",endbit,sep="") 8 from.prefix <- paste(from.prefix, chap, "-", sep="") 4 to.prefix <- paste(to.prefix, chap, "-", sep="")  for(i in 1:length(from)){ 8 if (from[i]<=9) ltext <- paste("0",from[i],sep="") " else ltext <- paste(from[i]) 4 if (to[i]<=9) rtext <- paste("0",to[i],sep="") , else rtext <- paste(to[i]) 8 txt<-paste("mv ", from.prefix, ltext, ".eps", " ", 2 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)  }  }ў§­џ ~/r-book/ed2/Art/Ўџ ~/r-book/ed2/Art/ЏџVџ@4@ўАџVџ@5@ ўБџcџў џџWџXџўўџYџrџZџWџ /ch[џ ў?№ўўџВџ1џYџRџYџўўў џZџjџВџў@ўџ chapџВџўўџЖџџ 0Вџџ ўўўџ­џџ­џЖџ -џ ўўџЎџџЎџЖџ -џ ўўiџjџVџ?№RџЏџўў џ џзџ1џЏџjџў@"ўџ ltextџ 01џЏџjџўџ ўўџЗџџ1џЏџjџўўўў џзџ1џАџjџў@"ўџ rtextџ 01џАџjџўџ ўўџИџџ1џАџjџўўўўџ џџ mv ­џЗџ .eps  ЎџИџ .epsџ ўўџ backupџ cp ­џЗџ .eps   archiveџ ўў џБџ systemЙџўў џБџКџ џўў-џЙџў-џ џўўўў џџ4 =function(width=3.75, height=3.75, color=FALSE, trellis=FALSE, < device=c("","pdf","ps"), path=getwd(), file=NULL, D format=c("nn-nn", "name"), split="\\.", pointsize=c(8,4),  fonts=NULL,  horiz=FALSE, ...){ ) if(!trellis)pointsize <- pointsize[1]  funtxt <- sys.call(1) B nam <- strsplit(as.character(funtxt), "(", fixed=TRUE)[[1]][1] 3 suffix <- switch(device, ps=".eps", pdf=".pdf") N if(is.character(path) & nchar(path)>1 & substring(path, nchar(path))!="/") & path <- paste(path, "/", sep="") - if(is.null(file)) if(format[1]=="nn-nn"){ C if(!is.null(split))dotsplit <- strsplit(nam, split)[[1]] else  dotsplit <- nam 8 if(length(dotsplit)==1)dotsplit <- c("", dotsplit) E nn2 <- paste(if(nchar(dotsplit[2])==1)"0" else "", dotsplit[2],  sep="")  if(nchar(dotsplit[1])>0){ O numstart <- which(unlist(strsplit(dotsplit[1], "")) %in% paste(0:9))[1] / nn1 <- substring(dotsplit[1], numstart) G nn1 <- paste(if(nchar(nn1) == 1) "0" else "", nn1, "-", sep="")  } else nn1 <- "" % file <- paste(nn1, nn2, sep="")  } else file <- nam L 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] 6 dev.fun <- switch(dev.out, pdf=pdf, ps=postscript)  if(trellis){  library(lattice)  if(device=="ps") 1 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, D color = color, width=width, height=height, ...) R trellis.par.set(list(fontsize=list(text=pointsize[1], points=pointsize[2]))) } else  if (dev.out!=""){  print(c(width, height))  if(device=="ps") F dev.fun(file=file, paper="special", horiz=horiz, fonts=fonts, M width=width, height=height, pointsize=pointsize[1], ...) else 6 dev.fun(file=file, paper="special", fonts=fonts, L width=width, height=height, pointsize=pointsize[1], ...)  } D if(trellis)trellis.par.set(list(fontsize=list(text=pointsize[1], < points=pointsize[2])))  }ў§џ@џ@:џ cџ  џџ  pdf psўWџXџўuџўDџџ nn-nn nameўJџ \.џџ@ @ў fontsўЉџ <џћўёџ  )ђџђєџ native.encѕџAбрТŽіџ POSIXt POSIXctўїџ /tmp/johnm.hardcopy.Rјџ /Users/johnm/r/ch5ўўіџ srcfileўіџ srcrefў ђџМџіџ srcrefў   BђџМџіџ srcrefў   3ђџМџіџ srcrefў   &ђџМџіџ srcrefў  ђџМџіџ srcrefў ђџМџіџ srcrefў -ђџМџіџ srcrefў :ђџМџіџ srcrefў ђџМџіџ srcrefў 6ђџМџіџ srcrefў  2ђџМџіџ srcrefў 34<ђџМџіџ srcrefўђџМџў џ џ‚џcџўџџ1џџ?№ўўўџdџeџ?№ўўџHџ1џrџZџfџdџў ([џ ў?№ў?№ўўџmџвџ џnџ .epsoџ .pdfўў џжџжџ is.characterWџўGџjџWџў?№ўў џhџWџjџWџўў /ўўџWџџWџ /џ ўўў џƒџuџў џZџ1џDџ?№ў nn-nnўёџ ђџМџіџ srcrefў 8ђџМџіџ srcrefў ђџМџіџ srcrefў ђџМџіџ srcrefў %ђџМџіџ srcrefўђџМџў џ џ‚џƒџJџўўџgџrџZџHџJџў?№ўўџgџHџўў џZџRџgџў?№ўџgџџ gџўўўџ nn2џ џZџjџ1џgџ@ўў?№ў 0 ў1џgџ@ўџ ўў џGџjџ1џgџ?№ўўўёџ  OђџМџіџ srcrefў  /ђџМџіџ srcrefў  GђџМџіџ srcrefўђџМџў џџ numstart1џ/џ %in% unlistZџ1џgџ?№ў ўўџVџ@"ўўўў?№ўўџ nn1hџ1џgџ?№ўПџўўџТџџ џZџjџТџў?№ў 0 ўТџ -џ ўўўџТџ ўўџuџџТџОџџ ўўўџuџHџўўў џжџGџjџuџў@ўZџhџuџ7џBџjџuџўjџmџўў?№ўўmџўўџmџ ўўџuџџWџuџmџџ ўў-џџ Output will be directed to file:uџўўџpџ1џ џ?№ўўџqџвџpџoџoџnџrџўў џcџёџ !!ђџМџіџ srcrefў "'DђџМџіџ srcrefў ))RђџМџіџ srcrefўђџМџў џЩџsџў џZџ џ psўtџuџuџ џqџ:џ:џЉџЉџЛџЛџџџџџ<џўtџuџuџ џqџЛџЛџ:џ:џџџџџ<џўўvџ=џwџ=џAџ1џџ?№ўYџ1џџ@ўўўўў џ џpџ ўёџ ,,ђџМџіџ srcrefў -1FђџМџіџ srcrefўђџМџў џ-џџџџўў џZџ џ psўqџuџuџxџ specialЉџЉџЛџЛџџџџџџ1џџ?№ў<џўqџuџuџxџ specialЛџЛџџџџџџ1џџ?№ў<џўўўўў џcџvџ=џwџ=џAџ1џџ?№ўYџ1џџ@ўўўўўўў