# function without parentheses prints the code:
plotfunction (x, y, ...)
UseMethod("plot")
<bytecode: 0x7f7f782f3b58>
<environment: namespace:base>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:
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 codeLet’s see the code that actually does the work of plot.lm() and plot.merMod():
stats:::plot.lmfunction (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.merModfunction (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>