Code for ‘Data Analysis And Graphics Using R’, 3rd edn, CUP, 2010

Chapter 6: Multiple Linear Regression
Packages required: “DAAG”, “lattice”, “MASS”, “MCMCpack”, “compositions”

The script that follows is designed to be executed as it stands. It sets up
functions that create Ch 6 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

g6.1 <-
function(device="")
{
  if(device!="")hardcopy(width=2.25, height=2.25, device=device)
  oldpar<-par(mar=c(4.1,4.1,2.1,1.1), mgp=c(2.5,0.75,0), pty="s")
  on.exit(par(oldpar))
  xy<-par()$usr[c(1,4)]
  if(!require(DAAG))return("Package 'DAAG' is not installed -- cannot proceed")  
  xlim <- range(allbacks$volume)
  xlim <- xlim+c(-.075,.075)*diff(xlim)
  ## Plot of weight vs volume: data frame allbacks (DAAG)
  plot(weight ~ volume, data=allbacks, pch=c(16,1)[unclass(cover)], lwd=1.25,
       xlim=xlim)
  ## unclass(cover) gives the integer codes that identify levels
  ## As text() does not accept the parameter data, use with()
  ## to specify the data frame.
  with(allbacks,text(weight ~ volume, labels=paste(1:15), cex=0.75, offset=0.35,
                     pos=c(2,4)[unclass(cover)]))
  par(xpd=TRUE)
  legend(x=660, y=1700, pch=c(16,1), legend=c("hardback  ","softback"),
         horiz=T, bty="n", xjust=0.5, x.intersp=0.75)
  par(xpd=FALSE)
  if(device!="")dev.off()
}

g6.2 <-
function(device="", stats=F){
    if(!require(DAAG))return("Package 'DAAG' is not installed -- cannot proceed")  
    if(device!="")hardcopy(width=5.5, height=1.4, device=device)
    oldpar <- par(mfrow=c(1,4), mar=c(4.1,4.1,2.1,0.6),
                  mgp=c(2.0,.5,0),oma=c(0,0,0,1), pty="s")
    on.exit(par(oldpar))
    allbacks.lm<-lm(weight~-1+volume+area,data=allbacks)
    plot(allbacks.lm,caption=
         c("A: Resids vs Fitted", 
           "B: Normal Q-Q", "C: Scale-Location",
           "D: Resids vs leverage"))
    if(stats)
      print(summary(allbacks.lm))
    if(device!="")dev.off()
    invisible(allbacks.lm)
  }

g6.3 <-
function(device="", dframe=nihills[,c("dist", "climb", "time")]){
    if(device!="")hardcopy(width=5, height=2.9,
                           device=device, pointsize=c(6,4), trellis=TRUE)
    oldpar <- par(fig=c(0,0.5,0,1), mar=c(0.5,0.5,2.6,0.5))
    on.exit(par(oldpar))
    if(!require(DAAG))return("Package 'DAAG' is not installed -- cannot proceed")
    if(!require(lattice))return("Package 'lattice' is not installed -- cannot proceed")
    plot.new()
    mtext(side=3, line=1.5, "A: Untransformed scales:", adj=0, cex=0.7)
    par(fig=c(0.5,1,0,1), mar=c(0.5,0.5,2.6,0.5),new=T)
    mtext(side=3, line=1.5, "B: Logarithmic scales", adj=0, cex=0.7)    
    hills.splom <- splom(~dframe, cex.labels=1.2,
                         axis.line.tck=0.5,
                         varnames=c("dist\n(miles)","climb\n(feet)",
                           "time\n(hours)"), newpage=FALSE, xlab="")    
    print(hills.splom, position=c(0, 0, 0.5, 1), newpage=FALSE)
    lhills.splom <- splom(~log(dframe), cex.labels=1.2,
                          axis.line.tck=0.5,
                          varnames=c("dist\n(log miles)",
                            "climb\n(log feet)", "time\n(log hours)"), xlab="")   
    print(lhills.splom, position=c(0.5, 0, 1, 1), newpage=FALSE)
    if(device!="")dev.off()
    invisible()
  }

