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 w1 and w1.
#' 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.
#'
#'
#' @param m              The object generated by metadiag. The model object must be fitted with the options: re = "sm" and split.w = TRUE.
#'
#' @param group          An optional argument which is a variable name indicating a group factor.
#'                       If set, then the plot is colored by groups.
#' @param title          The title of the plot.
#'
#' @param group.colors   A character vector with two color names.
#'
#' @seealso \code{\link{metadiag}}.
#'
#' @keywords file
#'
#' @examples
#'
#' ## execute analysis
#' \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)
#' #Relationship between conflict and study design
#' plotw(m.ep, group = "design")
#'
#' }
#'
#' @import ggplot2
#' @import grid
#' @import gridExtra
#' @import R2jags
#' @import rjags
#' @importFrom stats median quantile

#' @export
plotw <- function(m, group=NULL,
                  title = "Posterior quantiles (25%, 50%, 75%)",
                  group.colors = c("blue", "red"))
{

  re <- m$re
  if(re!="sm")stop("This plot function is only for scale mixtures of normals.")


  split.w <- m$split.w
  if(split.w == FALSE)stop("This plot function is for split.w==TRUE")

  # Patch for the "note" no-visible-binding-for-global-variable
  w1.post <- m$BUGSoutput$sims.list$w1
  w2.post <- m$BUGSoutput$sims.list$w2

  nstudies <- dim(w1.post)[2]  #for the number of studies


  if(is.null(group)){
    dat.w <- data.frame(study.name = factor(1:nstudies),
                        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_string(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", size = 0.5, lty =2) +
      xlab("Study") +
      ylab("Weights") +
      ggtitle(title)
    } else if (!is.null(group))
      {
    dat.w <- data.frame(study.name = factor(1:nstudies),
                        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 = m$data[, group])

    theplot <- ggplot(dat.w, aes_string(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", size = 0.5, lty =2) +
      xlab("Study") +
      ylab("Weights") +
      ggtitle("Posteriors quantiles (25%, 50%, 75%)") +
      scale_color_manual(values = group.colors)

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

  return(theplot)
}