1# Modified Nov. 24, 2009 by S. Weisberg to use showLabels 2# rather than showExtremes 3# 11 & 20 January 2010: changed lty=3 to lty=1 for fitted curve. J. Fox 4# 14 April 2010: set id.n = 0. J. Fox 5# 15 April 2010; rewrite showLabels 6# 25 May 2010 added grid() to plots, S. Weisberg 7# 15 August 2010, fixed so col= works correctly with plot, but not Boxplot 8# 15 August 2010, deleted pch= argument, as it wasn't used 9# 17 January 2011, allow spline terms; plot against 10# predict(model, type="terms")[[term.name]] 11# 1 February 2011 default for AsIs changed to TRUE 12# 31 March 2011 tukeyNonaddTest updated to check that yhat^2 is not 13# a linear combination of other predictors (as in 1-way anova). 14# 6 April 2011 omit printing lack-of-fit if no lack-of-fit test is possible 15# 16 June 2011 allow layout=NA, in which case the layout is not set in this 16# function, so it is the responsibility of the user 17# 10 Feb 2013: adjusted colinearity check in tukeyNonaddTest 18# 21 March 2013: fixed nonconstant variance test with missing values for glms 19# 11 July 2013: wording changes 20# 11 July 2013: 'groups' arg for residualPlot and residualPlots. 21# 19 July 2014: type='rstudent' fixed 22# 7 October 2014: trapped error resulting from groups= when n<3 23# 25 April 2016: checks for na.action=na.exclude and changes it to na.omit for compatibility with Rcmdr. sw 24# 2017-02-13: consolidated id and smooth arguments. John 25# 2017-11-30: substitute carPalette() for palette(). J. Fox 26# 2019-11-14: change class(x) == "y" to inherits(x, "y") 27# 2018-08-06: enabled spread and var for smoothers. J. Fox 28 29residualPlots <- function(model, ...){UseMethod("residualPlots")} 30 31residualPlots.default <- function(model, terms= ~ . , 32 layout=NULL, ask, main="", 33 fitted=TRUE, AsIs=TRUE, plot=TRUE, tests=TRUE, groups, ...){ 34 mf <- if(!is.null(terms)) termsToMf(model, terms) else NULL 35 # Added for compatibility with Rcmdr 36 if(inherits(model$na.action, "exclude")) model <- update(model, na.action=na.omit) 37 # End addition 38 groups <- if (!missing(groups)) { 39 termsToMf(model, as.formula(paste("~", 40 deparse(substitute(groups)))))$mf.vars[, 2, drop=FALSE] 41 } else { 42 if(is.null(mf$mf.groups)) NULL else 43 mf$mf.groups[, 2, drop=FALSE] 44 } 45 mf <- mf$mf.vars 46 vform <- update(formula(model), attr(mf, "terms")) 47 if(any(is.na(match(all.vars(vform), all.vars(formula(model)))))) 48 stop("Only regressors in the formula can be plotted.") 49 terms <- attr(mf, "term.labels") # this is a list 50 vterms <- attr(terms(vform), "term.labels") 51 # drop interactions (order > 1) 52 vterms <- setdiff(vterms, terms[attr(mf, "order") > 1]) 53 # keep only terms that are numeric or integer or factors or poly 54 good <- NULL 55 for (term in vterms) if( 56 (AsIs == TRUE & inherits(model$model[[term]], "AsIs")) | 57 inherits(model$model[[term]], "numeric") | 58 inherits(model$model[[term]], "integer") | 59 (inherits(model$model[[term]], "factor") & is.null(groups)) | 60 inherits(model$model[[term]], "matrix") | 61 inherits(model$model[[term]], "poly")) good <- c(good, term) 62 nt <- length(good) + fitted 63 nr <- 0 64 if (nt == 0) stop("No plots specified") 65 if (nt > 1 & plot == TRUE & (is.null(layout) || is.numeric(layout))) { 66 if(is.null(layout)){ 67 layout <- switch(min(nt, 9), c(1, 1), c(1, 2), c(2, 2), c(2, 2), 68 c(3, 2), c(3, 2), c(3, 3), c(3, 3), c(3, 3)) 69 } 70 ask <- if(missing(ask) || is.null(ask)) prod(layout)<nt else ask 71 op <- par(mfrow=layout, ask=ask, no.readonly=TRUE, 72 oma=c(0, 0, 1.5, 0), mar=c(5, 4, 1, 2) + .1) 73 on.exit(par(op)) 74 } 75 ans <- NULL 76 if(!is.null(good)){ 77 for (term in good){ 78 nr <- nr + 1 79 qtest <- if(is.null(groups)) 80 residualPlot(model, term, plot=plot, ...) else 81 residualPlot(model, term, plot=plot, groups=groups, ...) 82 if(!is.null(qtest)){ 83 ans <- rbind(ans, qtest) 84 row.names(ans)[nr] <- term} 85 } } 86 # Tukey's test 87 if (fitted == TRUE){ 88 tuk <- if(is.null(groups)) 89 residualPlot(model, "fitted", plot=plot, ...) else 90 residualPlot(model, "fitted", plot=plot, groups=groups, ...) 91 if (!is.null(tuk) & class(model)[1] == "lm"){ 92 ans <- rbind(ans, tuk) 93 row.names(ans)[nr + 1] <- "Tukey test" 94 ans[nr + 1, 2] <- 2*pnorm(abs(ans[nr + 1, 1]), lower.tail=FALSE)}} 95 if(plot == TRUE) mtext(side=3, outer=TRUE, main, cex=1.2) 96 if(!is.null(ans)) { 97 dimnames(ans)[[2]] <- c("Test stat", "Pr(>|Test stat|)") 98 return(if(tests == FALSE | !is.null(groups)) invisible(ans) else 99 if(all(is.na(ans))) warning("No possible lack-of-fit tests") else 100 printCoefmat(ans, has.Pvalue=TRUE, na.print="")) } else 101 invisible(NULL) 102} 103 104residualPlots.lm <- function(model, ...) { 105 residualPlots.default(model, ...) 106} 107 108residualPlots.glm <- function(model, ...) { 109 residualPlots.default(model, ...) 110} 111 112residualPlot <- function(model, ...) UseMethod("residualPlot") 113 114residualPlot.default <- function(model, variable = "fitted", type = "pearson", 115 groups, 116 plot = TRUE, 117 linear = TRUE, 118 quadratic = if(missing(groups)) TRUE else FALSE, 119 smooth=FALSE, id=FALSE, 120 col = carPalette()[1], col.quad = carPalette()[2], 121 pch=1, 122 xlab, ylab, lwd = 1, lty = 1, 123 grid=TRUE, key=!missing(groups), ...) { 124 id <- applyDefaults(id, defaults=list(method="r", n=2, cex=1, col=carPalette()[1], location="lr"), type="id") 125 if (isFALSE(id)){ 126 id.n <- 0 127 id.method <- "none" 128 labels <- id.cex <- id.col <- id.location <- NULL 129 } 130 else{ 131 labels <- id$labels 132 if (is.null(labels)) labels <- names(na.omit(residuals(model))) 133 id.method <- id$method 134 id.n <- if ("identify" %in% id.method) Inf else id$n 135 id.cex <- id$cex 136 id.col <- id$col 137 id.location <- id$location 138 } 139 140 smoother.args <- applyDefaults(smooth, defaults=list(smoother=loessLine, span=2/3, 141 var=FALSE, col=carPalette()[3]), type="smooth") 142 if (!isFALSE(smoother.args)) { 143 smoother <- smoother.args$smoother 144 col.smooth <- smoother.args$col 145 smoother.args$smoother <- smoother.args$col <- NULL 146 if (is.null(smoother.args$spread)) smoother.args$spread <- smoother.args$var 147 } 148 else smoother <- NULL 149 string.capitalize <- function(string) { 150 paste(toupper(substring(string, 1, 1)), substring(string, 2), sep="")} 151 # if(missing(labels)) 152 # labels <- names(residuals(model)[!is.na(residuals(model))]) 153 ylab <- if(!missing(ylab)) ylab else 154 paste(string.capitalize(type), "residuals") 155 column <- match(variable, names(model$model)) 156 # Added for compatibility with Rcmdr 157 if(inherits(model$na.action, "exclude")) model <- update(model, na.action=na.omit) 158 # End addition 159 if(is.na(column) && variable != "fitted") 160 stop(paste(variable, "is not a regressor in the mean function")) 161 horiz <- if(variable == "fitted") predict(model) else model$model[[column]] 162 lab <- if(variable == "fitted") { 163 if(inherits(model, "glm")) 164 "Linear Predictor" else "Fitted values"} else variable 165 lab <- if(!missing(xlab)) xlab else lab 166 if(class(horiz)[1] == "ordered") horiz <- factor(horiz, ordered=FALSE) 167 ans <- 168 if(inherits(horiz, "poly")) { 169 horiz <- horiz[ , 1] 170 lab <- paste("Linear part of", lab) 171 c(NA, NA)} 172 else if (inherits(horiz, "matrix")) { 173 horiz <- try(predict(model, type="terms"), silent=TRUE) 174 if(inherits(horiz, "try-error")) 175 stop("Could not plot spline terms") 176 warning("Splines replaced by a fitted linear combination") 177 horiz <- horiz[ , variable] 178 c(NA, NA) 179 } 180 else if (inherits(horiz, "factor")) c(NA, NA) 181 else residCurvTest(model, variable) 182 # are there groups 183 if(!missing(groups)){ 184 if(is.data.frame(groups)){ 185 groups.name <- names(groups)[1] 186 groups <- groups[, 1, drop=TRUE] 187 } else 188 groups.name <- deparse(substitute(groups)) 189 groups <- if(class(groups)[1] == "factor") groups else factor(groups, ordered=FALSE) 190 if(key){ 191 mar3 <- 1.1 + length(levels(groups)) 192 op <- par(mar=c(5.1, 4.1, mar3, 2.1)) 193 on.exit(par(op)) 194 } 195 colors <- if(length(col) >=length(levels(groups))) col else carPalette() 196 col <- colors[as.numeric(groups)] 197 pchs <- if(length(pch) >= length(levels(groups))) pch else 1:length(levels(groups)) 198 pch <- pchs[as.numeric(groups)] 199 } 200 theResiduals <- switch(type, "rstudent"=rstudent(model), 201 "rstandard"=rstandard(model), residuals(model, type=type)) 202 if(plot==TRUE){ 203 if(inherits(horiz, "factor")) { 204 idm <- if(is.list(id.method)) { 205 lapply(id.method, function(x) if(x[1]=="xy") "y" else x)} else { 206 if(id.method[1] == "xy") "y"} 207 Boxplot(theResiduals, horiz, xlab=lab, ylab=ylab, labels=labels, 208 id.method=idm, id.n=id.n, id.cex=id.cex, 209 id.col=id.col, id.location=id.location, ...) 210 abline(h=0, lty=2) } else 211 { 212 plot(horiz, theResiduals, xlab=lab, ylab=ylab, type="n", ...) 213 if(grid){ 214 grid(lty=1, equilogs=FALSE) 215 box()} 216 points(horiz, theResiduals, col=col, pch=pch, ...) 217 if(linear){ 218 if(missing(groups)){abline(h=0, lty=2, lwd=2)} else { 219 for (g in 1:length(levels(groups))) 220 try(abline(lm(theResiduals ~ horiz, 221 subset=groups==levels(groups)[g]), lty=2, lwd=2, 222 col=colors[g]), silent=TRUE) 223 }} 224 if(quadratic){ 225 new <- seq(min(horiz), max(horiz), length=200) 226 if(missing(groups)){ 227 if(length(unique(horiz)) > 2){ 228 lm2 <- lm(theResiduals ~ poly(horiz, 2)) 229 lines(new, predict(lm2, list(horiz=new)), lty=1, lwd=2, col=col.quad) 230 }} else { 231 for (g in 1:length(levels(groups))){ 232 if(length(unique(horiz)) > 2){ 233 lm2 <- lm(theResiduals~poly(horiz, 2), 234 subset=groups==levels(groups)[g]) 235 lines(new, predict(lm2, list(horiz=new)), lty=1, lwd=1.5, col=colors[g]) 236 }}}} 237 if(is.function(smoother)) 238 if(missing(groups)){ 239 smoother(horiz, theResiduals, col.smooth, log.x=FALSE, log.y=FALSE, 240 spread=smoother.args$spread, smoother.args=smoother.args)} else 241 for (g in 1:length(levels(groups))){ 242 sel <- groups == levels(groups)[g] 243 smoother(horiz[sel], theResiduals[sel], colors[g], log.x=FALSE, log.y=FALSE, 244 spread=smoother.args$spread, smoother.args=smoother.args)} 245 if(key & !missing(groups)){ 246 items <- paste(groups.name, levels(groups), sep= " = ") 247 plotArrayLegend("top", items=items, col.items=colors, pch=pchs) 248 } 249 showLabels(horiz, theResiduals, labels=labels, 250 method=id.method, n=id.n, cex=id.cex, 251 col=id.col, location=id.location) 252 } 253 } 254 invisible(ans)} 255 256residCurvTest <- function(model, variable) {UseMethod("residCurvTest")} 257residCurvTest.lm <- function(model, variable) { 258 if(variable == "fitted") tukeyNonaddTest(model) else { 259 if(is.na(match(variable, attr(model$terms, "term.labels")))) 260 stop(paste(variable, "is not a term in the mean function")) else { 261 xsqres <- qr.resid(model$qr, model.frame(model)[[variable]]^2) 262 r <- residuals(model, type="pearson") 263 m1 <- lm(r ~ xsqres, weights=weights(model)) 264 df.correction <- sqrt((df.residual(model)-1) / df.residual(m1)) 265 test <- summary(m1)$coef[2, 3] * df.correction 266 c(Test=test, Pvalue=2 * pt(-abs(test), df.residual(model)-1)) 267 }}} 268 269residCurvTest.glm <- function(model, variable) { 270 if(variable == "fitted") c(NA, NA) else { 271 if(is.na(match(variable, attr(model$terms, "term.labels")))) 272 stop(paste(variable, "is not a term in the mean function")) else { 273 newmod <- paste(" ~ . + I(", variable, "^2)") 274 m2 <- update(model, newmod, start=NULL) 275 c(Test= test<-deviance(model)-deviance(m2), Pvalue=1-pchisq(test, 1)) 276 }}} 277 278residCurvTest.negbin <- function(model, variable) { 279 if(variable == "fitted") c(NA, NA) else { 280 if(is.na(match(variable, attr(model$terms, "term.labels")))) 281 stop(paste(variable, "is not a term in the mean function")) else { 282 newmod <- paste(" ~ . + I(", variable, "^2)") 283 m2 <- update(model, newmod, start=NULL) 284 c(Test= test<-m2$twologlik - model$twologlik, Pvalue=1-pchisq(test, 1)) 285 }}} 286 287tukeyNonaddTest <- function(model){ 288 tol <- model$qr$tol 289 qr <- model$qr 290 fitsq <- predict(model, type="response")^2 291 fitsq <- qr.resid(qr, fitsq/sqrt(sum(fitsq^2))) 292 if(sd(fitsq) < tol) { 293 return(c(Test=NA, Pvalue=NA)) 294 } else { 295 r <- residuals(model, type="pearson") 296 m1 <- lm(r ~ fitsq, weights=weights(model)) 297 df.correction <- sqrt((df.residual(model) - 1)/df.residual(m1)) 298 tukey <- summary(m1)$coef[2, 3] * df.correction 299 c(Test=tukey, Pvalue=2*pnorm(-abs(tukey))) 300 } 301} 302 303 304 305residualPlot.lm <- function(model, ...) { 306 residualPlot.default(model, ...) 307} 308 309residualPlot.glm <- function(model, variable = "fitted", type = "pearson", 310 plot = TRUE, quadratic = FALSE, 311 smooth=TRUE, ...){ 312 residualPlot.default(model, variable=variable, type=type, plot=plot, 313 quadratic=quadratic, smooth=smooth, ...) 314} 315 316