R/plotw.R

In bamdit: Bayesian Meta-Analysis of Diagnostic Test Data

#' Plot for the conflict of evidence parameters w1 and w2
#'
#' Conflict of evidence plot: this plot displays the posterior distribution of the study's weights.
#' If split.w = TRUE, two weights (w1 and w2) are shown. These weights indicate potential conflict
#' of evidence of the studies. The weight w1 indicates deviations with respect to the specificity
#' and w2 to the sensitivity. If split.w = FALSE, a single weight w is plotted.
#'
#' @param m              The object generated by metadiag. The model object must be fitted with the option re = "sm".
#'
#' @param group          An optional argument which is a variable name indicating a group factor.
#'                       If set, then the plot is colored by groups.
#'
#' @param study.names    An optional argument specifying the column name in m$data containing study names.
#'                       If provided, these names are used on the study axis instead of numeric study IDs.
#'
#' @param title          The title of the plot.
#'
#' @param group.colors   A character vector with two color names.
#'
#' @param x.lab          Text with the label of the x-axis.
#'                       Default is "Variance inflation factor (1/w)".
#'
#' @seealso \code{\link{metadiag}}.
#'
#' @keywords file
#'
#' @examples
#'
#' \dontrun{
#' data(ep)
#' ep$design = factor(ep$d1,labels = c("prospective", "retrospective"))
#' m.ep <- metadiag(ep, re = "sm", re.model = "SeSp",
#'                 split.w = TRUE,
#'                 df.estimate = TRUE)
#'
#' plotw(m.ep)
#'
#' # split.w = FALSE case
#' m.ep2 <- metadiag(ep, re = "sm", re.model = "SeSp",
#'                  split.w = FALSE)
#' plotw(m.ep2)
#' }
#'
#' @import ggplot2
#' @import grid
#' @import gridExtra
#' @import R2jags
#' @import rjags
#' @importFrom stats median quantile
#'
#' @export
plotw <- function(m, group=NULL, study.names=NULL,
                  title = "Posterior quantiles (25%, 50%, 75%)",
                  group.colors = c("blue", "red"),
                  x.lab = "Variance inflation factor (1/w)")
{

  study.name <- w.median <- w.lower <- w.upper <- NULL
  re <- m$re
  if(re!="sm")stop("This plot function is only for scale mixtures of normals.")

  split.w <- m$split.w

  # ==============================
  # CASE 1: split.w = FALSE
  # ==============================
  if(split.w == FALSE){

    w.post <- m$BUGSoutput$sims.list$w
    nstudies <- dim(w.post)[2]

    # Study labels
    if(is.null(study.names)){
      study.labels <- factor(1:nstudies)
    } else {
      study.labels <- factor(m$data[, study.names])
    }

    if(is.null(group)){

      dat.w <- data.frame(
        study.name = study.labels,
        w.median = apply(w.post, 2, median),
        w.lower  = apply(w.post, 2, quantile, prob = 0.25),
        w.upper  = apply(w.post, 2, quantile, prob = 0.75)
      )

      theplot <- ggplot(dat.w,
                        aes(x = study.name, y = w.median,
                            ymin = w.lower, ymax = w.upper)) +
        coord_flip() +
        geom_pointrange(lwd = 0.8) +
        geom_hline(yintercept = 1, colour = "black",
                   linewidth = 0.5, lty =2) +
        xlab("Study") +
        ylab(x.lab) +
        ggtitle(title)

    } else if (!is.null(group)) {

      dat.w <- data.frame(
        study.name = study.labels,
        w.median = apply(w.post, 2, median),
        w.lower  = apply(w.post, 2, quantile, prob = 0.25),
        w.upper  = apply(w.post, 2, quantile, prob = 0.75),
        group = m$data[, group]
      )

      theplot <- ggplot(dat.w,
                        aes(x = study.name, y = w.median,
                            ymin = w.lower, ymax = w.upper,
                            color = group)) +
        coord_flip() +
        geom_pointrange(lwd =0.8) +
        geom_hline(yintercept = 1, colour = "black",
                   linewidth = 0.5, lty =2) +
        xlab("Study") +
        ylab(x.lab) +
        ggtitle(title) +
        scale_color_manual(values = group.colors)

    } else {
      stop("Argument 'group' has to be a factor!")
    }

    return(theplot)
  }

  # ==============================
  # CASE 2: split.w = TRUE
  # ==============================

  w1.post <- m$BUGSoutput$sims.list$w1
  w2.post <- m$BUGSoutput$sims.list$w2

  nstudies <- dim(w1.post)[2]

  if(is.null(study.names)){
    study.labels <- factor(1:nstudies)
  } else {
    study.labels <- factor(m$data[, study.names])
  }

  if(is.null(group)){

    dat.w <- data.frame(study.name = rep(study.labels, 2),
                        w.median = c(apply(w1.post, 2, median),
                                     apply(w2.post, 2, median)),
                        w.lower = c(apply(w1.post, 2, quantile, prob = 0.25),
                                    apply(w2.post, 2, quantile, prob = 0.25)),
                        w.upper = c(apply(w1.post, 2, quantile, prob = 0.75),
                                    apply(w2.post, 2, quantile, prob = 0.75)),
                        component = gl(2, nstudies, labels = c("TPR", "FPR")))

    theplot <- ggplot(dat.w,
                      aes(x = study.name, y = w.median, ymin = w.lower,
                          ymax = w.upper)) +
      coord_flip() +
      geom_pointrange(lwd = 0.8) +
      facet_grid(. ~ component) +
      geom_hline(yintercept = 1, colour = "black", linewidth = 0.5, lty =2) +
      xlab("Study") +
      ylab(x.lab) +
      ggtitle(title)

  } else if (!is.null(group)) {

    dat.w <- data.frame(study.name = rep(study.labels, 2),
                        w.median = c(apply(w1.post, 2, median),
                                     apply(w2.post, 2, median)),
                        w.lower = c(apply(w1.post, 2, quantile, prob = 0.25),
                                    apply(w2.post, 2, quantile, prob = 0.25)),
                        w.upper = c(apply(w1.post, 2, quantile, prob = 0.75),
                                    apply(w2.post, 2, quantile, prob = 0.75)),
                        component = gl(2, nstudies, labels = c("TPR", "FPR")),
                        group = rep(m$data[, group], 2))

    theplot <- ggplot(dat.w, aes(x = study.name, y = w.median, ymin = w.lower,
                                 ymax = w.upper, color = group)) +
      coord_flip() +
      geom_pointrange(lwd =0.8) +
      facet_grid(. ~ component) +
      geom_hline(yintercept = 1, colour = "black", linewidth = 0.5, lty =2) +
      xlab("Study") +
      ylab(x.lab) +
      ggtitle(title) +
      scale_color_manual(values = group.colors)

  } else {
    stop("Argument 'group' has to be a factor!")
  }

  return(theplot)
}