R/plot.hmr.R

In jarbes: Just a Rather Bayesian Evidence Synthesis

#' Generic plot function for hmr object in jarbes.
#'
#' @param x The object generated by the hmr 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 AD.colour Colour of the location of the baseline risk of the aggregated data AD
#' @param IPD.colour Colour of the location of the baseline risk of the individual participant data (IPD) data
#' @param Study.Types Vector of text for the label of the study types
#'
#' @param ... \dots
#'
#' @export
#'
plot.hmr = function(x,
                    x.lim = c(-5, 2.8),
                    y.lim = c(-2, 1.0),
                    x.lab = "Event rate of The Control Group (logit scale)",
                    y.lab = "No improvement <- Effectiveness -> Improvement",
                    title.plot = "HMR: Effectiveness Against Baseline Risk",
                    AD.colour = "red",
                    IPD.colour = "blue",
                    Study.Types = c("AD-RCTs", "IPD-RWD"),
                    ...) {

  x.axis=x.end=Design=upper.hat=lower.hat=logitPc=TE=seTE=NULL
  x.line=median.hat=y.axis=y.end=NULL
  object = x

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

  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)

    mu.phi = object$BUGSoutput$sims.list$mu.phi
      mu.1 = object$BUGSoutput$sims.list$mu.1

  X.baseline = c(mean(mu.1), mean(mu.1 + mu.phi))

  vlines = data.frame(X.baseline = X.baseline,
                      Study.Types = Study.Types)

  p = ggplot(dat.lines, aes(x = x.line, y = median.hat)) +
    # Posterior intervals for the line  between TE and Baseline risk ...
        geom_line(aes(x = x.line, y = median.hat), colour = "black", size = 1.5) +
        geom_line(aes(x = x.line, y = upper.hat), colour = "black", size = 1.5) +
        geom_line(aes(x = x.line, y = lower.hat), colour = "black", size = 1.5) +
        scale_x_continuous(name = x.lab, limits = x.lim) +
        scale_y_continuous(name = y.lab, limits = y.lim)+
    # Reference horizontal line at no effect ...
        geom_hline(yintercept = 0, colour = "black", size = 0.5, lty = 2) +
    # Location of baselines AD and IPD
        geom_vline(data = vlines,
               aes(xintercept = X.baseline,
                   colour = Study.Types,
                   lty = Study.Types),
                   size = 1.5) +
    scale_color_manual(values= c(AD.colour, IPD.colour))+
    # 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.25) +
        ggtitle(title.plot)+ theme_bw()


  return(suppressWarnings(print(p)))
}