R/plot.copas.R

In metasens: Statistical Methods for Sensitivity Analysis in Meta-Analysis

#' 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)
}