g13.1 <-
function(device=""){
    oldpal <- palette(gray(c(0,.6)))
    if(device!="") hardcopy(width=2.75, height=2.75,
                            device=device)
    oldpar <- par(mar = par()$mar + c(-0.5, 0, -2, 0))
    on.exit(par(oldpar))
    require(DAAG); data(socsupport)
    not.na <- apply(socsupport[,9:19],1,function(x)!any(is.na(x)))
    not.na[36] <- FALSE
    ss.pcp1 <- princomp(socsupport[not.na, 9:19], cor=TRUE)$scores[,1]
    attach(socsupport)
    plot(BDI[not.na] ~ ss.pcp1, col=as.numeric(gender[not.na]),
         pch=as.numeric(gender[not.na]),
         xlab ="First principal component", ylab="BDI")
    topleft <- par()$usr[c(1,4)]
    legend(topleft[1], topleft[2], col=1:2, pch=1:2,
           legend=levels(gender), cex=0.7)
    detach(socsupport)
    if(device!="")dev.off()
    palette(oldpal)
    invisible()
  }

g13.2 <-
function(df=nsw74psid1, device="", lty=c(1,1),
           colr=c("black","gray"), lwd=1:2){
    if(device!="")hardcopy(width=5.5, height=1.6,
                           device=device, pointsize=9)
    oldpar <- par(mfrow=c(1,5), mar=c(4.1,2.6,2.1,0.1), oma=c(0,2.1,2.1,1.1),
                  mgp=c(2.25,.5,0), pty="s")
    on.exit(par(oldpar))
    attach(df)
    k <- 0
    for(i in c("educ","age","re74","re75","re78")){
      k <- k+1
      y <- df[,i]
      d0 <- density(y[trt==0])
      d1 <- density(y[trt==1])
      ylim <- range(c(d0$y, d1$y))
##      if(k %in% c(1,4))
##        plot(d0$x, d0$y, ylim=ylim, main="", xlab=i, type="l", ylab="Density",
##             lwd=lwd[1], col=colr[1])
##      else
        plot(d0$x, d0$y, ylim=ylim, main="", xlab=i, type="l", ylab="",
             lwd=lwd[1], col=colr[1], cex.lab=1.5)
      lines(d1$x, d1$y, lty=lty[2], lwd=lwd[2], col=colr[2])
    }
#    plot(0:1, 0:1, type="n", axes=FALSE, xlab="", ylab="")
    mtext(side=2, line=0.1, "Density", outer=TRUE)
    par(fig=c(0,1,0,1), mar=rep(0,4), oma=rep(0,4), new=TRUE)
    xpos <- mean(par()$usr[1:2])
    ypos <- par()$usr[4]
    par(xpd=TRUE)
    legend(x=xpos, y=ypos, lty=lty, col=colr, lwd=lwd, bty="n",
           legend=c("control", "treated"), horiz=T, inset=-3,
           xjust=0.75, yjust=0)
    if(device!="")dev.off()
    detach(df)
  }

g13.3 <-
function(device="", df=nsw74psid1, maxf=30, bw=TRUE)
{
  library(splines)
  if(device!="")hardcopy(width=2.25, height=2.25, device=device,
                         trellis=TRUE, pointsize=c(8,5))
  oldpar <- par(mar = c(4.1, 4.1, 1.1, 1.1),
                mgp = c(2.25, 0.5, 0) )
  on.exit(par(oldpar))
  common <- multilap(maxf=maxf)
  offset <- mean(sapply(nsw74psid1[,8:10],
                        function(x)unique(sort(x))[2]))
  if(bw)trellis.par.set(superpose.symbol=list(col="gray", pch=c(1,3)),
                        superpose.line=list(col=c("black", "gray"),
                          lwd=c(1,2)), lty=c(1,1))
  A.lm <- lm(log(re78+offset) ~ trt + 
             (ns(age,4) + educ + log(re74+offset) +
              log(re75+offset)) + I(re74>0) + I(re75>0)+
             (black + hisp + marr + nodeg), subset=common,
             data = df)
  print(summary(A.lm)$coef)
  trt <- factor(df$trt[common], labels=c("ctl","trt"))
  print(xyplot(resid(A.lm) ~ fitted(A.lm), groups=trt,
               auto.key=list(columns=2, lines=TRUE),
               type=c("p","smooth")))
  if(device!="")dev.off()
  invisible()
}

