R/MCMCplot.R
In caseyyoungflesh/MCMCvis: Tools to Visualize, Manipulate, and Summarize MCMC Output
Defines functions
MCMCplot
Documented in
MCMCplot
#' Caterpillar plots of posterior distributions from MCMC output
#'
#' Visualize posterior distributions from MCMC output for specific parameters of interest using caterpillar plots. Color of median dot represents the overlap of the posterior distribution with 0 (or other specified value).
#'
#'
#' @param object Object containing MCMC output. See DETAILS below.
#'
#' @param object2 Optional second object containing MCMC output. If specified, parameter estimates from each model will be displayed in a paired manner. Parameter names for \code{'object'} and \code{'object2'} must be identical. See DETAILS below.
#'
#' @param params Character string (or vector of character strings) denoting parameters to be plotted.
#'
#' Default \code{'all'} plots posteriors for all parameters. See VALUE below.
#'
#' @param HPD Logical specifying whether to calculate equal-tailed credible intervals (\code{HPD = FALSE}) or highest posterior density intervals (\code{HPD = TRUE}) for the selected parameters. Default is \code{HPD = FALSE}.
#'
#' @param ci Numeric vector of length 2, where each element is (0,100] and represents the width of an equal-tailed (\code{HPD = FALSE}) or highest posterior density (\code{HPD = TRUE}) credible interval. The first element of this vector corresponds to the thicker (narrower) credible interval displayed on the plot (default is 0.5) and the second element of this vector corresponds to the thinner (wider) credible interval (default is 0.95). The first credible interval width (\code{ci[1]}) must be less than or equal to the width of the second credible interval (\code{ci[2]}).
#'
#' @param excl Character string (or vector of character strings) denoting parameters to exclude. Used in conjunction with \code{params} argument to select parameters of interest.
#'
#' @param ISB Ignore Square Brackets (ISB). Logical specifying whether square brackets should be ignored in the \code{params} and \code{excl} arguments. If \code{TRUE}, square brackets are ignored. If \code{FALSE}, square brackets are not ignored. This allows partial names to be used when specifying parameters of interest. Use \code{exact} argument to specify whether input from \code{params} and \code{excl} arguments should be matched exactly.
#'
#' @param exact Logical specifying whether input from \code{params} and \code{excl} arguments should be matched exactly (after ignoring square brackets if \code{ISB = FALSE}). If \code{TRUE}, input from \code{params} and \code{excl} are matched exactly (after taking \code{ISB} argument into account). If \code{FALSE}, input from \code{params} and \code{excl} are matched using regular expression format (after taking \code{ISB} argument into account).
#'
#' @param ref Value indicating where vertical reference line should be created and what value to use a reference for caterpillar median coloration.
#'
#' Default is \code{ref = 0}.
#'
#' Argument \code{NULL} will plot no reference line.
#'
#' @param ref_ovl Logical specifying whether the style/color of plotted median dots and CI should be changed based on whether the specified credible intervals (50 \% and 95 \% by default) overlap the reference line. See DETAILS for more information.
#'
#' @param col Character string (or vector of character strings) specifying which color to render estimates on plot. When \code{ref_ovl = TRUE}, this argument has no effect and colors plotted will be based on the credible intervals and reference line. Number of specified colors must equal the number of specified parameters or one.
#'
#' @param col2 Character string (or vector of character strings) specifying which color to render estimates on plot for \code{object2} (if specified). Number of specified colors must equal the number of specified parameters or one. Red by default.
#'
#' @param offset Value indicating how much to offset plotted posteriors when \code{object2} is specified (i.e., control the amount of space between the two plotted posteriors for each parameter). The distance from one set of parameters to another corresponds to a value of 1.
#'
#' @param rank Logical specifying whether output should be ranked. If \code{TRUE} posteriors will be ranked in decreasing order (based on specified measure of centrality) from top down.
#'
#' @param horiz Logical specifying orientation of plot. If \code{TRUE} posteriors will be plotted running horizontally (parallel to the x-axis). If \code{FALSE} posteriors will be plotted running vertically (perpendicular to the x-axis).
#'
#' @param xlim Numerical vector of length 2, indicating range of x-axis. Only applicable if \code{horiz = TRUE}.
#'
#' @param ylim Numerical vector of length 2, indicating range of y-axis. Only applicable if \code{horiz = FALSE}.
#'
#' @param xlab Character string labeling x-axis. Only applicable if \code{horiz = TRUE}.
#'
#' Default label is 'Parameter Estimate'. Option \code{NULL} will return plot with no label on x-axis.
#' @param ylab Character string labeling y-axis. Only applicable if \code{horiz = FALSE}.
#'
#' Default label is 'Parameter Estimate'. Option \code{NULL} will return plot with no label on y-axis.
#' @param main Character string indicating title of plot.
#'
#' @param labels Character string (or vector of character strings if plotting > 1 parameter) labeling parameter estimates along y-axis (if \code{horiz = FALSE}) or x-axis (if \code{horiz = TRUE}).
#'
#' Default option will use parameter names from \code{object}.
#'
#' Option \code{NULL} will return plot with no labels on axis.
#'
#' @param guide_lines Logical specifying whether to plot reference lines for each parameter in order to better visualize which parameter names correspond to each posterior.
#'
#' @param guide_axis Logical specifying whether a second axis should be plotted (x-axis if \code{HORIZ = TRUE}, y-axis if \code{HORIZ = FALSE}) to help interpret values on plot.
#'
#' @param sz_labels Number specifying size of text for parameter labels on axis.
#'
#' @param sz_med Number specifying size of points represents posterior medians.
#'
#' @param sz_thick Number specifying thickness of 50 percent CI line (thicker line).
#'
#' @param sz_thin Number specifying thickness of 95 percent CI line (thinner line).
#'
#' @param sz_ax Number specifying thickness of axis and ticks.
#'
#' @param sz_ax_txt Number specifying size of text for axis label.
#'
#' @param sz_tick_txt Number specifying size of text for tick labels on axis.
#'
#' @param sz_main_txt Number specifying size of text for main title.
#'
#' @param pos_tick Numeric vector specifying where ticks on axis should be placed.
#'
#' @param mar Numerical vector of length 4 specifying plot margins - (BOTTOM, LEFT, TOP, RIGHT). Changes to the margin should be made within the function rather than using the \code{par} call.
#'
#' Default is c(5.1, 4.1, 4.1, 2.1) - the R plot default.
#'
#' @section Details:
#' Points represent posterior medians. Parameters where the smaller specified credible intervals (first value in vector provided to \code{ci} argument; 50\% by default) overlap 0 (or other value specified by the \code{ref} argument) are indicated by 'open' circles and a lighter color than that specified (e.g., gray will be displayed with the default black is specified for \code{col}). Parameters where the specific smaller credible intervals DO NOT overlap 0 (or other specified value) AND the larger specified credible intervals (second value in vector provided to \code{ci} argument; 95\% by default) DO overlap 0 (or other specified value) are indicated by 'closed' circles and a lighter color that that specified. Parameters where the larger specified credible intervals DO NOT overlap 0 (or other specified value) are indicated by 'closed' circles and the color specified (black by default). Thick lines represent the smaller specified credible intervals percent credible intervals (50\% by default) while thin lines represent the larger specified credible intervals (95 \% by default). \code{ref_ovl = TRUE} can be used to enable this feature. When two model objects are supplied to the function (with \code{object} and \code{object2}) and no argument is supplied to \code{col} or \code{col2}, light red lines can be interpreted as analogous to gray lines. When \code{ref_ovl = TRUE} and a color (or colors) other than the default are specified (with the \code{col} and \code{col2} arguments) lighter versions of the color specified are used in black of the light gray and/or light red lines.
#'
#' When \code{object2} is specified, paired caterpillar plots of each parameter are produced. For this reason, parameter names of \code{object} and \code{object2} specified with the \code{params} argument must be identical (to be used for comparing posterior estimates of similar models). \code{col} and \code{col2} arguments can be specified to change the color of output from \code{object} and \code{object2}, respectively. By default, output from \code{object} is plotted in black and \code{object2} is plotted in red. The \code{ref_ovl} argument can also be specified.
#'
#' \code{object} argument can be a \code{stanfit} object (\code{rstan} package), a \code{CmdStanMCMC} object (\code{cmdstanr} package), a \code{stanreg} object (\code{rstanarm} package), a \code{brmsfit} object (\code{brms} package), an \code{mcmc.list} object (\code{coda} and \code{rjags} packages), \code{mcmc} object (\code{coda} and \code{nimble} packages), \code{list} object (\code{nimble} package), an \code{R2jags} model object (\code{R2jags} package), a \code{jagsUI} model object (\code{jagsUI} package), or a matrix containing MCMC chains (each column representing MCMC output for a single parameter, rows representing iterations in the chain). The function automatically detects the object type and proceeds accordingly.
#'
#' @section Notes:
#'
#' When specifying \code{rank = TRUE} and specifying labels for \code{labels}, labels will be applied to parameters before they are ranked.
#'
#' Thanks to Cinner et al. 2016, whose Fig. 1 inspired this plot.
#'
#' @section References:
#'
#' Cinner, J. E., C. Huchery, M. A. MacNeil, N. A. J. Graham, T. R. McClanahan, J. Maina, E. Maire, J. N. Kittinger, C. C. Hicks, C. Mora, E. H. Allison, S. D'Agata, A. Hoey, D. A. Feary, L. Crowder, I. D. Williams, M. Kulbicki, L. Vigliola, L. Wantiez, G. Edgar, R. D. Stuart-Smith, S. A. Sandin, A. L. Green, M. J. Hardt, M. Beger, A. Friedlander, S. J. Campbell, K. E. Holmes, S. K. Wilson, E. Brokovich, A. J. Brooks, J. J. Cruz-Motta, D. J. Booth, P. Chabanet, C. Gough, M. Tupper, S. C. A. Ferse, U. R. Sumaila, and D. Mouillot. 2016. Bright spots among the world's coral reefs. Nature 535:416-419.
#'
#' @examples
#' #Load data
#' data(MCMC_data)
#'
#' #Plot MCMC output
#' MCMCplot(MCMC_data, labels = NULL)
#'
#' #Just 'beta' parameters
#' MCMCplot(MCMC_data, params = 'beta')
#'
#' #Just 'beta' parameters using highest posterior density intervals
#' MCMCplot(MCMC_data, params = 'beta', HPD = TRUE)
#
#' #Just 'beta[1]', 'beta[4]', and 'alpha[3]'
#' MCMCplot(MCMC_data, params = c('beta[1]', 'beta[4]', 'alpha[3]'), ISB = FALSE, exact = TRUE)
#'
#' #Just 'beta[1]', 'beta[4]', and 'alpha[3]' and change the credible interval widths
#' MCMCplot(MCMC_data, ci = c(50, 89), params = c('beta[1]', 'beta[4]', 'alpha[3]'),
#' ISB = FALSE, exact = TRUE)
#'
#' #Rank parameters by posterior mean
#' MCMCplot(MCMC_data, params = 'beta', rank = TRUE)
#'
#' #Create vertical plot
#' MCMCplot(MCMC_data, params = 'beta', horiz = FALSE)
#'
#' @export
#'
MCMCplot <- function(object,
object2 = NULL,
params = 'all',
HPD = FALSE,
ci = c(50, 95),
excl = NULL,
ISB = TRUE,
exact = TRUE,
ref = 0,
ref_ovl = FALSE,
col = 'black',
col2 = 'red',
offset = 0.1,
rank = FALSE,
horiz = TRUE,
xlim,
ylim,
xlab,
ylab,
main,
labels,
guide_lines = FALSE,
guide_axis = TRUE,
sz_labels = 1.2,
sz_med = 1.5,
sz_thick = 5,
sz_thin = 2,
sz_ax = 3,
sz_ax_txt = 1.3,
sz_tick_txt = 1.2,
sz_main_txt = 1.2,
pos_tick,
mar = c(5.1, 4.1, 4.1, 2.1)) {
####################################################################################################
# load object(s)
data <- MCMCchains(object, params = params, excl = excl, ISB = ISB, exact = exact)
if (!is.null(object2)) {
data2 <- MCMCchains(object2, params = params, excl = excl, ISB = ISB, exact = exact)
if (!identical(colnames(data), colnames(data2))) {
stop('`object2` must have identical column names to `object`.')
}
}
####################################################################################################
# user and non-user defined plotting parameters
if (missing(pos_tick)) pos_tick = NULL
LA <- 0.4 #amount to lighten colors when ref_ovl = TRUE
ref_col <- 'gray60' # color used for 0 line
guide_col <- 'gray80'
PL_SC <- 0.3 # how much whitespace flanks plotted estimates
if (length(ci) == 2 & typeof(ci) == 'double' & ci[1] <= ci[2] & ci[1] > 0 &
ci[1] <= 100 & ci[2] > 0 & ci[2] <= 100) {
thin <- ci[2] # CI for thin line
thick <- ci[1] # CI for thick line
} else {
stop("ci must be a numeric vector of length 2 where each element is (0, 100] and ci[1] <= ci[2]")
}
# Colors
# number of parameters
np <- NCOL(data)
if (np > 1) {
if (length(col) == 1) {
COL <- rep(col, np)
} else {
if (length(col) != np) {
stop('Number of specified colors must equal number of plotted parameters (or one).')
} else {
COL <- col
}
} # closes if (length(col) == 1)
if (!is.null(object2)) {
if (length(col2) == 1) {
COL2 <- rep(col2, np)
} else {
if (length(col2) != np) {
stop('Number of specified colors must equal number of plotted parameters (or one).')
} else {
COL2 <- col2
}
} # closes if (length(col2) == 1)
} # closes if (!is.null(object2))
} # closes f (np > 1)
if (np == 1) {
if (length(col) == 1) {
COL <- col
} else {
stop('Number of specified colors must equal number of plotted parameters (or one).')
}
if (!is.null(object2)) {
if (length(col2) == 1) {
COL2 <- col2
} else {
stop('Number of specified colors must equal number of plotted parameters (or one).')
}
} # closes if (!is.null(object2))
} # closes if (np == 1)
####################################################################################################
# Process data
pro_fun <- function(input, ...) {
chains <- as.data.frame(input)
tsrt <- apply(chains, 2, stats::median)
if (np > 1) {
len <- np
if (rank == TRUE) idx <- order(tsrt, decreasing = TRUE)
if (rank == FALSE) idx <- len:1
}
if (np == 1) {
len <- 1
idx <- 1
}
medians <- tsrt[idx]
if (HPD == FALSE) {
thick_ci <- c((100 - ((100 - thick) / 2)), ((100 - thick) / 2)) * 0.01
thin_ci <- c((100 - ((100 - thin) / 2)), ((100 - thin) / 2)) * 0.01
if (np == 1)
{
thick_q <- as.matrix(apply(chains, 2, function(x) stats::quantile(x, probs = thick_ci, na.rm = TRUE))[, idx], nrow = 2)
thin_q <- as.matrix(apply(chains, 2, function(x) stats::quantile(x, probs = thin_ci, na.rm = TRUE))[, idx], nrow = 2)
} else {
thick_q <- apply(chains, 2, function(x) stats::quantile(x, probs = thick_ci, na.rm = TRUE))[, idx]
thin_q <- apply(chains, 2, function(x) stats::quantile(x, probs = thin_ci, na.rm = TRUE))[, idx]
}
} else {
if (np == 1)
{
thick_q <- as.matrix(coda::HPDinterval(coda::as.mcmc(chains), prob = ci[1] / 100)[idx, ])
thin_q <- as.matrix(coda::HPDinterval(coda::as.mcmc(chains), prob = ci[2] / 100)[idx, ])
}
else {
thick_q <- t(coda::HPDinterval(coda::as.mcmc(chains), prob = ci[1] / 100)[idx, ])
thin_q <- t(coda::HPDinterval(coda::as.mcmc(chains), prob = ci[2] / 100)[idx, ])
}
}
return(list(len, idx, thick_q, thin_q, medians))
} # closes pro_fun
pf_out <- pro_fun(input = data)
# assign list elements to objects
len <- pf_out[[1]]
idx <- pf_out[[2]]
thick_q <- pf_out[[3]]
thin_q <- pf_out[[4]]
medians <- pf_out[[5]]
if (!is.null(object2)) {
pf_out2 <- pro_fun(input = data2)
thick_q2 <- pf_out2[[3]]
thin_q2 <- pf_out2[[4]]
medians2 <- pf_out2[[5]]
} else {
thin_q2 <- NULL
}
####################################################################################################
# Plot prep
# reference line
if(!is.null(ref)) {
marker <- ref
} else {
marker <- 0
}
# position for coloring with regards to reference line for credible intervals
if (ref_ovl == TRUE) {
# Determine which params have CI that overlap 0 (or ref line if specified)
# kept object names to avoid changing variable names - SHOULD CHANGE
black_cl <- c() #95% CI (default) does not overlap 0
gray_cl <- c() #50% CI (default) does not overlap 0
white_cl <- c() #Both 50% and 95% CI overlap 0
for (i in 1:len) {
if ((thin_q[1, i] > marker & thin_q[2, i] > marker) |
(thin_q[1, i] < marker & thin_q[2, i] < marker)) {
black_cl <- c(black_cl, i)
} else {
if ((thick_q[1, i] > marker & thick_q[2, i] > marker) |
(thick_q[1, i] < marker & thick_q[2, i] < marker)) {
gray_cl <- c(gray_cl, i)
} else {
white_cl <- c(white_cl, i)
}
}
} # closes for (i in 1:len)
if (horiz == FALSE) {
v_black_cl <- rev(black_cl)
v_gray_cl <- rev(gray_cl)
v_white_cl <- rev(white_cl)
}
if (!is.null(object2)) {
black_cl2 <- c()
gray_cl2 <- c()
white_cl2 <- c()
for (i in 1:len) {
if ((thin_q2[1, i] > marker & thin_q2[2, i] > marker) |
(thin_q2[1, i] < marker & thin_q2[2, i] < marker)) {
black_cl2 <- c(black_cl2, i)
} else {
if ((thick_q2[1, i] > marker & thick_q2[2, i] > marker) |
(thick_q2[1, i] < marker & thick_q2[2, i] < marker)) {
gray_cl2 <- c(gray_cl2, i)
} else {
white_cl2 <- c(white_cl2, i)
}
}
} # closes for (i in 1:len)
if (horiz == FALSE) {
v_black_cl2 <- rev(black_cl2)
v_gray_cl2 <- rev(gray_cl2)
v_white_cl2 <- rev(white_cl2)
}
} # closes (!is.null(object2))
} # closes (ref_ovl == TRUE)
####################################################################################################
# Horizontal plot - CI lines parallel to x-axis
if (horiz == TRUE) {
if (missing(xlim)) {
thin_comb <- c(thin_q, thin_q2)
rn <- diff(range(thin_comb, na.rm = TRUE)) * PL_SC
XLIM <- c((min(thin_comb, na.rm = TRUE) - rn), (max(thin_comb, na.rm = TRUE) + rn))
} else {
XLIM <- xlim
}
if (len > 20) {
tt <- 1.5
} else {
alpha <- 0.667
beta <- 0.0412
tt <- alpha + beta * len
} # closes if (len > 20)
YLIM <- c((1 - tt), (len + tt))
if (missing(xlab)) xlab <- 'Parameter Estimate'
if (is.null(xlab)) xlab <- ''
if (missing(main)) main <- ''
if (missing(labels)) {
labs <- names(medians)
} else {
if (!missing(labels)) {
if (is.null(labels)) {
labs <- rep('', len)
}
if (!is.null(labels)) {
if (length(labels) == len) {
labs <- labels[idx]
} else {
stop('Labels length not equal to number of parameters')
}
}
}
} # closes if (missing(labels))
# 0.2 inches per line - mar measured in lines
m_char <- max(sapply(labs, function(x){graphics::strwidth(x, cex = sz_labels, units = 'in')})) / 0.2
graphics::par(mar = c(mar[1], (m_char + (mar[2] - 3)), mar[3], mar[4] - 1))
#plot blank plot
graphics::plot(
medians,
(1:len),
xlim = XLIM,
ylim = YLIM,
type = "n",
ann = TRUE,
xaxt = 'n',
yaxt = "n",
bty = "n",
ylab = NA,
xlab = xlab,
cex.lab = sz_ax_txt,
yaxs = 'i') #cex.lab is axis label
# title
graphics::title(main, cex.main = sz_main_txt)
if (guide_axis == TRUE) {
# top axis params
graphics::axis(3, lwd.ticks = sz_ax, labels = FALSE, at = pos_tick, lwd = sz_ax)
}
# bottom axis params
graphics::axis(1, lwd.ticks = sz_ax, labels = TRUE, at = pos_tick, lwd = sz_ax, cex.axis = sz_tick_txt)
# left axis params (labels)
# las - 0 parallel to axis, 1 horiz, 2 perp to axis, 3 vert
graphics::axis(
2,
at = ((1:len)),
tick = FALSE,
labels = labs,
las = 1,
adj = 0,
line = -1,
cex.axis = sz_labels)
if (guide_lines == TRUE) {
# limits
# determine where ticks were placed - from axTicks source code:
xaxp <- graphics::par("xaxp")
XLIM2 <- c(xaxp[1], xaxp[2])
xs <- matrix(rep(XLIM2, len), nrow = 2)
ys <- rbind(1:len, 1:len)
graphics::matlines(xs, ys, lty = 1, col = guide_col)
}
# ref line
if(!is.null(ref)) graphics::abline(v=ref, lty = 2, lwd = 3, col = ref_col)
if (ref_ovl == TRUE) {
# positions to plot CI
blk_bnd <- rbind(black_cl, black_cl)
gry_bnd <- rbind(gray_cl, gray_cl)
wht_bnd <- rbind(white_cl, white_cl)
#lightened colors
lcol <- colorspace::lighten(COL, LA)
bl_col <- COL[black_cl]
gr_col <- lcol[gray_cl]
wh_col <- lcol[white_cl]
#if any colors are 'black', change them to 'grey60' (for consistency with older pkg versions)
bl_idx <- which(COL == 'black')
gr_col[which(gray_cl %in% bl_idx)] <- 'grey60'
wh_col[which(white_cl %in% bl_idx)] <- 'grey60'
# Black CI
if (!is.null(black_cl)) {
#only one object
if (is.null(object2)) {
# Thick
graphics::matlines(
thick_q[,black_cl],
blk_bnd,
type = 'l',
lty = 1,
lwd = sz_thick,
col = bl_col)
# Thin
graphics::matlines(
thin_q[,black_cl],
blk_bnd,
type = 'l',
lty = 1,
lwd = sz_thin,
col = bl_col)
} else {
#two objects
blk_bnd2 <- rbind(black_cl2, black_cl2)
bl_col2 <- COL2[black_cl2]
# object
graphics::matlines(
thick_q[,black_cl],
(blk_bnd + offset),
type = 'l',
lty = 1,
lwd = sz_thick,
col = bl_col)
graphics::matlines(
thin_q[,black_cl],
(blk_bnd + offset),
type = 'l',
lty = 1,
lwd = sz_thin,
col = bl_col)
# object2
graphics::matlines(
thick_q2[, black_cl2],
(blk_bnd2 - offset),
type = 'l',
lty = 1,
lwd = sz_thick,
col = bl_col2)
graphics::matlines(
thin_q2[, black_cl2],
(blk_bnd2 - offset),
type = 'l',
lty = 1,
lwd = sz_thin,
col = bl_col2)
} # closes if (is.null(object2))
} # closes if (!is.null(black_cl))
# Gray CI
if (!is.null(gray_cl)) {
#only one object
if (is.null(object2)) {
# Thick
graphics::matlines(
thick_q[,gray_cl],
gry_bnd,
type = 'l',
lty = 1,
lwd = sz_thick,
col = gr_col)
# Thin
graphics::matlines(
thin_q[,gray_cl],
gry_bnd,
type = 'l',
lty = 1,
lwd = sz_thin,
col = gr_col)
} else {
#two objects
gry_bnd2 <- rbind(gray_cl2, gray_cl2)
#lightened colors for obj 2
lcol2 <- colorspace::lighten(COL2, LA)
gr_col2 <- lcol2[gray_cl2]
#if any colors are 'black', change them to 'grey60' (for consistency with older pkg versions)
bl_idx2 <- which(COL2 == 'black')
gr_col2[which(gray_cl2 %in% bl_idx2)] <- 'grey60'
# object
graphics::matlines(
thick_q[,gray_cl],
(gry_bnd + offset),
type = 'l',
lty = 1,
lwd = sz_thick,
col = gr_col)
graphics::matlines(
thin_q[,gray_cl],
(gry_bnd + offset),
type = 'l',
lty = 1,
lwd = sz_thin,
col = gr_col)
#object2
graphics::matlines(
thick_q2[,gray_cl2],
(gry_bnd2 - offset),
type = 'l',
lty = 1,
lwd = sz_thick,
col = gr_col2)
graphics::matlines(
thin_q2[,gray_cl2],
(gry_bnd2 - offset),
type = 'l',
lty = 1,
lwd = sz_thin,
col = gr_col2)
} # closes if (is.null(object2))
} # closes if (!is.null(gray_cl))
# White CI
if (!is.null(white_cl)) {
#only one object
if (is.null(object2)) {
graphics::matlines(
thick_q[,white_cl],
wht_bnd,
type = 'l',
lty = 1,
lwd = sz_thick,
col = wh_col) # white (gray)
graphics::matlines(
thin_q[,white_cl],
wht_bnd,
type = 'l',
lty = 1,
lwd = sz_thin,
col = wh_col) # white (gray)
} else {
#two objects
wht_bnd2 <- rbind(white_cl2, white_cl2)
#lightened colors for obj 2
lcol2 <- colorspace::lighten(COL2, LA)
wh_col2 <- lcol2[white_cl2]
#if any colors are 'black', change them to 'grey60' (for consistency with older pkg versions)
bl_idx2 <- which(COL2 == 'black')
wh_col2[which(white_cl2 %in% bl_idx2)] <- 'grey60'
# object
graphics::matlines(
thick_q[,white_cl],
(wht_bnd + offset),
type = 'l',
lty = 1,
lwd = sz_thick,
col = wh_col) # white (gray)
graphics::matlines(
thin_q[,white_cl],
(wht_bnd + offset),
type = 'l',
lty = 1,
lwd = sz_thin,
col = wh_col) # white (gray)
# object2
graphics::matlines(
thick_q2[,white_cl2],
(wht_bnd2 - offset),
type = 'l',
lty = 1,
lwd = sz_thick,
col = wh_col2) # white (gray)
graphics::matlines(
thin_q2[,white_cl2],
(wht_bnd2 - offset),
type = 'l',
lty = 1,
lwd = sz_thin,
col = wh_col2) #white (gray)
} # closes if (is.null(object2))
} # closes if (!is.null(white_cl))
# Medians
if (is.null(object2)) {
#only one object
graphics::points(medians, 1:len, pch = 16, col = 'white', cex = sz_med)
graphics::points(medians[black_cl], black_cl, pch = 16, col = bl_col, cex = sz_med)
graphics::points(medians[gray_cl], gray_cl, pch = 16, col = gr_col, cex = sz_med)
graphics::points(medians[white_cl], white_cl, pch = 21, col = wh_col, cex = sz_med, lwd = 2)
} else {
# two objects
graphics::points(medians, (1:len + offset), pch = 16, col = 'white', cex = sz_med)
graphics::points(medians[black_cl], (black_cl + offset), pch = 16, col = bl_col, cex = sz_med)
graphics::points(medians[gray_cl], (gray_cl + offset), pch = 16, col = gr_col, cex = sz_med)
graphics::points(medians[white_cl], (white_cl + offset), pch = 21, col = wh_col, cex = sz_med, lwd = 2)
# object2
graphics::points(medians2, (1:len - offset), pch = 16, col = 'white', cex = sz_med)
graphics::points(medians2[black_cl2], (black_cl2 - offset), pch = 16, col = bl_col2, cex = sz_med)
graphics::points(medians2[gray_cl2], (gray_cl2 - offset), pch = 16, col = gr_col2, cex = sz_med)
graphics::points(medians2[white_cl2], (white_cl2 - offset), pch = 21, col = wh_col2, cex = sz_med, lwd = 2)
}
} else {
if (is.null(object2))
{
graphics::matlines(thick_q[,1:len], rbind(1:len, 1:len), type = 'l', lty = 1, lwd = sz_thick, col = COL[idx])
graphics::matlines(thin_q[,1:len], rbind(1:len, 1:len), type = 'l', lty = 1, lwd = sz_thin, col = COL[idx])
graphics::points(medians[1:len], 1:len, pch = 16, col = COL[idx], cex = sz_med)
} else {
# object
graphics::matlines(thick_q[,1:len], (rbind(1:len, 1:len) + offset), type = 'l',
lty = 1, lwd = sz_thick, col = COL[idx])
graphics::matlines(thin_q[,1:len], (rbind(1:len, 1:len) + offset), type = 'l',
lty = 1, lwd = sz_thin, col = COL[idx])
graphics::points(medians[1:len], (1:len + offset), pch = 16, col = COL[idx], cex = sz_med)
# object2
graphics::matlines(thick_q2[,1:len], (rbind(1:len, 1:len) - offset), type = 'l',
lty = 1, lwd = sz_thick, col = COL2[idx])
graphics::matlines(thin_q2[,1:len], (rbind(1:len, 1:len) - offset), type = 'l',
lty = 1, lwd = sz_thin, col = COL2[idx])
graphics::points(medians2[1:len], (1:len - offset), pch = 16, col = COL2[idx], cex = sz_med)
}
}
} # closes if (horiz == TRUE)
####################################################################################################
# Vertical plot - CI lines parallel to x-axis
if (horiz == FALSE) {
if (missing(ylim)) {
thin_comb <- c(thin_q, thin_q2)
rn <- diff(range(thin_comb, na.rm = TRUE)) * PL_SC
YLIM <- c((min(thin_comb, na.rm = TRUE) - rn), (max(thin_comb, na.rm = TRUE) + rn))
} else {
YLIM <- ylim
}
if (len > 20) {
tt <- 1.5
} else {
alpha <- 0.667
beta <- 0.0412
tt <- alpha + beta * len
}
XLIM = c((1 - tt), (len + tt))
if (missing(ylab)) ylab = 'Parameter Estimate'
if (is.null(ylab)) ylab = ''
if (missing(main)) main = ''
if (missing(labels)) {
labs = names(medians)
} else {
if (!missing(labels)) {
if (is.null(labels)) {
labs <- rep('', len)
}
if (!is.null(labels)) {
if (length(labels) == len) {
labs <- labels[idx]
} else {
stop('Labels length not equal to number of parameters')
}
}
}
} # closes if (missing(labels))
# 0.2 inches per line - mar measured in lines
m_char <- (max(sapply(labs, function(x){graphics::strwidth(x, cex = sz_labels, units = 'i')})) / 0.2)
# blank plot - do not display
grDevices::pdf(file = NULL)
graphics::plot(
(len:1),
medians,
xlim = XLIM,
ylim = YLIM,
type = "n",
ann = TRUE,
xaxt = 'n',
yaxt = "n",
bty = "n",
ylab = NA,
xlab = NA,
cex.lab = sz_ax_txt,
xaxs = 'i')
# create invisible ticks to determine where to put y-axis label
tickp <- graphics::axis(
2,
lwd.ticks = sz_ax,
labels = FALSE,
at = pos_tick,
lwd = sz_ax,
cex.axis = sz_tick_txt,
las = 1,
col = 'white',
col.ticks = 'white')
invisible(grDevices::dev.off())
# determine how long labels are
ml_tickp <- max(graphics::strwidth(tickp, cex = sz_tick_txt, units = 'in'))
# 5 lines/inch
ll <- 1.8 + 5 * ml_tickp
# set plot margins according to labels
graphics::par(mar = c((m_char + (mar[1] - 4)), mar[2] + 1 + ll - 3, mar[3] - 1.5, mar[4]))
# new blank plot
graphics::plot(
(len:1),
medians,
xlim = XLIM,
ylim = YLIM,
type = "n",
ann = TRUE,
xaxt = 'n',
yaxt = "n",
bty = "n",
ylab = NA,
xlab = NA,
cex.lab = sz_ax_txt,
xaxs = 'i')
# ticks
graphics::axis(
2,
lwd.ticks = sz_ax,
labels = TRUE,
at = pos_tick,
lwd = sz_ax,
cex.axis = sz_tick_txt,
las = 1)
# y-axis label
graphics::title(ylab = ylab, cex.lab = sz_ax_txt, line = ll)
# title
graphics::title(main, cex.main = sz_main_txt, line = 0.5)
if (guide_axis == TRUE) {
# right axis params
graphics::axis(4, lwd.ticks = sz_ax, labels = FALSE, at = pos_tick, lwd = sz_ax)
}
# bottom axis params (labels)
# las - 0 parallel to axis, 1 horizontal, 2 perpendicular to axis, 3 vertical
graphics::axis(
1,
at = (len:1) + 0.013,
tick = FALSE,
labels = labs,
las = 2,
adj = 0,
line = -1,
cex.axis = sz_labels)
if (guide_lines == TRUE) {
# limits
yaxp <- graphics::par("yaxp")
YLIM2 <- c(yaxp[1], yaxp[2])
ys <- matrix(rep(YLIM2, len), nrow = 2)
xs <- rbind(1:len, 1:len)
graphics::matlines(xs, ys, lty = 1, col = guide_col)
}
# ref line
if (!is.null(ref)) {
graphics::abline(h=ref, lty = 2, lwd = 3, col = ref_col)
}
if (ref_ovl == TRUE) {
# positions to plot CI
v_blk_bnd <- matrix(rep(rev(len + 1 - black_cl), 2), nrow = 2, byrow = TRUE)
v_gry_bnd <- matrix(rep(rev(len + 1 - gray_cl), 2), nrow = 2, byrow = TRUE)
v_wht_bnd <- matrix(rep(rev(len + 1 - white_cl), 2), nrow = 2, byrow = TRUE)
# Black CI
if (!is.null(black_cl)) {
if (is.null(object2)) {
# Thick
graphics::matlines(v_blk_bnd, thick_q[,v_black_cl], type = 'l', lty = 1, lwd = sz_thick, col = 'black')
# Thin
graphics::matlines(v_blk_bnd, thin_q[,v_black_cl], type = 'l', lty = 1, lwd = sz_thin, col = 'black')
} else {
v_blk_bnd2 <- matrix(rep(rev(len + 1 - black_cl2), 2), nrow = 2, byrow = TRUE)
# object
graphics::matlines((v_blk_bnd - offset), thick_q[,v_black_cl], type = 'l', lty = 1, lwd = sz_thick, col = 'black')
graphics::matlines((v_blk_bnd - offset), thin_q[,v_black_cl], type = 'l', lty = 1, lwd = sz_thin, col = 'black')
# object2
graphics::matlines((v_blk_bnd2 + offset), thick_q2[,v_black_cl2], type = 'l', lty = 1, lwd = sz_thick, col = 'black')
graphics::matlines((v_blk_bnd2 + offset), thin_q2[,v_black_cl2], type = 'l', lty = 1, lwd = sz_thin, col = 'black')
} # closes if (is.null(object2))
} # closes if (!is.null(black_cl))
# Gray CI
if (!is.null(gray_cl)) {
if (is.null(object2)) {
# Thick
graphics::matlines(v_gry_bnd, thick_q[,v_gray_cl], type = 'l', lty = 1, lwd = sz_thick, col = gr_col)
# Thin
graphics::matlines(v_gry_bnd, thin_q[,v_gray_cl], type = 'l', lty = 1, lwd = sz_thin, col = gr_col)
} else {
v_gry_bnd2 <- matrix(rep(rev(len + 1 - gray_cl2), 2), nrow = 2, byrow = TRUE)
# object
graphics::matlines((v_gry_bnd - offset), thick_q[,v_gray_cl], type = 'l', lty = 1, lwd = sz_thick, col = gr_col)
graphics::matlines((v_gry_bnd - offset), thin_q[,v_gray_cl], type = 'l', lty = 1, lwd = sz_thin, col = gr_col)
#object2
graphics::matlines((v_gry_bnd2 + offset), thick_q2[,v_gray_cl2], type = 'l', lty = 1, lwd = sz_thick, col = gr_col)
graphics::matlines((v_gry_bnd2 + offset), thin_q2[,v_gray_cl2], type = 'l', lty = 1, lwd = sz_thin, col = gr_col)
}
}
# White CI
if (!is.null(white_cl)) {
if (is.null(object2)) {
graphics::matlines(v_wht_bnd, thick_q[,v_white_cl], type = 'l', lty = 1, lwd = sz_thick, col = gr_col) #white (gray)
graphics::matlines(v_wht_bnd, thin_q[,v_white_cl], type = 'l', lty = 1, lwd = sz_thin, col = gr_col) #white (gray)
} else {
v_wht_bnd2 <- matrix(rep(rev(len + 1 - white_cl2), 2), nrow = 2, byrow = TRUE)
# object
graphics::matlines((v_wht_bnd - offset), thick_q[,v_white_cl], type = 'l', lty = 1, lwd = sz_thick, col = gr_col)
graphics::matlines((v_wht_bnd - offset), thin_q[,v_white_cl], type = 'l', lty = 1, lwd = sz_thin, col = gr_col)
# object2
graphics::matlines((v_wht_bnd2 + offset), thick_q2[,v_white_cl2], type = 'l', lty = 1, lwd = sz_thick, col = gr_col)
graphics::matlines((v_wht_bnd2 + offset), thin_q2[,v_white_cl2], type = 'l', lty = 1, lwd = sz_thin, col = gr_col)
}
}
# Medians
if (is.null(object2)) {
graphics::points(len:1, medians, pch = 16, col = 'white', cex = sz_med)
graphics::points(v_blk_bnd[1,], medians[v_black_cl], pch = 16, col = 'black', cex = sz_med)
graphics::points(v_gry_bnd[1,], medians[v_gray_cl], pch = 16, col = gr_col, cex = sz_med)
graphics::points(v_wht_bnd[1,], medians[v_white_cl], pch = 21, col = gr_col, cex = sz_med, lwd = 2)
} else {
# object
graphics::points((len:1 - offset), medians, pch = 16, col = 'white', cex = sz_med)
graphics::points((v_blk_bnd[1,] - offset), medians[v_black_cl], pch = 16, col = 'black', cex = sz_med)
graphics::points((v_gry_bnd[1,] - offset), medians[v_gray_cl], pch = 16, col = gr_col, cex = sz_med)
graphics::points((v_wht_bnd[1,] - offset), medians[v_white_cl], pch = 21, col = gr_col, cex = sz_med, lwd = 2)
# object2
graphics::points((len:1 + offset), medians2, pch = 16, col = 'white', cex = sz_med)
graphics::points((v_blk_bnd2[1,] + offset), medians2[v_black_cl2], pch = 16, col = 'black', cex = sz_med)
graphics::points((v_gry_bnd2[1,] + offset), medians2[v_gray_cl2], pch = 16, col = gr_col, cex = sz_med)
graphics::points((v_wht_bnd2[1,] + offset), medians2[v_white_cl2], pch = 21, col = gr_col, cex = sz_med, lwd = 2)
}
} else {
if (is.null(object2)) {
graphics::matlines(rbind(1:len, 1:len), thick_q[,len:1], type = 'l', lty = 1, lwd = sz_thick, col = COL[idx])
graphics::matlines(rbind(1:len, 1:len), thin_q[,len:1], type = 'l', lty = 1, lwd = sz_thin, col = COL[idx])
graphics::points(1:len, medians[len:1], pch = 16, col = COL[idx], cex = sz_med)
} else {
# object
graphics::matlines((rbind(1:len, 1:len) - offset), thick_q[,len:1], type = 'l', lty = 1, lwd = sz_thick, col = COL[idx])
graphics::matlines((rbind(1:len, 1:len) - offset), thin_q[,len:1], type = 'l', lty = 1, lwd = sz_thin, col = COL[idx])
graphics::points((1:len - offset), medians[len:1], pch = 16, col = COL[idx], cex = sz_med)
# object2
graphics::matlines((rbind(1:len, 1:len) + offset), thick_q2[,len:1], type = 'l', lty = 1, lwd = sz_thick, col = COL2[idx])
graphics::matlines((rbind(1:len, 1:len) + offset), thin_q2[,len:1], type = 'l', lty = 1, lwd = sz_thin, col = COL2[idx])
graphics::points((1:len + offset), medians2[len:1], pch = 16, col = COL2[idx], cex = sz_med)
} # closes if (is.null(object2))
} # closes if (ref_ovl == TRUE)
} # closes if (horiz == FALSE)
# restore graphics
graphics::par(mar = c(5, 4, 4, 2) + 0.1)
}
caseyyoungflesh/MCMCvis documentation built on March 18, 2024, 9:39 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions MCMCplot
Documented in MCMCplot
#' Caterpillar plots of posterior distributions from MCMC output
#'
#' Visualize posterior distributions from MCMC output for specific parameters of interest using caterpillar plots. Color of median dot represents the overlap of the posterior distribution with 0 (or other specified value).
#'
#'
#' @param object Object containing MCMC output. See DETAILS below.
#'
#' @param object2 Optional second object containing MCMC output. If specified, parameter estimates from each model will be displayed in a paired manner. Parameter names for \code{'object'} and \code{'object2'} must be identical. See DETAILS below.
#'
#' @param params Character string (or vector of character strings) denoting parameters to be plotted.
#'
#' Default \code{'all'} plots posteriors for all parameters. See VALUE below.
#'
#' @param HPD Logical specifying whether to calculate equal-tailed credible intervals (\code{HPD = FALSE}) or highest posterior density intervals (\code{HPD = TRUE}) for the selected parameters. Default is \code{HPD = FALSE}.
#'
#' @param ci Numeric vector of length 2, where each element is (0,100] and represents the width of an equal-tailed (\code{HPD = FALSE}) or highest posterior density (\code{HPD = TRUE}) credible interval. The first element of this vector corresponds to the thicker (narrower) credible interval displayed on the plot (default is 0.5) and the second element of this vector corresponds to the thinner (wider) credible interval (default is 0.95). The first credible interval width (\code{ci[1]}) must be less than or equal to the width of the second credible interval (\code{ci[2]}).
#'
#' @param excl Character string (or vector of character strings) denoting parameters to exclude. Used in conjunction with \code{params} argument to select parameters of interest.
#'
#' @param ISB Ignore Square Brackets (ISB). Logical specifying whether square brackets should be ignored in the \code{params} and \code{excl} arguments. If \code{TRUE}, square brackets are ignored. If \code{FALSE}, square brackets are not ignored. This allows partial names to be used when specifying parameters of interest. Use \code{exact} argument to specify whether input from \code{params} and \code{excl} arguments should be matched exactly.
#'
#' @param exact Logical specifying whether input from \code{params} and \code{excl} arguments should be matched exactly (after ignoring square brackets if \code{ISB = FALSE}). If \code{TRUE}, input from \code{params} and \code{excl} are matched exactly (after taking \code{ISB} argument into account). If \code{FALSE}, input from \code{params} and \code{excl} are matched using regular expression format (after taking \code{ISB} argument into account).
#'
#' @param ref Value indicating where vertical reference line should be created and what value to use a reference for caterpillar median coloration.
#'
#' Default is \code{ref = 0}.
#'
#' Argument \code{NULL} will plot no reference line.
#'
#' @param ref_ovl Logical specifying whether the style/color of plotted median dots and CI should be changed based on whether the specified credible intervals (50 \% and 95 \% by default) overlap the reference line. See DETAILS for more information.
#'
#' @param col Character string (or vector of character strings) specifying which color to render estimates on plot. When \code{ref_ovl = TRUE}, this argument has no effect and colors plotted will be based on the credible intervals and reference line. Number of specified colors must equal the number of specified parameters or one.
#'
#' @param col2 Character string (or vector of character strings) specifying which color to render estimates on plot for \code{object2} (if specified). Number of specified colors must equal the number of specified parameters or one. Red by default.
#'
#' @param offset Value indicating how much to offset plotted posteriors when \code{object2} is specified (i.e., control the amount of space between the two plotted posteriors for each parameter). The distance from one set of parameters to another corresponds to a value of 1.
#'
#' @param rank Logical specifying whether output should be ranked. If \code{TRUE} posteriors will be ranked in decreasing order (based on specified measure of centrality) from top down.
#'
#' @param horiz Logical specifying orientation of plot. If \code{TRUE} posteriors will be plotted running horizontally (parallel to the x-axis). If \code{FALSE} posteriors will be plotted running vertically (perpendicular to the x-axis).
#'
#' @param xlim Numerical vector of length 2, indicating range of x-axis. Only applicable if \code{horiz = TRUE}.
#'
#' @param ylim Numerical vector of length 2, indicating range of y-axis. Only applicable if \code{horiz = FALSE}.
#'
#' @param xlab Character string labeling x-axis. Only applicable if \code{horiz = TRUE}.
#'
#' Default label is 'Parameter Estimate'. Option \code{NULL} will return plot with no label on x-axis.
#' @param ylab Character string labeling y-axis. Only applicable if \code{horiz = FALSE}.
#'
#' Default label is 'Parameter Estimate'. Option \code{NULL} will return plot with no label on y-axis.
#' @param main Character string indicating title of plot.
#'
#' @param labels Character string (or vector of character strings if plotting > 1 parameter) labeling parameter estimates along y-axis (if \code{horiz = FALSE}) or x-axis (if \code{horiz = TRUE}).
#'
#' Default option will use parameter names from \code{object}.
#'
#' Option \code{NULL} will return plot with no labels on axis.
#'
#' @param guide_lines Logical specifying whether to plot reference lines for each parameter in order to better visualize which parameter names correspond to each posterior.
#'
#' @param guide_axis Logical specifying whether a second axis should be plotted (x-axis if \code{HORIZ = TRUE}, y-axis if \code{HORIZ = FALSE}) to help interpret values on plot.
#'
#' @param sz_labels Number specifying size of text for parameter labels on axis.
#'
#' @param sz_med Number specifying size of points represents posterior medians.
#'
#' @param sz_thick Number specifying thickness of 50 percent CI line (thicker line).
#'
#' @param sz_thin Number specifying thickness of 95 percent CI line (thinner line).
#'
#' @param sz_ax Number specifying thickness of axis and ticks.
#'
#' @param sz_ax_txt Number specifying size of text for axis label.
#'
#' @param sz_tick_txt Number specifying size of text for tick labels on axis.
#'
#' @param sz_main_txt Number specifying size of text for main title.
#'
#' @param pos_tick Numeric vector specifying where ticks on axis should be placed.
#'
#' @param mar Numerical vector of length 4 specifying plot margins - (BOTTOM, LEFT, TOP, RIGHT). Changes to the margin should be made within the function rather than using the \code{par} call.
#'
#' Default is c(5.1, 4.1, 4.1, 2.1) - the R plot default.
#'
#' @section Details:
#' Points represent posterior medians. Parameters where the smaller specified credible intervals (first value in vector provided to \code{ci} argument; 50\% by default) overlap 0 (or other value specified by the \code{ref} argument) are indicated by 'open' circles and a lighter color than that specified (e.g., gray will be displayed with the default black is specified for \code{col}). Parameters where the specific smaller credible intervals DO NOT overlap 0 (or other specified value) AND the larger specified credible intervals (second value in vector provided to \code{ci} argument; 95\% by default) DO overlap 0 (or other specified value) are indicated by 'closed' circles and a lighter color that that specified. Parameters where the larger specified credible intervals DO NOT overlap 0 (or other specified value) are indicated by 'closed' circles and the color specified (black by default). Thick lines represent the smaller specified credible intervals percent credible intervals (50\% by default) while thin lines represent the larger specified credible intervals (95 \% by default). \code{ref_ovl = TRUE} can be used to enable this feature. When two model objects are supplied to the function (with \code{object} and \code{object2}) and no argument is supplied to \code{col} or \code{col2}, light red lines can be interpreted as analogous to gray lines. When \code{ref_ovl = TRUE} and a color (or colors) other than the default are specified (with the \code{col} and \code{col2} arguments) lighter versions of the color specified are used in black of the light gray and/or light red lines.
#'
#' When \code{object2} is specified, paired caterpillar plots of each parameter are produced. For this reason, parameter names of \code{object} and \code{object2} specified with the \code{params} argument must be identical (to be used for comparing posterior estimates of similar models). \code{col} and \code{col2} arguments can be specified to change the color of output from \code{object} and \code{object2}, respectively. By default, output from \code{object} is plotted in black and \code{object2} is plotted in red. The \code{ref_ovl} argument can also be specified.
#'
#' \code{object} argument can be a \code{stanfit} object (\code{rstan} package), a \code{CmdStanMCMC} object (\code{cmdstanr} package), a \code{stanreg} object (\code{rstanarm} package), a \code{brmsfit} object (\code{brms} package), an \code{mcmc.list} object (\code{coda} and \code{rjags} packages), \code{mcmc} object (\code{coda} and \code{nimble} packages), \code{list} object (\code{nimble} package), an \code{R2jags} model object (\code{R2jags} package), a \code{jagsUI} model object (\code{jagsUI} package), or a matrix containing MCMC chains (each column representing MCMC output for a single parameter, rows representing iterations in the chain). The function automatically detects the object type and proceeds accordingly.
#'
#' @section Notes:
#'
#' When specifying \code{rank = TRUE} and specifying labels for \code{labels}, labels will be applied to parameters before they are ranked.
#'
#' Thanks to Cinner et al. 2016, whose Fig. 1 inspired this plot.
#'
#' @section References:
#'
#' Cinner, J. E., C. Huchery, M. A. MacNeil, N. A. J. Graham, T. R. McClanahan, J. Maina, E. Maire, J. N. Kittinger, C. C. Hicks, C. Mora, E. H. Allison, S. D'Agata, A. Hoey, D. A. Feary, L. Crowder, I. D. Williams, M. Kulbicki, L. Vigliola, L. Wantiez, G. Edgar, R. D. Stuart-Smith, S. A. Sandin, A. L. Green, M. J. Hardt, M. Beger, A. Friedlander, S. J. Campbell, K. E. Holmes, S. K. Wilson, E. Brokovich, A. J. Brooks, J. J. Cruz-Motta, D. J. Booth, P. Chabanet, C. Gough, M. Tupper, S. C. A. Ferse, U. R. Sumaila, and D. Mouillot. 2016. Bright spots among the world's coral reefs. Nature 535:416-419.
#'
#' @examples
#' #Load data
#' data(MCMC_data)
#'
#' #Plot MCMC output
#' MCMCplot(MCMC_data, labels = NULL)
#'
#' #Just 'beta' parameters
#' MCMCplot(MCMC_data, params = 'beta')
#'
#' #Just 'beta' parameters using highest posterior density intervals
#' MCMCplot(MCMC_data, params = 'beta', HPD = TRUE)
#
#' #Just 'beta[1]', 'beta[4]', and 'alpha[3]'
#' MCMCplot(MCMC_data, params = c('beta[1]', 'beta[4]', 'alpha[3]'), ISB = FALSE, exact = TRUE)
#'
#' #Just 'beta[1]', 'beta[4]', and 'alpha[3]' and change the credible interval widths
#' MCMCplot(MCMC_data, ci = c(50, 89), params = c('beta[1]', 'beta[4]', 'alpha[3]'),
#' ISB = FALSE, exact = TRUE)
#'
#' #Rank parameters by posterior mean
#' MCMCplot(MCMC_data, params = 'beta', rank = TRUE)
#'
#' #Create vertical plot
#' MCMCplot(MCMC_data, params = 'beta', horiz = FALSE)
#'
#' @export
#'
MCMCplot <- function(object,
object2 = NULL,
params = 'all',
HPD = FALSE,
ci = c(50, 95),
excl = NULL,
ISB = TRUE,
exact = TRUE,
ref = 0,
ref_ovl = FALSE,
col = 'black',
col2 = 'red',
offset = 0.1,
rank = FALSE,
horiz = TRUE,
xlim,
ylim,
xlab,
ylab,
main,
labels,
guide_lines = FALSE,
guide_axis = TRUE,
sz_labels = 1.2,
sz_med = 1.5,
sz_thick = 5,
sz_thin = 2,
sz_ax = 3,
sz_ax_txt = 1.3,
sz_tick_txt = 1.2,
sz_main_txt = 1.2,
pos_tick,
mar = c(5.1, 4.1, 4.1, 2.1)) {
####################################################################################################
# load object(s)
data <- MCMCchains(object, params = params, excl = excl, ISB = ISB, exact = exact)
if (!is.null(object2)) {
data2 <- MCMCchains(object2, params = params, excl = excl, ISB = ISB, exact = exact)
if (!identical(colnames(data), colnames(data2))) {
stop('`object2` must have identical column names to `object`.')
}
}
####################################################################################################
# user and non-user defined plotting parameters
if (missing(pos_tick)) pos_tick = NULL
LA <- 0.4 #amount to lighten colors when ref_ovl = TRUE
ref_col <- 'gray60' # color used for 0 line
guide_col <- 'gray80'
PL_SC <- 0.3 # how much whitespace flanks plotted estimates
if (length(ci) == 2 & typeof(ci) == 'double' & ci[1] <= ci[2] & ci[1] > 0 &
ci[1] <= 100 & ci[2] > 0 & ci[2] <= 100) {
thin <- ci[2] # CI for thin line
thick <- ci[1] # CI for thick line
} else {
stop("ci must be a numeric vector of length 2 where each element is (0, 100] and ci[1] <= ci[2]")
}
# Colors
# number of parameters
np <- NCOL(data)
if (np > 1) {
if (length(col) == 1) {
COL <- rep(col, np)
} else {
if (length(col) != np) {
stop('Number of specified colors must equal number of plotted parameters (or one).')
} else {
COL <- col
}
} # closes if (length(col) == 1)
if (!is.null(object2)) {
if (length(col2) == 1) {
COL2 <- rep(col2, np)
} else {
if (length(col2) != np) {
stop('Number of specified colors must equal number of plotted parameters (or one).')
} else {
COL2 <- col2
}
} # closes if (length(col2) == 1)
} # closes if (!is.null(object2))
} # closes f (np > 1)
if (np == 1) {
if (length(col) == 1) {
COL <- col
} else {
stop('Number of specified colors must equal number of plotted parameters (or one).')
}
if (!is.null(object2)) {
if (length(col2) == 1) {
COL2 <- col2
} else {
stop('Number of specified colors must equal number of plotted parameters (or one).')
}
} # closes if (!is.null(object2))
} # closes if (np == 1)
####################################################################################################
# Process data
pro_fun <- function(input, ...) {
chains <- as.data.frame(input)
tsrt <- apply(chains, 2, stats::median)
if (np > 1) {
len <- np
if (rank == TRUE) idx <- order(tsrt, decreasing = TRUE)
if (rank == FALSE) idx <- len:1
}
if (np == 1) {
len <- 1
idx <- 1
}
medians <- tsrt[idx]
if (HPD == FALSE) {
thick_ci <- c((100 - ((100 - thick) / 2)), ((100 - thick) / 2)) * 0.01
thin_ci <- c((100 - ((100 - thin) / 2)), ((100 - thin) / 2)) * 0.01
if (np == 1)
{
thick_q <- as.matrix(apply(chains, 2, function(x) stats::quantile(x, probs = thick_ci, na.rm = TRUE))[, idx], nrow = 2)
thin_q <- as.matrix(apply(chains, 2, function(x) stats::quantile(x, probs = thin_ci, na.rm = TRUE))[, idx], nrow = 2)
} else {
thick_q <- apply(chains, 2, function(x) stats::quantile(x, probs = thick_ci, na.rm = TRUE))[, idx]
thin_q <- apply(chains, 2, function(x) stats::quantile(x, probs = thin_ci, na.rm = TRUE))[, idx]
}
} else {
if (np == 1)
{
thick_q <- as.matrix(coda::HPDinterval(coda::as.mcmc(chains), prob = ci[1] / 100)[idx, ])
thin_q <- as.matrix(coda::HPDinterval(coda::as.mcmc(chains), prob = ci[2] / 100)[idx, ])
}
else {
thick_q <- t(coda::HPDinterval(coda::as.mcmc(chains), prob = ci[1] / 100)[idx, ])
thin_q <- t(coda::HPDinterval(coda::as.mcmc(chains), prob = ci[2] / 100)[idx, ])
}
}
return(list(len, idx, thick_q, thin_q, medians))
} # closes pro_fun
pf_out <- pro_fun(input = data)
# assign list elements to objects
len <- pf_out[[1]]
idx <- pf_out[[2]]
thick_q <- pf_out[[3]]
thin_q <- pf_out[[4]]
medians <- pf_out[[5]]
if (!is.null(object2)) {
pf_out2 <- pro_fun(input = data2)
thick_q2 <- pf_out2[[3]]
thin_q2 <- pf_out2[[4]]
medians2 <- pf_out2[[5]]
} else {
thin_q2 <- NULL
}
####################################################################################################
# Plot prep
# reference line
if(!is.null(ref)) {
marker <- ref
} else {
marker <- 0
}
# position for coloring with regards to reference line for credible intervals
if (ref_ovl == TRUE) {
# Determine which params have CI that overlap 0 (or ref line if specified)
# kept object names to avoid changing variable names - SHOULD CHANGE
black_cl <- c() #95% CI (default) does not overlap 0
gray_cl <- c() #50% CI (default) does not overlap 0
white_cl <- c() #Both 50% and 95% CI overlap 0
for (i in 1:len) {
if ((thin_q[1, i] > marker & thin_q[2, i] > marker) |
(thin_q[1, i] < marker & thin_q[2, i] < marker)) {
black_cl <- c(black_cl, i)
} else {
if ((thick_q[1, i] > marker & thick_q[2, i] > marker) |
(thick_q[1, i] < marker & thick_q[2, i] < marker)) {
gray_cl <- c(gray_cl, i)
} else {
white_cl <- c(white_cl, i)
}
}
} # closes for (i in 1:len)
if (horiz == FALSE) {
v_black_cl <- rev(black_cl)
v_gray_cl <- rev(gray_cl)
v_white_cl <- rev(white_cl)
}
if (!is.null(object2)) {
black_cl2 <- c()
gray_cl2 <- c()
white_cl2 <- c()
for (i in 1:len) {
if ((thin_q2[1, i] > marker & thin_q2[2, i] > marker) |
(thin_q2[1, i] < marker & thin_q2[2, i] < marker)) {
black_cl2 <- c(black_cl2, i)
} else {
if ((thick_q2[1, i] > marker & thick_q2[2, i] > marker) |
(thick_q2[1, i] < marker & thick_q2[2, i] < marker)) {
gray_cl2 <- c(gray_cl2, i)
} else {
white_cl2 <- c(white_cl2, i)
}
}
} # closes for (i in 1:len)
if (horiz == FALSE) {
v_black_cl2 <- rev(black_cl2)
v_gray_cl2 <- rev(gray_cl2)
v_white_cl2 <- rev(white_cl2)
}
} # closes (!is.null(object2))
} # closes (ref_ovl == TRUE)
####################################################################################################
# Horizontal plot - CI lines parallel to x-axis
if (horiz == TRUE) {
if (missing(xlim)) {
thin_comb <- c(thin_q, thin_q2)
rn <- diff(range(thin_comb, na.rm = TRUE)) * PL_SC
XLIM <- c((min(thin_comb, na.rm = TRUE) - rn), (max(thin_comb, na.rm = TRUE) + rn))
} else {
XLIM <- xlim
}
if (len > 20) {
tt <- 1.5
} else {
alpha <- 0.667
beta <- 0.0412
tt <- alpha + beta * len
} # closes if (len > 20)
YLIM <- c((1 - tt), (len + tt))
if (missing(xlab)) xlab <- 'Parameter Estimate'
if (is.null(xlab)) xlab <- ''
if (missing(main)) main <- ''
if (missing(labels)) {
labs <- names(medians)
} else {
if (!missing(labels)) {
if (is.null(labels)) {
labs <- rep('', len)
}
if (!is.null(labels)) {
if (length(labels) == len) {
labs <- labels[idx]
} else {
stop('Labels length not equal to number of parameters')
}
}
}
} # closes if (missing(labels))
# 0.2 inches per line - mar measured in lines
m_char <- max(sapply(labs, function(x){graphics::strwidth(x, cex = sz_labels, units = 'in')})) / 0.2
graphics::par(mar = c(mar[1], (m_char + (mar[2] - 3)), mar[3], mar[4] - 1))
#plot blank plot
graphics::plot(
medians,
(1:len),
xlim = XLIM,
ylim = YLIM,
type = "n",
ann = TRUE,
xaxt = 'n',
yaxt = "n",
bty = "n",
ylab = NA,
xlab = xlab,
cex.lab = sz_ax_txt,
yaxs = 'i') #cex.lab is axis label
# title
graphics::title(main, cex.main = sz_main_txt)
if (guide_axis == TRUE) {
# top axis params
graphics::axis(3, lwd.ticks = sz_ax, labels = FALSE, at = pos_tick, lwd = sz_ax)
}
# bottom axis params
graphics::axis(1, lwd.ticks = sz_ax, labels = TRUE, at = pos_tick, lwd = sz_ax, cex.axis = sz_tick_txt)
# left axis params (labels)
# las - 0 parallel to axis, 1 horiz, 2 perp to axis, 3 vert
graphics::axis(
2,
at = ((1:len)),
tick = FALSE,
labels = labs,
las = 1,
adj = 0,
line = -1,
cex.axis = sz_labels)
if (guide_lines == TRUE) {
# limits
# determine where ticks were placed - from axTicks source code:
xaxp <- graphics::par("xaxp")
XLIM2 <- c(xaxp[1], xaxp[2])
xs <- matrix(rep(XLIM2, len), nrow = 2)
ys <- rbind(1:len, 1:len)
graphics::matlines(xs, ys, lty = 1, col = guide_col)
}
# ref line
if(!is.null(ref)) graphics::abline(v=ref, lty = 2, lwd = 3, col = ref_col)
if (ref_ovl == TRUE) {
# positions to plot CI
blk_bnd <- rbind(black_cl, black_cl)
gry_bnd <- rbind(gray_cl, gray_cl)
wht_bnd <- rbind(white_cl, white_cl)
#lightened colors
lcol <- colorspace::lighten(COL, LA)
bl_col <- COL[black_cl]
gr_col <- lcol[gray_cl]
wh_col <- lcol[white_cl]
#if any colors are 'black', change them to 'grey60' (for consistency with older pkg versions)
bl_idx <- which(COL == 'black')
gr_col[which(gray_cl %in% bl_idx)] <- 'grey60'
wh_col[which(white_cl %in% bl_idx)] <- 'grey60'
# Black CI
if (!is.null(black_cl)) {
#only one object
if (is.null(object2)) {
# Thick
graphics::matlines(
thick_q[,black_cl],
blk_bnd,
type = 'l',
lty = 1,
lwd = sz_thick,
col = bl_col)
# Thin
graphics::matlines(
thin_q[,black_cl],
blk_bnd,
type = 'l',
lty = 1,
lwd = sz_thin,
col = bl_col)
} else {
#two objects
blk_bnd2 <- rbind(black_cl2, black_cl2)
bl_col2 <- COL2[black_cl2]
# object
graphics::matlines(
thick_q[,black_cl],
(blk_bnd + offset),
type = 'l',
lty = 1,
lwd = sz_thick,
col = bl_col)
graphics::matlines(
thin_q[,black_cl],
(blk_bnd + offset),
type = 'l',
lty = 1,
lwd = sz_thin,
col = bl_col)
# object2
graphics::matlines(
thick_q2[, black_cl2],
(blk_bnd2 - offset),
type = 'l',
lty = 1,
lwd = sz_thick,
col = bl_col2)
graphics::matlines(
thin_q2[, black_cl2],
(blk_bnd2 - offset),
type = 'l',
lty = 1,
lwd = sz_thin,
col = bl_col2)
} # closes if (is.null(object2))
} # closes if (!is.null(black_cl))
# Gray CI
if (!is.null(gray_cl)) {
#only one object
if (is.null(object2)) {
# Thick
graphics::matlines(
thick_q[,gray_cl],
gry_bnd,
type = 'l',
lty = 1,
lwd = sz_thick,
col = gr_col)
# Thin
graphics::matlines(
thin_q[,gray_cl],
gry_bnd,
type = 'l',
lty = 1,
lwd = sz_thin,
col = gr_col)
} else {
#two objects
gry_bnd2 <- rbind(gray_cl2, gray_cl2)
#lightened colors for obj 2
lcol2 <- colorspace::lighten(COL2, LA)
gr_col2 <- lcol2[gray_cl2]
#if any colors are 'black', change them to 'grey60' (for consistency with older pkg versions)
bl_idx2 <- which(COL2 == 'black')
gr_col2[which(gray_cl2 %in% bl_idx2)] <- 'grey60'
# object
graphics::matlines(
thick_q[,gray_cl],
(gry_bnd + offset),
type = 'l',
lty = 1,
lwd = sz_thick,
col = gr_col)
graphics::matlines(
thin_q[,gray_cl],
(gry_bnd + offset),
type = 'l',
lty = 1,
lwd = sz_thin,
col = gr_col)
#object2
graphics::matlines(
thick_q2[,gray_cl2],
(gry_bnd2 - offset),
type = 'l',
lty = 1,
lwd = sz_thick,
col = gr_col2)
graphics::matlines(
thin_q2[,gray_cl2],
(gry_bnd2 - offset),
type = 'l',
lty = 1,
lwd = sz_thin,
col = gr_col2)
} # closes if (is.null(object2))
} # closes if (!is.null(gray_cl))
# White CI
if (!is.null(white_cl)) {
#only one object
if (is.null(object2)) {
graphics::matlines(
thick_q[,white_cl],
wht_bnd,
type = 'l',
lty = 1,
lwd = sz_thick,
col = wh_col) # white (gray)
graphics::matlines(
thin_q[,white_cl],
wht_bnd,
type = 'l',
lty = 1,
lwd = sz_thin,
col = wh_col) # white (gray)
} else {
#two objects
wht_bnd2 <- rbind(white_cl2, white_cl2)
#lightened colors for obj 2
lcol2 <- colorspace::lighten(COL2, LA)
wh_col2 <- lcol2[white_cl2]
#if any colors are 'black', change them to 'grey60' (for consistency with older pkg versions)
bl_idx2 <- which(COL2 == 'black')
wh_col2[which(white_cl2 %in% bl_idx2)] <- 'grey60'
# object
graphics::matlines(
thick_q[,white_cl],
(wht_bnd + offset),
type = 'l',
lty = 1,
lwd = sz_thick,
col = wh_col) # white (gray)
graphics::matlines(
thin_q[,white_cl],
(wht_bnd + offset),
type = 'l',
lty = 1,
lwd = sz_thin,
col = wh_col) # white (gray)
# object2
graphics::matlines(
thick_q2[,white_cl2],
(wht_bnd2 - offset),
type = 'l',
lty = 1,
lwd = sz_thick,
col = wh_col2) # white (gray)
graphics::matlines(
thin_q2[,white_cl2],
(wht_bnd2 - offset),
type = 'l',
lty = 1,
lwd = sz_thin,
col = wh_col2) #white (gray)
} # closes if (is.null(object2))
} # closes if (!is.null(white_cl))
# Medians
if (is.null(object2)) {
#only one object
graphics::points(medians, 1:len, pch = 16, col = 'white', cex = sz_med)
graphics::points(medians[black_cl], black_cl, pch = 16, col = bl_col, cex = sz_med)
graphics::points(medians[gray_cl], gray_cl, pch = 16, col = gr_col, cex = sz_med)
graphics::points(medians[white_cl], white_cl, pch = 21, col = wh_col, cex = sz_med, lwd = 2)
} else {
# two objects
graphics::points(medians, (1:len + offset), pch = 16, col = 'white', cex = sz_med)
graphics::points(medians[black_cl], (black_cl + offset), pch = 16, col = bl_col, cex = sz_med)
graphics::points(medians[gray_cl], (gray_cl + offset), pch = 16, col = gr_col, cex = sz_med)
graphics::points(medians[white_cl], (white_cl + offset), pch = 21, col = wh_col, cex = sz_med, lwd = 2)
# object2
graphics::points(medians2, (1:len - offset), pch = 16, col = 'white', cex = sz_med)
graphics::points(medians2[black_cl2], (black_cl2 - offset), pch = 16, col = bl_col2, cex = sz_med)
graphics::points(medians2[gray_cl2], (gray_cl2 - offset), pch = 16, col = gr_col2, cex = sz_med)
graphics::points(medians2[white_cl2], (white_cl2 - offset), pch = 21, col = wh_col2, cex = sz_med, lwd = 2)
}
} else {
if (is.null(object2))
{
graphics::matlines(thick_q[,1:len], rbind(1:len, 1:len), type = 'l', lty = 1, lwd = sz_thick, col = COL[idx])
graphics::matlines(thin_q[,1:len], rbind(1:len, 1:len), type = 'l', lty = 1, lwd = sz_thin, col = COL[idx])
graphics::points(medians[1:len], 1:len, pch = 16, col = COL[idx], cex = sz_med)
} else {
# object
graphics::matlines(thick_q[,1:len], (rbind(1:len, 1:len) + offset), type = 'l',
lty = 1, lwd = sz_thick, col = COL[idx])
graphics::matlines(thin_q[,1:len], (rbind(1:len, 1:len) + offset), type = 'l',
lty = 1, lwd = sz_thin, col = COL[idx])
graphics::points(medians[1:len], (1:len + offset), pch = 16, col = COL[idx], cex = sz_med)
# object2
graphics::matlines(thick_q2[,1:len], (rbind(1:len, 1:len) - offset), type = 'l',
lty = 1, lwd = sz_thick, col = COL2[idx])
graphics::matlines(thin_q2[,1:len], (rbind(1:len, 1:len) - offset), type = 'l',
lty = 1, lwd = sz_thin, col = COL2[idx])
graphics::points(medians2[1:len], (1:len - offset), pch = 16, col = COL2[idx], cex = sz_med)
}
}
} # closes if (horiz == TRUE)
####################################################################################################
# Vertical plot - CI lines parallel to x-axis
if (horiz == FALSE) {
if (missing(ylim)) {
thin_comb <- c(thin_q, thin_q2)
rn <- diff(range(thin_comb, na.rm = TRUE)) * PL_SC
YLIM <- c((min(thin_comb, na.rm = TRUE) - rn), (max(thin_comb, na.rm = TRUE) + rn))
} else {
YLIM <- ylim
}
if (len > 20) {
tt <- 1.5
} else {
alpha <- 0.667
beta <- 0.0412
tt <- alpha + beta * len
}
XLIM = c((1 - tt), (len + tt))
if (missing(ylab)) ylab = 'Parameter Estimate'
if (is.null(ylab)) ylab = ''
if (missing(main)) main = ''
if (missing(labels)) {
labs = names(medians)
} else {
if (!missing(labels)) {
if (is.null(labels)) {
labs <- rep('', len)
}
if (!is.null(labels)) {
if (length(labels) == len) {
labs <- labels[idx]
} else {
stop('Labels length not equal to number of parameters')
}
}
}
} # closes if (missing(labels))
# 0.2 inches per line - mar measured in lines
m_char <- (max(sapply(labs, function(x){graphics::strwidth(x, cex = sz_labels, units = 'i')})) / 0.2)
# blank plot - do not display
grDevices::pdf(file = NULL)
graphics::plot(
(len:1),
medians,
xlim = XLIM,
ylim = YLIM,
type = "n",
ann = TRUE,
xaxt = 'n',
yaxt = "n",
bty = "n",
ylab = NA,
xlab = NA,
cex.lab = sz_ax_txt,
xaxs = 'i')
# create invisible ticks to determine where to put y-axis label
tickp <- graphics::axis(
2,
lwd.ticks = sz_ax,
labels = FALSE,
at = pos_tick,
lwd = sz_ax,
cex.axis = sz_tick_txt,
las = 1,
col = 'white',
col.ticks = 'white')
invisible(grDevices::dev.off())
# determine how long labels are
ml_tickp <- max(graphics::strwidth(tickp, cex = sz_tick_txt, units = 'in'))
# 5 lines/inch
ll <- 1.8 + 5 * ml_tickp
# set plot margins according to labels
graphics::par(mar = c((m_char + (mar[1] - 4)), mar[2] + 1 + ll - 3, mar[3] - 1.5, mar[4]))
# new blank plot
graphics::plot(
(len:1),
medians,
xlim = XLIM,
ylim = YLIM,
type = "n",
ann = TRUE,
xaxt = 'n',
yaxt = "n",
bty = "n",
ylab = NA,
xlab = NA,
cex.lab = sz_ax_txt,
xaxs = 'i')
# ticks
graphics::axis(
2,
lwd.ticks = sz_ax,
labels = TRUE,
at = pos_tick,
lwd = sz_ax,
cex.axis = sz_tick_txt,
las = 1)
# y-axis label
graphics::title(ylab = ylab, cex.lab = sz_ax_txt, line = ll)
# title
graphics::title(main, cex.main = sz_main_txt, line = 0.5)
if (guide_axis == TRUE) {
# right axis params
graphics::axis(4, lwd.ticks = sz_ax, labels = FALSE, at = pos_tick, lwd = sz_ax)
}
# bottom axis params (labels)
# las - 0 parallel to axis, 1 horizontal, 2 perpendicular to axis, 3 vertical
graphics::axis(
1,
at = (len:1) + 0.013,
tick = FALSE,
labels = labs,
las = 2,
adj = 0,
line = -1,
cex.axis = sz_labels)
if (guide_lines == TRUE) {
# limits
yaxp <- graphics::par("yaxp")
YLIM2 <- c(yaxp[1], yaxp[2])
ys <- matrix(rep(YLIM2, len), nrow = 2)
xs <- rbind(1:len, 1:len)
graphics::matlines(xs, ys, lty = 1, col = guide_col)
}
# ref line
if (!is.null(ref)) {
graphics::abline(h=ref, lty = 2, lwd = 3, col = ref_col)
}
if (ref_ovl == TRUE) {
# positions to plot CI
v_blk_bnd <- matrix(rep(rev(len + 1 - black_cl), 2), nrow = 2, byrow = TRUE)
v_gry_bnd <- matrix(rep(rev(len + 1 - gray_cl), 2), nrow = 2, byrow = TRUE)
v_wht_bnd <- matrix(rep(rev(len + 1 - white_cl), 2), nrow = 2, byrow = TRUE)
# Black CI
if (!is.null(black_cl)) {
if (is.null(object2)) {
# Thick
graphics::matlines(v_blk_bnd, thick_q[,v_black_cl], type = 'l', lty = 1, lwd = sz_thick, col = 'black')
# Thin
graphics::matlines(v_blk_bnd, thin_q[,v_black_cl], type = 'l', lty = 1, lwd = sz_thin, col = 'black')
} else {
v_blk_bnd2 <- matrix(rep(rev(len + 1 - black_cl2), 2), nrow = 2, byrow = TRUE)
# object
graphics::matlines((v_blk_bnd - offset), thick_q[,v_black_cl], type = 'l', lty = 1, lwd = sz_thick, col = 'black')
graphics::matlines((v_blk_bnd - offset), thin_q[,v_black_cl], type = 'l', lty = 1, lwd = sz_thin, col = 'black')
# object2
graphics::matlines((v_blk_bnd2 + offset), thick_q2[,v_black_cl2], type = 'l', lty = 1, lwd = sz_thick, col = 'black')
graphics::matlines((v_blk_bnd2 + offset), thin_q2[,v_black_cl2], type = 'l', lty = 1, lwd = sz_thin, col = 'black')
} # closes if (is.null(object2))
} # closes if (!is.null(black_cl))
# Gray CI
if (!is.null(gray_cl)) {
if (is.null(object2)) {
# Thick
graphics::matlines(v_gry_bnd, thick_q[,v_gray_cl], type = 'l', lty = 1, lwd = sz_thick, col = gr_col)
# Thin
graphics::matlines(v_gry_bnd, thin_q[,v_gray_cl], type = 'l', lty = 1, lwd = sz_thin, col = gr_col)
} else {
v_gry_bnd2 <- matrix(rep(rev(len + 1 - gray_cl2), 2), nrow = 2, byrow = TRUE)
# object
graphics::matlines((v_gry_bnd - offset), thick_q[,v_gray_cl], type = 'l', lty = 1, lwd = sz_thick, col = gr_col)
graphics::matlines((v_gry_bnd - offset), thin_q[,v_gray_cl], type = 'l', lty = 1, lwd = sz_thin, col = gr_col)
#object2
graphics::matlines((v_gry_bnd2 + offset), thick_q2[,v_gray_cl2], type = 'l', lty = 1, lwd = sz_thick, col = gr_col)
graphics::matlines((v_gry_bnd2 + offset), thin_q2[,v_gray_cl2], type = 'l', lty = 1, lwd = sz_thin, col = gr_col)
}
}
# White CI
if (!is.null(white_cl)) {
if (is.null(object2)) {
graphics::matlines(v_wht_bnd, thick_q[,v_white_cl], type = 'l', lty = 1, lwd = sz_thick, col = gr_col) #white (gray)
graphics::matlines(v_wht_bnd, thin_q[,v_white_cl], type = 'l', lty = 1, lwd = sz_thin, col = gr_col) #white (gray)
} else {
v_wht_bnd2 <- matrix(rep(rev(len + 1 - white_cl2), 2), nrow = 2, byrow = TRUE)
# object
graphics::matlines((v_wht_bnd - offset), thick_q[,v_white_cl], type = 'l', lty = 1, lwd = sz_thick, col = gr_col)
graphics::matlines((v_wht_bnd - offset), thin_q[,v_white_cl], type = 'l', lty = 1, lwd = sz_thin, col = gr_col)
# object2
graphics::matlines((v_wht_bnd2 + offset), thick_q2[,v_white_cl2], type = 'l', lty = 1, lwd = sz_thick, col = gr_col)
graphics::matlines((v_wht_bnd2 + offset), thin_q2[,v_white_cl2], type = 'l', lty = 1, lwd = sz_thin, col = gr_col)
}
}
# Medians
if (is.null(object2)) {
graphics::points(len:1, medians, pch = 16, col = 'white', cex = sz_med)
graphics::points(v_blk_bnd[1,], medians[v_black_cl], pch = 16, col = 'black', cex = sz_med)
graphics::points(v_gry_bnd[1,], medians[v_gray_cl], pch = 16, col = gr_col, cex = sz_med)
graphics::points(v_wht_bnd[1,], medians[v_white_cl], pch = 21, col = gr_col, cex = sz_med, lwd = 2)
} else {
# object
graphics::points((len:1 - offset), medians, pch = 16, col = 'white', cex = sz_med)
graphics::points((v_blk_bnd[1,] - offset), medians[v_black_cl], pch = 16, col = 'black', cex = sz_med)
graphics::points((v_gry_bnd[1,] - offset), medians[v_gray_cl], pch = 16, col = gr_col, cex = sz_med)
graphics::points((v_wht_bnd[1,] - offset), medians[v_white_cl], pch = 21, col = gr_col, cex = sz_med, lwd = 2)
# object2
graphics::points((len:1 + offset), medians2, pch = 16, col = 'white', cex = sz_med)
graphics::points((v_blk_bnd2[1,] + offset), medians2[v_black_cl2], pch = 16, col = 'black', cex = sz_med)
graphics::points((v_gry_bnd2[1,] + offset), medians2[v_gray_cl2], pch = 16, col = gr_col, cex = sz_med)
graphics::points((v_wht_bnd2[1,] + offset), medians2[v_white_cl2], pch = 21, col = gr_col, cex = sz_med, lwd = 2)
}
} else {
if (is.null(object2)) {
graphics::matlines(rbind(1:len, 1:len), thick_q[,len:1], type = 'l', lty = 1, lwd = sz_thick, col = COL[idx])
graphics::matlines(rbind(1:len, 1:len), thin_q[,len:1], type = 'l', lty = 1, lwd = sz_thin, col = COL[idx])
graphics::points(1:len, medians[len:1], pch = 16, col = COL[idx], cex = sz_med)
} else {
# object
graphics::matlines((rbind(1:len, 1:len) - offset), thick_q[,len:1], type = 'l', lty = 1, lwd = sz_thick, col = COL[idx])
graphics::matlines((rbind(1:len, 1:len) - offset), thin_q[,len:1], type = 'l', lty = 1, lwd = sz_thin, col = COL[idx])
graphics::points((1:len - offset), medians[len:1], pch = 16, col = COL[idx], cex = sz_med)
# object2
graphics::matlines((rbind(1:len, 1:len) + offset), thick_q2[,len:1], type = 'l', lty = 1, lwd = sz_thick, col = COL2[idx])
graphics::matlines((rbind(1:len, 1:len) + offset), thin_q2[,len:1], type = 'l', lty = 1, lwd = sz_thin, col = COL2[idx])
graphics::points((1:len + offset), medians2[len:1], pch = 16, col = COL2[idx], cex = sz_med)
} # closes if (is.null(object2))
} # closes if (ref_ovl == TRUE)
} # closes if (horiz == FALSE)
# restore graphics
graphics::par(mar = c(5, 4, 4, 2) + 0.1)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.