g6.4 <-
function(device="", col=1){
    if(device!="")hardcopy(width=5.5, height=1.4, device=device)
    oldpar <- par(mfrow=c(1,4), mar=c(4.1,4.1,2.1,0.6),
                  mgp=c(2.0,.5,0),oma=c(0,0,0,1), pty="s")
    on.exit(par(oldpar))
    nihills.lm <- lm(log(time) ~ log(climb) + log(dist), data = nihills)
    plot(nihills.lm, caption=
         c("A: Resids vs Fitted", 
           "B: Normal Q-Q", "C: Scale-Location", 
           "D: Resids vs leverage"), col=col)
    figtxt<-paste("Diagnostic plots for the regression of log(time)",
                  " on log(climb) and log(dist).", sep="")
    cat(figtxt, "\n")
    if(device!="")dev.off()
    invisible()
  }

g6.5 <-
function(device="")
{
  if(device!="")hardcopy(width=4, height=2, device=device, pointsize=8)
  oldpar <- par(mfrow = c(1, 2), mar = c(3.6, 4.1, 0.6, 2.1),
                mgp = c(2.25, 0.5, 0), pty="s", tcl=-0.35)
  par(mex = 1, cex = 1)
  on.exit(par(oldpar))
  lognihills <- log(nihills)
  names(lognihills) <- paste("log", names(nihills), sep="")
  nihills.lm <- lm(logtime ~ logdist + logclimb, data = lognihills)
  termplot(nihills.lm, partial.resid=TRUE, smooth=panel.smooth, 
           col.res="gray30")  
  if(device!="")dev.off()
  invisible()
}

g6.6 <-
function(device="", smooth=TRUE){
    if(device!="")hardcopy(width=2.75, height=2.75, device=device,
                           pointsize=c(5,3),
                           trellis=TRUE)
    print(splom(~litters, axis.line.tck=0.5,
                varnames=c("lsize\n(litter size)","bodywt\n(Body Weight)",
                  "brainwt\n(Brain Weight)"), axis.text.cex=0.8,
                varname.lineheight=2.5))
    figtxt<-paste("Scatterplot matrix for the litters data set.", sep="")
    print(figtxt)
    if(device!="")dev.off()
    invisible()
  }

g6.7 <-
function(device=""){
    if(!require(lattice))return("Package 'lattice' is not installed -- cannot proceed")  
    if(device!="")hardcopy(width=5.2, height=2.9, trellis=TRUE,
                           device=device, pointsize=c(7,5))
    logoddbooks <- log(oddbooks)
    plt1 <- splom(~log(oddbooks), varnames=c("log(thick)",
                                    "log(breadth)",
                           "log(height)", "log(weight)"),
                  between=list(x=2, y=2),
                   pscales=0, xlab="")
    figtxt<-paste("Scatterplot matrix for the oddbooks data.")
    print(plt1, pos=c(0,0,0.5,1))
    oddothers <- with(oddbooks,
                   data.frame(density = weight/(breadth*height*thick), 
                              area = breadth*height, thick=thick,
                              weight=weight))
    trellis.focus(name="toplevel", highlight=FALSE) 
    panel.text("A", x=0.05, y=0.975) 
    trellis.unfocus()
    plt2 <- splom(~log(oddothers),
                 varnames=c("log(density)", "log(area)", "log(thick)",
                  "log(weight)"),
                  between=list(x=2, y=2),
                 pscales=0, xlab="")
   print(plt2, pos=c(0.5,0,1,1), newpage=FALSE)
   trellis.focus(name="toplevel", highlight=FALSE) 
    panel.text("B", x=0.05, y=0.975) 
    trellis.unfocus()
    if(device!="")dev.off()
    invisible()
  }

