R/influencePlot.R

In car: Companion to Applied Regression

# influencePlot.R

# changed point marking, 25 November 2009 by S. Weisberg
#  deleted the cutoff for Cook's D, and the coloring of the circles
#  inserted default labeling of the id.n largest Cook D.
# 13 January 2009: changed to label points by all of hatvalues,
#  studentized residuals, and Cook's Ds. J. Fox
# 14 April 2010: set id.n = 0. J. Fox
# 23 April 2010: rewrote point marking, S. Weisberg
# 10 May 2010: fixed computation of n
# 2014-04-19: use labels for returned table rownames. J. Fox
# 2015-11-06: now returns Cook's distance, not its square root.  S. Weisberg
# 2017-02-12: consolidated id argument. J. Fox
# 2017-11-30: substitute carPalette() for palette(). J. Fox
# 2019-01-02: added lmerMod method. J. Fox
# 2022-09-21: Fill the bubble points by default. M. Friendly & J. Fox

# moved from Rcmdr 5 December 2006

colscale <- function(x, y, colors, min, max){
  n <- length(colors)
  polygon(x=c(x[1], x[2], x[2], x[1]), y=c(y[1], y[1], y[2], y[2]), xpd=TRUE)
  xincrement <- (x[2] - x[1])/n
  for (i in 1:n){
    xx <- c(x[1] + (i - 1)*xincrement, x[1] + i*xincrement)
    polygon(x=c(xx[1], xx[2], xx[2], xx[1]), y=c(y[1], y[1], y[2], y[2]),
            border=NA, col=colors[i], xpd=TRUE)
  }
  text(x[1], mean(y), labels=min, pos=2, xpd=TRUE)
  text(x[2], mean(y), labels=max, pos=4, xpd=TRUE)
}

influencePlot <- function(model, ...){
    UseMethod("influencePlot")
}

influencePlot.lm <- function(model, scale=10,  
                             xlab="Hat-Values", ylab="Studentized Residuals",
                             id=TRUE, 
                             fill=TRUE, fill.col=carPalette()[2], fill.alpha=0.5,
                             ...){
    id <- applyDefaults(id, defaults=list(method="noteworthy", n=2, cex=1, col=carPalette()[1], location="lr"), type="id")
    if (isFALSE(id)){
        id.n <- 0
        id.method <- "none"
        labels <- id.cex <- id.col <- id.location <- NULL
    }
    else{
        labels <- id$labels
        if (is.null(labels)) labels <- names(na.omit(residuals(model)))
        id.method <- id$method
        id.n <- if ("identify" %in% id.method) Inf else id$n
        id.cex <- id$cex
        id.col <- id$col
        id.location <- id$location
    }
    
    hatval <- hatvalues(model)
    rstud <- rstudent(model)
    if (missing(labels)) labels <- names(rstud)
    cook <- sqrt(cooks.distance(model))
    scale <- scale/max(cook, na.rm=TRUE)
    p <- length(coef(model))
    n <- sum(!is.na(rstud))
    plot(hatval, rstud, xlab=xlab, ylab=ylab, type="n", ...) 
    abline(v=c(2, 3)*p/n, lty=2)
    abline(h=c(-2, 0, 2), lty=2)
    points(hatval, rstud, cex=scale*cook, ...) 
    if (fill) {
      cols <- scales::alpha(fill.col, alpha=fill.alpha*(cook^2/max(cook)^2))
      points(hatval, rstud, cex=scale*cook, col=cols, pch=16)  
      usr <- par("usr")
      left <- usr[1] + 0.2*(usr[2] - usr[1])
      right <- usr[1] + 0.8*(usr[2] - usr[1])
      bot <- usr[4] + strheight("a")
      top <- bot + 2*strheight("A")
      colors <- scales::alpha(fill.col, alpha=seq(0, 1, length=100))
      colscale(c(left, right), c(bot, top), colors, "Cook's D: 0", signif(max(cook)^2, 3))
      }
    if(id.method == "noteworthy"){
        which.rstud <- order(abs(rstud), decreasing=TRUE)[1:id.n]
        which.cook <- order(cook, decreasing=TRUE)[1:id.n]
        which.hatval <- order(hatval, decreasing=TRUE)[1:id.n]
        which.all <- union(which.rstud, union(which.cook, which.hatval))
        id.method <- which.all
    }
    noteworthy <- if (!isFALSE(id)) showLabels(hatval, rstud, labels=labels, method=id.method, 
                             n=id.n, cex=id.cex, col=id.col, location = id.location)
    else NULL
    if (length(noteworthy > 0)){
        result <- data.frame(StudRes=rstud[noteworthy], Hat=hatval[noteworthy],
                             CookD=cook[noteworthy]^2)
        if (is.numeric(noteworthy)) rownames(result) <- labels[noteworthy]
        return(result)
    }
    else return(invisible(NULL))
}

influencePlot.lmerMod <- function(model, ...){
  influencePlot.lm(model, ...)
}