R/diagnostic.metarisk.R

In jarbes: Just a Rather Bayesian Evidence Synthesis

#' @title Diagnostic function for metarisk object in jarbes
#'
#' @description This function performers a specially designed diagnostic for a metarisk object
#'

#' @param object The object generated by the function hmr.
#' @param median.w Change color if median of a weight > median.w. The default value is 1.5.
#' @param study.names Character vector containing names of the studies used.
#' @param size.forest Size of the center symbol mark in the forest-plot lines
#' @param lwd.forest  Thickness of the lines in the forest-plot
#' @param shape.forest Type of symbol for the center mark in the forest-plot lines
#' @param ... \dots
#'

#'
#' @import ggplot2
#'
#' @export

diagnostic.metarisk = function(object,
                               median.w = 1.5,
                               study.names,
                               size.forest = 0.4,
                               lwd.forest = 0.2,
                               shape.forest = 23,
                               ...) {
  if (object$re != "sm")
    stop("This plot function is only for scale mixtures of normals.")

  x=y=ylo=yhi=NULL

  # Posteriors for the studies' weights
  # Forest plot for weights ...
  if (object$split.w == FALSE) {
    w = object$BUGSoutput$sims.list$w
    w.s = apply(w, 2, median)
    w.u = apply(w, 2, quantile, prob = 0.75)
    w.l = apply(w, 2, quantile, prob = 0.25)

    n.studies = length(w.s)

    w.col = ifelse(w.s<median.w, "red", "blue")

    if (missing(study.names)) {
      study.names = 1:n.studies
    }

    dat.weights = data.frame(x = as.character(study.names),
                             y = w.s,
                             ylo  = w.l,
                             yhi  = w.u)
  } else {
    w1 = object$BUGSoutput$sims.list$w1
    w2 = object$BUGSoutput$sims.list$w2
    w1.s = apply(w1, 2, median)
    w2.s = apply(w2, 2, median)
    w1.u = apply(w1, 2, quantile, prob = 0.75)
    w2.u = apply(w2, 2, quantile, prob = 0.75)
    w1.l = apply(w1, 2, quantile, prob = 0.25)
    w2.l = apply(w2, 2, quantile, prob = 0.25)

    n.studies = length(w1.s)

    w1.col = ifelse(w1.s<median.w, "red", "blue")
    w2.col = ifelse(w2.s<median.w, "red", "blue")
    w.col = c(w1.col, w2.col)

    if (missing(study.names)) {
      study.names = 1:n.studies
    }

    dat.weights = data.frame(x = as.character(study.names),
                             y = c(w1.s, w2.s),
                             ylo = c(w1.l, w2.l),
                             yhi = c(w1.u, w2.u),
                             component =
                               gl(2, n.studies,
                                  labels = c("Baseline Risk",
                                             "Treatment Effect")))
  }

  p = ggplot(dat.weights, aes(x = x, y = y,
                              ymin = ylo, ymax = yhi,
                              size = size.forest      # Point size
                              )) +
        geom_pointrange(colour = w.col,
                      lwd = lwd.forest,      # Thickness of the lines
                    shape = shape.forest)+
    coord_flip() +
    geom_hline(yintercept = 1, lty = 2) +
    xlab("Study") +
    ylab("Outliers detection weights") +
    ggtitle("Weights") +
    theme_bw()

  if (object$split.w == TRUE) p = p + facet_grid(. ~ component)

  return(p)
}