R/forest.netmeta.R

In netmeta: Network Meta-Analysis using Frequentist Methods

#' Forest plot for network meta-analysis
#' 
#' @description
#' Draws a forest plot in the active graphics window (using grid
#' graphics system).
#'
#' @aliases forest.netmeta plot.netmeta
#' 
#' @param x An object of class \code{netmeta}.
#' @param pooled A character string indicating whether results for the
#'   common (\code{"common"}) or random effects model
#'   (\code{"random"}) should be plotted. Can be abbreviated.
#' @param reference.group Reference treatment(s).
#' @param baseline.reference A logical indicating whether results
#'   should be expressed as comparisons of other treatments versus the
#'   reference treatment (default) or vice versa.
#' @param labels An optional vector with treatment labels.
#' @param equal.size A logical indicating whether all squares should
#'   be of equal size. Otherwise, the square size is proportional to
#'   the precision of estimates.
#' @param leftcols A character vector specifying columns to be plotted
#'   on the left side of the forest plot or a logical value (see
#'   Details).
#' @param leftlabs A character vector specifying labels for
#'   (additional) columns on left side of the forest plot (see
#'   Details).
#' @param rightcols A character vector specifying columns to be
#'   plotted on the right side of the forest plot or a logical value
#'   (see Details).
#' @param rightlabs A character vector specifying labels for
#'   (additional) columns on right side of the forest plot (see
#'   Details).
#' @param digits Minimal number of significant digits for treatment
#'   effects and confidence intervals, see \code{print.default}.
#' @param small.values A character string specifying whether small
#'   treatment effects indicate a beneficial (\code{"desirable"}) or
#'   harmful (\code{"undesirable"}) effect, can be abbreviated; see
#'   \code{\link{netrank}}.
#' @param nsim Number of simulations to calculate SUCRAs.
#' @param digits.prop Minimal number of significant digits for
#'   P-scores, SUCRAs and direct evidence proportions, see
#'   \code{\link{print.default}} and \code{\link{netrank}}.
#' @param smlab A label printed at top of figure. By default, text
#'   indicating either common or random effects model is printed.
#' @param sortvar An optional vector used to sort treatments (must be
#'   of same length as the total number of treatments).
#' @param backtransf A logical indicating whether results should be
#'   back transformed in forest plots. If \code{backtransf = TRUE},
#'   results for \code{sm = "OR"} are presented as odds ratios rather
#'   than log odds ratios, for example.
#' @param lab.NA A character string to label missing values.
#' @param add.data An optional data frame with additional columns to
#'   print in forest plot (see Details).
#' @param drop.reference.group A logical indicating whether the
#'   reference group should be printed in the forest plot.
#' @param col.by The colour to print information on subgroups.
#' @param print.subgroup.name A logical indicating whether the name of
#'   the grouping variable should be printed in front of the group
#'   labels.
#' @param \dots Additional arguments for \code{\link{forest.meta}}
#'   function.
#' 
#' @details
#' A forest plot, also called confidence interval plot, is drawn in
#' the active graphics window.
#' 
#' Argument \code{sortvar} can be either a numeric or character vector
#' with length of number of treatments. If \code{sortvar} is numeric
#' the \code{\link[base]{order}} function is utilised internally to
#' determine the order of values. If \code{sortvar} is character it
#' must be a permutation of the treatment names. It is also possible
#' to provide either \code{sortvar = Pscore}, \code{sortvar =
#' "Pscore"}, \code{sortvar = -Pscore}, or \code{sortvar = "-Pscore"}
#' in order to sort treatments according to the ranking generated by
#' \code{\link{netrank}} which is called internally. It is also
#' possible to use "SUCRA" instead of "Pscore". Similar expressions
#' are possible to sort by treatment comparisons (\code{sortvar = TE},
#' etc.), standard error (\code{sortvar = seTE}), number of studies
#' with direct treatment comparisons (\code{sortvar = k}), and direct
#' evidence proportion (\code{sortvar = prop.direct}, see also
#' \code{\link{netmeasures}}).
#'
#' The arguments \code{leftcols} and \code{rightcols} can be used to
#' specify columns which are plotted on the left and right side of the
#' forest plot, respectively. The following columns are available:
#' \tabular{ll}{
#' \bold{Name} \tab \bold{Definition} \cr
#' \code{"studlab"} \tab Treatments \cr
#' \code{"TE"} \tab Network estimates (either from common or random
#'   effects model) \cr
#' \code{"seTE"} \tab Corresponding standard errors \cr
#' \code{"Pscore"} \tab P-scores (see \code{\link{netrank}}) \cr
#' \code{"SUCRA"} \tab SUCRAs (see \code{\link{netrank}}) \cr
#' \code{"n.trts"} \tab Number of participants per treatment arm \cr
#' \code{"k"} \tab Number of studies in pairwise comparisons \cr
#' \code{"prop.direct"} \tab Direct evidence proportions (see
#'   \code{\link{netmeasures}}) \cr
#' \code{"effect"} \tab (Back-transformed) network estimates \cr
#' \code{"ci"} \tab Confidence intervals \cr
#' \code{"effect.ci"} \tab (Back-transformed) network estimates and
#'   confidence intervals
#' }
#'
#' As a sidenote, the rather odd column name \code{"studlab"} to
#' describe the treatment comparisons comes from internally calling
#' \code{\link{forest.meta}} which uses study labels as the essential
#' information.
#' 
#' Argument \code{add.data} can be used to add additional columns to
#' the forest plot. This argument must be a data frame with row names
#' equal to the treatment names in R object \code{x}, i.e.,
#' \code{x$trts}.
#' 
#' See help page of \code{\link{forest.meta}} for more information on
#' the generation of forest plots and additional arguments.
#' 
#' @author Guido Schwarzer \email{guido.schwarzer@@uniklinik-freiburg.de}
#' 
#' @seealso \code{\link{forest.meta}}
#' 
#' @keywords hplot
#' 
#' @examples
#' data(smokingcessation)
#' 
#' # Transform data from arm-based format to contrast-based format
#' #
#' p1 <- pairwise(list(treat1, treat2, treat3),
#'   event = list(event1, event2, event3), n = list(n1, n2, n3),
#'   data = smokingcessation, sm = "OR")
#' 
#' # Conduct random effects network meta-analysis
#' #
#' net1 <- netmeta(p1, common = FALSE)
#'
#' forest(net1)
#' 
#' \dontrun{
#' data(Senn2013)
#' 
#' # Conduct network meta-analysis
#' #
#' net2 <- netmeta(TE, seTE, treat1, treat2, studlab,
#'   data = Senn2013, sm = "MD")
#' 
#' forest(net2, ref = "plac")
#' forest(net2, xlim = c(-1.5, 1), ref = "plac",
#'   xlab = "HbA1c difference", rightcols = FALSE)
#' 
#' # Random effects effect model
#' #
#' net3 <- netmeta(TE, seTE, treat1, treat2, studlab,
#'   data = Senn2013, sm = "MD", common = FALSE)
#' 
#' forest(net3, xlim = c(-1.5, 1), ref = "plac",
#'   xlab = "HbA1c difference")
#' 
#' # Add column with P-Scores on right side of forest plot
#' #
#' forest(net3, xlim = c(-1.5, 1), ref = "plac",
#'   xlab = "HbA1c difference",
#'   rightcols = c("effect", "ci", "Pscore"),
#'   just.addcols = "right")
#' 
#' # Add column with P-Scores on left side of forest plot
#' #
#' forest(net3, xlim = c(-1.5, 1), ref = "plac",
#'   xlab = "HbA1c difference",
#'   leftcols = c("studlab", "Pscore"),
#'   just.addcols = "right")
#' 
#' # Sort forest plot by descending P-Score
#' #
#' forest(net3, xlim = c(-1.5, 1), ref = "plac",
#'   xlab = "HbA1c difference",
#'   rightcols = c("effect", "ci", "Pscore"),
#'   just.addcols = "right",
#'   sortvar = -Pscore)
#' 
#' # Drop reference group and sort by and print number of studies with
#' # direct treatment comparisons
#' #
#' forest(net3, xlim = c(-1.5, 1), ref = "plac",
#'   xlab = "HbA1c difference",
#'   leftcols = c("studlab", "k"),
#'   leftlabs = c("Contrast\nto Placebo", "Direct\nComparisons"),
#'   sortvar = -k,
#'   drop = TRUE,
#'   smlab = "Random Effects Model")
#' }
#' 
#' @method forest netmeta
#' @export


