Possible extensions to the default diagnostic plots for lm (in R and in general)?

Question

I started digging a bit into the plot.lm function, this function gives six plots for lm, they are:

1. a plot of residuals against fitted values
2. a Scale-Location plot of sqrt(| residuals |) against fitted values
3. a Normal Q-Q plot, a plot of Cook's distances versus row labels
4. a plot of residuals against leverages
5. a plot of Cook's distances against leverage/(1-leverage)

And I am wondering what other common/useful extensions of current plots exists for linear models, and how can they be done in R? (links to articles of packages are also welcomed)

So the boxcox function (from {MASS}) is an example of another useful diagnostic plot (and such an answer would be great), however, I am more curious about variations/extensions on existing default diagnostic plots for lm in R (although general other remarks on the topic are always welcomed).

Here are some simple examples of what I mean:

#Some example code for all of us to refer to
set.seed(2542)
x1 <- rnorm(100)
x2 <- runif(100, -2,2)
eps <- rnorm(100,0,2)
y <- 1 + 2*x1 + 3*x2 + eps
y[1:4] <- 14 # adding some contaminated points
fit <- lm(y~x1+x2)

#plot(y~x1+x2)
#summary(fit)

To plot the residuals vs each of the potential x

plot(resid(fit)~x1); abline (h = 0)
plot(resid(fit)~x2); abline (h = 0)
# plot(resid(fit)~x1+x2) # you can also use this, but then you wouldn't be able to use the abline on any plot but the last one

To add the the 0-1 line (how is this line called in English?!) to the qqplot so to see how much the qqline deviates from it

plot(fit, which = 2); abline(0,1, col = "green")

To plot the qq-plot using externally studentized residuals

# plot(fit, which = 2); abline(0,1, col = "green") # The next command is just like this one
qqnorm(rstandard(fit), ylim = c(-2.2,4.2)); qqline(rstudent(fit), lty = 2) ;abline(0,1, col = "green")
qqnorm(rstudent(fit), ylim = c(-2.2,4.2)); qqline(rstudent(fit), lty = 2) ;abline(0,1, col = "green")
# We can note how the "bad" points are more extreme when using the rstudent

Answer 1

Package car has quite a lot of useful functions for diagnostic plots of linear and generalized linear models. Compared to vanilla R plots, they are often enhanced with additional information. I recommend you try example("<function>") on the following functions to see what the plots look like. All plots are described in detail in chapter 6 of Fox & Weisberg. 2011. An R Companion to Applied Regression. 2nd ed.

• residualPlots() plots Pearson residuals against each predictor (scatterplots for numeric variables including a Lowess fit, boxplots for factors)
• marginalModelPlots() displays scatterplots of the response variable against each numeric predictor, inluding a Lowess fit
• avPlots() displays partial-regression plots: for each predictor, this is a scatterplot of a) the residuals from the regression of the response variable on all other predictors against b) the residuals from the regression of the predictor against all other predictors
• qqPlot() for a quantile-quantile plot which includes a confidence envelope
• influenceIndexPlot() displays each value for Cook's distance, hat-value, p-value for outlier test, and studentized residual in a spike-plot against the observation index
• influencePlot() gives a bubble-plot of studentized residuals against hat-values, with the size of the bubble corresponding to Cook's distance, also see dfbetaPlots() and leveragePlots()
• boxCox() displays a profile of the log-likelihood for the transformation parameter $\lambda$ in a Box-Cox power-transform
• crPlots() is for component + residual plots, a variant of which are CERES plots (Combining conditional Expectations and RESiduals), provided by ceresPlots()
• spreadLevelPlot() is for assessing non-constant error variance and displays absolute studentized residuals against fitted values
• scatterplot() provides much-enhanced scatterplots inluding boxplots along the axes, confidence ellipses for the bivariate distribution, and prediction lines with confidence bands
• scatter3d() is based on package rgl and displays interactive 3D-scatterplots including wire-mesh confidence ellipsoids and prediction planes, make sure to run example("scatter3d")

In addition, have a look at bplot() from package rms for another approach to illustrating the common distribution of three variables.

Answer 2

This answer focus on what's available in base R, rather than external packages, although I agree that Fox's package is worth to adopt.

The function influence() (or its wrapper, influence.measures()) returns most of what we need for model diagnostic, including jacknifed statistics. As stated in Chambers and Hastie's Statistical Models in S (Wadsworth & Brooks, 1992), it can be used in combination to summary.lm(). One of the example provided in the so-called "white book" (pp. 130-131) allows to compute standardized (residuals with equal variance) and studentized (the same with a different estimate for SE) residuals, DFBETAS (change in the coefficients scaled by the SE for the regression coefficients), DFFIT (change in the fitted value when observation is dropped), and DFFITS (the same, with unit variance) measures without much difficulty.

Based on your example, and defining the following objects:

lms <- summary(fit)
lmi <- influence(fit)
e <- residuals(fit)
s <- lms$sigma
xxi <- diag(lms$cov.unscaled)
si <- lmi$sigma
h <- lmi$hat
bi <- coef(fit) - coef(lmi)

we can compute the above quantities as follows:

std. residuals    e / (s * (1-h)^.5
stud. residuals   e / (si * (1-h)^.5
dfbetas           bi / (si %o% xxi^.5 
dffit             h * e / (1-h)
dffits            h^.5 * e / (si * (1-h))

(This is Table 4.1, p. 131.)

Chambers and Hastie give the following S/R code for computing DFBETAS:

dfbetas <- function(fit, lms = summary(fit), lmi = lm.influence(fit)) {
  xxi <- diag(lms$cov.unscaled)
  si <- lmi$sigma
  bi <- coef(fit) - coef(lmi)
  bi / (si %o% xxi^0.5)
}

Why do I mention that approach? Because, first, I find this is interesting from a pedagogical perspective (that's what I am using when teaching introductory statistics courses) as it allows to illustrate what can be computed from the output of a fitted linear model fitted in R (but the same would apply with any other statistical package). Second, as the above quantities will be returned as simple vectors or matrices in R, that also means that we can choose the graphics device we want---lattice or ggplot--- to display those statistics, or use them to enhance an existing plot (e.g., highlight DFFITS values in a scatterplot by varying point size cex).