R/plot.metarisk.R

In jarbes: Just a Rather Bayesian Evidence Synthesis

#' Generic plot function for metarisk object in jarbes.
#'
#' @param x The object generated by the metarisk 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 Text for setting a title in the plot.
#'
#' @param ... \dots
#'
#' @import ggplot2
#'
#' @export

plot.metarisk = function(x,
                         x.lim = c(-5, 2.8),
                         y.lim = c(-2, 1.0),
                         x.lab = "Rate of The Control Group (logit scale)",
                         y.lab = "No improvement <- Treatment effect -> Improvement",
                         title.plot = "Treatment Effect Against Baseline Risk",
                         ...) {
  object = x

  x.line=median.hat=upper.hat=lower.hat=logitPc=TE=seTE=NULL

  # External validity
  a0.f = object$BUGSoutput$sims.list$beta.0 # Intercept centered at the mean baseline risk mu.1.
  b0.f = object$BUGSoutput$sims.list$beta.1 # Slope between baseline risk and relative treatment effect.
  mu.1 = object$BUGSoutput$sims.list$mu.1

  x = seq(x.lim[1], x.lim[2], length = 50)

  B = length(a0.f)

  y.f = rep(0, length(x)*B)
  dim(y.f) = c(length(x), B)

  y.f2 = rep(0, length(x)*B)
  dim(y.f2) = c(length(x), B)


  for(i in 1:length(x))
  {
    y.f[i, ] = a0.f + b0.f*x[i]
  }


  ###########################################################################
  dat.lines = data.frame(x.line = x,
                         median.hat = apply(y.f, 1, median),
                         upper.hat  = apply(y.f, 1, quantile, prob=0.95),
                         lower.hat  = apply(y.f, 1, quantile, prob=0.05))
  #############################################################################
  dat = object$data
  incr.e = 1
  incr.c = 1

  if (object$two.by.two == FALSE) {
    dat$yc = dat[,1]
    dat$nc = dat[,2]
    dat$yt = dat[,3]
    dat$nt = dat[,4]
  }

  n11 = dat$yt
  n12 = dat$yc

  n21 = dat$nt - dat$yt
  n22 = dat$nc - dat$yc


  dat$TE = log(((n11 + incr.e) * (n22 + incr.c)) /
                 ((n12 + incr.e) * (n21 + incr.c)))
  dat$seTE = sqrt((1/(n11 + incr.e) + 1/(n12 + incr.e) +
                     1/(n21 + incr.c) + 1/(n22 + incr.c)))

  Pc = ((n12 + incr.c)/(dat$nc + 2*incr.c))

  dat$logitPc = log(Pc/(1-Pc))

  rm(list=c("Pc", "n11","n12","n21","n22","incr.c", "incr.e"))

  dat.points = data.frame(TE = dat$TE,
                          seTE = dat$seTE,
                          logitPc = dat$logitPc,
                          N.total = dat$nt + dat$nc)


  p = ggplot(dat.lines, aes(x = x.line, y = median.hat)) +
    scale_x_continuous(name = x.lab, limits = x.lim) +
    scale_y_continuous(name = y.lab, limits = y.lim)+
    # Posterior intervals for the line  between TE and Baseline risk ...
    geom_line(aes(x = x.line, y = median.hat), colour = "red", size = 1.5)+
    geom_line(aes(x = x.line, y = upper.hat), colour = "blue", size = 1.5) +
    geom_line(aes(x = x.line, y = lower.hat), colour = "blue", size = 1.5) +

    # Reference line at no effect
    geom_hline(yintercept = 0, colour = "black", size = 0.5, lty = 2) +
    scale_size_area(max_size = 10) +

    # Forest plot for studies' results
    geom_pointrange(data = dat.points,
                    aes(x = logitPc,
                        y = TE,
                        ymin = TE - seTE,
                        ymax = TE + seTE),
                    lwd = 0.8, alpha = 0.5) +
    geom_vline(aes(xintercept = mean(mu.1)),size = 1, lty = 2) +
    ggtitle(title.plot)+
    theme_bw()

  return(suppressWarnings(print(p)))
  }