g6.8 <-
function(device=""){
    if(device!="")hardcopy(width=4.8, height=1.6, device=device, pointsize=9)
    logoddbooks <- log(oddbooks)
    oldpar <- par(mfrow=c(1,3), mar=c(4.1,4.1,1.6,2.1),
                  mgp=c(2.5,.75,0), oma=c(0,0,0,0), pty="s")
    on.exit(par(oldpar))
    attach(logoddbooks)
    u1 <- lm(weight~I(thick+breadth+height))
    plot(weight~I(thick+breadth+height), pch=16,
         xlab="log(thick)+log(breadth)+log(height)", ylab="log(weight)")
    abline(u1)
    mtext(side=3,line=0.5, "A", adj=0, cex=0.9)
    u2 <- lm(weight~thick+I(breadth+height))
    plot(weight~fitted(u2), xlab="Predicted log(weight)", pch=16,
         ylab="log(weight)")
    abline(0,1)
    mtext(side=3,line=0.5, "B", adj=0, cex=0.9)
    u3 <- lm(weight~thick+breadth+height)
    plot(weight~fitted(u2), xlab="Predicted log(weight)", pch=16,
         ylab="log(weight)")
    abline(0,1)
    mtext(side=3,line=0.5, "C", adj=0, cex=0.9)
    if(device!="")dev.off()
  }

g6.9 <-
function(device="")
{
  if(!require(DAAG))return("Package 'DAAG' is not installed -- cannot proceed")
  if(device!="")hardcopy(width=2.5, height=2.5, device=device,
                         pointsize=7)
  oldpar<-par(mar=c(4.1,4.1,2.1,1.1), mgp=c(2.5,0.75,0), pty="s")
  on.exit(par(oldpar))
  nihills.lm<- lm(log(time) ~ log(dist)+ log(climb), data= nihills) 
  plot(nihills.lm, which=5, add.smooth=FALSE, sub.caption="")
  if(device!="")dev.off()
}

g6.10 <-
function(device="")
{
  if(device!="")hardcopy(width=4.25, height=1.5, device=device, fonts="mono")
  oldpar<-par(mar=c(5.1,3.6,2.1,1.1), mgp=c(2.25,0.75,0), las=1)
  on.exit(par(oldpar))
  allbacks.lm0 <- lm(weight ~ -1+volume+area, data=allbacks)
  z <- dfbetas(allbacks.lm0)
  plot(range(z), c(0.25,1), xlab="Standardized change in coefficient", ylab="",
       axes=F, type="n")
  chh <- par()$cxy[2]
  par(xpd=T)
  axis(1)
  abline(h=par()$usr[3])
  m <- dim(z)[1]
  points(z[,1], rep(0.5,m), pch="|")
  points(z[,2], rep(1,m), pch="|")
  idrows <- (1:m)[apply(z,1,function(x)any(abs(x)>2))]
  if(length(idrows)>0)text(z[idrows,],rep(c(.5,1)-0.25*chh,rep(length(idrows),2)),
                           rep(paste(idrows), 2), pos=1, cex=0.8)
  par(family="mono")
  text(rep(par()$usr[1],2), c(.5,1), c("volume", "area"), pos=2)
  title(sub="dfbetas(allbacks.lm0)")
  par(family="sans")
  if(device!="")dev.off()
  
}

