4  How S3 and S4 methods work

Both S3 and S4 accomplish the goals of encapsulation and polymorphism. In fact, the polymorphism is accomplished in nearly the same way.

In object-oriented programming, a method is a function that belongs to an object. R’s handling of methods is an extremely important way that the language is flexible and intuitive. When you call the plot() function, you call a generic method that decides what specific method should be used based on the object’s class name. First, let’s recall that you can see the code for a function by writing the function name without parentheses:

# function without parentheses prints the code:
plot
function (x, y, ...) 
UseMethod("plot")
<bytecode: 0x7f7f782f3b58>
<environment: namespace:base>

The UseMethod("plot") part causes R to select from all the functions named like plot.<classname>(), which you can list this way:

methods(plot)
 [1] plot.acf*           plot.data.frame*    plot.decomposed.ts*
 [4] plot.default        plot.dendrogram*    plot.density*      
 [7] plot.ecdf           plot.factor*        plot.formula*      
[10] plot.function       plot.hclust*        plot.histogram*    
[13] plot.HoltWinters*   plot.isoreg*        plot.lm*           
[16] plot.medpolish*     plot.mlm*           plot.ppr*          
[19] plot.prcomp*        plot.princomp*      plot.profile.nls*  
[22] plot.raster*        plot.spec*          plot.stepfun       
[25] plot.stl*           plot.table*         plot.ts            
[28] plot.tskernel*      plot.TukeyHSD*     
see '?methods' for accessing help and source code

Let’s see the code that actually does the work of plot.lm() and plot.merMod():

