################################################################
#
    null.residuals <- function(model.lm) {  
      if(exists("y.null", envir=parent.frame())) gaga <- y.null else gaga <- NULL
      y.null <<- rnorm(length(residuals(model.lm)))
      null.lm <- update(model.lm, y.null ~ .)
      if(is.null(gaga)) rm(y.null, envir=parent.frame()) else y.null <<- gaga
      return(residuals(null.lm) * sqrt(deviance(model.lm) / deviance(null.lm)))
    }
# The relevant lines above begin with "y.null <<- ..." and "null.lm <- ...",
# whereas the two lines "if(..." are silly complications of R semantics
# that I don't want to explain.

    res.vs.fits <- function(model.lm, ident=F, new=T, pch=16, cex=.5, ...) {
      if(new) windows(width=9, height=6)
      par(mfcol=c(2,3), mgp=c(1.5,0.5,0), mar=c(2.5,2.5,1,1))
      fit <- fitted.values(model.lm)
      res <- residuals(model.lm)
      rg <- c(-1,1)*max(abs(res))
      for(i in 1:5) {
        null.res <<- null.residuals(model.lm)
        plot(x=fit, y=null.res,
             xlab="Fitted", ylab="Null Resid", ylim=rg,
             pch=pch, cex=cex, ...)
        abline(h=0, col="gray", lwd=.5)
        lines(lowess(null.res ~ fit), col="red", lwd=.5)
      }
      plot(x=fit, y=res, xlab="Fitted", ylab="Actual Resid", ylim=rg,
           pch=pch, cex=cex, ...)
      abline(h=0, col="gray", lwd=.5)
      lines(lowess(res ~ fit), col="red", lwd=.5)
      if(ident) {
        cat("Identify odd points in the last plot; abort when done...\n")
        identify(x=fit, y=res, labels=names(res))
      }
    }
    res.vs.qnorm <- function(model.lm, ident=F, new=T, pch=16, cex=.5, ...) {
      if(new) windows(width=9, height=6)
      par(mfcol=c(2,3), mgp=c(1.5,0.5,0), mar=c(2.5,2.5,1,1))
      res.sort <- sort(residuals(model.lm))
      n <- length(res.sort)
      quant <- qnorm(1:n/(n+1))
      rg <- c(-1,1)*max(abs(res.sort))
      sd <- sqrt(deviance(model.lm) / df.residual(model.lm));  
      for(i in 1:5) {
        plot(quant, sort(null.residuals(model.lm)),
             xlab="Normal Quantiles", ylab="Null Resid", ylim=rg,
             pch=pch, cex=cex, ...)
        abline(a=0, b=sd, col="gray", lwd=.5)
      }
      plot(quant, res.sort, xlab="Norm. Quantiles", ylab="Actual Resid", ylim=rg,
           pch=pch, cex=cex, ...)
      abline(a=0, b=sd, col="gray", lwd=.5)
      if(ident) {
        cat("Identify odd points in the last plot; abort when done...\n")
        identify(x=quant, y=res.sort, labels=names(res.sort))
      }
    }
    res.vs.order <- function(model.lm, ident=F, new=T, pch=16, cex=.5, col="red", lwd=2, f=0.1) {
      if(new) windows(width=9, height=6)
      par(mfcol=c(2,3), mgp=c(1.5,0.5,0), mar=c(2.5,2.5,1,1))
      r.obs <- residuals(model.lm)
      ord <- 1:length(r.obs)
      rg <- c(-1,1)*max(abs(r.obs))
      for(i in 1:5) {
        r.null <- null.residuals(model.lm)
        plot(ord, r.null, xlab="Order", ylab="Null Resid", ylim=rg, type="n")
        lines(lowess(x=ord, y=r.null, f), col=col, lwd=lwd)
        points(ord, r.null, pch=pch, cex=cex)
      }
      plot(ord, r.obs, xlab="Order", ylab="Actual Res", ylim=rg, type="n")
      lines(lowess(x=ord, y=r.obs, f), col=col, lwd=lwd)
      points(ord, r.obs, pch=pch, cex=cex)
      if(ident) {
        cat("Identify odd points in the last plot; abort when done...\n")
        identify(x=ord, y=r.obs, labels=names(r.obs))
      }
    }
    res.vs.lag <- function(model.lm, ident=F, new=T, pch=16, cex=.5, col="red", lwd=2) {
      if(new) windows(width=9, height=6)
      par(mfcol=c(2,3), mgp=c(1.5,0.5,0), mar=c(2.5,2.5,1,1))
      r.obs <- residuals(model.lm)
      rg <- c(-1,1)*max(abs(r.obs))
      N <- length(r.obs)
      for(i in 1:5) {
        r.null <- null.residuals(model.lm)
        plot(r.null[-N], r.null[-1], xlab="Lag", ylab="Null Resid", xlim=rg, ylim=rg, type="n")
        abline(lm(r.null[-1] ~ r.null[-N]), col=col, lwd=lwd)
        points(r.null[-N], r.null[-1], pch=pch, cex=cex)
      }
      plot(r.obs[-N], r.obs[-1], xlab="Lag", ylab="Null Resid", xlim=rg, ylim=rg, type="n")
      abline(lm(r.obs[-1] ~ r.obs[-N]), col=col, lwd=lwd)
      points(r.obs[-N], r.obs[-1], pch=pch, cex=cex)
      if(ident) {
        cat("Identify odd points in the last plot; abort when done...\n")
        identify(x=r.obs[-N], y=r.obs[-1], labels=names(r.obs))
      }
    }

    stdize <- function(x) (x-min(x))/(max(x)-min(x))
    res.vs.predictors <- function(model.lm,
                                  preds=names(model.lm$model)[-1],
                                  new=T, pch=16, cex=.5, f=0.5, ...) {
      r.obs <- stdize(residuals(model.lm))*.9+1.1
      r.null <- stdize(null.residuals(model.lm))*.9
      if(new) windows(width=length(preds)*3, height=2*3)
      par(mfcol=c(1,length(preds)), mar=c(.5,.5,.5,.5))
      for(pred in preds) {
        x <- model.lm$model[,pred];  x.range <- range(x)
        plot(x.range, c(0,2), type="n", xaxt="n", yaxt="n", xlab="", ylab="", ...)
        text(mean(x.range), 1, pred)
        points(x, r.obs, pch=pch, cex=cex, xlab=pred, ylab="Actual Res")
        lines(lowess(x,r.obs))
        points(x, r.null, pch=pch, cex=cex, xlab=pred, ylab="Null Res")
        lines(lowess(x,r.null))
      }
    }  
    partial.res.plot <- function(x, y, new=T, cex=0.5, pch=16, ...) {
      p <- ncol(x);  x <- as.matrix(x)   # Why can't lm take data frames here?
      if(new) windows(width=ceiling(p/4)*3, height=4*3)
      par(mfcol=c(4,ceiling(p/4)), mgp=c(1.5,0.5,0), mar=c(2.5,2.5,0.5,0.5))
      for(i in 1:p)
        plot(x=resid(lm(x[,i] ~ x[,-i])),
             y=resid(lm(y ~ x[,-i])),
             cex=cex, pch=pch, xlab=colnames(x)[i], ylab="y...", ...)
    }