# 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(1:lenx, sx, n = leny)$y
    if (leny > lenx) 
        sy <- approx(1: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))
}
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)
}