g6.11 <-
function(df=nihills, pred = T, device="")
{
  if(device!="")hardcopy(width=4.25, height=1.85, device=device, pointsize=7)
  obs <- df$time
  hills.logm <- lm(log(time) ~ log(dist) + log(climb), 
                   data = nihills)
  hills.ci<-predict(hills.logm, interval="confidence")
  lognihills<- log(nihills) 
  names(lognihills)<- paste("log",names(nihills),sep="") 
  hills2.lm <- lm(logtime ~ poly(logdist,2) + logclimb, data = lognihills)
  citimes2 <- predict(hills2.lm, interval="confidence") 
  hat <- hills.ci[,"fit"]
  hat2 <- citimes2[,"fit"]
  low <- hills.ci[,"lwr"]-hat
  hi <- hills.ci[,"upr"]-hat
  low2 <- citimes2[,"lwr"]-hat2
  hi2 <- citimes2[,"upr"]-hat2
  ylim <- range(c(low, hi, low2, hi2))
  xlim <- range(c(hat, hat2))
  oldpar <- par(mar = c(4.1,4.1,3.1,1.1), mgp = c(2.5, 0.75, 0), mfrow=c(1,2),
                pty="s", tcl=-0.35)
  on.exit(par(oldpar))
  r <- cor(hat, obs)
  plot(hat, low, type="n", xlab = "Observed time", ylab = "Residual/(Fitted Value)", 
       xlim=xlim, ylim = ylim, xaxt="n", yaxt="n")
  mtext(side=3, line=0.5, adj=0, "A: Logarithmic scales")
  xlab <- pretty(exp(hat))
  axis(1, at=log(xlab), labels=paste(xlab))
  ylab <- pretty(exp(ylim))
  axis(2, at=log(ylab), labels=paste(ylab))
  points(hat, obs-exp(hat), pch=16, cex=0.65)
  points(hat2, obs-exp(hat2), pch=16, col="gray", cex=0.65)
  ord <- order(hat)
  lines(hat[ord], low[ord], col = "black")
  lines(hat[ord], hi[ord], col = "black")
  ord <- order(hat2)
  lines(hat2[ord], low2[ord], col = "grey60", lwd=1.5)
  lines(hat2[ord], hi2[ord], col = "grey60", lwd=1.5)
  ## Unlogged
  ulow <- exp(low+hat)-exp(hat)
  ulow2 <- exp(low2+hat2)-exp(hat2)
  uhi <- exp(hi+hat)-exp(hat)
  uhi2 <- exp(hi2+hat2)-exp(hat2)
  ylim <- range(c(ulow,ulow2,uhi, uhi2, obs-exp(hat)))
  xlim <- range(c(exp(hat),exp(hat2)))
  plot(exp(hat), ulow, type="n", xlab = "Observed time", ylab = "Residual", 
       xlim=xlim, ylim = ylim)
  mtext(side=3, line=0.5, adj=0, "B: Untransformed scales (time)")
  points(exp(hat), obs-exp(hat), pch=16, cex=0.65)
  points(exp(hat2), obs-exp(hat2), pch=16, col="gray", cex=0.65)
  ord <- order(hat)
  lines(exp(hat[ord]), ulow[ord], col = "black")
  lines(exp(hat[ord]), uhi[ord], col = "black")
  ord <- order(hat2)
  lines(exp(hat2[ord]), ulow2[ord], col = "grey60", lwd=1.5)
  lines(exp(hat2[ord]), uhi2[ord], col = "grey60", lwd=1.5)
  topleft <- par()$usr[c(1, 4)]
  chw <- par()$cxy[1]
  chh <- par()$cxy[2]
  if(device!="")dev.off()
  invisible(hills.logm)
}

g6.12 <-
function(device="", dframe=hills2000[,c("dist", "climb", "time")]){
    if(device!="")hardcopy(width=2.75, height=2.75,
                           device=device, pointsize=c(7,4), trellis=TRUE,
                           color=FALSE)
    if(!require(lattice))return("Package 'lattice' is not installed -- cannot proceed")   
    set1 <- simpleTheme(alpha=0.8)
    lhills.splom <- splom(~log(dframe), cex.labels=1.2, axis.line.tck=0.5,
                          axis.text.cex=0.75,
                          par.settings=set1, varnames=c("dist\n(log miles)",
                                               "climb\n(log feet)", "time\n(log hours)"), xlab="")   
    print(lhills.splom)
    if(device!="")dev.off()
    invisible()
  }

g6.13 <-
function(device="", dframe=hills2000[, c("dist", "climb", "time")])
{
  if(device!="")hardcopy(width=4, height=2.15, device=device)
  oldpar <- par(mfrow=c(1,2), mar=c(4.1,4.1,2.1,0.6),
                mgp=c(2.0,.5,0),oma=c(0,0,0,1), pty="s", tcl=-0.35)
  on.exit(par(oldpar))
  lhills2k.lm <- lm(log(time) ~ log(climb) + log(dist), data = dframe)
  plot(lhills2k.lm, caption="", which=1)
  mtext(side=3, line=0.5, "A: lm fit", adj=0)
  if(!require(MASS))return("Package 'MASS' is not installed -- cannot proceed")    
  lhills2k.lqs <- lqs(log(time) ~ log(climb) + log(dist), data = dframe)
  reres <- residuals(lhills2k.lqs)
  refit <- fitted(lhills2k.lqs)
  big3 <- which(abs(reres) >= sort(abs(reres), decreasing=TRUE)[3])
  plot(reres ~ refit, xlab="Fitted values (resistant fit)",
       ylab="Residuals (resistant fit)")
  lines(lowess(reres ~ refit), col=2)
  mtext(side=3, line=0.5, "B: Resistant (lqs) fit", adj=0)
  text(reres[big3] ~ refit[big3], labels=rownames(dframe)[big3],
       pos=4-2*(refit[big3] > mean(refit)), cex=0.8)
  if(device!="")dev.off()
  invisible()
}

