R/effect.hmr.R

In jarbes: Just a Rather Bayesian Evidence Synthesis

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

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


#' @title Posterior distribution of Effectiveness for a subgroup of patients
#'
#' @description This function estimates the posterior distribution for a subgroup of patients identified with the function hmr (Hierarchical Meta-Regression).
#'
#' @param object The object generated by the function hmr.
#'
#' @param B.lower Lower limit of bias correction. The default is 0 meaning no bias correction.
#' @param B.upper Upper limit of bias correction. The default is 3 meaning three times bias correction.
#' @param k       Covariable number indicating the subgroup.
#'
#' @param level Vector with the probability levels of the contour plot. The default values are: 0.5, 0.75, and 0.95.
#' @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 bias plot.

#' @param kde2d.n The number of grid points in each direction for the non-parametric density estimation. The default is 25.
#' @param marginals If TRUE the marginal histograms of the posteriors are added to the plot.
#'
#' @param bin.hist The number of bins in for the histograms. The default value is 30.
#' @param color.line The color of the contour lines. The default is "black.
#' @param color.hist The color of the histogram bars. The default is "white".
#' @param color.data.points The color of the data points. The default is "black".
#' @param alpha.data.points Transparency of the data points.
#'
#' @param S The number of sample values from the joint posterior distribution used to approximate the contours. The default is S=5000.
#' @param display.probability Logical, if TRUE the figure display probabilities.
#' @param line.no.effect Horizontal line used as reference for no effect.
#' @param font.size.title Font size of the title.
#'
#' @param ... \dots
#'
#'
#' @import ggplot2
#' @import ggExtra
#' @import MASS
#'
#' @export

effect.hmr = function(object,
                      # Range of the bias correction ...
                             B.lower = 0,
                             B.upper = 3,
                             k = 1,
                             # Parameters for the contour plot...
                             level = c(0.5, 0.75, 0.95),
                             x.lim = c(-9, 5),
                             y.lim = c(-1, 5),
                             x.lab = "Baseline risk",
                             y.lab = "Effectiveness",
                             title.plot = paste("Posterior Effectiveness for a subgroup (50%, 75% and 95%)"),
                             kde2d.n = 25,
                             marginals = TRUE,
                             bin.hist = 30,
                             color.line = "black",
                             color.hist = "white",
                             color.data.points = "black",
                             alpha.data.points = 0.1,
                             S = 5000,
                             display.probability = FALSE,
                             line.no.effect = 0,
                             font.size.title = 20,
                             ...) {

  x=y=ylo=yhi=kde2d=pi.bias=bias=dens.z=baseline=NULL

  # Data preparation
  alpha.0 = object$BUGSoutput$sims.list$alpha.0
  alpha.1 = object$BUGSoutput$sims.list$alpha.1
  mu.phi = object$BUGSoutput$sims.list$mu.phi
  mu.1 = object$BUGSoutput$sims.list$mu.1
  beta.K = object$BUGSoutput$sims.list$beta.IPD[,k]
  n.B = length(mu.1)
  B = runif(n.B, B.lower, B.upper)

  theta.1.B.star = mu.1 +(1-B)* mu.phi + beta.K
  theta.2.B.star = alpha.0 + alpha.1*(theta.1.B.star - mu.1)

  if(display.probability==TRUE){
    dat.post = data.frame(x = 1/(1+exp(-1*theta.1.B.star)),
                          y = exp(theta.2.B.star),
                          Bias.Correction = B)
  }
  else{
    dat.post = data.frame(x = theta.1.B.star,
                          y = theta.2.B.star,
                          Bias.Correction = B)
  }

  dat.post = dat.post[sample(1:S), ]
  cut.point = line.no.effect

  # Base plot ..................................................................
  # baseplot = ggplot(dat.post, aes(x = x, y = y, color = Bias.Correction)) +
  baseplot = ggplot(dat.post, aes(x = x, y = y)) +
    geom_point(size=0.01, color = color.data.points, alpha = alpha.data.points)+
    scale_x_continuous(name = x.lab, limits = x.lim) +
    scale_y_continuous(name = y.lab, limits = y.lim) +
    ggtitle(title.plot) +
    geom_hline(yintercept=cut.point, linetype="dashed", color = "black")+
    theme_bw()+
    theme(plot.title=element_text(size = font.size.title))

  # Non-parametric .............................................................
  # Estimation of the nonparametric density ...
  x.nopar = dat.post$x
  y.nopar = dat.post$y

  dens = kde2d(x.nopar, y.nopar, n = kde2d.n)
  dx = diff(dens$x[1:2])
  dy = diff(dens$y[1:2])
  sz = sort(dens$z)
  c1 = cumsum(sz) * dx * dy
  Levels.nonpar = approx(c1, sz, xout = 1 - level)$y

  densdf = data.frame(expand.grid(baseline = dens$x,
                                  effect = dens$y),
                      dens.z = as.vector(dens$z))

  # Non-parametric .............................................................
  finalplot = baseplot + geom_contour(data = densdf,
                                      aes(baseline, effect, z = dens.z),
                                      colour = color.line,
                                      breaks = Levels.nonpar,
                                      lwd =1)

  if(marginals==TRUE){
    p.effect = ggMarginal(finalplot,
                        type = "histogram",
                        fill = color.hist,
                        bins = bin.hist)}
  else{p.effect = finalplot}

  #.............................................................................
  return(p.effect)
}