`g4.1` <- function(mu = 10, sigma = 1, n = c(4, 9), device="") { if(device!="")hardcopy(width=4, height=2, pointsize=8, device=device) titl <- paste(": Sampling Distribution of the Mean, for Samples", "of sizes n =", paste(n, collapse = " & "), "\nfrom a Normal Population with Mean = 10 and SD = 1.", sep = "") n2 <- length(n) z <- range(pretty(mu + (c(-3.4, 3.4) * sigma), 50)) ylim<-c(0,dnorm(0)*sqrt(max(n)))*1.08 xy <- list(one=list(n=n[1],xlim=z,ylim=ylim,atx=pretty(z), aty=pretty(ylim)), two=list(n=n[2],xlim=z,ylim=ylim,atx=pretty(z),aty=NULL)) panel.dist <- function(data, mu = 10, sigma = 1, ...) { vline <- function(x, y, lty = 1, col = 1) { if(length(y) == 1) y <- c(0, y) lines(c(x, x), y, lty = lty, col = col) } nn<-data$n z <- pretty(mu + (c(-3, 3) * sigma), 50) qz <- pnorm((z - mu)/sigma) dz <- dnorm((z - mu)/sigma) lines(z, dz, yaxs = "s") chw <- strwidth("O") chh <- strheight("O")*1.25 z2 <- min(z[z > mu + 2 * sigma]) nz <- match(z2, z) text(z2 + chw/3.5, dz[nz] + chh/2, "n = 1", adj = 0, cex = 0.75) vline(mu - sigma, y = c( - chh/4, 0) + dnorm(-1)) vline(mu + sigma, y = c( - chh/4, 0) + dnorm(1)) vline(mu - sigma, y = c(chh/4, 0)) vline(mu + sigma, y = c(chh/4, 0)) vline(mu, y = dnorm(0), lty = 1) zn <- mu + (z - mu)/sqrt(nn) dzn <- dnorm(((zn - mu) * sqrt(nn))/sigma) * sqrt(nn) abline(h = 0) lines(zn, dzn, col = 2, lwd=2) vline(mu - sigma/sqrt(nn), y = dnorm(1) * sqrt(nn), col = 2) vline(mu + sigma/sqrt(nn), y = dnorm(1) * sqrt(nn), col = 2) text(mu + chw*6.5/nn, dnorm(0) * sqrt(nn)-chh/2, paste("Sampling \ndistribution \nof means of", "\nsamples of \nsize n =", nn), col = 2, adj = 0, cex = 0.7) atx <- data$atx aty <- data$aty if(!is.null(atx))text(atx,rep(-chh*0.4,length(atx)), paste(atx),cex=0.75) if(!is.null(aty))axis(2,at=aty) } panelplot(xy,panel=panel.dist,totrows=1,totcols=2,oma=c(1,4,1,2)) titl <- paste(titl, "\nThe dark curve shows the population distribution,", "\ni.e., the curve is for n=1.") mtext(side = 2, line = 3.25, "Probability density") if(device!="")dev.off() cat(titl,"\n") } `g4.10` <- function(x1=two65$ambient, x2=two65$heated, nsim=2000, plotit=T, fignum="4.10", device=""){ if(device!="")hardcopy(width=2.25, height=2.25, device=device) oldpar<-par(mar=par()$mar-c(1,0,1,0)) on.exit(par(oldpar)) library(DAAG); data(two65) n1 <- length(x1) n2<-length(x2) n<-n1+n2 x<-c(x1,x2) dbar <- mean(x2)-mean(x1) z <- array(,nsim) for(i in 1:nsim){ mn <- sample(n,n2,replace=FALSE) dbardash <- mean(x[mn]) - mean(x[-mn]) z[i] <- dbardash } pval <- (sum(z>abs(dbar)) + sum (z < -abs(dbar)))/nsim if(plotit){plot(density(z), xlab="", main="", yaxs="i", ylim=c(0,0.08), cex.axis=0.8) abline(v=dbar) abline(v=-dbar, lty=2) mtext(side=3,line=0.5, text=expression(bar(italic(x))[2]-bar(italic(x))[1]), at=dbar) mtext(side=3,line=0.5, text=expression(-(bar(italic(x))[2]- bar(italic(x))[1])), at=-dbar)} print(signif(pval,3)) ## Second try for(i in 1:nsim){ mn <- sample(n,n2,replace=FALSE) dbardash <- mean(x[mn]) - mean(x[-mn]) z[i] <- dbardash } pval <- (sum(z>abs(dbar)) + sum (z < -abs(dbar)))/nsim if(plotit){lines(density(z),lty=5,lwd=2) print(signif(pval,3)) } if(device!="")dev.off() } `g4.2` <- function(device=""){ if(device!="")hardcopy(width=4.5, pointsize=8, device=device, height=2.2) oldpar <- par(mfrow = c(1,2), mar=par()$mar-c(.5, 0, 1.5, 0), mgp=c(2.5, .75, 0), xpd=T) on.exit(par(oldpar)) x <- seq(from=-4.2, to = 4.2, length.out = 50) ylim <- c(0, dnorm(0)) ylim[2] <- ylim[2]+0.1*diff(ylim) h1 <- dnorm(x) h3 <- dt(x, 3) h8 <- dt(x,8) plot(x, h1, type="l", xlab = "", xaxs="i", ylab = "", yaxs="i", bty="L", ylim=ylim) mtext(side=3,line=0.75, "A", adj=-0.2) lines(x, h8, col="grey60") mtext(side=1, line=2.5, "No. of SEMs from mean") mtext(side=2, line=2.5, "Probability density") chh <- par()$cxy[2] topleft <- par()$usr[c(1,4)] + c(0, 0.4*chh) legend(topleft[1], topleft[2], col=c("black","grey60"), lty=c(1,1), legend=c("Normal","t (8 d.f.)"), bty="n", cex=0.8) plot(x, h1, type="l", xlab = "", xaxs="i", ylab = "", yaxs="i", bty="L", ylim=ylim) mtext(side=3,line=0.75, "B", adj=-0.2) lines(x, h3, col="grey60") mtext(side=1, line=2.5, "No. of SEMs from mean") mtext(side=2, line=2.5, "Probability density") legend(topleft[1], topleft[2], col=c("black","grey60"), lty=c(1,1), legend=c("Normal","t (3 d.f.)"), bty="n", cex=0.8) if(device!="")dev.off() } `g4.2a` <- function(fignum="4.2a", pscript=FALSE){ if(pscript)psfile(fignum=fignum, width=4.5, height=2.2) oldpar <- par(mfrow = c(1,2), mar=par()$mar-c(.5, 0, 1.5, 0), mgp=c(2.5, .75, 0)) on.exit(par(oldpar)) x <- seq(from=-4.2, to = 4.2, length.out = 50) ylim <- c(0, dnorm(0)) h1 <- dnorm(x) h3 <- dt(x, 3) h8 <- dt(x,8) plot(x, h1, type="l", xlab = "", xaxs="i", ylab = "", yaxs="i", bty="L", ylim=ylim) mtext(side=3,line=1, "A", adj=-0.2) lines(x, h8, col="grey60") mtext(side=1, line=2.5, "No. of SEMs from mean") mtext(side=2, line=2.5, "Probability density") chh <- par()$cxy[2] topleft <- par()$usr[c(1,4)] + c(0, 0.75*chh) legend(topleft[1], topleft[2], col=c("black","grey60"), lty=c(1,1), legend=c("Normal","t (8 d.f.)"), bty="n", cex=0.8) plot(x, h1, type="l", xlab = "", xaxs="i", ylab = "", yaxs="i", bty="L", ylim=ylim) mtext(side=3,line=1, "B", adj=-0.2) lines(x, h3, col="grey60") mtext(side=1, line=2.5, "No. of SEMs from mean") mtext(side=2, line=2.5, "Probability density") legend(topleft[1], topleft[2], col=c("black","grey60"), lty=c(1,1), legend=c("Normal","t (3 d.f.)"), bty="n", cex=0.8) if(pscript)dev.off() } `g4.3` <- function(cump=0.975, device=""){ if(device!="")hardcopy(width=4.5, height=2.2, pointsize=8, device=device) oldpar <- par(mfrow = c(1,2), mar=par()$mar-c(.5, 0, 1.5, 0), mgp=c(2.5, .75, 0)) on.exit(par(oldpar)) x <- seq(from=-3.9, to = 3.9, length.out = 50) ylim <- c(0, dnorm(0)) plotfun <- function(cump, dfun = dnorm, qfun=qnorm, ytxt = "Normal probability density", txt1="qnorm", txt2="") { h <- dfun(x) plot(x, h, type="l", xlab = "", xaxs="i", xaxt="n", ylab = ytxt, yaxs="i", bty="L", ylim=ylim) axis(1, at=c(-2, 0), cex=0.8) axis(1, at=c((-3):3), labels=F) tailp <- 1-cump z <- qfun(cump) ztail <- pretty(c(z,4),20) htail <- dfun(ztail) polygon(c(z,z,ztail,max(ztail)), c(0,dfun(z),htail,0), col="gray") text(0, 0.5*dfun(z)+0.08*dfun(0), paste(round(tailp, 3), " + ", round(1-2*tailp,3), "\n= ", round(cump, 3), sep=""), cex=0.8) lines(rep(z, 2), c(0, dfun(z))) lines(rep(-z, 2), c(0, dfun(z)), col="gray60") chh <- par()$cxy[2] arrows(z, -2*chh,z,-0.1*chh, length=0.1, xpd=T) text(z, -3*chh, paste(txt1, "(", cump, txt2, ")", "\n= ", round(z,2), sep=""), xpd=T) x1 <- z + .3 y1 <- dfun(x1)*0.35 y0 <- dfun(0)*0.2 arrows(-2.75, y0, -x1, y1, length=0.1, col="gray60") arrows(2.75, y0, x1, y1, length=0.1) text(-2.75, y0+0.5*chh, tailp, col="gray60") text(2.75, y0+0.5*chh, tailp) } ytxt <- "t probability density (8 d.f.)" plotfun(cump=cump) mtext(side=3, line=1.5, "A", adj=-0.2) ytxt <- "t probability density (8 d.f.)" plotfun(cump=cump, dfun=function(x)dt(x, 8), qfun=function(x)qt(x, 8), ytxt=ytxt, txt1="qt", txt2=", 8") mtext(side=3, line=1.5, "B", adj=-0.2) if(device!="")dev.off() } `g4.4` <- function(device="") { if(device!="")hardcopy(width=4.75, height=2.5, pointsize=8, device=device) library(DAAG); data(pair65) titl <- paste( "Second versus first member, for each pair. The first", "\npanel is for the elastic band data. The second (from", "\nDarwin) is for plants of the species Reseda lutea") oldpar <- par(mfrow = c(1, 2), mar = par()$mar - c(.5, -0.5, .5, 0), pty="s", oma=c(0,0,0,.5),mgp = c(2, 0.5, 0)) on.exit(par(oldpar)) onesamp(dset = pair65, x = "ambient", y = "heated", xlab = "Amount of stretch (ambient)", ylab = "Amount of stretch (heated)") ## Data set mignonette holds Darwin's data on Reseda lutea, from ## Darwin, Charles 1877. The Effects of Cross and Self Fertilisation ## in the Vegetable Kingdom. Appleton and Company, New York, p.118. ## Data were in five pots, holding 5,5,5,5,4 pairs of plants ## respectively. onesamp(dset = mignonette, x = "self", y = "cross", xlab = "Height of Self-fertilised Plant", ylab = "Height of Cross-Fertilised Plant", dubious = 0) cat("\n", titl, "\n") if(device!="")dev.off() invisible() } `g4.5` <- function(device="") { if(device!="")hardcopy(width=2.4, height=2.4, pointsize=8, device=device) oldpar<-par(mar=c(4.1,4.1,1.1,1.1), mfrow=c(1,1), mgp=c(2.5,0.75,0)) on.exit(par(oldpar)) library(DAAG) x2 <- chisq.test(rareplants) x2E <- stack(data.frame(t(x2$expected))) habitat <- rep(c(1,2,3), 4) plot(x2E$values ~ habitat, xlab="habitat", ylab="Expected Number of Species",axes=FALSE, xlim=c(0.5, 3.5), pch=16) text(x2E$values ~ habitat, labels=x2E$ind, pos=rep(c(4,4,2,2),3), cex=0.85) axis(1, at=seq(1,3), labels=c("D", "W", "WD")) axis(2) box() if(device!="")dev.off() invisible() } `g4.6` <- function(device="") { if(device!="")hardcopy(width=3.0, height=1.75, device=device, trellis=TRUE) oldpar <- par(mar = c(4.1,4.1,1.1,1.1), mgp=c(2.5,0.75,0)) on.exit(par(oldpar)) tomato <- data.frame(weight= c(1.5, 1.9, 1.3, 1.5, 2.4, 1.5, # water 1.5, 1.2, 1.2, 2.1, 2.9, 1.6, # Nutrient 1.9, 1.6, 0.8, 1.15, 0.9, 1.6), # Nutrient+24D trt = rep(c("water", "Nutrient", "Nutrient+24D"), c(6, 6, 6))) tomato$trt <- relevel(tomato$trt, ref="water") gph <- with(tomato, stripplot(trt~weight, aspect=0.6)) print(gph) figtxt <- paste("Weights of tomato plants" ) cat(figtxt,"\n") if(device!="")dev.off() invisible() } `g4.7` <- function(device=""){ if(device!="")hardcopy(width=4.2, height=2.0, device=device) tomato <- data.frame(weight= c(1.5, 1.9, 1.3, 1.5, 2.4, 1.5, # water 1.5, 1.2, 1.2, 2.1, 2.9, 1.6, # Nutrient 1.9, 1.6, 0.8, 1.15, 0.9, 1.6), # Nutrient+24D trt = rep(c("water", "Nutrient", "Nutrient+24D"), c(6, 6, 6))) tomato$trt <- relevel(tomato$trt, ref="water") tomato.aov <- aov(weight~trt, data=tomato) onewayPlot(obj=tomato.aov) if(device!="")dev.off() } `g4.8` <- function(figtxt = "Interaction plot of rice data", mc = F, device="") { if(device!="")hardcopy(width=5, height=2.35, pointsize=8, device=device) oldpar <- par(mgp=c(2.5,0.75,0), pty="s", las=1) on.exit(par(oldpar)) attach(rice) par(fig=c(0,0.555,0,1), mar=c(4.1,9.1,1.1,1.1), las=1) stripchart(ShootDryMass ~ trt, pch=1, cex=1, col="gray40", xlim=c(0,140), ylim=c(.5, 6.25), xlab="Shoot dry mass (g)") box(col="gray75") av <- tapply(ShootDryMass, INDEX=trt, FUN=mean) stripchart(av ~ factor(names(av), levels=names(av)), pch=3, cex=1.5, add=TRUE) par(fig=c(0.555,1,0,1), mar = c(4.1,5.1,1.1,1.6), new=TRUE) ## Interaction plot for ShootDryMass with factors fert and variety interaction.plot(fert, variety, ShootDryMass, legend = FALSE, xlab="Fertilizer treatment", ylab="Mean of Shoot dry mass") detach(rice) box(col="gray75") if(device!="")dev.off() invisible() } `g4.9` <- function(mc = F, type="box", device="") { if(device!="")hardcopy(width=2, height=2, pointsize=8, device=device) figtxt <- "Distance travelled, relative to distance up a 20 degree ramp." oldpar <- par(mar = c(4.1,4.1,1.1,1.1), pty="s", mgp=c(2.5,0.75,0)) on.exit(par(oldpar)) library(DAAG); data(modelcars) attach(modelcars) stripchart(distance.traveled~starting.point, vertical=T, pch=16, cex=0.75, xlab = "Distance up ramp", ylab="Distance travelled") detach(modelcars) if(device!="")dev.off() invisible() } `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="") 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]))) }