g6.14 <-
function(df = hills2000[, c("dist", "climb", "time")], device="")
{
  if(device!="")hardcopy(width=4, height=2.15, device=device, pointsize=8)
  oldpar <- par(mfrow = c(1, 2), mar = c(4.6, 4.1, 2.1, 2.1),
                mgp = c(2.5, 0.75, 0), pty="s", tcl=-0.35)
  if(!require(DAAG))return("Package 'DAAG' is not installed -- cannot proceed")  
  lhills2k <- log(df)
  names(lhills2k) <- paste("log", names(df), sep="")
  par(mex = 1, cex = 1)
  on.exit(par(oldpar))
  u <- lm(logtime ~ logdist+logclimb, data=lhills2k)
  termplot(u, partial.resid=TRUE, smooth=panel.smooth, col.res="gray30")
  if(device!="")dev.off()
  invisible()
}

g6.15 <-
function(device=""){
    if(device!="")hardcopy(width=5.5, height=1.4, device=device)
    oldpar <- par(mfrow=c(1,4), mar=c(4.1,4.1,2.1,0.6),
                  mgp=c(2.0,.5,0),oma=c(0,0,0,1), pty="s")
    on.exit(par(oldpar))
    lhills2k <- log(hills2000[, c("dist", "climb", "time")])
    names(lhills2k) <- paste("log", names(lhills2k), sep="")    
    lhills.lm2 <- lm(logtime ~ poly(logdist,2)+logclimb, data=lhills2k[-42, ])
    plot(lhills.lm2, caption=
         c("A: Resids vs Fitted", 
           "B: Normal Q-Q", "C: Scale-Location", 
           "D: Cook's distance"))
    if(device!="")dev.off()
    invisible()
  }

g6.16 <-
function(device=""){
    if(device!="")hardcopy(width=4.25, height=4.25,
                           device=device, pointsize=7)                
    if(!exists("coxite")){
      if(!require(compositions))return("Package 'compositions' is not installed -- cannot proceed")      
      data(Coxite)         # compositions pkg does not have lazy loading
      coxite <- as.data.frame(Coxite)
    }
    pairs(coxite, oma=rep(2,4))
    if(device!="")dev.off()
  }

g6.17 <-
function(mu=20, n=100,a=15, b=2.5, sigma=12.5, device=""){
    if(device!="")hardcopy(width=2.25, height=2.25,
                           device=device, pointsize=7)
    if(!require(compositions))return("Package 'compositions' is not installed -- cannot proceed")     
    oldpar <- par(mar=c(4.1,4.6, 1.1, 1.6), mgp=c(2.5,0.5,0),
                  pty="s")
    on.exit(par(oldpar))
    data(Coxite)         # compositions pkg does not have lazy loading
    coxite <- as.data.frame(Coxite)     # Convert matrix to data frame
    coxiteAll.lm <- lm(porosity ~ A+B+C+D+E+depth, data=coxite)
    coxite.hat <- predict(coxiteAll.lm, interval="confidence")
    hat <- coxite.hat[,"fit"]
    plot(porosity ~ hat, data=coxite,
         xlab="Fitted values",
         ylab="Fitted values, with 95% CIs\n(Points are observed porosities)",
         type="n", pty="s", tcl=-0.35)
    with(coxite, points(porosity ~ hat, col="gray45"))
    lines(hat, hat)
    ord <- order(hat)
    sebar <- function(x, y1, y2, eps=0.15){
      lines(rep(x,2), c(y1,y2))
      lines(c(x-eps,x+eps), rep(y1,2))
      lines(c(x-eps,x+eps), rep(y2,2))
    }
    q <- ord[round(quantile(1:length(hat), (1:9)/10))]
    for(i in q)sebar(hat[i], coxite.hat[i,"lwr"], coxite.hat[i,"upr"])
    if(device!="") dev.off()
  }

