################################################################ # 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...", ...) }