Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

R | combine plots that use par(mfrow = ...) internally

Tags:

r

graphics

par

Take plot.acf as an example. Both acf and pacf call this function internally. How can i plot them side by side?

Example:

TS <- ts.union(mdeaths, fdeaths)
acf(TS)
pacf(TS)

I tried to use par(mfrow = c(2,4)) and layout to combine them, but stats:::plot.acf overwrites this. The expected output would be:

enter image description here

like image 257
Rentrop Avatar asked Dec 29 '25 18:12

Rentrop


1 Answers

This isn't an ideal solution, but you can redefine what it means to plot an ACF/PACF by defining plot.acf().

First store the existing version.

old.plot.acf <- plot.acf

Now you can use stats:::plot.acf to get the source and copy/paste into the editor. Remove the part that resets mfrow.

plot.acf <- function(x, ci = 0.95, type = "h", xlab = "Lag", ylab = NULL,
                     ylim = NULL, main = NULL, ci.col = "blue",
                     ci.type = c("white", "ma"), max.mfrow = 6,
                     ask = Npgs > 1 && dev.interactive(), 
                     mar = if (nser > 2) c(3, 2, 2, 0.8) else par("mar"),
                     oma = if (nser > 2) c(1, 1.2, 1, 1) else par("oma"),
                     mgp = if (nser > 2) c(1.5, 0.6, 0) else par("mgp"),
                     xpd = par("xpd"), cex.main = if (nser > 2) 1 else
                     par("cex.main"), verbose = getOption("verbose"), ...) 
{
    ci.type <- match.arg(ci.type)
    if ((nser <- ncol(x$lag)) < 1L) 
        stop("x$lag must have at least 1 column")
    if (is.null(ylab)) 
        ylab <- switch(x$type, correlation = "ACF", covariance = "ACF (cov)", 
                       partial = "Partial ACF")
    if (is.null(snames <- x$snames)) 
        snames <- paste("Series ", if (nser == 1L) 
            x$series
            else 1L:nser)
    with.ci <- ci > 0 && x$type != "covariance"
    with.ci.ma <- with.ci && ci.type == "ma" && x$type == "correlation"
    if (with.ci.ma && x$lag[1L, 1L, 1L] != 0L) {
        warning("can use ci.type=\"ma\" only if first lag is 0")
        with.ci.ma <- FALSE
    }
    clim0 <- if (with.ci) 
        qnorm((1 + ci)/2)/sqrt(x$n.used)
    else c(0, 0)
    Npgs <- 1L
    nr <- nser
    if (nser > 1L) {
        sn.abbr <- if (nser > 2L) 
            abbreviate(snames)
        else snames
        if (nser > max.mfrow) {
            Npgs <- ceiling(nser/max.mfrow)
            nr <- ceiling(nser/Npgs)
        }

        ### NOT INCLUDED: mfrow = rep(nr, 2L)

        opar <- par(mar = mar, oma = oma, 
                    mgp = mgp, ask = ask, xpd = xpd, cex.main = cex.main)
        on.exit(par(opar))
        if (verbose) {
            message("par(*) : ", appendLF = FALSE, domain = NA)
            str(par("mfrow", "cex", "cex.main", "cex.axis", "cex.lab", 
                    "cex.sub"))
        }
    }
    if (is.null(ylim)) {
        ylim <- range(x$acf[, 1L:nser, 1L:nser], na.rm = TRUE)
        if (with.ci) 
            ylim <- range(c(-clim0, clim0, ylim))
        if (with.ci.ma) {
            for (i in 1L:nser) {
                clim <- clim0 * sqrt(cumsum(c(1, 2 * x$acf[-1, 
                                                           i, i]^2)))
                ylim <- range(c(-clim, clim, ylim))
            }
        }
    }
    for (I in 1L:Npgs) for (J in 1L:Npgs) {
        dev.hold()
        iind <- (I - 1) * nr + 1L:nr
        jind <- (J - 1) * nr + 1L:nr
        if (verbose) 
            message("Page [", I, ",", J, "]: i =", paste(iind, 
                                                         collapse = ","), "; j =", paste(jind, collapse = ","), 
                    domain = NA)
        for (i in iind) for (j in jind) if (max(i, j) > nser) {
            frame()
            box(col = "light gray")
        }
        else {
            clim <- if (with.ci.ma && i == j) 
                clim0 * sqrt(cumsum(c(1, 2 * x$acf[-1, i, j]^2)))
            else clim0
            plot(x$lag[, i, j], x$acf[, i, j], type = type, xlab = xlab, 
                 ylab = if (j == 1) 
                     ylab
                 else "", ylim = ylim, ...)
            abline(h = 0)
            if (with.ci && ci.type == "white") 
                abline(h = c(clim, -clim), col = ci.col, lty = 2)
            else if (with.ci.ma && i == j) {
                clim <- clim[-length(clim)]
                lines(x$lag[-1, i, j], clim, col = ci.col, lty = 2)
                lines(x$lag[-1, i, j], -clim, col = ci.col, lty = 2)
            }
            title(if (!is.null(main)) 
                main
                else if (i == j) 
                    snames[i]
                else paste(sn.abbr[i], "&", sn.abbr[j]), line = if (nser > 
                                                                        2) 
                    1
                else 2)
        }
        if (Npgs > 1) {
            mtext(paste("[", I, ",", J, "]"), side = 1, line = -0.2, 
                  adj = 1, col = "dark gray", cex = 1, outer = TRUE)
        }
        dev.flush()
    }
    invisible()
}

Now that this is defined locally, you can set mfrow as needed, do your plotting, then reset the function or clear it from the namespace.

plot.acf <- old.plot.acf

To avoid also having to change plot.pacf() as well, you can just use acf(..., type="partial"), which gets the PACF.

like image 190
Alex A. Avatar answered Dec 31 '25 11:12

Alex A.



Donate For Us

If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!