R/diagnostic.hmr.R

In jarbes: Just a Rather Bayesian Evidence Synthesis

#' Generic diagnostic function.
#' @param object The object generated by the function hmr.
#' @param ... \dots
#'
#'
#' @export

diagnostic = function(object, ...) UseMethod("diagnostic", object)



#' @title Diagnostic function for hmr object in jarbes
#'
#' @description This function performers a specially designed diagnostic for a hmr object
#'
#' @param object The object generated by the function hmr.
#' @param median.w Change colour 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 mu.phi     Prior-to-posterior sensitivity analysis of mu.phi. Default value is TRUE.
#' @param mu.phi.x.lim.low  Lower limit of the prior to posterior plot for mu.phi
#' @param mu.phi.x.lim.up  Upper limit of the prior to posterior plot for mu.phi

#' @param colour.hist.mu.phi colour of the posterior mu.phi histogram
#' @param colour.prior.mu.phi colour of the prior of mu.phi
#' @param colour.posterior.mu.phi colour of the posterior of mu.phi

#' @param title.plot.mu.phi  Text for the title in the mu phi plot.
#' @param title.plot.weights Text for the title of the posterior weights.
#' @param ... \dots
#'
#'
#' @import ggplot2
#' @import gridExtra
#'
#' @export
#'
diagnostic.hmr = function(object,
                          # Forest plot
                          median.w = 1.5,
                          study.names,
                          size.forest = 0.4,    # Point size
                          lwd.forest = 0.2,     # Thickness of the lines
                          shape.forest = 23,
                          mu.phi = TRUE,
                          mu.phi.x.lim.low = -10,
                          mu.phi.x.lim.up = 10,
                          colour.hist.mu.phi = "royalblue",
                          colour.prior.mu.phi = "black",
                          colour.posterior.mu.phi = "blue",
                          title.plot.mu.phi = "Prior-to-Posterior Sensitivity",
                          title.plot.weights = "Outlier Detection",
                          ...) {
  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.95)
    w.l = apply(w, 2, quantile, prob = 0.05)

    n.studies = length(w.s)

    w.col = ifelse(w.s<median.w, "red", "blue")
    w.col[length(w.s)] = "green"

    if (missing(study.names)) {
      study.names = 1:(n.studies-1)
    }
    study.names = c(as.character(study.names), "Cohort (individual data)")

    dat.weights = data.frame(x = 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")
    w1.col[n.studies] = "green"
    w2.col = ifelse(w2.s<median.w, "red", "blue")
    w2.col[n.studies] = "green"
    w.col = c(w1.col, w2.col)

    if (missing(study.names)) {
      study.names = 1:(n.studies-1)
    }
    study.names = c(as.character(study.names), "Cohort (individual data)")

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

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

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

  if (mu.phi == TRUE) {
    mu.phi = object$BUGSoutput$sims.list$mu.phi
    mean.mu.phi = object$prior$mean.mu.phi
      sd.mu.phi = object$prior$sd.mu.phi

    # Posterior distribution of mu.phi.........................................
    df.mu.phi = data.frame(x = mu.phi)

    mycolours = c("Prior" = colour.prior.mu.phi, "Posterior" = colour.posterior.mu.phi)

    diagnostic.plot1 = ggplot(df.mu.phi,
                       aes(x=x, colour = "Posterior"))+
      geom_density()+
      geom_histogram(aes(y = stat(density)),
                 fill = colour.hist.mu.phi, colour = "black",
                 alpha = 0.3, bins=60) +
      geom_vline(xintercept = mean.mu.phi,
                 colour = colour.prior.mu.phi, size = 1, lty = 2) +
      geom_vline(xintercept = mean(mu.phi),
                 colour = "blue", size = 1, lty = 2) +
      stat_function(fun = dlogis,
                 n = 101,
                 args = list(location = mean.mu.phi, scale = sd.mu.phi),
                 size = 1, colour = colour.prior.mu.phi) +
      scale_colour_manual(values = mycolours)+
      labs(x=expression(mu[phi]), y="Density", colour = "Density") +
      ggtitle(title.plot.mu.phi)+
      xlim(c(mu.phi.x.lim.low, mu.phi.x.lim.up)) +
      theme_bw()

    # densities mu AG vs ID...
    return(grid.arrange(diagnostic.plot1, diagnostic.plot2, ncol = 2, nrow = 1))
  } else {
    return(diagnostic.plot2)
  }

}