Nothing
## Note by Valentijn: I can't seem to make it pass all the checks for forest.numeric or forest.default
# (by Thomas, below) work without making a generic method. But that masks the generic from metafor.
# The issue is that i was trying to unmask the generic from metafor...
#' Forest plot
#'
#' Generate a forest plot by specifying the various effect sizes, confidence intervals and summary estimate.
#'
#' @author Thomas Debray <thomas.debray@gmail.com>
#' @author Valentijn de Jong <Valentijn.M.T.de.Jong@gmail.com>
#'
#' @param \dots Additional arguments, which are currently ignored.
#'
#' @details This is a generic function. See \link{forest.default} for making forest plots of summary statistics,
#' \link{forest.metapred} for plotting \link{metapred} objects, and \link{forest.mp.cv.val} for plotting
#' \code{mp.cv.val} objects.
#'
#' @export forest
forest <- function(...)
UseMethod("forest")
#' Forest plot
#'
#' Generate a forest plot by specifying the various effect sizes, confidence intervals and summary estimate.
#' @param theta Numeric vector with effect size for each study
#' @param theta.ci.lb Numeric vector specifying the lower bound of the confidence interval of the effect sizes
#' @param theta.ci.ub Numeric vector specifying the upper bound of the confidence interval of the effect sizes
#' @param theta.slab Character vector specifying the study labels
#' @param theta.summary Meta-analysis summary estimate of the effect sizes
#' @param theta.summary.ci.lb Lower bound of the confidence (or credibility) interval of the summary estimate
#' @param theta.summary.ci.ub Upper bound of the confidence (or credibility) interval of the summary estimate
#' @param theta.summary.pi.lb Lower bound of the (approximate) prediction interval of the summary estimate.
#' @param theta.summary.pi.ub Upper bound of the (approximate) prediction interval of the summary estimate.
#' @param title Title of the forest plot
#' @param sort By default, studies are sorted by ascending effect size (\code{sort="asc"}). Set to \code{"desc"} for
#' sorting in reverse order, or any other value to ignore sorting.
#' @param theme Theme to generate the forest plot. By default, the classic dark-on-light ggplot2 theme is used.
#' See \link[ggplot2]{ggtheme} for more information.
#' @param predint.linetype The linetype of the prediction interval
#' @param xlim The \code{x} limits \code{(x1, x2)} of the forest plot
#' @param xlab Optional character string specifying the X label
#' @param refline Optional numeric specifying a reference line
#' @param label.summary Optional character string specifying the label for the summary estimate
#' @param label.predint Optional character string specifying the label for the (approximate) prediction interval
#' @param nrows.before.summary How many empty rows should be introduced between the study results and the summary estimates
#' @param study.digits How many significant digits should be used to print the stuy results
#' @param study.shape Plotting symbol to use for the study results. By default, a filled square is used.
#' @param col.diamond The filling color for the diamond representing the summary estimate.
#' E.g. "red", "blue", or hex color code ("#2e8aff")
#' @param col.predint Line color for the prediction interval. E.g. "red", "blue", or hex color code ("#2e8aff")
#' @param size.study Line width for the study results in mm
#' @param size.predint Line width for the prediction interval in mm
#' @param lty.ref Line type for the reference line
#' @param \dots Additional arguments, which are currently ignored.
#'
#' @author Thomas Debray <thomas.debray@gmail.com>
#'
#' @method forest default
#'
#' @return An object of class \code{ggplot}
#'
#' @import metafor
#' @import ggplot2
#' @importFrom stats reorder
#'
#' @export
# This could also be named forest.numeric
forest.default <- function(theta,
theta.ci.lb,
theta.ci.ub,
theta.slab,
theta.summary,
theta.summary.ci.lb, theta.summary.ci.ub,
theta.summary.pi.lb, theta.summary.pi.ub,
title,
sort = "asc",
theme = theme_bw(),
predint.linetype = 1,
xlim,
xlab = "",
refline = 0,
label.summary = "Summary Estimate",
label.predint = "Prediction Interval",
nrows.before.summary = 1,
study.digits = 2,
study.shape = 15,
col.diamond = "white",
col.predint = "black",
size.study = 0.5,
size.predint = 1,
lty.ref = "dotted",
...) {
requireNamespace("ggplot2")
if (missing(theta)) stop("Study effect sizes are missing!")
if (missing(theta.ci.lb) | missing(theta.ci.ub)) stop("Confidence intervals of effect sizes missing!")
if (missing(theta.slab)) stop("Study labels are missing!")
num.studies <- unique(c(length(theta), length(theta.ci.lb), length(theta.ci.ub), length(theta.slab)))
if (length(num.studies) > 1) stop(paste("Too few studies for a forest plot!"))
#Extract data
yi <- theta
k <- length(theta)
if (missing(theta.slab)) {
if (!is.null(attr(theta, "slab"))) {
slab <- attr(theta, "slab")
} else {
slab <- paste("Study", seq_len(k))
}
} else {
if (length(theta.slab) == 1 && is.na(theta.slab))
slab <- rep("", k)
else
slab <- theta.slab
}
add.predint <- TRUE # Add prediction interval by default
if (missing(theta.summary.pi.lb) | missing(theta.summary.pi.ub)) {
theta.summary.pi.lb <- theta.summary.pi.ub <- NA
}
if (NA %in% c(theta.summary.pi.lb, theta.summary.pi.ub)) {
add.predint <- FALSE
}
# Determine ordering
i.index <- 1:length(yi)
if (sort == "asc") {
i.index <- order(yi, decreasing = FALSE)
} else if (sort == "desc") {
i.index <- order(yi, decreasing = TRUE)
}
# Order the study results
scat <- rep(1, num.studies) #indicator variable for data points
slab <- slab[i.index]
yi <- yi[i.index]
ci.lb <- theta.ci.lb[i.index]
ci.ub <- theta.ci.ub[i.index]
if (!missing(theta.summary)) {
if (!add.predint) {
scat <- c(scat, 0)
slab <- c(slab, label.summary)
yi <- c(yi, theta.summary)
ci.lb <- c(ci.lb, theta.summary.ci.lb)
ci.ub <- c(ci.ub, theta.summary.ci.ub)
} else {
scat <- c(scat, 0, 0)
slab <- c(slab, label.summary, label.predint)
yi <- c(yi, theta.summary, theta.summary)
ci.lb <- c(ci.lb, theta.summary.ci.lb, theta.summary.pi.lb)
ci.ub <- c(ci.ub, theta.summary.ci.ub, theta.summary.pi.ub)
}
}
ALL <- data.frame("study" = slab,
"mean" = yi,
"m.lower" = ci.lb,
"m.upper" = ci.ub,
"order" = length(yi):1,
"scat" = scat)
# Add extra space between the study results and the summary estimates
rows.summaries <- which(ALL$study %in% c(label.summary, label.predint))
ALL$order[rows.summaries] <- ALL$order[rows.summaries] - nrows.before.summary
# reorder factor levels based on another variable (by yi)
ALL$study.ES_order <- reorder(ALL$study, ALL$order, mean)
# Information for secondary axis
labels_axis2 <- paste(format(ALL$mean, digits = study.digits, nsmall=2), " [",
format(ALL$m.lower, digits = study.digits, nsmall=2), " ; ",
format(ALL$m.upper, digits = study.digits, nsmall=2), "]", sep = "")
p <- with(ALL, ggplot(ALL[!is.na(ALL$mean), ],
aes(x = order, y = mean, ymin = m.lower, ymax = m.upper)) +
theme +
theme(axis.ticks.y.right = element_blank(),
panel.grid.minor.y = element_blank()) +
geom_pointrange(data = subset(ALL, scat == 1), shape = study.shape, size = size.study) +
scale_x_continuous(
breaks = order,
labels = study,
sec.axis = ggplot2::sec_axis(
~ .,
breaks = order,
labels = labels_axis2)
) +
coord_flip() +
ylab(xlab) +
xlab(""))
if (!missing(xlim)) {
p <- p + ylim(xlim)
}
# Add title
if (!missing(title)) {
p <- p + ggtitle(title)
}
# Add refline
if (!missing(refline)) {
if (is.numeric(refline)) {
p <- p + geom_hline(yintercept = refline, linetype = lty.ref)
}
}
# Add meta-analysis summary
if (!missing(theta.summary)) {
g2 <- with(ALL, subset(ALL, study == label.summary))
g2$ci.upper <- theta.summary.ci.ub
g2$ci.lower <- theta.summary.ci.lb
# Add confidence interval of the summary estimate
p <- p + with(g2, geom_errorbar(data = g2, aes(ymin = ci.lower, ymax = ci.upper, x = order), width = 0.5, size=1.0))
# Add summary estimate
p <- p + with(g2, geom_point(data = g2, aes(x = order, y = mean), shape=23, size=3, fill = col.diamond))
# Add diaomond
#df_diamond <- data.frame(row = c(g2$order+0.5, g2$order, g2$order-0.5, g2$order), ci = c(0, 1, 0, 0))
#print(df_diamond)
#p <- p + geom_polygon(data = df_diamond, aes(x = row, y = ci))
# Add (approximate) prediction interval
if (add.predint) {
g3 <- with(ALL, subset(ALL, study == label.predint))
g3$pi.upper <- theta.summary.pi.ub
g3$pi.lower <- theta.summary.pi.lb
p <- p + with(g3, geom_errorbar(data = g3,
aes(ymin = pi.lower, ymax = pi.upper, x = order),
size = size.predint,
width = 0.5,
color = col.predint,
linetype = predint.linetype))
}
}
p
}
#' Forest Plots
#'
#' Function to create forest plots for objects of class \code{"mm_perf"}.
#'
#' @param x An object of class \code{"mm_perf"}
#' @param \ldots Additional arguments which are passed to \link{forest}.
#'
#' @details The forest plot shows the performance estimates of each study with corresponding confidence
#' intervals.
#'
#' @references
#' Lewis S, Clarke M. Forest plots: trying to see the wood and the trees. \emph{BMJ}. 2001; 322(7300):1479--80.
#'
#' @examples
#' data(EuroSCORE)
#'
#' # Calculate the c-statistic and its standard error
#' est1 <- ccalc(cstat = c.index, cstat.se = se.c.index, cstat.cilb = c.index.95CIl,
#' cstat.ciub = c.index.95CIu, N = n, O = n.events, data = EuroSCORE, slab = Study)
#' plot(est1)
#'
#' # Calculate the total O:E ratio and its standard error
#' est2 <- oecalc(O = n.events, E = e.events, N = n, data = EuroSCORE, slab = Study)
#' plot(est2)
#'
#' @keywords forest
#'
#' @author Thomas Debray <thomas.debray@gmail.com>
#'
#' @return An object of class \code{ggplot}
#'
#' @method plot mm_perf
#' @export
plot.mm_perf <- function(x, ...) {
forest(theta = x$theta, theta.ci.lb = x$theta.cilb, theta.ci.ub = x$theta.ciub, theta.slab = rownames(x), ...)
}
# Multiple plot function
#
# ggplot objects can be passed in ..., or to plotlist (as a list of ggplot objects)
# - cols: Number of columns in layout
# - layout: A matrix specifying the layout. If present, 'cols' is ignored.
#
# If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE),
# then plot 1 will go in the upper left, 2 will go in the upper right, and
# 3 will go all the way across the bottom.
#
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
#library(grid)
# Make a list from the ... arguments and plotlist
plots <- c(list(...), plotlist)
numPlots = length(plots)
# If layout is NULL, then use 'cols' to determine layout
if (is.null(layout)) {
# Make the panel
# ncol: Number of columns of plots
# nrow: Number of rows needed, calculated from # of cols
layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
ncol = cols, nrow = ceiling(numPlots/cols))
}
if (numPlots==1) {
print(plots[[1]])
} else {
# Set up the page
grid.newpage()
pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
# Make each plot, in the correct location
for (i in 1:numPlots) {
# Get the i,j matrix positions of the regions that contain this subplot
matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
layout.pos.col = matchidx$col))
}
}
}
#' Posterior distribution of estimated model parameters
#'
#' Generate a plot of the posterior distribution
#'
#' @author Thomas Debray <thomas.debray@gmail.com>
#' @param \dots Additional arguments, which are currently ignored.
#' @return A \code{ggplot} object.
#'
#' @details This is a generic function.
#'
#' @export dplot
dplot <- function(...)
UseMethod("dplot")
#' Gelman-Rubin-Brooks plot
#'
#' This plot shows the evolution of Gelman and Rubin's shrink factor as the number of iterations increases. The code is adapted from
#' the R package coda.
#'
#' @param \ldots Additional arguments which are currently not used
#' @return A \code{ggplot} object.
#'
#' @author Thomas Debray <thomas.debray@gmail.com>
#'
#' @return An object of class \code{ggplot}
#'
#' @export gelmanplot
gelmanplot <- function(...)
UseMethod("gelmanplot")
#' Plot the autocorrelation of a Bayesian meta-analysis model
#'
#' Function to display autocorrelation of a fitted Bayesian meta-analysis model.
#'
#' @author Thomas Debray <thomas.debray@gmail.com>
#' @param \dots Additional arguments, which are currently ignored.
#' @return A \code{ggplot} object.
#'
#' @details This is a generic function.
#'
#' @export acplot
acplot <- function(...)
UseMethod("acplot")
#' Plot the running means of a Bayesian meta-analysis model
#'
#' Function to display running means of a fitted Bayesian meta-analysis model.
#'
#' @author Thomas Debray <thomas.debray@gmail.com>
#' @param \dots Additional arguments, which are currently ignored.
#' @return A \code{ggplot} object.
#'
#' @details This is a generic function.
#'
#' @export rmplot
rmplot <- function(...)
UseMethod("rmplot")
#' Plot the autocorrelation of a Bayesian meta-analysis model
#'
#' Function to display autocorrelation of a fitted Bayesian meta-analysis model.
#'
#' @param x An object of class \code{"mcmc.list"}
#' @param P Optional dataframe describing the parameters to plot and their respective names
#' @param greek Logical value indicating whether parameter labels have to be parsed to get Greek letters. Defaults to FALSE.
#' @param nLags Integer indicating the number of lags of the autocorrelation plot.
#' @param \ldots Additional arguments which are passed to ggs_autocorrelation
#' @return A \code{ggplot} object.
#'
#' @keywords meta-analysis autocorrelation
#'
#' @author Thomas Debray <thomas.debray@gmail.com>
#'
#' @export
#' @importFrom dplyr group_by do %>%
acplot.mcmc.list <- function(x, P, nLags = 50, greek = FALSE, ...) {
requireNamespace("ggmcmc")
if (!missing(P)) {
S <- ggmcmc::ggs(x, par_labels = P, sort = FALSE)
D <- subset(S, S$ParameterOriginal %in% P$Parameter)
} else {
D <- ggmcmc::ggs(x, sort = FALSE)
}
g <- ggmcmc::ggs_autocorrelation(D = D, nLags = nLags, greek = greek, ...)
g <- g + theme(strip.text.y = element_text(angle = 0))
return(g)
}
#' Plot the running means of a Bayesian meta-analysis model
#'
#' Function to display running means of a fitted Bayesian meta-analysis model.
#'
#' @param x An object of class \code{"mcmc.list"}
#' @param P Optional dataframe describing the parameters to plot and their respective names
#' @param greek Logical value indicating whether parameter labels have to be parsed to get Greek letters. Defaults to FALSE.
#' @param \ldots Additional arguments which are passed to ggs_running
#' @return A \code{ggplot} object.
#'
#'
#' @author Thomas Debray <thomas.debray@gmail.com>
#'
#' @export
rmplot.mcmc.list <- function(x, P, greek = FALSE, ...) {
requireNamespace("ggmcmc")
if (!missing(P)) {
S <- ggmcmc::ggs(x, par_labels = P, sort = FALSE)
D <- subset(S, S$ParameterOriginal %in% P$Parameter)
} else {
D <- ggmcmc::ggs(x, sort = FALSE)
}
g <- ggmcmc::ggs_running(D = D, greek = greek, ...)
g <- g + theme(strip.text.y = element_text(angle = 0))
return(g)
}
#' Posterior distribution of estimated model parameters
#'
#' Generate a plot of the posterior distribution
#'
#' @param x An object of class \code{"mcmc.list"}
#' @param P Optional dataframe describing the parameters to plot and their respective names
#' @param plot_type Optional character string to specify whether a density plot (\code{"dens"}) or
#' histogram (\code{"hist"}) should be displayed.
#' @param \ldots Additional arguments which are currently not used
#' @return A \code{ggplot} object.
#'
#'
#' @keywords meta-analysis density distribution
#'
#' @author Thomas Debray <thomas.debray@gmail.com>
#'
#' @export
dplot.mcmc.list <- function(x, P, plot_type = "dens", ...) {
requireNamespace("ggmcmc")
if (!missing(P)) {
S <- ggmcmc::ggs(x, par_labels = P, sort = FALSE)
S <- subset(S, S$ParameterOriginal %in% P$Parameter)
} else {
S <- ggmcmc::ggs(x, sort = FALSE)
}
if (plot_type == "dens") {
ggmcmc::ggs_density(S)
} else if (plot_type == "hist") {
ggmcmc::ggs_histogram(S)
} else {
stop("Invalid plot type")
}
}
#' Gelman-Rubin-Brooks plot
#'
#' This plot shows the evolution of Gelman and Rubin's shrink factor as the number of iterations increases. The code is adapted from
#' the R package coda.
#'
#' @param x An mcmc object
#' @param P Optional dataframe describing the parameters to plot and their respective names
#' @param confidence The coverage probability of the confidence interval for the potential scale reduction factor
#' @param max.bins Maximum number of bins, excluding the last one.
#' @param autoburnin Logical flag indicating whether only the second half of the series should be used in the computation.
#' If set to TRUE (default) and start(x) is less than end(x)/2 then start of series will be adjusted so that only second half of series is used.
#' @param greek Logical value indicating whether parameter labels have to be parsed to get Greek letters. Defaults to false.
#' @param \ldots Additional arguments which are currently not used
#' @return A \code{ggplot} object.
#'
#'
#' @author Thomas Debray <thomas.debray@gmail.com>
#'
#' @return An object of class \code{ggplot}
#'
#' @export
gelmanplot.mcmc.list <- function(x, P, confidence = 0.95, max.bins = 50, autoburnin = TRUE, greek = FALSE, ...) {
requireNamespace("coda")
nbin <- min(floor((coda::niter(x) - 50)/coda::thin(x)), max.bins)
if (nbin < 1) {
stop("Insufficient iterations to produce Gelman-Rubin plot")
}
binw <- floor((coda::niter(x) - 50)/nbin)
last.iter <- c(seq(from = start(x) + 50 * coda::thin(x), by = binw *
coda::thin(x), length = nbin), end(x))
shrink <- array(dim = c(nbin + 1, coda::nvar(x), 2))
dimnames(shrink) <- list(last.iter, coda::varnames(x),
c("median", paste(50 * (confidence + 1), "%", sep = "")))
for (i in 1:(nbin + 1)) {
shrink[i, , ] <- coda::gelman.diag(window(x, end = last.iter[i]),
confidence = confidence,
transform = FALSE,
autoburnin = autoburnin,
multivariate = FALSE)$psrf
}
all.na <- apply(is.na(shrink[, , 1, drop = FALSE]), 2, all)
if (any(all.na)) {
cat("\n******* Error: *******\n")
cat("Cannot compute Gelman & Rubin's diagnostic for any chain \n")
cat("segments for variables", coda::varnames(x)[all.na], "\n")
cat("This indicates convergence failure\n")
}
labels_gelmandiag <- names(shrink[1,1,])
ggdatMed <- data.frame(cbind(shrinkfactor = as.vector(shrink[,,labels_gelmandiag[1]]),
parameter = rep(colnames(shrink[,,1]), each = nrow(shrink[,,1])),
type = labels_gelmandiag[1],
last.iter = rep(last.iter, ncol(shrink[,,1]))))
ggdatBnd <- data.frame(cbind(shrinkfactor = as.vector(shrink[,,labels_gelmandiag[2]]),
parameter = rep(colnames(shrink[,,2]), each = nrow(shrink[,,2])),
type = labels_gelmandiag[2],
last.iter = rep(last.iter, ncol(shrink[,,2]))))
ggdat <- rbind(ggdatMed, ggdatBnd)
ggdat$last.iter <- as.numeric(as.character(ggdat$last.iter))
ggdat$shrinkfactor <- as.numeric(as.character(ggdat$shrinkfactor))
if (!missing(P)) {
ggdat <- subset(ggdat, ggdat$parameter %in% P$Parameter)
ggdat$parameter <- factor(ggdat$parameter, levels = P$Parameter, labels = P$Label)
}
if (greek) {
g <- with(ggdat, ggplot(ggdat, aes(last.iter, shrinkfactor, color = type)) +
geom_line() + facet_wrap(~ parameter, labeller = label_parsed) +
xlab("Last iteration in chain") +
ylab("Shrink factor"))
} else {
g <- with(ggdat, ggplot(ggdat, aes(last.iter, shrinkfactor, color = type)) +
geom_line() + facet_wrap(~ parameter) +
xlab("Last iteration in chain") +
ylab("Shrink factor"))
}
return(g)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.