Scatterplot matrix with linear model fit and correlation (base graphics)

pairs.cor.R

#' @param formula,data,subset,na.action arguments to `model.frame`. If the first argument is a data frame, the rest is ignored. 
#' @param method argument to `cor.test`. A character string indicating which correlation coefficient is to be used for the test. One of "pearson", "kendall", or "spearman", can be abbreviated.
#' @param n integer; the number of x values at which to evaluate.
#' @param lty,col,cex graphical parameters fit line type, line colour, character expansion for labels.
#' @param \dots optional, additional arguments passed to `points`
#' @example
#' data <- data.frame(x = runif(100), x1 = rlnorm(100), x2 = rnorm(100))
#' pairs.cor(~ x + log(x1) + x2, data)
pairs.cor <-
function(formula, data = NULL, subset = NULL, na.action = na.fail,
    method = c("pearson", "kendall", "spearman"), n = 25,
    lty = c(1L, 2L, 2L), col = par("col"), cex = 1,
        main = switch(method, "pearson" = "Pearson's product-moment correlation",
               "kendall" = "Kendall's rank correlation tau",
               "spearman" = "Spearman's rank correlation rho"
              ), ...) {
         
    
    data <- if(inherits(formula, "formula")) {
        cl <- match.call()
       names(formals(model.frame.default))
        cl <- cl[c(TRUE, names(cl)[-1] %in% names(formals(model.frame.default)))]
        cl[[1]] <- as.name("model.frame")
        eval.parent(cl) 
    } else formula

    
    col <- rep(col, length.out = 3L)
    lty <- rep(lty, length.out = 3L)
    method <- match.arg(method)
    stat <- switch(method, "pearson" = quote(r[p]),
               "kendall" = quote(tau),
               "spearman" = quote(rho)
              )
    pairs(data, upper.panel = function(x, y, ...) {
        points(x, y, ...)
        if(length(unique(x)) > 1) {
            i <- is.finite(x) & is.finite(y)
            x <- x[i]
            y <- y[i]
            fm <- lm(y ~ x)
            nx <- seq(min(x), max(x), length.out = n)
            ny <- predict(fm, list(x = nx), interval = "confidence")
            matplot(nx, ny, type = "l", lty = lty, col = col, add = TRUE)
        }
    }, lower.panel = function(x, y) {
        test <- cor.test(x, y, method = method)
        text(grconvertX(.5, "npc", "user"), grconvertY(.5, "npc", "user"),
            bquote(atop(
                .(stat) == .(round(test$estimate, 2)),
                p == .(round(test$p.value, digits = 3))
            )), cex = cex)
    }, main = main)
}

Example:

dat <- data.frame(x = runif(100), x1 = rlnorm(100), x2 = rnorm(100))
pairs.cor(~ x + log(x1) + x2, dat)