RDX2 X   g3.1 source function(fignum=3.1, device="") { + if(device!="")hardcopy(width=2, height=2, ' device=device) H oldpar <- par(mar = c(4.1,4.1,1.1,1.1), mgp = c(2.5, 0.5, 0), pty="s")  on.exit(par(oldpar)) x <- 0:5  y <- 0.5 * 9.8 * x^2  par(col="gray25", pty="s") 7 curve(0.5*9.8*x^2, from=0, to=5, xaxs="i", yaxs="i", 8 xlab = "Time (sec)", ylab = "Distance fallen (m)")  chh <- par()$cxy[2] ' text(rep(0.2,2), 107.5-c(0,1.75*chh), N expression("Distance " == " "*frac(1,2)*phantom(i)*g*t^2, F "where "*g == 9.8*phantom(i)*m*"."*"sec"^{-2}), adj=0, cex=0.8) < figtxt <- paste("Distance of drop of stone in free fall, " , ,"versus time.", sep = "")  cat("\n",figtxt,"\n")  if(device!="")dev.off() } fignum@ device  { if !=  hardcopy width@ height@ <- oldpar par mar c@ffffff@ffffff?񙙙?񙙙 mgp@? pty s on.exit    x :@  y *?@# ^@  col gray25 s curve?@#@ from to@ xaxs i yaxs i xlab Time (sec) ylab Distance fallen (m)  chh [ $  cxy@ text rep?ə@ -@Z?  expression == Distance    frac?@ phantom i g t@( where ,@#*+ m . sec&@ adj cex?陙  figtxt paste (Distance of drop of stone in free fall,  versus time. sep  cat  1    dev.off g3.2 function(fignum=3.2, device="") { : if(device!="")hardcopy(width=2, height=2, device=device) = oldpar <- par(mar = c(4.1,4.1,1.1,1.1), xaxs="i", yaxs="i", . mgp = c(2.5, 0.5, 0), pty="s")  on.exit(par(oldpar))  library(DAAG); data(roller) ( plot(depression~weight, data = roller, " xlim=c(0,13), ylim=c(0,32), A xlab = "Weight of roller (t)", ylab = "Depression (mm)",  pch = 16)  abline(0, 2.5) H figtxt <- paste("Depression (mm), versus weight of roller.", sep = "")  cat("\n",figtxt,"\n")  if(device!="")dev.off() }@    @ @   @ffffff@ffffff?񙙙?񙙙 i i@? s   library DAAG data roller plot ~ depression weight9: xlim@* ylim@@ Weight of roller (t) Depression (mm) pch@0 abline@ 12 )Depression (mm), versus weight of roller.3 4  1   5 g3.3> 2function(y = roller$depression, x = roller$weight,  device="") { @ if(device!="")hardcopy(width=4.25, height=2.25, device=device) oldpar <- par(mfrow = c(1, 2), ) mar = c(4.1,4.6,1.6,1.1), > mgp = c( 2.5, 0.75, 0),oma=c(0,0,0,0),pty="s")  on.exit(par(oldpar)) # pretext <- c(reg = "A", lo = "B") for(curve in c("reg", "lo")) { 2 plot(x, y, xlab = "Roller weight (t)", ylab = , "Depression in lawn (mm)", pch = 4) 9 mtext(side = 3, line = 0.25, pretext[curve], adj = 0) ! topleft <- par()$usr[c(1, 4)]  chw <- strwidth("O")  chh <- strheight("O") @ points(topleft[1]+rep(0.75,2)*chw,topleft[2]-c(0.5,1.3)*chh,  pch=c(4,1)) ? text(topleft[1]+rep(1.25,2)*chw, topleft[2]-c(0.5,1.3)*chh, : c("Data values", "Fitted values"),adj=0,cex=0.65) H text(topleft[1]+1.25*chw, topleft[2]-2*chh,"(smooth)",adj=0,cex=.65)  if(curve[1] == "reg") {  u <- lm(y ~ -1 + x)  abline(0, u$coef[1])  yhat <- predict(u)  } else {  lawn.lm<-lm(y~x+I(x^2))  yhat<-predict(lawn.lm)  xnew<-pretty(x,20)  b<-lawn.lm$coef & ynew<-b[1]+b[2]*xnew+b[3]*xnew^2  lines(xnew,ynew)  }  here <- y < yhat F yyhat <- as.vector(rbind(y[here], yhat[here], rep(NA, sum(here)))) @ xx <- as.vector(rbind(x[here], x[here], rep(NA, sum(here)))) ) lines(xx, yyhat, lty = 2, col="gray")  here <- y > yhat F yyhat <- as.vector(rbind(y[here], yhat[here], rep(NA, sum(here)))) @ xx <- as.vector(rbind(x[here], x[here], rep(NA, sum(here)))) ) lines(xx, yyhat, lty = 1, col="gray")  n <- length(y) 4 ns <- min((1:n)[y - yhat >= 0.75*max(y - yhat)]) $ ypos <- 0.5 * (y[ns] + yhat[ns])  chw <- par()$cxy[1] P text(x[ns] - 0.25*chw, ypos, "+ve residual", adj = 1,cex=0.65, col="gray30")  points(x, yhat, pch = 1) - ns <- (1:n)[y - yhat == min(y - yhat)][1] $ ypos <- 0.5 * (y[ns] + yhat[ns]) N text(x[ns] + 0.4*chw, ypos, "-ve residual", adj = 0,cex=0.65,col="gray30")  } : titl <- "Lawn depression (mm) versus roller weight (t)."  titl <- paste(titl, D "\nIn (a) a line has been fitted, while in (b) the", C "\nloess method was used to fit a smoothed curve.", G "\nResiduals (the `rough') appear as vertical lines. ", H "\nPositive residuals are black lines, while negative ", * "\nresiduals are dotted.")  if(device!="")dev.off() invisible() }":=":>   @ @    mfrow?@@ffffff@ffffff??񙙙@? oma s    pretext reg A lo B for reg lo; Roller weight (t) Depression in lawn (mm)A@ mtext side@ line?!F/  topleft!"  usr?@  chw strwidth O   strheight O points +!M?%?@O&!M@?? A@?$S!M?%?@O&!M@??  Data values Fitted values/0?$S!M??O&!M@@  (smooth)/0?(!? reg  u lm<S&?B!"T coef?  yhat predictT  lawn.lmU<S I@ WXY  xnew pretty@4  b"YV  ynewSS!]?!]@[!]@[@ lines[^  here <W  yyhat as.vector rbind!`!W`%  sum`  xxcd!`!`% e`_fb lty@ gray ` >W bcd!`!W`% e` fcd!`!`% e`_fbg? gray  n length  ns min! (?i >=&W? max&W  ypos?mS!k!Wk O!" #?$&!k?Op +ve residual/?0? gray30RWA? k!!m?i(&Wl&W? p?mS!k!Wk$S!k?ٙOp -ve residual/0? gray30  titl .Lawn depression (mm) versus roller weight (t). q2q 0 In (a) a line has been fitted, while in (b) the / loess method was used to fit a smoothed curve. 3 Residuals (the `rough') appear as vertical lines.  4 Positive residuals are black lines, while negative   residuals are dotted. 5 invisible g3.4 function(device="") { : if(device!="")hardcopy(width=2, height=2, device=device) J oldpar <- par(mar = c(4.1,4.1,1.1,1.1), mgp=c(2.5,0.5,0), cex.axis=0.85)  on.exit(par(oldpar))  z <- seq(-3,3,length=101) G plot(z, dnorm(z), type="l", ylab="normal density", yaxs="i", bty="L", . xlab="Distance, in SDs, from the mean") J polygon(c(z[z <= 1.0],1.0),c(dnorm(z[z <= 1.0]), dnorm(-3)), col="grey")  chh <- par()$cxy[2] 4 arrows(-1.8, 0.32, -0.25, 0.2, length=0.07, xpd=T)  cump <- round(pnorm(1), 3) K text(-1.8, 0.32+0.75*chh, paste("pnorm(1)\n", "=", cump), xpd=T, cex=0.8)  if(device!="")dev.off() invisible() }$   @ @   @ffffff@ffffff?񙙙?񙙙@? cex.axis?333333    z seq&@@j@Y@;u dnormu type l normal density i bty L Distance, in SDs, from the mean polygon!u <=u??w!u{u?w&@ grey  !"  cxy@ arrows&??zG{&??əj?Q xpd T  cump round pnorm?@$&?S?zG{? 2 pnorm(1)  =}~0?陙 5r g3.5D @function(mu = 10, sigma = 1, n = 1, m = 50, four = 4, nrep = 5, B seed=21, fignum=3.5, pscript=F, totrows=1, device="") { F if(device!="")psfile(fignum, width=5.2, height=1.4, device=device, % pointsize=8) $ if(!is.null(seed))set.seed(seed) D titl <- paste("\nEach panel shows a simulated distribution of ", B m, " values ", "from a normal \ndistribution ", 0 "with mean = 10 and sd = 1.", M " The underlying theoretical \nnormal curve is overlaid", . " on the lower left panel.") B oldpar <- par(mgp = c(2.25, 0.5, 0),mar=par()$mar-c(1,0,1.5,0)  )  on.exit(par(oldpar))  if(is.null(totrows)) $ totrows <- floor(sqrt(nrep)) $ totcols <- ceiling(nrep/totrows) 7 z <- range(pretty(mu + (c(-3.4, 3.4) * sigma), 50)) ? xy <- data.frame(x=rep(0,nrep),y=rep(0,nrep),n=rep(n,nrep), 8 mm=rep(m,nrep),four=rep(four,nrep)) 0 fac <- factor( paste("Simulation", 1:nrep), 6 lev <- paste("Simulation", 1:nrep) ) xlim<-z " ## ylim<-c(0,dnorm(0)*sqrt(n))  ylim <- c(0,1)  xy <- split(xy,fac) K xy<-lapply(1:length(xy),function(i){c(as.list(xy[[i]]), list(xlim=xlim, N ylim=ylim))}) = panel.mean <- function(data, mu = 10, sigma = 1, n2 = 1, 7 mm = 100, nrows, ncols, ...)  { 1 vline <- function(x, y, lty = 1, col = 1) { 9 lines(c(x, x), c(0, y), lty = lty, col = col) }  n2<-data$n[1]  mm<-data$mm[1]  our<-data$four[1] 5 ## Four characters in each unit interval of x  nmid <- round(four * 4) ( nn <- array(0, 2 * nmid + 1) )######################################### > z <- mu+seq(from=-3.4*sigma, to=3.4*sigma, length=mm)  atx<-pretty(z) # qz <- pnorm((z - mu)/sigma) # dz <- dnorm((z - mu)/sigma)  chw <- sigma/four " chh <- strheight("O")*0.75 htfac <- (mm * chh)/four  if(nrows==1&&ncols==1) lines(z, dz * htfac) " if(nrows==1)axis(1,at=atx) * y <- rnorm(mm, mu, sigma/sqrt(n2)) + pos <- round((y - mu)/sigma * four)  for(i in 1:mm) { 6 nn[nmid + pos[i]] <- nn[nmid + pos[i]] + 1 xpos <- chw * pos[i] A text(mu + xpos, nn[nmid + pos[i]] * chh - chh/4, "x") }  } B panelplot(xy,panel=panel.mean,totrows=totrows,totcols=totcols,  oma=rep(2.5,4))  if(device!="")dev.off()  cat(titl, "\n") } mu@$ sigma?i?.@I four@ nrep@ seed@5@  pscript F totrows?  srcref  srcfile encoding native.enc timestampAҳC@ class POSIXt POSIXct filename /tmp/johnm.g3.5.R wd /Users/johnm/r/ch3 srcfile srcref %% srcref $$ srcref  .. srcref   srcref  srcref $$ srcref $$ srcref 77 srcref 88 srcref 66 srcref    srcref  srcref  srcref NN srcref  @ srcref AB srcref CC srcref DD srcref  psfile @ ?ffffff pointsize@  ! is.null set.seed q2 . Each panel shows a simulated distribution of .  values  from a normal distribution  with mean = 10 and sd = 1. 6 The underlying theoretical normal curve is overlaid  on the lower left panel.   @?&" ??    floor sqrt  totcols ceiling / u range\Sm&@ 333333@ 333333@I  xy data.frame%%i%i mm%.%  fac factor2 Simulation?  lev2 Simulation? ?u @?  split  lapply?j function+ (((( srcref )L)L srcref as.list [[+ list??@@ /function(i){c(as.list(xy[[i]]), list(xlim=xlim, M ylim=ylim))}  panel.mean9@$? n2?@Y nrows ncols ... "" srcref # &  srcref ' '  srcref ( (  srcref ) )  srcref + +  srcref , ,$ $ srcref . .= = srcref / /  srcref 0 0# # srcref 1 1# # srcref 2 2  srcref 3 3" " srcref 4 4  srcref 5 6  srcref 7 7" " srcref 8 8* * srcref 9 9+ + srcref ; ?  srcref  vlineg?? $ $  srcref % %9 9 srcref_gg function(x, y, lty = 1, col = 1) { 9 lines(c(x, x), c(0, y), lty = lty, col = col) } !"9i? !"9?  our!"9?  nmid@  nn arrayS@? uSv&@ 333333@ 333333j  atx\u  qzm&u  dzwm&u O  Q O?  htfacm  &&(?(?_u(? axis? at  rnorm  posm&I+? ;; srcref < <6 6 srcref = =  srcref > >A A srcref !S!+S!S!+?  xposO!+$S&!S!+  @ x! +function(data, mu = 10, sigma = 1, n2 = 1, 7 mm = 100, nrows, ncols, ...)  { 1 vline <- function(x, y, lty = 1, col = 1) { 9 lines(c(x, x), c(0, y), lty = lty, col = col) }  n2<-data$n[1]  mm<-data$mm[1]  our<-data$four[1] 5 ## Four characters in each unit interval of x  nmid <- round(four * 4) ( nn <- array(0, 2 * nmid + 1) )######################################### > z <- mu+seq(from=-3.4*sigma, to=3.4*sigma, length=mm)  atx<-pretty(z) # qz <- pnorm((z - mu)/sigma) # dz <- dnorm((z - mu)/sigma)  chw <- sigma/four " chh <- strheight("O")*0.75 htfac <- (mm * chh)/four  if(nrows==1&&ncols==1) lines(z, dz * htfac) " if(nrows==1)axis(1,at=atx) * y <- rnorm(mm, mu, sigma/sqrt(n2)) + pos <- round((y - mu)/sigma * four)  for(i in 1:mm) { 6 nn[nmid + pos[i]] <- nn[nmid + pos[i]] + 1 xpos <- chw * pos[i] A text(mu + xpos, nn[nmid + pos[i]] * chh - chh/4, "x") }  } panelplot panelE%@@ 54q   g3.6 @function(mu = 10, sigma = 1, n = 1, m = 50, four = 4, nrep = 5, 5 seed=21, fignum=3.6, totrows=1, device="") { : if(device!="")hardcopy(width=5.2, height=1.6, trellis=T, 7 device=device, pointsize=c(8,6)) " if(!is.null(seed))set.seed(seed)  if(is.null(totrows)) totrows <- floor(sqrt(nrep)) " totcols <- ceiling(nrep/totrows) 5 z <- range(pretty(mu + (c(-3.4, 3.4) * sigma), 50)) xlim<-z D xy <- data.frame(y = rnorm(50*nrep, 10, 1), fac=factor(rep(1:nrep, @ rep(50, nrep)))) J qq <- qqmath(~y|fac, data=xy, par.strip.text=list(cex=0), layout=c(5,1), 2 xlab="",ylab="", aspect=1, cex=0.5) print(qq)  if(device!="")dev.off() }@$?i?.@I@@@5@ ?    native.encAҳC8 POSIXt POSIXct /tmp/johnm.g3.6.R /Users/johnm/r/ch3 srcfile srcref 77 srcref "" srcref    srcref   "" srcref   55 srcref     srcref  @@ srcref 22 srcref    srcref  srcref  @ ? trellis~@ @   u\Sm&@ 333333@ 333333@I ?u @I@$?%?%@I  qq qqmath< |9 par.strip.text0 layout@?   aspect?0? print 5 g3.7 function(device=""){ < if(device!="")hardcopy(width=4.5, height=2.6, trellis=T, ) device=device)  library(DAAG); data(pair65)  attach(pair65) ? qreference(test = heated-ambient, nrep=8, nrows=2, xlab="")  detach(pair65)  if(device!="")dev.off()  }$   @ @~789 pair65 attach qreference test& heated ambient@ @  detach 5 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  path getwd  pathtag strsplit /ch fixed ? 2!j .R3  23 /  objnames objects pattern ^g envir sys.frame hardcopy4  Dump to file:   dump 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()  } 7 lattice7 grid  @ @  color  trellis.par.set layout.heights key.top? axis.top?  hp1.plt1 xyplot< o2 wattsPerKg groups id9 humanpower1 subscripts TU<  hat fittedT panel.superposex 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' Oxygen intake ( ml.min&? .kg&? ) legend top fun textGrob args label A position?Q? T lme< random<9  hp1.plt2<9 TU< Tx 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  B T?\(?? newpage  trellis.focus panel row? column?  arglist trellis.panelArgs"""x lg@ trellis.unfocus 5 gsave 7function(fnam=NULL, prefix="~/r-book/ed2/figures/figs", @ splitchar="/ch", xtras=c("renum.fun","renum.files")){  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)  } ~/r-book/ed2/figures/figs /ch xtras renum.fun renum.files   /ch ? 2!j .RData3  23 /  ^g4  Dump to file:   save file 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()  } trellis.device postscript  test.eps?  zz.d dotplot< variety yield9 barley grid.text ABC pushViewport viewport grid.layout@? layout.pos.row?  upViewport@  popViewport@5 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@4@@5@  doit   /ch ?  endbit!j 23  23 I+?j  txt2 .!+  <-  .!+3  eval parse$ renum.files Nfunction(from.prefix="~/r-book/second/Art/", to.prefix="~/r-book/second/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(txt)  }  } ~/r-book/second/Art/ ~/r-book/second/Art/@4@@5@    /ch ? !j( nchar@  chap2 $2 03  2$ -3  2$ -3 I+?j{!+@"  ltext2 0!+3  %2!+{!+@"  rtext2 0!+3  &2!+ 2 mv % .eps  & .eps3   backup2 cp % .eps   archive3  system'(