g6.18 <-
function(device=""){
    if(device!="")hardcopy(width=4.5, height=3.5, trellis=TRUE,
                           device=device, pointsize=c(8,4))
    if(device=="ps")parset <- simpleTheme(col=c("gray40","black"),
         col.line=c("gray40","black")) else
    parset <- simpleTheme(alpha=0.75, col=c("gray40","black"), 
                          col.line=c("gray40","black"), lty=1:2)    
    errorsINx(gpdiff=0, parset=parset, plotit=TRUE)
    if(device!="") dev.off()
  }

g6.19 <-
function(device=""){
    if(device!="")hardcopy(width=5.25, height=2.25, trellis=TRUE,
                           device=device, pointsize=c(7,4))
    if(device=="ps")parset <- simpleTheme(col=c("gray40","black"),
         col.line=c("gray40","black")) else
    parset <- simpleTheme(alpha=0.75, col=c("gray40","black"), pch=c(0,1),
                          col.line=c("gray40","black"), lty=1:2)
    errorsINx(gpdiff=1.5, timesSDx=(1:2)*0.8, layout=c(3,1), parset=parset)
    if(device!="") dev.off()
  }

g6.20 <-
function(min=30, device=""){
    if(device!="")hardcopy(width=4.5, height=1.75, trellis=TRUE, color=FALSE,
                           device=device, pointsize=c(8,5))
    ## gabalong <- stack(gaba[c(2,5,8),], select=-min)
    ## gabalong$min <- rep(c(30,90,150), 6)
    gabalong <- stack(gaba[paste(min),], select=-min)
    gabalong$sex <- factor(rep(c("male", "female","all"), rep(2,3)),
                           levels=c("female","male","all"))
    gabalong$treatment <- factor(rep(c("Baclofen","No baclofen"), 3),
                                 levels=c("No baclofen","Baclofen"))
    print(
          stripplot(sex~values, groups=treatment, data=gabalong,
                    par.settings=list(superpose.symbol=list(cex=c(1.35,1.35),
                                        pch=c(1,16))),
                    xlab=list("Average reduction: 30 min vs 0 min",
                      cex=1.0),
                    scales=list(cex=1.0, tck=0.5),
                    panel=function(x,y,...){
                      panel.stripplot(x,y,...)
                      ltext(x,y,paste(c(3,9,15,7,22,12)), pos=1, cex=0.8)
                    }, auto.key=list(space="right", points=TRUE, cex=1.0))
          )
    ## plot(rep(c(30, 90, 150),6), unlist(gaba[c(2,5,8), 2:7]),
    ##     pch=rep(c(16,15,1,0,16,15),rep(3,6)),
    ##     col=rep(c("gray","black"), c(12,6)),
    ##     xlab="Reduction in VAS pain rating")
    ## text(rep(c(30, 90, 150),6), unlist(gaba[c(2,5,8), 2:7]),
    ##     labels=rep(c(paste(c(3,15,9,7)),12,22), rep(3,6)), pos=4)
    if(device!="")dev.off()
    gabalong
  }

g6.21 <-
function(min=30, device=""){
    if(!require(MCMCpack))return("The package MCMCpack must be installed.")
    if(device!="")hardcopy(width=5, height=5, color=FALSE,
                           device=device, pointsize=c(8,5))
    opar <- par(mgp=c(2,.75,0), mar=c(4.1,3.6,2.6,1.1), oma=rep(0,4),
                no.readonly=TRUE)   
  mat <- matrix(c(1:8), byrow=TRUE, ncol=2)
  layout(mat, widths=rep(c(2,1),4), heights=rep(1,8))    
posterior <- MCMCregress(log(time) ~ log(climb/dist) + log(dist),
                          data=nihills)
    plot(posterior, auto.layout=FALSE, ask=FALSE, col="gray40")
     if (device != "") dev.off() 
  }

gdump <-
function(fnam=NULL, prefix="~/dbeta/figures/figs",
           xtras=c("hardcopy","renum.fun","renum.files"), 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)), xtras)
    cat("\nDump to file:", fnam, "\n")
    print(objnames)
    dump(objnames, fnam)
  }

gsave <-
function(fnam=NULL, prefix="~/dbeta/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)
  }

