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:
A different approach than my other answer: Plot the ACF using ggplot2.
ggacf <- function(x, ci=0.95, type="correlation", xlab="Lag", ylab=NULL,
ylim=NULL, main=NULL, ci.col="blue", lag.max=NULL) {
x <- as.data.frame(x)
x.acf <- acf(x, plot=F, lag.max=lag.max, type=type)
ci.line <- qnorm((1 - ci) / 2) / sqrt(x.acf$n.used)
d.acf <- data.frame(lag=x.acf$lag, acf=x.acf$acf)
g <- ggplot(d.acf, aes(x=lag, y=acf)) +
geom_hline(yintercept=0) +
geom_segment(aes(xend=lag, yend=0)) +
geom_hline(yintercept=ci.line, color=ci.col, linetype="dashed") +
geom_hline(yintercept=-ci.line, color=ci.col, linetype="dashed") +
theme_bw() +
xlab("Lag") +
ggtitle(ifelse(is.null(main), "", main)) +
if (is.null(ylab))
ylab(ifelse(type=="partial", "PACF", "ACF"))
else
ylab(ylab)
g
}
This seeks to create a similar interface to plot.acf(). Then you can use all of the great features available to ggplot2 plots from the gridExtra package.
library(ggplot2)
library(gridExtra)
grid.arrange(ggacf(lh), ggacf(lh, type="partial"), ncol=2)
Then you get this:
Unfortunately grid.arrange() doesn't work with base graphics, hence the ggplot2 suggestion.
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.
You can use the PerformanceAnalytics package:
library(PerformanceAnalytics)
chart.ACFplus(TS)