stats:::plot.lm
function (x, which = c(1, 2, 3, 5), caption = list("Residuals vs Fitted", 
    "Normal Q-Q", "Scale-Location", "Cook's distance", "Residuals vs Leverage", 
    expression("Cook's dist vs Leverage* " * h[ii]/(1 - h[ii]))), 
    panel = if (add.smooth) function(x, y, ...) panel.smooth(x, 
        y, iter = iter.smooth, ...) else points, sub.caption = NULL, 
    main = "", ask = prod(par("mfcol")) < length(which) && dev.interactive(), 
    ..., id.n = 3, labels.id = names(residuals(x)), cex.id = 0.75, 
    qqline = TRUE, cook.levels = c(0.5, 1), cook.col = 8, cook.lty = 2, 
    cook.legendChanges = list(), add.smooth = getOption("add.smooth"), 
    iter.smooth = if (isGlm) 0 else 3, label.pos = c(4, 2), cex.caption = 1, 
    cex.oma.main = 1.25, extend.ylim.f = 0.08) 
{
    dropInf <- function(x, h) {
        if (any(isInf <- h >= 1)) {
            warning(gettextf("not plotting observations with leverage one:\n  %s", 
                paste(which(isInf), collapse = ", ")), call. = FALSE, 
                domain = NA)
            x[isInf] <- NaN
        }
        x
    }
    if (!inherits(x, "lm")) 
        stop("use only with \"lm\" objects")
    if (!is.numeric(which) || any(which < 1) || any(which > 6)) 
        stop("'which' must be in 1:6")
    if ((isGlm <- inherits(x, "glm"))) 
        binomialLike <- family(x)$family == "binomial"
    show <- rep(FALSE, 6)
    show[which] <- TRUE
    r <- if (isGlm) 
        residuals(x, type = "pearson")
    else residuals(x)
    yh <- predict(x)
    w <- weights(x)
    if (!is.null(w)) {
        wind <- w != 0
        r <- r[wind]
        yh <- yh[wind]
        w <- w[wind]
        labels.id <- labels.id[wind]
    }
    n <- length(r)
    if (any(show[2L:6L])) {
        s <- if (inherits(x, "rlm")) 
            x$s
        else if (isGlm) 
            sqrt(summary(x)$dispersion)
        else sqrt(deviance(x)/df.residual(x))
        hii <- (infl <- influence(x, do.coef = FALSE))$hat
        if (any(show[4L:6L])) {
            cook <- cooks.distance(x, infl)
        }
    }
    if (any(show[c(2L, 3L, 5L)])) {
        ylab5 <- ylab23 <- if (isGlm) 
            "Std. Pearson resid."
        else "Standardized residuals"
        rs <- dropInf(if (isGlm) 
            rstandard(x, type = "pearson")
        else (if (is.null(w)) 
            r
        else sqrt(w) * r)/(s * sqrt(1 - hii)), hii)
    }
    if (any(show[5L:6L])) {
        r.hat <- range(hii, na.rm = TRUE)
        isConst.hat <- all(r.hat == 0) || diff(r.hat) < 1e-10 * 
            mean(hii, na.rm = TRUE)
    }
    if (any(show[c(1L, 3L)])) 
        l.fit <- if (isGlm) 
            "Predicted values"
        else "Fitted values"
    if (is.null(id.n)) 
        id.n <- 0L
    else {
        id.n <- as.integer(id.n)
        if (id.n < 0L || id.n > n) 
            stop(gettextf("'id.n' must be in {1,..,%d}", n), 
                domain = NA)
    }
    if (id.n > 0L) {
        if (is.null(labels.id)) 
            labels.id <- paste(1L:n)
        iid <- 1L:id.n
        show.r <- sort.list(abs(r), decreasing = TRUE)[iid]
        if (any(show[2L:3L])) 
            show.rs <- sort.list(abs(rs), decreasing = TRUE)[iid]
        text.id <- function(x, y, ind, adj.x = TRUE, usr = par("usr")) {
            labpos <- if (adj.x) 
                label.pos[(x > mean(usr[1:2])) + 1L]
            else 3
            text(x, y, labels.id[ind], cex = cex.id, xpd = TRUE, 
                pos = labpos, offset = 0.25)
        }
    }
    getCaption <- function(k) if (length(caption) < k) 
        NA_character_
    else as.graphicsAnnot(caption[[k]])
    if (is.null(sub.caption)) {
        cal <- x$call
        if (!is.na(m.f <- match("formula", names(cal)))) {
            cal <- cal[c(1, m.f)]
            names(cal)[2L] <- ""
        }
        cc <- deparse(cal, 80)
        nc <- nchar(cc[1L], "c")
        abbr <- length(cc) > 1 || nc > 75
        sub.caption <- if (abbr) 
            paste(substr(cc[1L], 1L, min(75L, nc)), "...")
        else cc[1L]
    }
    one.fig <- prod(par("mfcol")) == 1
    if (ask) {
        oask <- devAskNewPage(TRUE)
        on.exit(devAskNewPage(oask))
    }
    if (show[1L]) {
        ylim <- range(r, na.rm = TRUE)
        if (id.n > 0) 
            ylim <- extendrange(r = ylim, f = extend.ylim.f)
        dev.hold()
        plot(yh, r, xlab = l.fit, ylab = "Residuals", main = main, 
            ylim = ylim, type = "n", ...)
        panel(yh, r, ...)
        if (one.fig) 
            title(sub = sub.caption, ...)
        mtext(getCaption(1), 3, 0.25, cex = cex.caption)
        if (id.n > 0) {
            y.id <- r[show.r]
            y.id[y.id < 0] <- y.id[y.id < 0] - strheight(" ")/3
            text.id(yh[show.r], y.id, show.r)
        }
        abline(h = 0, lty = 3, col = "gray")
        dev.flush()
    }
    if (show[2L]) {
        ylim <- range(rs, na.rm = TRUE)
        ylim[2L] <- ylim[2L] + diff(ylim) * 0.075
        dev.hold()
        qq <- qqnorm(rs, main = main, ylab = ylab23, ylim = ylim, 
            ...)
        if (qqline) 
            qqline(rs, lty = 3, col = "gray50")
        if (one.fig) 
            title(sub = sub.caption, ...)
        mtext(getCaption(2), 3, 0.25, cex = cex.caption)
        if (id.n > 0) 
            text.id(qq$x[show.rs], qq$y[show.rs], show.rs)
        dev.flush()
    }
    if (show[3L]) {
        sqrtabsr <- sqrt(abs(rs))
        ylim <- c(0, max(sqrtabsr, na.rm = TRUE))
        yl <- as.expression(substitute(sqrt(abs(YL)), list(YL = as.name(ylab23))))
        yhn0 <- if (is.null(w)) 
            yh
        else yh[w != 0]
        dev.hold()
        plot(yhn0, sqrtabsr, xlab = l.fit, ylab = yl, main = main, 
            ylim = ylim, type = "n", ...)
        panel(yhn0, sqrtabsr, ...)
        if (one.fig) 
            title(sub = sub.caption, ...)
        mtext(getCaption(3), 3, 0.25, cex = cex.caption)
        if (id.n > 0) 
            text.id(yhn0[show.rs], sqrtabsr[show.rs], show.rs)
        dev.flush()
    }
    if (show[4L]) {
        if (id.n > 0) {
            show.r <- order(-cook)[iid]
            ymx <- cook[show.r[1L]] * 1.075
        }
        else ymx <- max(cook, na.rm = TRUE)
        dev.hold()
        plot(cook, type = "h", ylim = c(0, ymx), main = main, 
            xlab = "Obs. number", ylab = "Cook's distance", ...)
        if (one.fig) 
            title(sub = sub.caption, ...)
        mtext(getCaption(4), 3, 0.25, cex = cex.caption)
        if (id.n > 0) 
            text.id(show.r, cook[show.r], show.r, adj.x = FALSE)
        dev.flush()
    }
    if (show[5L]) {
        ylim <- range(rs, na.rm = TRUE)
        if (id.n > 0) {
            ylim <- extendrange(r = ylim, f = extend.ylim.f)
            show.rs <- order(-cook)[iid]
        }
        do.plot <- TRUE
        if (isConst.hat) {
            if (missing(caption)) 
                caption[[5L]] <- "Constant Leverage:\n Residuals vs Factor Levels"
            if (nf <- length(xlev <- x$xlevels)) {
                facvars <- names(xlev)
                mf <- model.frame(x)[facvars]
                dm <- data.matrix(mf)
                nlev <- unname(lengths(xlev))
                ff <- if (nf == 1) 
                  1
                else rev(cumprod(c(1, nlev[nf:2])))
                facval <- (dm - 1) %*% ff
                xx <- facval
                dev.hold()
                plot(facval, rs, xlim = c(-1/2, sum((nlev - 1) * 
                  ff) + 1/2), ylim = ylim, xaxt = "n", main = main, 
                  xlab = "Factor Level Combinations", ylab = ylab5, 
                  type = "n", ...)
                axis(1, at = ff[1L] * (1L:nlev[1L] - 1/2) - 1/2, 
                  labels = xlev[[1L]])
                mtext(paste(facvars[1L], ":"), side = 1, line = 0.25, 
                  adj = -0.05)
                abline(v = ff[1L] * (0:nlev[1L]) - 1/2, col = "gray", 
                  lty = "F4")
                panel(facval, rs, ...)
                abline(h = 0, lty = 3, col = "gray")
                dev.flush()
            }
            else {
                message(gettextf("hat values (leverages) are all = %s\n and there are no factor predictors; no plot no. 5", 
                  format(mean(r.hat))), domain = NA)
                frame()
                do.plot <- FALSE
            }
        }
        else {
            xx <- hii
            xx[xx >= 1] <- NA
            dev.hold()
            plot(xx, rs, xlim = c(0, max(xx, na.rm = TRUE)), 
                ylim = ylim, main = main, xlab = "Leverage", 
                ylab = ylab5, type = "n", ...)
            panel(xx, rs, ...)
            abline(h = 0, v = 0, lty = 3, col = "gray")
            if (one.fig) 
                title(sub = sub.caption, ...)
            if (length(cook.levels)) {
                p <- x$rank
                usr2 <- par("usr")[2L]
                hh <- seq.int(min(r.hat[1L], r.hat[2L]/100), 
                  usr2, length.out = 101)
                for (crit in cook.levels) {
                  cl.h <- sqrt(crit * p * (1 - hh)/hh)
                  lines(hh, cl.h, lty = cook.lty, col = cook.col)
                  lines(hh, -cl.h, lty = cook.lty, col = cook.col)
                }
                if (!is.null(cook.legendChanges)) 
                  do.call(legend, modifyList(list(x = "bottomleft", 
                    legend = "Cook's distance", lty = cook.lty, 
                    col = cook.col, text.col = cook.col, bty = "n", 
                    x.intersp = 1/4, y.intersp = 1/8), cook.legendChanges))
                xmax <- min(0.99, usr2)
                ymult <- sqrt(p * (1 - xmax)/xmax)
                aty <- sqrt(cook.levels) * ymult
                axis(4, at = c(-rev(aty), aty), labels = paste(c(rev(cook.levels), 
                  cook.levels)), mgp = c(0.25, 0.25, 0), las = 2, 
                  tck = 0, cex.axis = cex.id, col.axis = cook.col)
            }
            dev.flush()
        }
        if (do.plot) {
            mtext(getCaption(5), 3, 0.25, cex = cex.caption)
            if (id.n > 0) {
                y.id <- rs[show.rs]
                y.id[y.id < 0] <- y.id[y.id < 0] - strheight(" ")/3
                text.id(xx[show.rs], y.id, show.rs)
            }
        }
    }
    if (show[6L]) {
        g <- dropInf(hii/(1 - hii), hii)
        ymx <- max(cook, na.rm = TRUE) * 1.025
        dev.hold()
        plot(g, cook, xlim = c(0, max(g, na.rm = TRUE)), ylim = c(0, 
            ymx), main = main, ylab = "Cook's distance", xlab = expression("Leverage  " * 
            h[ii]), xaxt = "n", type = "n", ...)
        panel(g, cook, ...)
        athat <- pretty(hii)
        axis(1, at = athat/(1 - athat), labels = paste(athat))
        if (one.fig) 
            title(sub = sub.caption, ...)
        p <- x$rank
        bval <- pretty(sqrt(p * cook/g), 5)
        usr <- par("usr")
        xmax <- usr[2L]
        ymax <- usr[4L]
        for (i in seq_along(bval)) {
            bi2 <- bval[i]^2
            if (p * ymax > bi2 * xmax) {
                xi <- xmax + strwidth(" ")/3
                yi <- bi2 * xi/p
                abline(0, bi2, lty = cook.lty)
                text(xi, yi, paste(bval[i]), adj = 0, xpd = TRUE)
            }
            else {
                yi <- ymax - 1.5 * strheight(" ")
                xi <- p * yi/bi2
                lines(c(0, xi), c(0, yi), lty = cook.lty)
                text(xi, ymax - 0.8 * strheight(" "), paste(bval[i]), 
                  adj = 0.5, xpd = TRUE)
            }
        }
        mtext(getCaption(6), 3, 0.25, cex = cex.caption)
        if (id.n > 0) {
            show.r <- order(-cook)[iid]
            text.id(g[show.r], cook[show.r], show.r, usr = usr)
        }
        dev.flush()
    }
    if (!one.fig && par("oma")[3L] >= 1) 
        mtext(sub.caption, outer = TRUE, cex = cex.oma.main)
    invisible()
}
<bytecode: 0x7f7f829df350>
<environment: namespace:stats>
lme4:::plot.merMod
function (x, form = resid(., type = "pearson") ~ fitted(.), abline, 
    id = NULL, idLabels = NULL, grid, ...) 
{
    object <- x
    if (!inherits(form, "formula")) 
        stop("\"form\" must be a formula")
    allV <- all.vars(asOneFormula(form, id, idLabels))
    allV <- allV[is.na(match(allV, c("T", "F", "TRUE", "FALSE", 
        ".obs")))]
    if (length(allV) > 0) {
        data <- getData(object)
        if (is.null(data)) {
            alist <- lapply(as.list(allV), as.name)
            names(alist) <- allV
            alist <- c(list(as.name("data.frame")), alist)
            mode(alist) <- "call"
            data <- eval(alist, sys.parent(1))
        }
        else if (any(naV <- is.na(match(allV, names(data))))) 
            stop(allV[naV], " not found in data")
    }
    else data <- NULL
    dots <- list(...)
    args <- if (length(dots) > 0) 
        dots
    else list()
    data <- as.list(c(as.list(cbind(data, .obs = seq(nrow(data)))), 
        . = list(object)))
    covF <- getCovariateFormula(form)
    .x <- eval(covF[[2]], data)
    if (!is.numeric(.x)) {
        stop("Covariate must be numeric")
    }
    argForm <- ~.x
    argData <- data.frame(.x = .x, check.names = FALSE)
    if (is.null(args$xlab)) {
        if (is.null(xlab <- attr(.x, "label"))) 
            xlab <- deparse(covF[[2]])
        args$xlab <- xlab
    }
    respF <- getResponseFormula(form)
    if (!is.null(respF)) {
        .y <- eval(respF[[2]], data)
        if (is.null(args$ylab)) {
            if (is.null(ylab <- attr(.y, "label"))) 
                ylab <- deparse(respF[[2]])
            args$ylab <- ylab
        }
        argForm <- .y ~ .x
        argData[, ".y"] <- .y
    }
    grpsF <- getGroupsFormula(form)
    if (!is.null(grpsF)) {
        gr <- splitFormula(grpsF, sep = "*")
        for (i in seq_along(gr)) {
            auxGr <- all.vars(gr[[i]])
            for (j in auxGr) argData[[j]] <- eval(as.name(j), 
                data)
        }
        argForm <- as.formula(paste(if (length(argForm) == 2) 
            "~ .x |"
        else ".y ~ .x |", deparse(grpsF[[2]])))
    }
    args <- c(list(argForm, data = argData), args)
    if (is.null(args$strip)) {
        args$strip <- function(...) strip.default(..., style = 1)
    }
    if (is.null(args$cex)) 
        args$cex <- par("cex")
    if (is.null(args$adj)) 
        args$adj <- par("adj")
    if (!is.null(id)) {
        idResType <- "pearson"
        id <- switch(mode(id), numeric = {
            if (id <= 0 || id >= 1) stop(shQuote("id"), " must be between 0 and 1")
            abs(resid(object, type = idResType))/sigma(object) > 
                -qnorm(id/2)
        }, call = eval(asOneSidedFormula(id)[[2]], data), stop(shQuote("id"), 
            " can only be a formula or numeric."))
        if (is.null(idLabels)) {
            idLabels <- getIDLabels(object)
        }
        else {
            if (inherits(idLabels, "formula")) {
                idLabels <- getIDLabels(object, idLabels)
            }
            else if (is.vector(idLabels)) {
                if (length(idLabels <- unlist(idLabels)) != length(id)) {
                  stop("\"idLabels\" of incorrect length")
                }
            }
            else stop("\"idLabels\" can only be a formula or a vector")
        }
        idLabels <- as.character(idLabels)
    }
    if (missing(abline)) {
        abline <- if (missing(form)) 
            c(0, 0)
        else NULL
    }
    assign("abl", abline)
    if (length(argForm) == 3) {
        if (is.numeric(.y)) {
            plotFun <- "xyplot"
            if (is.null(args$panel)) {
                args <- c(args, panel = list(function(x, y, subscripts, 
                  ...) {
                  x <- as.numeric(x)
                  y <- as.numeric(y)
                  dots <- list(...)
                  if (grid) panel.grid()
                  panel.xyplot(x, y, ...)
                  if (any(ids <- id[subscripts])) {
                    ltext(x[ids], y[ids], idLabels[subscripts][ids], 
                      cex = dots$cex, adj = dots$adj)
                  }
                  if (!is.null(abl)) {
                    if (length(abl) == 2) panel.abline(a = abl, 
                      ...) else panel.abline(h = abl, ...)
                  }
                }))
            }
        }
        else {
            plotFun <- "bwplot"
            if (is.null(args$panel)) {
                args <- c(args, panel = list(function(x, y, ...) {
                  if (grid) panel.grid()
                  panel.bwplot(x, y, ...)
                  if (!is.null(abl)) {
                    panel.abline(v = abl[1], ...)
                  }
                }))
            }
        }
    }
    else {
        plotFun <- "histogram"
        if (is.null(args$panel)) {
            args <- c(args, panel = list(function(x, ...) {
                if (grid) panel.grid()
                panel.histogram(x, ...)
                if (!is.null(abl)) {
                  panel.abline(v = abl[1], ...)
                }
            }))
        }
    }
    if (missing(grid)) {
        grid <- (plotFun == "xyplot")
    }
    do.call(plotFun, as.list(args))
}
<bytecode: 0x7f7f85beb078>
<environment: namespace:lme4>