{{ message }}

Instantly share code, notes, and snippets.

# ellisp/dualplot.R

Last active Jul 20, 2020
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 dualplot <- function(x1, y1, y2, x2 = x1, col = c("#C54E6D", "#009380"), lwd = c(1, 1), colgrid = NULL, mar = c(3, 6, 3, 6) + 0.1, ylab1 = paste(substitute(y1), collapse = ""), ylab2 = paste(substitute(y2), collapse = ""), nxbreaks = 5, yleg1 = paste(gsub("\n\$", "", ylab1), "(left axis)"), yleg2 = paste(ylab2, "(right axis)"), ylim1 = NULL, ylim2 = NULL, ylim.ref = NULL, xlab = "", main = NULL, legx = "topleft", legy = NULL, silent = FALSE, bty = "n", ...){ # Base graphics function for drawing dual axis line plot. # Assumed to be two time series on a conceptually similar, non-identical scale # # Assumes data is in sequence of x1 and of x2 ie ordered by time # # Use with caution! # Please don't use to show growth rates and the original # series at the same time! # # Peter Ellis, 16-27 August 2016, GNU GPL-3 # most parameters should be obvious: # x1 and y1 are the x and y coordinates for first line # x2 and y2 are the x and y coordinates for second line. Often x2 will == x1, but can be overridden # ylim1 and ylim2 are the vertical limits of the 2 axes. Recommended NOT to use these, as # the default algorithm will set them in a way that makes the axes equivalent to using an index (for # positive data) or mean of each series +/- 3 standard deviations (if data include negatives) # ylim.ref the two numbers in the two series to use as the reference point for converting them to indices # when drawing on the page. If both elements are 1, both series will start together at the left of the plot. # nbreaks is number of breaks in horizontal axis # lwd and mar are graphics parameters (see ?par) # colgrid is colour of gridlines; if NULL there are no gridlines # ylab1 and ylab2 are the labels for the two y axes # yleg1 and yleg2 are the labels for the two series in the legend # xlab and main are for x label and main title as in plot() # legx and legy are x and y position fed through to legend() # ... is parameters to pass to legend() # Note that default colours were chosen by colorspace::rainbow_hcl(2, c = 80, l = 50) # strip excess attributes (eg xts etc) from the two vertical axis variables ylab1 <- as.character(ylab1) ylab2 <- as.character(ylab2) y1 <- as.numeric(y1) y2 <- as.numeric(y2) # is ylim.ref is NULL, calculate a good default if(is.null(ylim.ref)){ if (length(y1) == length(y2)){ ylim.ref <- c(1, 1) } else { if (min(x1) > min(x2)){ ylim.ref <- c(1, which(abs(x2 - min(x1)) == min(abs(x2 - min(x1))))) } else { ylim.ref <- c(which(abs(x1 - min(x2)) == min(abs(x1 - min(x2)))), 1) } } } oldpar <- par(mar = mar) xbreaks <- round(seq(from = min(c(x1, x2)), to = max(c(x1, x2)), length.out = nxbreaks)) # unless ylim1 or ylim2 were set, we set them to levels that make it equivalent # to a graphic drawn of indexed series (if all data positive), or to the mean # of each series +/- three standard deviations if some data are negative if(is.null(ylim1) & is.null(ylim2)){ if(min(c(y1, y2), na.rm = TRUE) < 0){ message("With negative values ylim1 or ylim2 need to be chosen by a method other than treating both series visually as though they are indexed. Defaulting to mean value +/- 3 times the standard deviations.") ylim1 <- c(-3, 3) * sd(y1, na.rm = TRUE) + mean(y1, na.rm = TRUE) ylim2 <- c(-3, 3) * sd(y2, na.rm = TRUE) + mean(y2, na.rm = TRUE) } if(ylim.ref > length(y1)){ stop("ylim.ref must be a number shorter than the length of the first series.") } if(ylim.ref > length(y2)){ stop("ylim.ref must be a number shorter than the length of the second series.") } if(!silent) message("The two series will be presented visually as though they had been converted to indexes.") # convert the variables to indexes (base value of 1 at the time specified by ylim.ref) ind1 <- as.numeric(y1) / y1[ylim.ref] ind2 <- as.numeric(y2) / y2[ylim.ref] # calculate y axis limits on the "index to 1" scale indlimits <- range(c(ind1, ind2), na.rm = TRUE) # convert that back to the original y axis scales ylim1 = indlimits * y1[ylim.ref] ylim2 = indlimits * y2[ylim.ref] } else { if(!silent) warning("You've chosen to set at least one of the vertical axes limits manually. Up to you, but it is often better to leave it to the defaults.") } # draw first series - with no axes. plot(x1, y1, type = "l", axes = FALSE, lwd = lwd, xlab = xlab, ylab = "", col = col, main = main, xlim = range(xbreaks), ylim = ylim1) # add in the gridlines if wanted: if(!is.null(colgrid)){ grid(lty = 1, nx = NA, ny = NULL, col = colgrid) abline(v = xbreaks, col = colgrid) } # add in the left hand vertical axis and its label axis(2, col = col, col.axis= col, las=1 ) ## las=1 makes horizontal labels mtext(paste0("\n", ylab1, "\n"), side = 2, col = col, line = 1.5) # Allow a second plot on the same graph par(new=TRUE) # Plot the second series: plot(x2, y2, xlab="", ylab="", axes = FALSE, type = "l", lwd = lwd, col = col, xlim = range(xbreaks), ylim = ylim2) ## add second vertical axis (on right) and its label mtext(paste0("\n", ylab2, "\n"), side = 4, col = col, line = 4.5) axis(4, col = col, col.axis = col, las=1) # Draw the horizontal time axis axis(1, at = xbreaks, labels = xbreaks) # Add Legend legend(x = legx, y = legy, legend=c(yleg1, yleg2), text.col = col, lty = c(1, 1), lwd = lwd, col = col, bty = bty, ...) par(oldpar) }

### ajtanskanen commented Apr 13, 2018

This code has an error if the first element of y1 is NaN. To fix this, add code such as
while ((is.na(y1[ylim.ref]))&(ylim.ref<length(y1))) {
ylim.ref<-ylim.ref+1
}
while ((is.na(y2[ylim.ref]))&(ylim.ref<length(y2))) {
ylim.ref<-ylim.ref+1
}
after line 50.