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 names Add IPD names to the plot.
#'
#' @param name.side Set the side of each name in the plot relative to the vertical line.
#'
#' @param ... \dots
#'
#' @export
#'
plot.hmr = 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",
                    names=NULL,
                    name.side=NULL, ...) {

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

  if (missing(names) & !missing(name.side))stop("name.side must have same length as names.")

  if (missing(names)) names = c()
  if (missing(name.side)) {
    name.side = c()
    for (n in names) {
      name.side = c(name.side, "right")
    }
  }

  if ((length(names) != length(name.side)))stop("names and name.side must have the same length.")


  object = x

  # 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.

  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
  beta.IPD = object$BUGSoutput$sims.list$beta.IPD
  colnames(beta.IPD) = object$beta.names

  x.axis = c(mean(mu.1) + 0.7,
             mean(mu.1 + mu.phi) + 0.7)
  x.end = c(mean(mu.1),
            mean(mu.1 + mu.phi))
  for (i in 1:length(names)) {
    n = paste0("beta.", names[i])
    if (n %in% colnames(beta.IPD)) {
      if (name.side[i] == "left") offset = -1
      else offset = 0.7
      x.axis = c(x.axis, mean(mu.1 - mu.phi + beta.IPD[,n] + offset))
      x.end = c(x.end, mean(mu.1 - mu.phi + beta.IPD[,n]))
    } else {
      warning(paste("name", names[i], "not found in IPD names"))
    }
  }

  vlines = data.frame(x.axis = x.axis,
                      x.end = x.end,
                      y.axis = -1.3,
                      y.end = -1.8,
                      text = c("RCT", "OS", names),
                      Design = c("RCT", "OS", rep("RCT+OS", length(names))))


  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)+
        geom_text(data = vlines,           ## Text for baseline risk groups
              aes(label = text,
                  x = x.axis ,
                  y = y.axis )) +
        geom_segment(data = vlines,
                 aes(x = x.axis ,       ## Start of arrow on x axis
                     xend = x.end  ,    ## End of arrow on x axis
                     y = y.axis - 0.1,  ## Start of arrow on y axis
                     yend = y.end ),    ## End of arrow on y axis
                 size = 1,
                 arrow = arrow(length = unit(0.3, "cm")) ) +
    # Location of  groups' baseline
    geom_vline(data = vlines, aes(xintercept = x.end,
                                  colour = Design, lty = Design),
               size = 1.5) +
    scale_color_manual(values=c("red", "blue","black")) +
    # 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) +
    # Refrence 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) +
    ggtitle(title.plot)+
    theme_bw()

  return(suppressWarnings(print(p)))
}