renum.fun <-
function(from.prefix="g", to.prefix="g",from=20:7, to=21:8, doit=F){
          path <- getwd()
    pathtag <- strsplit(path, "/ch", fixed=TRUE)[[1]]
    endbit <- pathtag[length(pathtag)]
    from.prefix <- paste(from.prefix, endbit, sep="")
    to.prefix <- paste(to.prefix, endbit, sep="")
    for(i in 1:length(to))
      {txt<-paste(to.prefix,".",to[i]," <- ", from.prefix,".",from[i],sep="")
       if(doit)eval(parse(text=txt),envir=sys.frame(0))
       print(txt)
     }
  }

pkgs <- c("DAAG", "lattice", "MASS", "MCMCpack", "compositions")
z <- sapply(pkgs, require, character.only=TRUE, warn.conflicts=FALSE, quietly=TRUE)

## ##
## ## Markov Chain Monte Carlo Package (MCMCpack)
## ## Copyright (C) 2003-2014 Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park
## ##
## ## Support provided by the U.S. National Science Foundation
## ## (Grants SES-0350646 and SES-0350613)
## ##
## 
## Attaching package: 'tensorA'
## 
## The following object is masked from 'package:base':
## 
##     norm
## 
## 
## Attaching package: 'robustbase'
## 
## The following object is masked from 'package:DAAG':
## 
##     milk
## 
## 
## Attaching package: 'bayesm'
## 
## The following object is masked from 'package:MCMCpack':
## 
##     rdirichlet
## 
## Welcome to compositions, a package for compositional data analysis.
## Find an intro with "? compositions"

if(any(!z)){
  notAvail <- paste(names(z)[!z], collapse=", ")
  warning(paste("The following packages should be installed:", notAvail))
}

g6.1()

g6.2()

g6.3()

g6.4()

## Diagnostic plots for the regression of log(time) on log(climb) and log(dist).

g6.5()

g6.6()

## [1] "Scatterplot matrix for the litters data set."

g6.7()

g6.8()

g6.9()

g6.10()

g6.11()

g6.12()

g6.13()

g6.14()

g6.15()

g6.16()

g6.17()

g6.18()

##               Intercept  Slope
## No error in x     14.96 1.4858
## 0.4sx             17.76 1.2648
## 0.8sx             22.46 0.8821
## 1.2sx             25.60 0.6565
## 1.6sx             28.28 0.4130
## 2sx               30.88 0.2220

g6.19()

##               Intercept:ctl Offset:trt  Slope
## No error in x         15.20    0.09713 1.4931
## 0.8sx                 22.06    1.08829 0.9519
## 1.6sx                 28.55    1.45975 0.4451

g6.20()

##   values    ind    sex   treatment
## 1  1.311   mbac   male    Baclofen
## 2  1.647    mpl   male No baclofen
## 3  3.479   fbac female    Baclofen
## 4  4.151    fpl female No baclofen
## 5  3.118  avbac    all    Baclofen
## 6  2.743 avplac    all No baclofen

g6.21()

renum.files <-
function(from.prefix="~/dbeta/Art/", to.prefix="~/dbeta/Art/",
           from=20:7, to=21:8, doit=F){
    path <- getwd()
    pathtag <- strsplit(path, "/ch", fixed=TRUE)[[1]]
    endbit <- pathtag[length(pathtag)]
    if(nchar(endbit)==2)chap <- paste(endbit) else
    chap <- paste("0",endbit,sep="")
    from.prefix <- paste(from.prefix, chap, "-", sep="")
    to.prefix <- paste(to.prefix, chap, "-", sep="")
    for(i in 1:length(from)){
      if (from[i]<=9) ltext <- paste("0",from[i],sep="")
      else ltext <- paste(from[i])
      if (to[i]<=9) rtext <- paste("0",to[i],sep="")
      else rtext <- paste(to[i])            
      txt<-paste("mv ", from.prefix, ltext, ".pdf", " ",
                 to.prefix, rtext, ".pdf", sep="")
      backup<-paste("cp ", from.prefix, ltext, ".pdf", " ",
                    "archive", sep="")
      if(doit)system(backup)
      if(doit)system(txt)
      print(backup)
      print(txt)
    }
  }