g13.4 <-
function(df=nsw74psid1, device="", maxf=30, maxf.score=30, bw=TRUE){
    if(device!="")hardcopy(width=4.25, height=2.35, trellis=TRUE, 
                           device=device, pointsize=c(8,5), color=!bw)
    oldpar <- par(mar=c(2.1,4.1,2.6,2.1), cex=0.7,
                  mgp=c(2.25,.5,0), pty="s")
    on.exit(par(oldpar))
    common <- multilap(maxf=maxf)
    offset <- mean(sapply(nsw74psid1[,8:10],
                          function(x)unique(sort(x))[2]))
    disc.glm <- glm(trt ~ (ns(age,4) + educ + log(re74+offset) +
                           log(re75+offset)) + I(re74>0) + I(re75>0)+
                    (black + hisp + marr + nodeg), subset=common,
                    data = df, family=binomial)    
    Pscores <- predict(disc.glm)
    nsw74psidB <- nsw74psid1[common, ]
    trt <-  nsw74psidB$trt
    par(fig=c(0, 0.5, 0, 1))
    if(bw)colr <- c("black","gray") else
    colr <- trellis.par.get()$superpose.symbol$col
    par(xpd=TRUE)
    xchop <- overlapDensity(Pscores[trt==0], Pscores[trt==1], lty=c(1,1),
                            compare.numbers=FALSE, col=colr,
                            ratio=c(1/maxf.score, maxf.score))
    mtext(side=3, line=2.85, adj=-0.2, "A", cex=0.8)
    par(xpd=FALSE)
    par(fig=c(0.5, 1, 0, 1))
    par(new=TRUE)
    plot(0:1,0:1, xlab="", ylab="",axes=FALSE, type="n")
    mtext(side=3, line=2.85, adj=-0.35, "B", cex=0.8)
    overlap <- Pscores > xchop[1] & Pscores < xchop[2]
    print(with(subset(nsw74psidB, overlap), table(re78>0, trt)))
    AS.glm <- glm(I(re78>0) ~ trt + Pscores, data=nsw74psidB, 
                  subset =  overlap, family=binomial)    
    AS.lm <- lm(log(re78+47) ~ trt + Pscores, data=nsw74psidB, 
                subset =  overlap&nsw74psidB$re78>0)
    print(summary(AS.glm)$coef)
    print(summary(AS.lm)$coef)
    print(sum(overlap))
    trt <- factor(trt, labels=c("ctl", "trt"))
    linecols <- c(rgb(0,0.5,0.5), rgb(0.65,0,0.65))

    if(bw) trellis.par.set(superpose.symbol=list(col=c("gray40", "gray40"),
                             pch=c(1,3)),
                           superpose.line=list(col=c("black", "gray"),
                             lty=c(1,1), lwd=c(1,2)),
                           layout.heights=list(axis.top=0.6))
    else
      trellis.par.set(superpose.symbol=list(pch=c(1,3)),
                      superpose.line=list(lty=c(1,1), col=linecols,
                        lwd=c(1.5,1.5)))
    print(trellis.par.get()$superpose.line$lty)
    gph <- xyplot(resid(AS.lm) ~ fitted(AS.lm),
                  groups=trt[overlap&nsw74psidB$re78>0], 
                  type=c("p","smooth"),
                  auto.key=list(columns=2, lines=TRUE), aspect=1)
    print(gph, position=c(0.5, 0, 1, 1), newpage=FALSE)
    if(device!="")dev.off()
  }

gdump <-
function(fnam=NULL, prefix="~/r-book/ed2/figures/figs",
           xtras=c("hardcopy","renum.obj","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="~/r-book/ed2/figures/figs",
           splitchar="/ch", xtras=c("hardcopy","renum.obj","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)
  }

hardcopy <-
function(width=3.75, height=3.75, color=FALSE, trellis=FALSE,
           device=c("","pdf","ps"), path="~/r-book/ed2/Art/", 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])))
  }

renum.obj <-
function(from.prefix="g", to.prefix="g",from=2:5, to=1:4,
           prenum=paste(rev(strsplit(path, "/ch", fixed=TRUE)[[1]])[1], ".",
             sep=""), details=FALSE, doit=FALSE, archive=TRUE){
    if(length(from) != length(to))
      stop(paste("to vector has length ", to, "; from vector has length ",
                 from,sep=""))
    path <- getwd()
#    pathtag <- strsplit(path, "/ch", fixed=TRUE)[[1]]
#    endbit <- pathtag[length(pathtag)]
    from.prefix <- paste(from.prefix, prenum, sep="")
    to.prefix <- paste(to.prefix, prenum, sep="")
    namsfrom <- paste(from.prefix, from, sep="")
    namsto <- paste(to.prefix, to, sep="")
    tmpfrom <- paste("tmp", from, sep="")
    from2to <- paste(namsto, " <- ", namsfrom, sep="")
    cat(paste(from2to, collapse="\n"))
    cat("\n")
    if(doit)
      for(i in 1:length(to)){
         if(!exists(namsfrom[i]))
           stop(paste("From file", namsfrom[i], "was not found"))
        if(!exists(namsto[i]))
           stop(paste("To file", namsto[i], "was not found"))
         from2tmp <- paste(tmpfrom[i], " <- ", namsfrom[i], sep="")
         if(details)print(from2tmp)
         eval(parse(text=from2tmp), envir=sys.frame(0))
       }
      if(doit)
        for(i in 1:length(to)){
          tmp2to <- paste(namsto[i], " <- ", tmpfrom[i], sep="")
          rmtxt <- paste("rm(", tmpfrom[i], ")", sep="")
          if(details){
            print(tmp2to)
            print(rmtxt)
          }
          eval(parse(text=tmp2to), envir=sys.frame(0))
          eval(parse(text=rmtxt), envir=sys.frame(0))
        }
    if(archive & doit){
      save(list=namsfrom, file="objfrom.RData")
      print("From objects have been archived to 'objfrom.RData'", quote=FALSE)
    }
  }

renum.files <-
function(from.prefix="~/r-book/ed2/Art/", to.prefix="~/r-book/ed2/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, ".eps", " ",
                 to.prefix, rtext, ".eps", sep="")
      backup<-paste("cp ", from.prefix, ltext, ".eps", " ",
                    "archive", sep="")
      if(doit)system(backup)
      if(doit)system(txt)
      print(backup)
      print(txt)
    }
  }