Chapter 2: Styles of Data Analysis

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()

plot of chunk unnamed-chunk-1

g2.3()

plot of chunk unnamed-chunk-1

g2.4()

plot of chunk unnamed-chunk-1

g2.5()

plot of chunk unnamed-chunk-1

g2.6()

plot of chunk unnamed-chunk-1

g2.7()

plot of chunk unnamed-chunk-1

g2.8()

plot of chunk unnamed-chunk-1

g2.9()

plot of chunk unnamed-chunk-1

g2.10()

plot of chunk unnamed-chunk-1

g2.11()

plot of chunk unnamed-chunk-1

g2.12()

plot of chunk unnamed-chunk-1

g2.13()

plot of chunk unnamed-chunk-1

g2.14()

plot of chunk unnamed-chunk-1

g2.15()
## [1] 0.878 0.928

plot of chunk unnamed-chunk-1