# Redefine qqplot so it adds, instead of setting up a new plot:
qqplot <- 
function (x, y, plot.it = TRUE, xlab = deparse(substitute(x)), 
    ylab = deparse(substitute(y)), add.to.plot = FALSE, ...) 
{
    sx <- sort(x)
    sy <- sort(y)
    lenx <- length(sx)
    leny <- length(sy)
    if (leny < lenx) 
        sx <- approx(1L:lenx, sx, n = leny)$y
    if (leny > lenx) 
        sy <- approx(1L:leny, sy, n = lenx)$y
    if (plot.it)
      if(add.to.plot)
        points(sx, sy, ...)
      else
        plot(sx, sy, xlab = xlab, ylab = ylab, ...)
    invisible(list(x = sx, y = sy))
}
# Two-sample comparison with empirical qq-plot, with permutation null band:
qqplot.null <- function(x, y, xlab="x", ylab="y",
                        jitter=0, nrep=200, rg=range(c(x,y), na.rm=T), 
                        pch=16, cex=0.8, lwd=3,
                        mgp=c(2,0.5,0), mar=c(3,3,1,1)+0.1, ...) {
  par(mgp=mgp, mar=mar)
  nx <- length(x);  ny <- length(y);  xy <- c(x,y)
  if(jitter>0) xy <- xy + runif(nx+ny, -jitter, jitter)
  qqplot(rg, rg, xlab=xlab, ylab=ylab, xlim=rg, ylim=rg, type="n", ...)
  if(nrep>0) for(irep in 1:nrep) {
    xynull <- sample(xy);  xnull <- xynull[1:nx]
    ynull <- xynull[(nx+1):(nx+ny)]
    qqplot(xnull, ynull, type="l", col="gray50", add=T, ...)
  }
  lines(rg,rg)
  xx <- xy[1:nx];  yy <- xy[(nx+1):(nx+ny)]
  qqplot(xx, yy, type="l", lwd=lwd, add=T)
  qqplot(xx, yy, type="p", pch=pch, cex=cex, add=T)
}

box.cox <- function(y, lambda=1) {
  if(lambda==0) log(y) else (y^lambda-1)/lambda }

#
qqnorm.plus <- 
qqnorm.null.conf <-
  function(dat, box.cox.par=1, nsim=100,
           ylab="Standardized Empirical Quantiles",
           new.window=T, null=T, conf=T) {
    if(new.window) if(null & conf) { windows(12, 6.5); par(mfrow=c(1,2)) } else windows(6.5,6.5)
    par(mgp=c(2,.8,0), mar=c(4,4,4,2))
    y <- sort(box.cox(dat, box.cox.par))
    y <- (y - mean(y))/sd(y)                   # standardize
    n <- length(y)
    x <- qnorm((1:n)/(n+1))                      # normal quantiles
    r <- range(c(x,y))
    # Null band:
    if(null) {
      plot(x=r, y=r, pch=16, cex=.4, type="n",
           xlab="Normal Quantiles", ylab=ylab)
      for(i in 1:nsim) {
        z <- sort(rnorm(length(x)))
        z <- (z - mean(z))/sd(z)                 # standardize
        points(x=x, y=z, pch=16, cex=.3, col="gray60")  # overplot
      }
      lines(x=x, y=y, lwd=2) #
      lines(range(c(x,y)), range(c(x,y)))
      title("Normal Null Band")
    }
    # Bootstrap confidence band:
    if(conf) {
      plot(x=r, y=r, pch=16, cex=.4, type="n",
           xlab="Normal Quantiles", ylab=ylab)
      for(i in 1:nsim) {
        z <- sort(sample(y, size=length(y), replace=T))
        z <- (z - mean(z))/sd(z)                 # standardize
        points(x=x, y=z, pch=16, cex=.3, col="gray60")  # overplot
      }
      lines(x=x, y=y, lwd=2) #
      lines(range(c(x,y)), range(c(x,y)))
      title("Bootstrap Confidence Band")
    }
  } # end of:  qqnorm.null.conf <-