R/diagnostic.bmeta.R

In jarbes: Just a Rather Bayesian Evidence Synthesis

#' @title Diagnostic function for bmeta object in jarbes
#'
#' @description This function performers an approximated Bayesian cross-validation for a b3lmeta object
#'
#' @param object The object generated by the function bmeta.
#' @param post.p.value.cut Posterior p-value cut point to assess outliers.
#' @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.bmeta = function(object,
                               post.p.value.cut = 0.05,
                               median.w = 1.5,
                               study.names = NULL,
                               size.forest = 0.4,
                               lwd.forest = 0.2,
                               shape.forest = 23,
                               ...) {

  x=y=ylo=yhi=NULL

  # Data preparation
  y.ghost = object$BUGSoutput$sims.list$y.ghost
     g.m = apply(y.ghost, 2, median)
     g.u = apply(y.ghost, 2, quantile, prob = 0.95)
     g.l = apply(y.ghost, 2, quantile, prob = 0.025)

  n.studies = length(g.m)

  TE = object$data$TE

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

  # Posterior p-values to detect outliers...
  p.vec = NULL
  for(i in 1:n.studies)
  {
    p1 = sum(y.ghost[,i]<TE[i])/length(y.ghost[,i])
    p2 = sum(y.ghost[,i]>TE[i])/length(y.ghost[,i])
    p.val = min(p1, p2)
    p.vec = c(p.vec, p.val)
    }

  p.col = ifelse(p.vec < post.p.value.cut, "red", "blue")

  data.plot = data.frame(
             x = study.names,
            TE = TE,
            g.m = g.m,
            ylo  = g.l,
            yhi  = g.u,
            p.vec = p.vec,
            p.col = p.col)


  p = ggplot(data.plot, aes(x = x, y = TE,
                            ymin = ylo, ymax = yhi,
                            size = size.forest      # Point size
                              )) +
        geom_pointrange(colour = p.col,
                      lwd = lwd.forest,             # Thickness of the lines
                    shape = shape.forest)+
    # geom_pointrange(aes(x =x, y = TE,ymin = ylo, ymax = yhi),
    #                 colour = p.col)+
    coord_flip() +
    xlab("Study") +
    ylab("Posterior Predictive observation") +
    ggtitle("Bayesian Cross-Valdiation") +
    theme_bw()


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

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