R/plot.bcmeta.R

In jarbes: Just a Rather Bayesian Evidence Synthesis

#' Generic plot function for bcmeta object in jarbes.
#'
#' @param x The object generated by the bcmeta function.
#'
#' @param x.lim   Numeric vector of length 2 specifying the x-axis limits.
#' @param y.lim   Numeric vector of length 2 specifying the y-axis limits.
#' @param x.lab Text with the label of the x-axis.
#' @param y.lab Text with the label of the y-axis.
#' @param title.plot.1 Text for the posterior means by component (biased and bias corrected).
#' @param title.plot.2 Text for the posterior mean (pooled and predictive).
#' @param ... \dots
#'
#' @import ggplot2
#' @importFrom tidyr pivot_longer
#' @export

plot.bcmeta = function(x,
                       x.lim = c(-3, 3),
                       y.lim = c(0, 2),
                       x.lab = "Treatment Effect: log(OR)",
                       y.lab = "Posterior",
                       title.plot.1 = "Model Components",
                       title.plot.2 = "Bias Corrected Meta-Analysis",
                       ...) {
  object = x

  Component=Type=mu=NULL

  # Posterior: components .........................................

  mu = object$BUGSoutput$sims.list$mu
  mu = data.frame(mu)

  mu_names = c("Bias Corrected (BC)", "Biased")
  names(mu) = mu_names
  mu.groups = pivot_longer(data=mu, cols=mu_names)

  names(mu.groups) = c("Component", "mu")

  plot1 = ggplot(mu.groups, aes(x = mu, group = Component,
                                fill = Component, lty = Component))+
            geom_density(alpha=.2, lwd = 1) +
            scale_x_continuous(limits = x.lim ) +
            scale_y_continuous(limits = y.lim )+
            xlab(x.lab)+
            ylab(y.lab)+
            geom_vline(aes(xintercept=0), lty = 2)+
            ggtitle(title.plot.1)+
            theme_bw()

  # Posterior: pooled mean and predictive ..........................

  mu.wide = data.frame(mu = object$BUGSoutput$sims.list$mu[,1],
                       mu.new = object$BUGSoutput$sims.list$mu.new)

  mu.names = c("Pooled BC", "Predictive BC")
  names(mu.wide) = mu.names
  mu.long = pivot_longer(mu.wide, mu.names)
  names(mu.long) = c("Type", "mu")

  plot2 = ggplot(mu.long, aes(x = mu, group = Type,
                                fill = Type, lty = Type))+
    geom_density(alpha=.2, lwd = 1) +
    scale_x_continuous(limits = x.lim ) +
    scale_y_continuous(limits = y.lim )+
    xlab(x.lab)+
    ylab(y.lab)+
    geom_vline(aes(xintercept=0), lty = 2)+
    ggtitle(title.plot.2)+
    theme_bw()

  return(suppressWarnings(
            grid.arrange(plot1, plot2, ncol = 2, nrow = 1)
            ) )

  # mu = data.frame(mu.bc = mu[ , 1] , mu.bias = mu[ , 2])
  # p = ggplot(post.res, aes(x = mu.bc))+
  #     geom_histogram(aes(y=..density..), binwidth=.08,
  #                  colour="black", fill="white") +
  #     scale_x_continuous(limits = x.lim ) +
  #     scale_y_continuous(limits = y.lim )+
  #     xlab(x.lab)+
  #     ylab(y.lab)+
  #     geom_density(alpha=.2, lwd = 1) +
  #     geom_density(aes(x=mu.bias),  alpha=.2, lwd = 1, lty = 2) +
  #     geom_vline(aes(xintercept=0), lty = 2)+
  #     ggtitle(title.plot)+
  #   # geom_text(data = df.text, aes(x=x, y=y, label = study.type), size = 5)+
  #   # annotate(geom = "segment",
  #   #          x = -2.3, y = 1.3,
  #   #          xend = -1.5, yend = 0.85,
  #   #          arrow = arrow(length = unit(2, "mm")))+
  #   # annotate(geom = "segment",
  #   #          x = -1.8, y = 2,
  #   #          xend = -0.9, yend = 1.5,
  #   #          arrow = arrow(length = unit(2, "mm")))+
  #   theme_bw()
  #
  # return(suppressWarnings(print(p)))
  #
}