#' Display results of Copas selection modelling
#'
#' Four plots (selectable by 'which') are currently available: (1)
#' funnel plot, (2) contour plot, (3) treatment effect plot, (4)
#' p-value for residual publication bias plot. By default, all plots
#' are provided.
#'
#' Takes an object created by the \code{copas} function and draws up
#' to four plots to display the results of the Copas selection
#' modelling.
#'
#' The argument \code{which} specifies the plots to be drawn; plot
#' numbers below will be produced by setting \code{which=1}, etc.
#'
#' Plot 1: Funnel plot of studies in meta-analysis. Vertical grey line
#' is usual random effects estimate (DerSimonian-Laird method);
#' vertical broken line is common effects estimate.
#'
#' Plot 2: Plot of contours of treatment effect (estimated by the
#' Copas model) as the selection probability varies (the selection
#' probability is a function of gamma0 and gamma1 - see
#' \code{help(copas)} or the reference below).
#'
#' Plot 3: Assuming the contours of treatment effect in Plot 2 are
#' locally parallel, the results can be summarised in terms of the
#' probability of publishing the study with the largest standard
#' error. This plot displays the results of doing this, showing how
#' the estimated treatment effect (and \code{100*level}\% confidence
#' interval) vary as the probability of publishing the study with the
#' largest standard error decreases.
#'
#' The three horizontal grey lines are the usual random effects
#' treatment estimate (center) +/- the \code{100*level}\% confidence
#' interval (upper/lower grey lines).
#'
#' Plot 4: For any degree of selection (i.e. probability of the study
#' with largest SE being published), we can calculate a p-value for
#' the hypothesis that no further selection remains unexplained in the
#' data. These plot displays these p-values against the probability
#' that the study with the largest SE is published.
#'
#' Under the copas selection model, probabilities of the smallest
#' study being published which correspond to p-values for residual
#' selection bias that are larger than 0.1 are more plausible. The
#' corresponding treatment effect in plot 3 is thus the most plausible
#' under the copas selection model.
#'
#' \bold{Note}
#'
#' In the current version, fine control of the graphics parameters for
#' the individual panels is not possible. However, all the data used
#' to create the plots can be extracted manually from the object
#' created by the \code{copas} function (see attributes list for
#' \code{copas}) and used to create tailor-made plots.
#'
#' @param x An object of class \code{copas}, generated by the
#' \code{copas} function
#' @param which Specify plots required: 1:4 produces all plots
#' (default); 3 produces plot 3 etc; c(1,3) produces plots 1 and 3,
#' and so on.
#' @param main Specify plot captions. Must be of same length as
#' argument \code{which}.
#' @param xlim.pp A vector of x-axis limits for plots 3 and 4,
#' i.e. for the probability of publishing the study with largest
#' standard deviation. E.g. to specify limits between 0.3 and 0.1
#' set \code{xlim.pp=c(0.3,0.1)}.
#' @param orthogonal.line A logical indicating whether the orthogonal
#' line should be displayed in plot 2 (contour plot).
#' @param lines (Diagnostic use only) A logical indicating whether
#' regression lines should be plotted in contour plot. These
#' regression lines attempt to summarise each contour of constant
#' treatment effect by a straight line, prior to calculating the
#' orthogonal line. Regression lines with a positive adjusted
#' \code{R^2} will be printed in green color, others will be printed
#' in red color.
#' @param warn A number setting the handling of warning messages. It
#' is not uncommon for numerical problems to be encountered during
#' estimation over the grid of (gamma0, gamma1) values. Usually this
#' does not indicate a serious problem. This option specifies what
#' to do with warning messages. \code{warn=-1}: ignore all
#' warnings; \code{warn=0} (the default): store warnings till
#' function finishes; if there are less than 10, print them,
#' otherwise print a message saying warning messages were generated;
#' \code{warn=1}: print warnings as they occur; \code{warn=2}: stop
#' the function when the first warning is generated. For further
#' details see \code{help(options)}.
#' @param \dots Other arguments (to check for deprecated argument
#' 'caption').
#'
#' @author James Carpenter \email{James.Carpenter@@lshtm.ac.uk}, Guido
#' Schwarzer \email{guido.schwarzer@@uniklinik-freiburg.de}
#'
#' @seealso \code{\link{copas}}, \code{\link{summary.copas}},
#' \code{\link[meta]{metabias}}, \code{\link[meta]{metagen}}
#'
#' @references
#'
#' Carpenter JR, Schwarzer G, Rücker G, Künstler R (2009):
#' Empirical evaluation showed that the Copas selection model provided
#' a useful summary in 80\% of meta-analyses.
#' \emph{Journal of Clinical Epidemiology},
#' \bold{62}, 624--31
#'
#' Schwarzer G, Carpenter J, Rücker G (2010):
#' Empirical evaluation suggests Copas selection model preferable to
#' trim-and-fill method for selection bias in meta-analysis.
#' \emph{Journal of Clinical Epidemiology},
#' \bold{63}, 282--8
#'
#' @keywords hplot
#'
#' @examples
#' data(Fleiss1993bin, package = "meta")
#'
#' # Perform meta-analysis (outcome measure is OR = odds ratio)
#' #
#' m1 <- metabin(d.asp, n.asp, d.plac, n.plac, data = Fleiss1993bin, sm = "OR")
#'
#' # Perform Copas analysis
#' #
#' cop1 <- copas(m1)
#'
#' # Plot results
#' #
#' plot(cop1)
#'
#' # Only show plots 1 and 2 (without orthogonal line)
#' #
#' plot(cop1, which = 1:2, orth = FALSE)
#'
#' # Another example showing use of more arguments
#' # Note the use of "\n" to create a new line in the caption
#' #
#' plot(cop1, which = 3, xlim.pp = c(1, 0.5),
#' main = "Variation in estimated treatment\n effect with selection")
#'
#' @method plot copas
#' @export
#' @export plot.copas
#'
#' @importFrom graphics abline axTicks axis box contour mtext par plot
#' title
plot.copas <- function(x,
which = 1:4,
main = c("Funnel plot", "Contour plot",
"Treatment effect plot",
"P-value for residual selection bias"),
xlim.pp,
orthogonal.line = TRUE,
lines = FALSE,
warn = -1,
...) {
##
## Check class and arguments
##
chkclass(x, "copas")
##
missing.which <- missing(which)
chknumeric(which, min = 1, max = 4)
if (any(duplicated(which)))
stop("Duplicated values in argument 'which'.", call. = FALSE)
if (any(diff(which) < 0))
stop("Values of argument 'which' must be increasing.", call. = FALSE)
##
missing.main <- missing(main)
##
missing.xlim.pp <- missing(xlim.pp)
if (!missing.xlim.pp)
chknumeric(xlim.pp, length = 2)
##
chklogical(orthogonal.line)
chklogical(lines)
chknumeric(warn, length = 1)
##
## Check for deprecated arguments in '...'
##
args <- list(...)
## Check whether first argument is a list. In this case only use
## this list as input.
if (length(args) > 0 && is.list(args[[1]]))
args <- args[[1]]
##
## Argument 'caption' replaced by 'main'
##
main <- deprecated(main, missing.main, args, "caption")
##
if (missing.main) {
if (chkdeprecated(names(args), old = "caption", warn = FALSE)) {
if (length(main) == 4 & !missing.which)
main <- main[which]
}
else {
main <- main[which]
}
}
##
if (length(main) != length(which))
stop(if (length(main) > length(which)) "More" else "Fewer",
" titles (argument 'main') provided than figures.",
call. = FALSE)
##
if (length(main) == 4)
caption <- main
else {
caption <- rep("", 4)
caption[which] <- main
}
oldwarn <- options()$warn
on.exit(options(warn = oldwarn))
options(warn = warn)
TE <- x$TE
seTE <- x$seTE
TE.random <- x$TE.random
seTE.random <- x$seTE.random
gamma0 <- x$gamma0
gamma1 <- x$gamma1
TE.contour <- x$TE.contour
levels <- x$regr$levels
slope <- x$slope
x.slope <- x$x.slope
y.slope <- x$y.slope
TE.slope <- x$TE.slope
seTE.slope <- x$seTE.slope
publprob <- x$publprob
pval.rsb <- x$pval.rsb
sm <- x$sm
##
ord <- order(publprob)
show <- rep(FALSE, 4)
show[which] <- TRUE
##
if (sum(show) == 1) oldpar <- par(pty = "s")
if (sum(show) == 2) oldpar <- par(mfrow = c(1, 2), pty = "s")
if (sum(show) > 2) oldpar <- par(mfrow = c(2, 2),
mar = c(4.0, 4.1, 1.5, 0.5),
pty = "s")
##
on.exit(par(oldpar), add = TRUE)
if (is.relative.effect(sm))
sm <- paste("log ", sm, sep = "")
## Plot 1: funnel plot
##
if (show[1]) {
funnel(metagen(TE, seTE, level = 0.95,
common = TRUE, random = FALSE, sm = sm),
xlab = "", ylab = "")
abline(v = TE.random, lty = 1, col = "darkgray")
##
mtext(sm, side = 1, line = 2)
mtext("Standard error", side = 2, line = 2)
box()
title(main = caption[1])
}
## Plot 2: contour plot
##
## NB this is a square plot, with axes transformed.
##
gamma0.min <- min(gamma0)
gamma0.max <- max(gamma0)
gamma1.min <- min(gamma1)
gamma1.max <- max(gamma1)
##
gamma0.rescale <- (gamma0 - gamma0.min) / (gamma0.max - gamma0.min)
gamma1.rescale <- (gamma1 - gamma1.min) / (gamma1.max - gamma1.min)
##
if (show[2]) {
contour(gamma0.rescale, gamma1.rescale, TE.contour,
xlim = c(0, 1), ylim = c(0, 1),
xlab = "", ylab = "",
labcex = 0.8, lty = 2, col = "black", axes = FALSE, cex = 2,
levels = levels)
##
axis(side = 1, at = axTicks(1),
labels = round(axTicks(1) * (gamma0.max - gamma0.min) + gamma0.min, 2))
axis(side = 2, at = axTicks(2),
labels = round(axTicks(2) * (gamma1.max - gamma1.min) + gamma1.min, 2))
##
mtext("Values of gamma0", side = 1, line = 2)
mtext("Values of gamma1", side = 2, line = 2)
##
box()
##
## text(1,1, round(TE.random, 2)) not particularly useful
##
## Add estimated orthogonal line
##
if (orthogonal.line) {
lines(gamma0.rescale, 1 + slope * (gamma0.rescale - 1), lty = 1)
points(x.slope, y.slope)
}
##
if (lines) {
if (is.null(x$min.r.squared))
min.r.squared <- -100
else
min.r.squared <- x$min.r.squared
##
for (i in seq(along = x$regr$intercepts)[x$regr$adj.r.squareds > 0])
abline(x$regr$intercepts[i], x$regr$slopes[i], col = "green")
for (i in seq(along = x$regr$intercepts)[x$regr$adj.r.squareds <=0 ])
abline(x$regr$intercepts[i], x$regr$slopes[i], col = "red", lty = 3)
}
##
title(main = caption[2])
}
## Plot 3:
## mean (plus level%-CI) against prob of publishing trial
## with largest SE down orthogonal line
##
xvalue <- publprob[ord]
yvalue <- TE.slope[ord]
##
ci.y <- ci(yvalue, seTE.slope[ord], level = x$level.ma)
ci.y$lower[is.infinite(ci.y$lower)] <- NA
ci.y$upper[is.infinite(ci.y$upper)] <- NA
##
if (show[3]) {
if (missing.xlim.pp)
xlim <- range(xvalue, na.rm = TRUE)[c(2, 1)]
else
xlim <- xlim.pp
##
xlim[is.infinite(xlim)] <- c(1, 0)[is.infinite(xlim)]
##
if (diff(xlim) == 0)
xlim <- xlim + c(0.0001, -0.0001)
##
##
plot(xvalue, yvalue,
type = "l",
xlim = xlim,
ylim = c(
min(c(0, ci.y$lower), na.rm = TRUE),
max(c(ci.y$upper, 0), na.rm = TRUE)
),
xlab = "", ylab = "", axes = FALSE)
##
axis(side = 1, axTicks(1), labels = round(axTicks(1), 2))
axis(side = 2, axTicks(2), labels = round(axTicks(2), 2))
##
mtext("Probability of publishing the trial with largest se",
side = 1, line = 2)
mtext(sm, side = 2, line = 2)
##
points(xvalue, yvalue)
##
abline(h = 0)
##
abline(h = TE.random, lty = 1, col = "darkgray")
abline(h = x$lower.random, lty = 1, col = "gray")
abline(h = x$upper.random, lty = 1, col = "gray")
##
lines(xvalue, ci.y$lower, lty = 1)
lines(xvalue, ci.y$upper, lty = 1)
box()
title(main = caption[3])
}
## Plot 4:
## goodness of fit (p-value) against prob of publishing
## trial with largest SE
##
yvalue <- pval.rsb[ord]
sel.y <- !is.na(yvalue)
xvalue <- xvalue[sel.y]
yvalue <- yvalue[sel.y]
##
if (show[4]) {
if (missing.xlim.pp)
xlim <- range(xvalue, na.rm = TRUE)[c(2, 1)]
else
xlim <- xlim.pp
##
xlim[is.infinite(xlim)] <- c(1, 0)[is.infinite(xlim)]
##
if (diff(xlim) == 0)
xlim <- xlim + c(0.0001, -0.0001)
##
##
if (length(xvalue) == 0) {
plot(xvalue, yvalue,
xlim = xlim, ylim = c(0, 1),
xlab = "", ylab = "", axes = FALSE)
}
else {
plot(loess(yvalue ~ xvalue),
type = "l",
xlim = xlim, ylim = c(0, 1),
xlab = "", ylab = "", axes = FALSE)
}
##
points(xvalue, yvalue)
##
axis(side = 1, axTicks(1), labels = round(axTicks(1), 2))
axis(side = 2, axTicks(2), labels = round(axTicks(2), 2))
##
mtext("Probability of publishing the trial with largest se",
side = 1, line = 2)
mtext("P-value for residual selection bias", side = 2, line = 2)
##
abline(h = x$sign.rsb, lty = 2)
box()
##
title(main = caption[4])
}
invisible(NULL)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.