g3.1 <- function(fignum=3.1, device="") { if(device!="")hardcopy(width=2, height=2, device=device) 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") curve(0.5*9.8*x^2, from=0, to=5, xaxs="i", yaxs="i", xlab = "Time (sec)", ylab = "Distance fallen (m)") chh <- par()$cxy[2] text(rep(0.2,2), 107.5-c(0,1.75*chh), expression("Distance " == " "*frac(1,2)*phantom(i)*g*t^2, "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() } 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), xlab = "Weight of roller (t)", ylab = "Depression (mm)", pch = 16) abline(0, 2.5) figtxt <- paste("Depression (mm), versus weight of roller.", sep = "") cat("\n",figtxt,"\n") if(device!="")dev.off() } g3.3 <- function(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")) { plot(x, y, xlab = "Roller weight (t)", ylab = "Depression in lawn (mm)", pch = 4) 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) 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 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 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) 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.65, col="gray30") 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.65,col="gray30") } titl <- "Lawn depression (mm) versus roller weight (t)." titl <- paste(titl, "\nIn (a) a line has been fitted, while in (b) the", "\nloess method was used to fit a smoothed curve.", "\nResiduals (the `rough') appear as vertical lines. ", "\nPositive residuals are black lines, while negative ", "\nresiduals are dotted.") if(device!="")dev.off() invisible() } g3.4 <- function(device="") { if(device!="")hardcopy(width=2, height=2, device=device) 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) plot(z, dnorm(z), type="l", ylab="normal density", yaxs="i", bty="L", xlab="Distance, in SDs, from the mean") polygon(c(z[z <= 1.0],1.0),c(dnorm(z[z <= 1.0]), dnorm(-3)), col="grey") chh <- par()$cxy[2] arrows(-1.8, 0.32, -0.25, 0.2, length=0.07, xpd=T) cump <- round(pnorm(1), 3) text(-1.8, 0.32+0.75*chh, paste("pnorm(1)\n", "=", cump), xpd=T, cex=0.8) if(device!="")dev.off() invisible() } g3.5 <- function(mu = 10, sigma = 1, n = 1, m = 50, four = 4, nrep = 5, seed=21, fignum=3.5, pscript=F, totrows=1, device="") { if(device!="")psfile(fignum, width=5.2, height=1.4, device=device, pointsize=8) if(!is.null(seed))set.seed(seed) titl <- paste("\nEach panel shows a simulated distribution of ", m, " values ", "from a normal \ndistribution ", "with mean = 10 and sd = 1.", " The underlying theoretical \nnormal curve is overlaid", " on the lower left panel.") 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) 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), mm=rep(m,nrep),four=rep(four,nrep)) fac <- factor( paste("Simulation", 1:nrep), lev <- paste("Simulation", 1:nrep) ) xlim<-z ## ylim<-c(0,dnorm(0)*sqrt(n)) ylim <- c(0,1) xy <- split(xy,fac) xy<-lapply(1:length(xy),function(i){c(as.list(xy[[i]]), list(xlim=xlim, ylim=ylim))}) panel.mean <- function(data, mu = 10, sigma = 1, n2 = 1, mm = 100, nrows, ncols, ...) { vline <- function(x, y, lty = 1, col = 1) { lines(c(x, x), c(0, y), lty = lty, col = col) } n2<-data$n[1] mm<-data$mm[1] our<-data$four[1] ## 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) { nn[nmid + pos[i]] <- nn[nmid + pos[i]] + 1 xpos <- chw * pos[i] text(mu + xpos, nn[nmid + pos[i]] * chh - chh/4, "x") } } panelplot(xy,panel=panel.mean,totrows=totrows,totcols=totcols, oma=rep(2.5,4)) if(device!="")dev.off() cat(titl, "\n") } g3.6 <- function(mu = 10, sigma = 1, n = 1, m = 50, four = 4, nrep = 5, seed=21, fignum=3.6, totrows=1, device="") { if(device!="")hardcopy(width=5.2, height=1.6, trellis=T, 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) z <- range(pretty(mu + (c(-3.4, 3.4) * sigma), 50)) xlim<-z xy <- data.frame(y = rnorm(50*nrep, 10, 1), fac=factor(rep(1:nrep, rep(50, nrep)))) qq <- qqmath(~y|fac, data=xy, par.strip.text=list(cex=0), layout=c(5,1), xlab="",ylab="", aspect=1, cex=0.5) print(qq) if(device!="")dev.off() } 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() } 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) } 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")){ 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]))) }