This suite of functions can be used to compute some of the regression (leave-one-out deletion) diagnostics for linear and generalized linear models discussed in Belsley, Kuh and Welsch (1980), Cook and Weisberg (1982), etc.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48
influence.measures(model) rstandard(model, ...) ## S3 method for class 'lm' rstandard(model, infl = lm.influence(model, do.coef = FALSE), sd = sqrt(deviance(model)/df.residual(model)), type = c("sd.1", "predictive"), ...) ## S3 method for class 'glm' rstandard(model, infl = influence(model, do.coef = FALSE), type = c("deviance", "pearson"), ...) rstudent(model, ...) ## S3 method for class 'lm' rstudent(model, infl = lm.influence(model, do.coef = FALSE), res = infl$wt.res, ...) ## S3 method for class 'glm' rstudent(model, infl = influence(model, do.coef = FALSE), ...) dffits(model, infl = , res = ) dfbeta(model, ...) ## S3 method for class 'lm' dfbeta(model, infl = lm.influence(model, do.coef = TRUE), ...) dfbetas(model, ...) ## S3 method for class 'lm' dfbetas(model, infl = lm.influence(model, do.coef = TRUE), ...) covratio(model, infl = lm.influence(model, do.coef = FALSE), res = weighted.residuals(model)) cooks.distance(model, ...) ## S3 method for class 'lm' cooks.distance(model, infl = lm.influence(model, do.coef = FALSE), res = weighted.residuals(model), sd = sqrt(deviance(model)/df.residual(model)), hat = infl$hat, ...) ## S3 method for class 'glm' cooks.distance(model, infl = influence(model, do.coef = FALSE), res = infl$pear.res, dispersion = summary(model)$dispersion, hat = infl$hat, ...) hatvalues(model, ...) ## S3 method for class 'lm' hatvalues(model, infl = lm.influence(model, do.coef = FALSE), ...) hat(x, intercept = TRUE)
an R object, typically returned by
influence structure as returned by
(possibly weighted) residuals, with proper default.
standard deviation to use, see default.
hat values H[i,i], see default.
type of residuals for
the X or design matrix.
should an intercept column be prepended to
further arguments passed to or from other methods.
The primary high-level function is
influence.measures which produces a
"infl" object tabular display showing the DFBETAS for
each model variable, DFFITS, covariance ratios, Cook's distances and
the diagonal elements of the hat matrix. Cases which are influential
with respect to any of these measures are marked with an asterisk.
cooks.distance provide direct access to the
corresponding diagnostic quantities. Functions
rstudent give the standardized and Studentized residuals
respectively. (These re-normalize the residuals to have unit variance,
using an overall and leave-one-out measure of the error variance
Values for generalized linear models are approximations, as described in Williams (1987) (except that Cook's distances are scaled as F rather than as chi-square values). The approximations can be poor when some cases have large influence.
sd arguments are there
to encourage the use of these direct access functions, in situations
where, e.g., the underlying basic influence measures (from
lm.influence or the generic
Note that cases with
weights == 0 are dropped from all
these functions, but that if a linear model has been fitted with
na.action = na.exclude, suitable values are filled in for the
cases excluded during fitting.
For linear models,
rstandard(*, type = "predictive") provides
leave-one-out cross validation residuals, and the “PRESS”
statistic (PREdictive Sum of Squares, the same as
the CV score) of model
hat() exists mainly for S (version 2)
compatibility; we recommend using
dfbetas, the method
for linear models also works for generalized linear models.
Several R core team members and John Fox, originally in his ‘car’ package.
Belsley, D. A., Kuh, E. and Welsch, R. E. (1980) Regression Diagnostics. New York: Wiley.
Cook, R. D. and Weisberg, S. (1982) Residuals and Influence in Regression. London: Chapman and Hall.
Williams, D. A. (1987) Generalized linear model diagnostics using the deviance and single case deletions. Applied Statistics 36, 181–191.
Fox, J. (1997) Applied Regression, Linear Models, and Related Methods. Sage.
Fox, J. (2002) An R and S-Plus Companion to Applied Regression. Sage Publ.
Fox, J. and Weisberg, S. (2011) An R Companion to Applied Regression, second edition. Sage Publ; https://socserv.mcmaster.ca/jfox/Books/Companion/index.html.
‘plotmath’ for the use of
hat in plot annotation.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55
require(graphics) ## Analysis of the life-cycle savings data ## given in Belsley, Kuh and Welsch. lm.SR <- lm(sr ~ pop15 + pop75 + dpi + ddpi, data = LifeCycleSavings) inflm.SR <- influence.measures(lm.SR) which(apply(inflm.SR$is.inf, 1, any)) # which observations 'are' influential summary(inflm.SR) # only these inflm.SR # all plot(rstudent(lm.SR) ~ hatvalues(lm.SR)) # recommended by some plot(lm.SR, which = 5) # an enhanced version of that via plot(<lm>) ## The 'infl' argument is not needed, but avoids recomputation: rs <- rstandard(lm.SR) iflSR <- influence(lm.SR) identical(rs, rstandard(lm.SR, infl = iflSR)) ## to "see" the larger values: 1000 * round(dfbetas(lm.SR, infl = iflSR), 3) cat("PRESS :"); (PRESS <- sum( rstandard(lm.SR, type = "predictive")^2 )) stopifnot(all.equal(PRESS, sum( (residuals(lm.SR) / (1 - iflSR$hat))^2))) ## Show that "PRE-residuals" == L.O.O. Crossvalidation (CV) errors: X <- model.matrix(lm.SR) y <- model.response(model.frame(lm.SR)) ## Leave-one-out CV least-squares prediction errors (relatively fast) rCV <- vapply(seq_len(nrow(X)), function(i) y[i] - X[i,] %*% .lm.fit(X[-i,], y[-i])$coef, numeric(1)) ## are the same as the *faster* rstandard(*, "pred") : stopifnot(all.equal(rCV, unname(rstandard(lm.SR, type = "predictive")))) ## Huber's data [Atkinson 1985] xh <- c(-4:0, 10) yh <- c(2.48, .73, -.04, -1.44, -1.32, 0) lmH <- lm(yh ~ xh) summary(lmH) im <- influence.measures(lmH) im plot(xh,yh, main = "Huber's data: L.S. line and influential obs.") abline(lmH); points(xh[im$is.inf], yh[im$is.inf], pch = 20, col = 2) ## Irwin's data [Williams 1987] xi <- 1:5 yi <- c(0,2,14,19,30) # number of mice responding to dose xi mi <- rep(40, 5) # number of mice exposed glmI <- glm(cbind(yi, mi -yi) ~ xi, family = binomial) summary(glmI) signif(cooks.distance(glmI), 3) # ~= Ci in Table 3, p.184 imI <- influence.measures(glmI) imI stopifnot(all.equal(imI$infmat[,"cook.d"], cooks.distance(glmI)))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.