Packages required: “DAAG”, “lattice”, “grid”, “MASS”, “vcd”
The script that follows is designed to be executed as it stands. It sets up
functions that create Ch 2 graphs, then runs those functions. If executed by
clicking on the RStudio “Compile Notebook …” command, it will be processed
through R Markdown, generating a document that includes code, output and
graphs generated by executing the functions. If some needed packages are
missing, warning messages will appear. See further:
http://rmarkdown.rstudio.com/r_notebook_format.html
g2.1 <-
function(device="", path="~/dbeta/Art/")
{
if(device!="")hardcopy(width=5.5, height=2, path=path,
device=device)
oldpar <- par(mar=c(4.1,3.6,3.1,0.4), tcl=-0.3,
mgp=c(2.5,0.75,0), pty="s", cex=0.65)
on.exit(par(oldpar))
if(!require(DAAG))return("Package 'DAAG' is not installed -- cannot proceed")
ftotlngth <- with(possum, totlngth[sex == "f"])
dens <- density(ftotlngth)
xlim <- range(dens$x)
ylim <- range(dens$y)
par(fig=c(0,0.26,0,1))
hist(ftotlngth, breaks = 72.5 + (0:5) * 5, xlim=xlim,
ylim = c(0, 22), xlab="Total length (cm)", main="")
mtext(side=3, line=2.1, "A", at=67.5, cex=0.75)
mtext(side=3, line=0.8, text ="Breaks at 72.5, 77.5, ...", cex=0.65)
par(fig=c(0.22,0.48,0,1), new=TRUE)
hist(ftotlngth, breaks = 75 + (0:5) * 5, xlim=xlim, ylab="",
ylim = c(0, 22), xlab="Total length (cm)", yaxt="n", main="")
mtext(side=3, line=2.1, "B", at=67.5, cex=0.75)
mtext(side=3, line=0.8, text="Breaks at 75, 80, ...", cex=0.65)
par(fig=c(0.52, 0.78,0,1), new=TRUE)
hist(ftotlngth, breaks = 72.5 + (0:5) * 5, freq=F,
probability = T, xlim
= xlim, ylim = ylim, xlab="Total length (cm)", main="")
mtext(side=3, line=2.1, "C", at=67.5, cex=0.75)
mtext(side=3, line=0.8, text="Breaks as in A", at=81, cex=0.65)
lines(dens)
par(fig=c(0.74,1,0,1), new=TRUE)
hist(ftotlngth, breaks = 75 + (0:5) * 5, freq=F,
probability = T, xlim = xlim, ylim = ylim, yaxt="n",
xlab="Total length (cm)", ylab="", main="")
mtext(side=3, line=2.1, "D", at=67.5, cex=0.75)
mtext(side=3, line=0.8, text="Breaks as in B", at=81, cex=0.65)
lines(dens)
par(fig=c(0,1,0,1))
if(device!="")dev.off()
invisible()
}
g2.3 <-
function(dset = possum, x = totlngth, here = possum$sex ==
"f", device="", ytex=0.85, cex.jm=0.825)
{
yglim <- c(0,1)
if(device!="")hardcopy(device=device, width=4.25, height=2)
else yglim <- c(0.2,0.675)
dname <- as.character(substitute(dset))
xnam <- as.character(substitute(x))
x <- dset[here, xnam]
n <- length(x)
if(dname == "possum") {
xlab <- switch(xnam,
totlngth = "Total length (cm)",
pes = "Length of foot (cm)")
}
else xlab <- xnam
oldpar <- par(mar = c(4.1, 0.6, 0.6, 0.6), mgp=c(2.25,0.5,0), cex=cex.jm)
on.exit(par(oldpar))
z <- boxplot(list(val = x), plot = F)
xlim <- range(c(z$stats, z$out))
xlim <- xlim + c(-0.025,0.05) * diff(xlim)
ylim <- c(.55,1.5)
par(fig=c(0, 0.665, yglim[1], yglim[2]))
plot.new()
plot.window(xlim, ylim)
top <- 0.7
bxp(z, at=top, boxwex = 0.15, xlab = "", xlim=xlim, ylim=ylim, horiz=TRUE,
add=TRUE)
chh <- par()$cxy[2]
chw <- par()$cxy[1]
text(z$stats[5], top+0.35*chh,
"Largest value \n(there are no outliers)", adj = 0,
srt=90)
text(z$stats[4], top+0.65*chh, "upper quartile", adj = 0, srt=90)
text(z$stats[3], top+0.65*chh, "median", adj = 0, srt=90)
text(z$stats[2], top+0.65*chh, "lower quartile", adj = 0, srt=90)
text(z$stats[1], top+0.35*chh, "Smallest value \n(outliers excepted)", adj = 0, srt=90)
if(!is.null(z$out)) text(z$out[1], top+0.35*chh, "Outlier", adj = 0, srt=90)
# lines(c(90, 90), z$stats[c(2, 4)])
av <- mean(z$stats[c(2, 4)])
q1 <- z$stats[2]
q3 <- z$stats[4]
axis(1, at = c(q1, q3), tck = 0.02, labels = F)
botm<-par()$usr[3]
text(c(q1, q3), rep(top-chh,2),
c(format(round(q1, 2)), format(round(q3, 1))),adj=0.5)
qtop <- q3 + 0.5 * chh
mtext(side=1,line=2.15, xlab, cex=cex.jm)
par(fig=c(0.675, 1, yglim[1], yglim[2]), new=T)
plot(0:1, 0:1, bty="n", axes=F, xlab="", ylab="", type="n", tcl=-0.3)
text(0, ytex, "Inter-quartile range", adj = 0)
text(0.15, ytex - 1.15 * chh, paste("= ", format(round(q3, 2)), "-",
format(round(q1, 2)), "\n= ",
format(round(q3 - q1, 2))),
adj = 0)
here <- !is.na(x)
sd <- sqrt(var(x[here]))
text(0, ytex - 5 * chh, paste("Compare\n", "0.75 x Inter-quartile range",
"\n =", format(round(0.75 * (q3 - q1), 1)),
"\nwith", "standard deviation\n =",
format(round(sd, 1))), adj = 0)
n <- sum(here)
if(device!="")dev.off()
par(fig=c(0,1,0,1))
par(oldpar)
invisible()
}
g2.4 <-
function(device="")
{
if(device!="")hardcopy(device=device, width=4.8, height=3.8)
oldpar <- par(mar=c(3.6,5.1,2.1,1.1), oma=c(0.5,0,0,0),
mgp=c(3.5,0.75,0), tcl=-0.4)
on.exit(par(oldpar))
par(fig=c(0, 1, .4, 1))
plot(log(measles,10), xlab="", ylim=c(0,log(5000*1000, 10)),
ylab=" Deaths; Population (log scale)", yaxt="n")
ytikpoints <- c(1, 10, 100,1000, 1000000, 5000000)
axis(2, at=log10(ytikpoints), labels=paste(ytikpoints), cex=.75, las=2)
londonpop <- ts(c(1088,1258,1504,1778,2073,2491,2921,3336,3881,4266,
4563,4541,4498,4408), start=1801, end=1931, deltat=10)
points(log(londonpop*1000,10), pch=16, cex=.5)
mtext(side=3, line=0.5, "A (1629-1939)", adj=0)
par(fig=c(0, 1, 0, .45), mar=c(2.1,5.1,3.1,1.1), new=TRUE)
plot(window(measles, start=1840, end=1882), xlab="Year", yaxt="n",
ylim=c(0,4600), ylab="Deaths; Population in 1000s")
points(londonpop, pch=16, cex=.5)
axis(2, at=(1:4)*1000, cex=.75, las=2)
mtext(side=3, line=0.5, "B (1841-1881)", adj=0)
par(fig=c(0,1,0,1))
if(device!="")dev.off()
invisible()
}
g2.5 <-
function(device=""){
if(device!="")hardcopy(device=device, width=2, height=2)
oldpar <- par(mar = c(4.1,4.1,1.1,1.1), mgp=c(2.5,0.75,0), pty="s")
on.exit(par(oldpar))
xyrange <- range(milk)
plot(four ~ one, data = milk, xlim = xyrange, ylim = xyrange, tcl=-0.3,
pch = 16, cex=.65)
rug(milk$one, ticksize=0.04, col="gray40")
rug(milk$four, side = 2, ticksize=0.04, col="gray40")
abline(0, 1)
if(device!="")dev.off()
}
g2.6 <-
function(df=fruitohms, device=""){
if(device!="")hardcopy(width=4, height=2,
device=device)
oldpar<-par(mar=c(4.1,3.6, 1.6, 0.6), mgp=c(2.25, 0.5,0), pty="s",
oma=c(0,1.1,0,1.1), tcl=-0.4
)
on.exit(par(oldpar))
par(mgp=c(2.75,1,0))
par(mfrow=c(1,2))
plot(ohms~juice, data=df, cex=0.8, xlab="Apparent juice content (%)",
ylab="Resistance (kOhm)", yaxt="n")
mtext(side=3, line=0.25, "A", adj=0)
axis(2, at=(1:5)*2000, labels=paste((1:5)*2))
plot(ohms~juice,data=df,cex=0.8,xlab="Apparent juice content (%)",
ylab="",yaxt="n")
mtext(side=3, line=0.25, "B", adj=0)
axis(2, at=(1:5)*2000, labels=paste((1:5)*2))
lines(lowess(fruitohms$juice,fruitohms$ohms), lwd=2, col="gray40")
if(device!="")dev.off()
}
g2.7 <-
function(dset = Animals, show = "lines", device="", color = F)
{
library(MASS, quietly=TRUE)
if(!require(MASS))return("Package 'MASS' is not installed -- cannot proceed")
if(device!="")hardcopy(width=4.5, height=2.25,
device=device)
oldpar <- par(mfcol = c(1, 2), mar = c(3.1,3.1,2.1,3.1),
oma=c(0,1.1,0,1.1),
mgp = c(2.25, 0.65, 0), tck=-0.03, pty="s")
on.exit(par(oldpar))
fig1txt <- paste("(a) Untransformed scale")
fig2txt <- paste("(b) Logarithmic scale, both axes")
xlab <- "Body weight (kg x 100)"
ylab <- "Brain weight (g)"
dset$body <- dset$body/100
plot(brain ~ body, data=dset, xlab = xlab, ylab = ylab, type = "n")
points(dset$body, dset$brain)
mtext(side=3, line=1, "A", adj=-0.1)
eqscplot(log10(dset$body), log10(dset$brain), pch = 1, axes = F, xlab =
xlab, ylab = ylab)
mtext(side=3, line=1, "B", adj=-0.1)
xpos <- sort(unique(round(log10(dset$body))))
ypos <- sort(unique(round(log10(c(0.1,dset$brain)))))
lab <- paste(10^xpos)
par(cex=0.85)
axis(1, at = xpos, label = lab,cex=.75)
axis(3, at = xpos)
axis(4, at = ypos)
par(mgp = c(2.5, 0.75, 0))
axis(2, at = ypos, label = paste(10^ypos, sep = ""), srt = 90)
mtext(side = 3, line = 2, "log10(Body weight)")
mtext(side = 4, line = 2, "log10(Brain weight)")
box()
if(device!="")dev.off()
}
g2.8 <-
function(dset = cuckoos, device="")
{
if(device!="")hardcopy(width=5.5, height=2,
device=device, trellis=TRUE, pointsize=c(8, 5))
## Two lattice graphs on one page: data frame cuckoos (DAAG)
if(!require(DAAG))return("Package 'DAAG' is not installed -- cannot proceed")
if(!require(grid))return("Package 'grid' is not installed -- cannot proceed")
if(!require(lattice))return("Package 'lattice' is not installed -- cannot proceed")
trellis.par.set(layout.heights=list(key.top=0.5, axis.top=0.6,
bottom.padding=0.25))
nam <- levels(cuckoos$species)
splitnam <- strsplit(nam,"\\.")
newnam <- sapply(splitnam, function(x)if(length(x)==1)x else
paste(x,collapse=" "))
cuckoos.strip <- stripplot(species ~ length, data=cuckoos,
factor.levels=splitnam, scales=list(tck=.5),
xlab=list("Length of egg (mm)", cex=0.9),
legend=list(top=list(fun=textGrob,
args=list(label="A", x=0,
gp=gpar(cex=0.9), just="left"))))
print(cuckoos.strip, position=c(0,0,0.575,1))
cuckoos.bw <- bwplot(species~length, factor.levels=splitnam,
xlab=list("Length of egg (mm)", cex=0.9),
scales=list(tck=0.5, y=list(alternating=0)),
data=cuckoos,
legend=list(top=list(fun=textGrob,
args=list(label="B", x=0,
gp=gpar(cex=0.9), just="left"))))
print(cuckoos.bw, newpage=FALSE, position=c(0.525,0,1,1))
if(device!="")dev.off()
invisible()
}
g2.9 <-
function(width=4.5, height=2, pointsize=c(7,4), device=""){
if(device!="") hardcopy(width=width, height=height, trellis=T,
color=TRUE, pointsize=pointsize, device=device)
if(!require(DAAG))return("Package 'DAAG' is not installed -- cannot proceed")
gph <- densityplot(~earconch | sex, groups=Pop, data=possum,
auto.key=list(space="right"), aspect=1,
scales=list(tck=0.5),
par.settings=simpleTheme(col=c("gray40","black")))
print(gph)
if(device!="")dev.off()
}
g2.10 <-
function(device=""){
if(device!="")hardcopy(width=5, height=4.5, color=FALSE, device=device,
trellis=TRUE, pointsize=7)
trellis.par.set(layout.heights=list(top.padding=0, bottom.padding=0,
sub=0, xlab=0))
simplejobsA.xyplot <-
xyplot(Ontario+Quebec+BC+Alberta+Prairies+Atlantic ~ Date,
outer=FALSE, data=jobs, type="b",
ylab="Number of workers", scales=list(y=list(log="e")),
auto.key=list(space="right", lines=TRUE, points=TRUE))
## Update trellis object to use improved x- and y-axis tick labels
datelabpos <- seq(from=95, by=0.5, length=5)
datelabs <- format(seq(from=as.Date("1Jan1995", format="%d%b%Y"),
by="6 month", length=5), "%b%y")
## Now create $y$-labels that have numbers, with log values underneath
ylabpos <- exp(pretty(log(unlist(jobs[,-7])), 5))
ylabels <- paste(round(ylabpos),"\n(", log(ylabpos), ")", sep="")
jobsA.xyplot <- update(simplejobsA.xyplot, xlab="",
scales=list(tck=0.5, x=list(at=datelabpos,
labels=datelabs),
y=list(at=ylabpos, labels=ylabels)))
jobsB.xyplot <-
xyplot(Ontario+Quebec+BC+Alberta+Prairies+Atlantic ~ Date,
data=jobs, type="b", layout=c(3,2), ylab="Number of jobs",
scales=list(y=list(relation="sliced", log=TRUE)),
outer=TRUE)
## Update trellis object to use improved x- and y-axis tick labels
datelabpos <- seq(from=95, by=0.5, length=5)
datelabs <- format(seq(from=as.Date("1Jan1995", format="%d%b%Y"),
by="6 month", length=5), "%b%y")
## Now create $y$-labels that have numbers, with log values underneath
ylabpos <- exp(pretty(log(unlist(jobs[,-7])), 100))
ylabels <- paste(round(ylabpos),"\n(", log(ylabpos), ")", sep="")
## Take names of provinces in order of number of Jobs in Dec96
jobsB.xyplot <-
update(jobsB.xyplot, xlab="", between=list(x=0.25, y=0.25),
scales=list(x=list(at=datelabpos, labels=datelabs),
y=list(at=ylabpos, labels=ylabels), tck=0.6) )
plot.new()
oldpar <- par(mar=c(2.6,3.6,2.6,1.6), fig=c(0,1,0.6,1), mgp=c(2.0,0.5,0))
on.exit(par(oldpar))
mtext(side=3, line=1, "A", adj=-0.15, cex=0.75)
par(fig=c(0,1,0,0.575), mar=c(2.6,2,2.6,1.6), mgp=c(1.1,0.5,0), new=TRUE)
mtext(side=3, line=1.5, "B", adj=-0.15, cex=0.75)
if(!require(DAAG))return("Package 'DAAG' is not installed -- cannot proceed")
print(jobsA.xyplot, position=c(0.1,0.55,0.9,1), newpage=FALSE)
print(jobsB.xyplot, position=c(0,0,1,0.575), newpage=FALSE)
if(device!="")dev.off()
}
g2.11 <-
function(device="", pointsize=8, fonts=c("sans","mono")){
if(!require(DAAG))return("Package 'DAAG' is not installed -- cannot proceed")
if(device!="")hardcopy(width=5.5, height=1.2,
pointsize=pointsize, fonts=fonts,
device=device)
oldpar <- par(mar=c(2.1, 0.6, 0.6 ,0.6), mgp=c(2, 0.5,0), cex=0.75)
on.exit(par(oldpar))
plot(c(1230,1860), c(0, 10), axes=FALSE, xlab="", ylab="", type="n", log="x")
xpoints <- c(1366, 1436, 1752, 1840)
axis(1, at=xpoints, labels=FALSE, tck=0.01, lty=1, col="gray40")
for(i in 1:4){
axis(1, at=xpoints[i], labels=substitute(italic(a), list(a=paste(xpoints[i]))), line=-2.25, lty=0,
cex=0.8)
lines(rep(xpoints[i],2), c(0, 0.15*par()$cxy[2]), lty=1)
}
axpos <- 1250*cumprod(c(1, rep(1.2,2)))
lab <- round(axpos)
axis(1, at=axpos, labels=lab)
lab2 <- lapply(round(log2(xpoints),3), function(x)substitute(2^a, list(a=x)))
axis(1, at=xpoints, labels=as.expression(lab2), line=-3.5, lty=0)
labe <- lapply(format(round(log(xpoints),3)), function(x)substitute(e^a, list(a=x)))
axis(1, at=xpoints, labels=as.expression(labe), line=-5, lty=0)
lab10 <- lapply(round(log10(xpoints),3), function(x)substitute(10^a, list(a=x)))
axis(1, at=xpoints, labels=as.expression(lab10), line=-6.5, lty=0)
## mtext(side=1, line=2.25, "Number of workers")
par(family="mono", xpd=TRUE)
axis(1, at=1220, labels="log=2", line=-3.5, lty=0, hadj=0)
axis(1, at=1220, labels='log="e"', line=-5, lty=0, hadj=0)
axis(1, at=1220, labels="log=10", line=-6.5, lty=0, hadj=0)
wid2 <- strwidth("log=2")
par(family="sans")
## axis(1, at=1240*10^wid2, labels=expression(italic(" (lattice)")), line=-3.5, lty=0, hadj=0)
if(device!="")dev.off()
}
g2.12 <-
function(device="", pointsize=c(8,4)){
if(!require(DAAG))return("Package 'DAAG' is not installed -- cannot proceed")
if(device!="")hardcopy(width=5.5, height=3.25,
pointsize=pointsize,
device=device, trellis=TRUE, color=FALSE)
par(fig=c(0, 0.515, 0, 1))
plot.new()
mtext(side=3, line=2.5, "A", adj=-0.225, cex=0.75)
par(fig=c(0.455, 1, 0, 1), new=TRUE)
mtext(side=3, line=2.5, "B", adj=-0.225, cex=0.75)
colr <- c("gray40","black")
plotchar <- c(1, 16)
targplot <- xyplot(csoa ~ it|sex*agegp, data=tinting, groups=target,
col=colr, pch=plotchar, aspect=1,
key=list(space="top", columns=2,
points=list(pch=plotchar, col=colr),
text=list(levels(tinting$target))),
scale=list(y=list(alternating=1),tck=0.5))
print(targplot, position=c(0, 0, 0.515, 1), newpage=FALSE)
colr <- c("skyblue1", "skyblue4")[c(2,1,2)]
plotchar <- c(1,16,16) # open, filled, filled
tintplot <- xyplot(csoa~it|sex*agegp, data=tinting, groups=tint,
col=colr, pch=plotchar, aspect=1,
type=c("p","smooth"), span=1.25,
key=list(space="top", columns=3,
points=list(pch=c(1,16,16), col=colr),
text=list(levels(tinting$tint), col=colr)),
scale=list(y=list(alternating=2), tck=0.5), ylab="")
print(tintplot, position=c(0.485, 0, 1, 1), newpage=FALSE)
if(device!="")dev.off()
}
g2.13 <-
function(device=""){
if(device!="")hardcopy(width=2.15, height=2.15,
device=device)
oldpar <- par(mar=c(2.1,2.1, 1.1, 1.1), xpd=TRUE, pty="s")
on.exit(par(oldpar))
stones <- array(c(81,6,234,36,192,71,55,25), dim=c(2,2,2),
dimnames=list(Sucess=c("yes","no"),
Method=c("open","ultrasound"),
Size=c("<2cm\n", ">=2cm\n")))
if(!require(vcd))return("Package 'vcd' is not installed -- cannot proceed")
mosaic(aperm(stones, 3:1), main=NULL,
labeling_args=list(gp_labels=gpar(fontsize=7),
gp_varnames=gpar(fontsize=8)),
legend_args=list(fontsize=7))
if(device!="")dev.off()
}
g2.14 <-
function(device="", pointsize=c(8,6))
{
if(device!=""){hardcopy(width=5, height=2.5, trellis=T,
color=F, pointsize=pointsize, device=device)
}
if(!require(DAAG))return("Package 'DAAG' is not installed -- cannot proceed")
kiwishade$block <- factor(kiwishade$block, levels=c("west","north","east"))
gph <- dotplot(shade~yield | block, data=kiwishade, col="gray", aspect=1,
panel=function(x,y,...){panel.dotplot(x,y,...)
av <- sapply(split(x,y),mean);
ypos <- unique(y)
lpoints(ypos~av, pch=3, col="black")},
key=list(space="top", columns=2,
text=list(c("Individual vine yields",
"Plot means (4 vines)")),
points=list(pch=c(3,3), col=c("gray","black"))),
layout=c(3,1))
print(gph)
if(device!="")dev.off()
invisible()
}
g2.15 <-
function(device="", seed=21)
{
par(pty="s")
if(!is.null(seed))set.seed(seed)
if(device!="")hardcopy(width=4.5, height=1.25,
device=device)
titl <- paste("Different relationships between y and x.", sep = "")
x1 <- x2 <- x3 <- (11:30)/5
y1 <- x1 + rnorm(20)/2
y2 <- 2 - 0.05 * x1 + 0.1 * ((x1 - 1.75))^4 + 1.25 * rnorm(20)
r <- round(cor(x1, y2), 3)
rho <- round(cor(rank(x1), rank(y2)), 3)
print(c(r, rho))
y3 <- (x1 - 3.85)^2 + 0.015 + rnorm(20)/4
theta <- ((2 * pi) * (1:20))/20
x4 <- 10 + 4 * cos(theta)
y4 <- 10 + 4 * sin(theta) + (0.5 * rnorm(20))
r1 <- cor(x1, y1)
xy <- data.frame(x = c(rep(x1, 3), x4), y = c(y1, y2, y3, y4),
gp = rep(1:4, rep(20, 4)))
xy<-split(xy,xy$gp)
xlimdf<-lapply(list(x1,x2,x3,x4),range)
ylimdf<-lapply(list(y1,y2,y3,y4),range)
xy<-lapply(1:4,function(i,u,v,w){list(xlim=v[[i]], ylim=w[[i]],
x=u[[i]]$x, y=u[[i]]$y)},
u=xy,v=xlimdf,w=ylimdf)
panel.corr<-function(data,...){
x<-data$x
y<-data$y
points(x, y, pch = 16)
chh <- par()$cxy[2]
x1 <- min(x)
y1 <- max(y) - chh/8
r1 <- cor(x, y)
text(x1, y1, paste(round(r1, 3)), cex = 1.0, adj = 0)
}
panelplot(xy,panel=panel.corr,totrows=1,totcols=4,
oma=rep(1,4))
titl <- paste(titl, " In the lower right \npanel, the ",
"Pearson correlation is ", r,
", while the Spearman rank \ncorrelation is ",
rho, ".", sep = "")
if(device!="")dev.off()
par(mfrow=c(1,1), pty="m")
invisible()
}
pkgs <- c("DAAG", "lattice", "grid", "MASS", "vcd")
z <- sapply(pkgs, require, character.only=TRUE, warn.conflicts=FALSE, quietly=TRUE)
if(any(!z)){
notAvail <- paste(names(z)[!z], collapse=", ")
warning(paste("The following packages should be installed:", notAvail))
}
g2.1()
g2.3()
g2.4()
g2.5()
g2.6()
g2.7()
g2.8()
g2.9()
g2.10()
g2.11()
g2.12()
g2.13()
g2.14()
g2.15()
## [1] 0.878 0.928