forest.netmeta <- function(x,
                           pooled = ifelse(x$random, "random", "common"),
                           reference.group = x$reference.group,
                           baseline.reference = x$baseline.reference,
                           labels = x$trts,
                           equal.size = TRUE,
                           leftcols = "studlab",
                           leftlabs,
                           rightcols = c("effect", "ci"),
                           rightlabs,
                           digits = gs("digits.forest"),
                           small.values = x$small.values,
                           nsim = 1000,
                           digits.prop = 2,
                           smlab = NULL,
                           sortvar = x$seq,
                           backtransf = x$backtransf,
                           lab.NA = ".",
                           add.data,
                           drop.reference.group = FALSE,
                           ##
                           col.by = "black",
                           print.subgroup.name = FALSE,
                           ##
                           ...) {
  
  
  ##
  ##
  ## (1) Check and set arguments
  ##
  ##
  chkclass(x, "netmeta")
  x <- updateversion(x)
  ##
  is.bin <- inherits(x, "netmetabin")
  ##
  pooled <- setchar(pooled, c("common", "random", "fixed"))
  pooled[pooled == "fixed"] <- "common"
  ##
  chklogical(equal.size)
  ##
  chknumeric(digits, min = 0, length = 1)
  ##
  small.values <- setsv(small.values)
  ##
  chknumeric(nsim, min = 1, length = 1)
  chknumeric(digits.prop, min = 0, length = 1)
  ##
  chklogical(baseline.reference)
  ##
  trts <- x$trts
  ##
  if (!missing(labels)) {
    ##
    labels <- catch("labels", match.call(), x, sys.frame(sys.parent()))
    ##
    if (is.null(labels))
      stop("Argument 'labels' must be not NULL.")
    ##
    if (length(labels) != length(trts))
      stop("Length of argument 'labels' must be equal to number of treatments.")
    ##
    names(labels) <- trts
  }
  ##
  chklogical(drop.reference.group)
  chklogical(print.subgroup.name)
  ##
  chklogical(backtransf)
  chkchar(lab.NA)
  ##
  stdlabs <- c("event.e", "n.e", "event.c", "n.c",
               "mean.e", "sd.e", "mean.c", "sd.c",
               "n", "time", "event",
               "TE", "seTE",
               "time.e", "time.c",
               "effect", "ci", "effect.ci",
               "w.common", "w.random")
  ##
  if (missing(leftlabs)) {
    leftlabs <- leftcols
    leftlabs[leftcols %in% stdlabs] <- NA
    ##
    if (length(reference.group) > 1)
      leftlabs[matchVar(leftcols, "studlab")] <- "Comparison"
    else
      leftlabs[matchVar(leftcols, "studlab")] <- "Treatment"
  }
  else if (length(leftcols) != length(leftlabs)) {
    if (length(reference.group) > 1)
      leftlabs[matchVar(leftcols, "studlab")] <- "Comparison"
    else
      leftlabs[matchVar(leftcols, "studlab")] <- "Treatment"
  }
  ##
  if (missing(rightlabs)) {
    rightlabs <- rightcols
    rightlabs[rightcols %in% stdlabs] <- NA
  }
  ##
  for (i in names(list(...))) {
    if (!is.null(setchar(i, "weight.study", stop.at.error = FALSE)))
      stop("Argument 'weight.study' set internally.", call. = TRUE)
    if (!is.null(setchar(i, "prediction", stop.at.error = FALSE)))
      stop("For prediction intervals see example in help file of ",
           "forest.netsplit().", call. = TRUE)
  }
  
  
  ##
  ##
  ## (2) Extract results for common and random effects model and
  ##     calculate P-scores and SUCRAs if calcSUCRA == TRUE
  ##
  ##
  one.rg <- length(reference.group) == 1
  ##
  sortvar.c <- deparse(substitute(sortvar))
  sortvar.c <- gsub("\"", "", sortvar.c)
  ##
  calcPscore <-
    anyCol(rightcols, "Pscore") || anyCol(leftcols, "Pscore") ||
    any(matchVar(sortvar.c, "Pscore")) || any(matchVar(sortvar.c, "-Pscore"))
  ##
  calcSUCRA <-
    anyCol(rightcols, "SUCRA") || anyCol(leftcols, "SUCRA") ||
    any(matchVar(sortvar.c, "SUCRA")) || any(matchVar(sortvar.c, "-SUCRA"))
  ##
  if (one.rg && reference.group == "") {
    warning("First treatment used as reference as argument ",
            "'reference.group' is unspecified.",
            call. = FALSE)
    reference.group <- trts[1]
  }
  ##
  reference.group <- setref(reference.group, trts, length = 0)
  ##
  if (pooled == "common") {
    TE   <- x$TE.common
    seTE <- x$seTE.common
    ##
    prop.direct <- x$P.common
    ##
    if (calcPscore)
      Pscore <- netrank(x, small.values = small.values,
                        method = "P-score")$ranking.common
    if (calcSUCRA) {
      x$common <- TRUE
      x$random <- FALSE
      SUCRA <- netrank(x, small.values = small.values,
                       method = "SUCRA", nsim = nsim)$ranking.common
    }
    ##
    text.pooled <- "Common Effects Model"
    ##
    if (x$method == "MH")
      text.pooled <- "Mantel-Haenszel Method"
    else if (x$method == "NCH")
      text.pooled <- "Non-Central Hypergeometric"
  }
  ##
  if (pooled == "random") {
    TE   <- x$TE.random
    seTE <- x$seTE.random
    ##
    prop.direct <- x$P.random
    ##
    ##
    if (calcPscore)
      Pscore <- netrank(x, small.values = small.values,
                        method = "P-score")$ranking.random
    if (calcSUCRA) {
      x$common <- FALSE
      x$random <- TRUE
      SUCRA <- netrank(x, small.values = small.values,
                       method = "SUCRA", nsim = nsim)$ranking.random
    }
    ##
    text.pooled <- "Random Effects Model"
  }
  ##
  if (is.null(smlab)) {
    if (one.rg) {
      if (baseline.reference)
        smlab <- paste0("Comparison: other vs '",
                        reference.group, "'\n(",
                        text.pooled,
                        ")")
      else
        smlab <- paste0("Comparison: '",
                        reference.group,
                        "' vs other \n(",
                        text.pooled,
                        ")")
    }
    else
      smlab  <- text.pooled
  }
  ##
  rightcols <- setCol(rightcols, "Pscore")
  rightlabs <- setLab(rightlabs, rightcols, "Pscore", "P-score")
  ##
  rightcols <- setCol(rightcols, "SUCRA")
  rightlabs <- setLab(rightlabs, rightcols, "SUCRA", "SUCRA")
  ##
  rightcols <- setCol(rightcols, "n.trts")
  rightlabs <- setLab(rightlabs, rightcols, "n.trts",
                      "Number of\nParticipants")
  ##
  rightcols <- setCol(rightcols, "k")
  rightlabs <- setLab(rightlabs, rightcols, "k", "Direct\nComparisons")
  ##
  rightcols <- setCol(rightcols, "prop.direct")
  rightlabs <- setLab(rightlabs, rightcols, "prop.direct",
                      "Direct Evidence\nProportion")
  ##
  leftcols <- setCol(leftcols, "Pscore")
  leftlabs <- setLab(leftlabs, leftcols, "Pscore", "P-score")
  ##
  leftcols <- setCol(leftcols, "SUCRA")
  leftlabs <- setLab(leftlabs, leftcols, "SUCRA", "SUCRA")
  ##
  leftcols <- setCol(leftcols, "n.trts")
  leftlabs <- setLab(leftlabs, leftcols, "n.trts", "Number of\nParticipants")
  ##
  leftcols <- setCol(leftcols, "k")
  leftlabs <- setLab(leftlabs, leftcols, "k", "Direct\nComparisons")
  ##
  leftcols <- setCol(leftcols, "prop.direct")
  leftlabs <- setLab(leftlabs, leftcols, "prop.direct",
                     "Direct Evidence\nProportion")
  
  
  ##
  ##
  ## (3) Extract comparisons with reference group
  ##
  ##
  dat <- data.frame(comparison = character(0),
                    treat = character(0),
                    TE = numeric(0), seTE = numeric(0),
                    Pscore = numeric(0),
                    SUCRA = numeric(0),
                    k = numeric(0),
                    prop.direct = numeric(0),
                    stringsAsFactors = FALSE)
  ##
  for (i in seq_along(reference.group)) {
    rg.i <- reference.group[i]
    ##
    if (baseline.reference)
      dat.i <- data.frame(comparison = rg.i,
                          treat = colnames(TE),
                          labels = labels,
                          TE = TE[, colnames(TE) == rg.i],
                          seTE = seTE[, colnames(seTE) == rg.i],
                          Pscore = if (calcPscore) Pscore else NA,
                          SUCRA = if (calcSUCRA) SUCRA else NA,
                          k = x$A.matrix[, colnames(TE) == rg.i],
                          prop.direct =
                            if (is.bin) prop.direct
                            else prop.direct[, colnames(TE) == rg.i],
                          stringsAsFactors = FALSE)
    else
      dat.i <- data.frame(comparison = rg.i,
                          treat = rownames(TE),
                          labels = labels,
                          TE = TE[rownames(TE) == rg.i, ],
                          seTE = seTE[rownames(seTE) == rg.i, ],
                          Pscore = if (calcPscore) Pscore else NA,
                          SUCRA = if (calcSUCRA) SUCRA else NA,
                          k = x$A.matrix[rownames(TE) == rg.i, ],
                          prop.direct =
                            if (is.bin) prop.direct
                            else prop.direct[rownames(TE) == rg.i, ],
                          stringsAsFactors = FALSE)
    ##
    if (!is.null(x$n.trts))
      dat.i$n.trts <- x$n.trts
    ##
    if (!missing(add.data)) {
      if (!is.data.frame(add.data))
        stop("Argument 'add.data' must be a data frame.",
             call. = FALSE)
      if (nrow(add.data) != length(trts))
        stop("Dataset 'add.data' must have ", nrow(dat.i),
             " rows (corresponding to number of treatments)",
             call. = FALSE)
      if (any(rownames(add.data) != trts))
        stop("Dataset 'add.data' must have the following row names:\n",
             paste(paste("'", trts, "'", sep = ""), collapse = " - "),
             call. = FALSE)
      ##
      dat.i <- cbind(dat.i, add.data)
    }
    ##
    ## Sort dataset according to argument sortvar
    ##
    if (any(matchVar(sortvar.c, "Pscore")))
      sortvar <- Pscore
    else if (any(matchVar(sortvar.c, "-Pscore")))
      sortvar <- -Pscore
    else if (any(matchVar(sortvar.c, "SUCRA")))
      sortvar <- SUCRA
    else if (any(matchVar(sortvar.c, "-SUCRA")))
      sortvar <- -SUCRA
    else if (any(matchVar(sortvar.c, "TE")))
      sortvar <- dat.i$TE
    else if (any(matchVar(sortvar.c, "-TE")))
      sortvar <- -dat.i$TE
    else if (any(matchVar(sortvar.c, "seTE")))
      sortvar <- dat.i$seTE
    else if (any(matchVar(sortvar.c, "-seTE")))
      sortvar <- -dat.i$seTE
    else if (any(matchVar(sortvar.c, "k")))
      sortvar <- dat.i$k
    else if (any(matchVar(sortvar.c, "-k")))
      sortvar <- -dat.i$k
    else if (any(matchVar(sortvar.c, "n.trts")))
      sortvar <- dat.i$n.trts
    else if (any(matchVar(sortvar.c, "-n.trts")))
      sortvar <- -dat.i$n.trts
    else if (any(matchVar(sortvar.c, "prop.direct")))
      sortvar <- dat.i$prop.direct
    else if (any(matchVar(sortvar.c, "-prop.direct")))
      sortvar <- -dat.i$prop.direct
    ##
    if (!is.null(sortvar)) {
      if (is.character(sortvar))
        sort <- setseq(sortvar, trts)
      else
        sort <- order(sortvar)
      ##
      dat.i <- dat.i[sort, ]
    }
    ##
    if (drop.reference.group)
      dat.i <- subset(dat.i, treat != rg.i)
    ##
    if (baseline.reference)
      dat.i$comparison <- paste0("Other vs '", dat.i$comparison, "'")
    else
      dat.i$comparison <- paste0("'", dat.i$comparison, "' vs other")
    ##
    dat <- rbind(dat, dat.i)
  }
  ##
  dat.out <- dat
  ##
  if ("Pscore" %in% names(dat))
    dat$Pscore <- formatN(dat$Pscore, digits = digits.prop,
                          text.NA = lab.NA)
  ##
  if ("SUCRA" %in% names(dat))
    dat$SUCRA <- formatN(dat$SUCRA, digits = digits.prop,
                         text.NA = lab.NA)
  ##
  if ("prop.direct" %in% names(dat))
    dat$prop.direct <- formatN(dat$prop.direct,
                               digits = digits.prop, text.NA = lab.NA)
  ##
  rm(TE)
  rm(seTE)
  
  
  ##
  ##
  ## (5) Generate forest plot
  ##
  ##
  treat <- dat$treat
  ##
  if (one.rg)
    m1 <- suppressWarnings(metagen(TE, seTE, data = dat,
                                   sm = x$sm,
                                   studlab = labels, backtransf = backtransf,
                                   method.tau = "DL", method.tau.ci = "",
                                   warn = FALSE))
  else
    m1 <- suppressWarnings(metagen(TE, seTE, data = dat,
                                   subgroup = dat$comparison,
                                   sm = x$sm,
                                   studlab = labels, backtransf = backtransf,
                                   method.tau = "DL", method.tau.ci = "",
                                   warn = FALSE))
  ##
  forest(m1,
         digits = digits,
         overall = FALSE, common = FALSE, random = FALSE,
         hetstat = FALSE, test.subgroup = FALSE,
         leftcols = leftcols,
         leftlabs = leftlabs,
         rightcols = rightcols,
         rightlabs = rightlabs,
         smlab = smlab,
         lab.NA = lab.NA,
         ##
         col.by = col.by,
         print.subgroup.name = print.subgroup.name,
         ##
         weight.study = if (equal.size) "same" else pooled,
         ##
         ...)
  

  rownames(dat.out) <- seq_len(nrow(dat.out))
  ##
  attr(dat.out, "pooled") <- pooled
  attr(dat.out, "small.values") <- small.values
  ##
  invisible(dat.out)
}





#' @rdname forest.netmeta
#' @method plot netmeta
#' @export
#'

plot.netmeta <- function(x, ...)
  forest(x, ...)