R/ppc.R

Defines functions ppc

Documented in ppc

#' Posterior predictive checks for assessing the fit to the observed data of Bayesian models implemented in \code{JAGS} using the function \code{\link{selection}}, \code{\link{pattern}} or \code{\link{hurdle}}
#'
#' The focus is restricted to full Bayesian models in cost-effectiveness analyses based on the function \code{\link{selection}}, \code{\link{pattern}} and \code{\link{hurdle}}, 
#' with the fit to the observed data being assessed through graphical checks based on the posterior replications generated from the model. 
#' Examples include the comparison of histograms, density plots, intervals, test statistics, evaluated using both the observed and replicated data.
#' Different types of posterior predictive checks are implemented to assess model fit using functions contained in the package \strong{bayesplot}. 
#' Graphics and plots are managed using functions contained in the package \strong{ggplot2} and \strong{ggthemes}.
#' @keywords posterior predictive checks MCMC replications 
#' @param x An object of class "missingHE" containing the posterior results of a full Bayesian model implemented using the function \code{\link{selection}}, 
#' \code{\link{pattern}} or \code{\link{hurdle}}.
#' @param type Type of posterior predictive check to be plotted for assessing model fit. Available choices include: 'histogram', 'boxplot', 
#' 'freqpoly', 'dens', 'dens_overlay' and ecdf_overlay', which compare the empirical and repicated distributions of the data; 'stat' and 'stat_2d', which compare
#' the value of some statistics evaluated on the observed data with the replicated values for those statistics from the posterior predictions; 'error_hist', 
#' 'error_scatter', 'error_scatter_avg' and 'error_binned', which display the predictive errors of the model; 'intervals' and 'ribbon', which compare medians and
#' central interval estmates of the replications with the observed data overlaid; 'scatter' and 'scatter_avg', which display scatterplots of the observed and replicated data.   
#' @param outcome The outcome variables that should be displayed. Use the names 'effects_arm1' and effects_arm2' for the effectiveness in the control and intervention arm; 
#' use costs_arm1' or 'costs_arm2' for the costs; use "effects" or "costs" for the respective outcome in both arms; use "all" for all outcomes.
#' @param ndisplay Number of posterior replications to be displayed in the plots.
#' @param theme Type of ggplot theme among some pre-defined themes, mostly taken from the package \strong{ggthemes}. For a full list of available themes see details.
#' @param scheme_set Type of scheme sets among some pre-defined schemes, mostly taken from the package \strong{bayesplot}. For a full list of available themes see details.
#' @param legend Position of the legend: available choices are: "top", "left", "right", "bottom" and "none".
#' @param ... Additional parameters that can be provided to manage the output of \code{ppc}. For more details see \strong{bayesplot}.  
#' @return A \code{ggplot} object containing the plots specified in the argument \code{type}.
#' @seealso \code{\link{selection}} \code{\link{pattern}} \code{\link{hurdle}} \code{\link{diagnostic}}
#' @details The funciton produces different types of graphical posterior predictive checks using the estimates from a Bayesian cost-effectiveness model implemented
#' with the function \code{\link{selection}}, \code{\link{pattern}} or \code{\link{hurdle}}. The purpose of these checks is to visually compare the distribution (or some relevant quantity)
#' of the observed data with respect to that from the replicated data for both effectiveness and cost outcomes in each treatment arm. Since predictive checks are meaningful 
#' only with respect to the observed data, only the observed outcome values are used to assess the fit of the model.
#' The arguments \code{theme} and \code{scheme_set} allow to customise the graphical aspect of the plots generated by \code{ppc} and allow to choose among a set of possible 
#' pre-defined themes and scheme sets taken form the package \strong{ggtheme} and \code{bayesplot}. For a complete list of the available character names for each theme and scheme set, see \strong{ggthemes} and \strong{bayesplot}.
#' @author Andrea Gabrio
#' @references 
#' Gelman, A. Carlin, JB., Stern, HS. Rubin, DB.(2003). \emph{Bayesian Data Analysis, 2nd edition}, CRC Press.
#' @import ggplot2 bayesplot ggpubr
#' @importFrom stats rnorm rbeta rgamma rlnorm rweibull rnbinom rbinom rpois rlogis rexp complete.cases
#' @export 
#' @examples
#' # For examples see the function \code{\link{selection}}, 
#' # \code{\link{pattern}} or \code{\link{hurdle}}
#' #

ppc <- function(x, type = "histogram", outcome = "all", ndisplay = 15, theme = NULL, scheme_set = NULL, legend = "top", ...) {
  if(!isTRUE(requireNamespace("ggplot2")) | !isTRUE(requireNamespace("ggthemes")) | !isTRUE(requireNamespace("bayesplot")) | !isTRUE(requireNamespace("ggpubr"))) {
    stop("You need to install the R packages 'ggplot2' and 'bayesplot'. Please run in your R terminal:\n install.packages('ggplot2', 'ggthemes', 'bayesplot', 'ggpubr')")
  }
  if(x$model_output$ppc == FALSE) {
    stop("No posterior estimates that can be used to generate replications are detected. Try running the function 'selection', 'pattern' or 'hurdle' and set the argument ppc = TRUE")
  }
  if(is.null(theme) == TRUE) {theme = "default" }
  if(!theme %in% c("default", "base", "calc", "economist", "excel", "few", "538", "gdocs", "hc", "par", "pander", "solarized", "stata", "tufte", "wsj")) {
    stop("You must provide a valid theme type among those available")
  }
  if(is.null(scheme_set) == TRUE) {scheme_set = "blue" }
  if(!scheme_set %in% c("blue", "brightblue", "gray", "darkgray", "green", "pink", "purple", "red", "teal", "yellow", "viridis", "viridisA", "viridisB", "viridisC",
                        "viridisE", "mix-x-y")) {
    stop("You must provide a valid scheme sets among those available")
  }
  if(!type %in% c("histogram", "boxplot", "freqpoly", "dens", "dens_overlay", "ecdf_overlay", "stat", "stat_2d", "error_hist", 
                  "error_scatter", "error_scatter_avg", "error_binned", "intervals", "ribbon", "scatter", "scatter_avg")) {
    stop("You must provide a plot type among those available")
  }
  if(!outcome %in% c("effects_arm1", "costs_arm1", "effects_arm2", "costs_arm2", "effects", "costs", "all")) {
    stop("You must provide a plot outcome among those available")
  }
  if(is.numeric(ndisplay) == FALSE |  ndisplay < 1 | ndisplay > dim(x$model_output$`model summary`$BUGSoutput$sims.list$mu_e1)[1]) {
    stop("number of replications to display is not valid")
  }
  if(!legend %in% c("top", "right", "left", "bottom", "none")) {
    stop("Please provide a valid name for 'legend' to indicate the position of the legend")
  }
  exArgs <- list(...)
  if(class(x) != "missingHE") {
    stop("Only objects of class 'missingHE' can be used")
  }
  if(exists("corr", where = exArgs)) {corr = exArgs$corr} else {corr = NULL }
  if(!is.null(corr) & outcome != "all" | !is.null(corr) & !is.function(corr)) {
    stop("The argument 'corr' must be a function and is available only if 'all' outcomes selected ")
  }
  if(exists("stat", where = exArgs)) {stat = exArgs$stat} else {stat = NULL }
  if(exists("bpp", where = exArgs)) {bpp = exArgs$bpp} else {bpp = FALSE }
  if(bpp == TRUE & type == "stat_2d") { stop("Bayesian p-values not available for 'stat_2d' type")}
  if(type == "stat" & is.null(corr)) {
  if(is.null(stat) == TRUE | length(stat) !=1 | !is.function(stat)) {
    stop("Please use the argument 'stat' to provide a single function for computing statistics on replications")
   } 
  }
  if(type == "stat_2d" & is.null(stat) == TRUE | type == "stat_2d" & length(stat) !=2 | type == "stat_2d" & any(sapply(stat, is.function)) == FALSE) {
    stop("Please use the argument 'stat' to provide a vector of two functions for computing statistics on replications")
  }
  index_mis_e1 <- which(is.na(x$data_set$effects$Control))
  index_mis_e2 <- which(is.na(x$data_set$effects$Intervention))
  index_mis_c1 <- which(is.na(x$data_set$costs$Control))
  index_mis_c2 <- which(is.na(x$data_set$costs$Intervention))
  if(length(index_mis_e1) == 0 | length(index_mis_e2) == 0) {
    stop("Missing values not found for the effects in both arms")
  } 
  if(length(index_mis_c1) == 0 | length(index_mis_c2) == 0) {
    stop("Missing values not found for the costs in both arms")
  } 
  e1_obs <- as.numeric(x$data_set$effects$Control[-index_mis_e1])
  e2_obs <- as.numeric(x$data_set$effects$Intervention[-index_mis_e2])
  c1_obs <- as.numeric(x$data_set$costs$Control[-index_mis_c1])
  c2_obs <- as.numeric(x$data_set$costs$Intervention[-index_mis_c2])
  n1_e <- length(e1_obs)
  n1_c <- length(c1_obs)
  n2_e <- length(e2_obs)
  n2_c <- length(c2_obs)
  e1_rep <- matrix(NA, nrow = x$model_output$`model summary`$BUGSoutput$n.iter, n1_e)
  c1_rep <- matrix(NA, nrow = x$model_output$`model summary`$BUGSoutput$n.iter, n1_c)
  e2_rep <- matrix(NA, nrow = x$model_output$`model summary`$BUGSoutput$n.iter, n2_e)
  c2_rep <- matrix(NA, nrow = x$model_output$`model summary`$BUGSoutput$n.iter, n2_c)
   if(x$model_output$dist_e == "norm") {
   mu_e1 <- x$model_output$`model summary`$BUGSoutput$sims.list$mu_e1[, -index_mis_e1]
   mu_e2 <- x$model_output$`model summary`$BUGSoutput$sims.list$mu_e2[, -index_mis_e2]
   if(is.null(dim(mu_e1)) == TRUE | is.null(dim(mu_e2)) == TRUE){
     stop("Not possible to plot posterior predictive checks with a single observed effect value in one arm.")
   }
    if(grepl("SELECTION", x$model_output$type) == TRUE){
     sigma_e1 <- 1 / sqrt(x$model_output$`model summary`$BUGSoutput$sims.list$tau_e[, 1])
     sigma_e2 <- 1 / sqrt(x$model_output$`model summary`$BUGSoutput$sims.list$tau_e[, 2])
     for(i in 1 : nrow(e1_rep)){ e1_rep[i, ] <- rnorm(n1_e, mean = mu_e1[i, ], sd = sigma_e1 [i]) } 
     for(i in 1 : nrow(e2_rep)){ e2_rep[i, ] <- rnorm(n2_e, mean = mu_e2[i, ], sd = sigma_e2 [i]) }
    } else if(grepl("PATTERN", x$model_output$type) == TRUE) {
     sigma_e1 <- 1 / sqrt(x$model_output$`model summary`$BUGSoutput$sims.list$tau_e_p1)
     sigma_e2 <- 1 / sqrt(x$model_output$`model summary`$BUGSoutput$sims.list$tau_e_p2)
     for(i in 1 : nrow(e1_rep)){ e1_rep[i, ] <- rnorm(n1_e, mean = mu_e1[i, ], sd = sigma_e1 [i, ]) }
     for(i in 1 : nrow(e2_rep)){ e2_rep[i, ] <- rnorm(n2_e, mean = mu_e2[i, ], sd = sigma_e2 [i, ]) }
    } else if(grepl("HURDLE", x$model_output$type) == TRUE) {
     sigma_e1 <- 1 / sqrt(x$model_output$`model summary`$BUGSoutput$sims.list$tau_e1)
     sigma_e2 <- 1 / sqrt(x$model_output$`model summary`$BUGSoutput$sims.list$tau_e2)
     for(i in 1 : nrow(e1_rep)){ e1_rep[i, ] <- rnorm(n1_e, mean = mu_e1[i, ], sd = sigma_e1 [i, ]) }
     for(i in 1 : nrow(e2_rep)){ e2_rep[i, ] <- rnorm(n2_e, mean = mu_e2[i, ], sd = sigma_e2 [i, ]) }
    }
   } else if(x$model_output$dist_e == "logis") {
     mu_e1 <- x$model_output$`model summary`$BUGSoutput$sims.list$mu_e1[, -index_mis_e1]
     mu_e2 <- x$model_output$`model summary`$BUGSoutput$sims.list$mu_e2[, -index_mis_e2]
     if(is.null(dim(mu_e1)) == TRUE | is.null(dim(mu_e2)) == TRUE){
       stop("Not possible to plot posterior predictive checks with a single observed effect value in one arm.")
     }
      if(grepl("SELECTION", x$model_output$type) == TRUE){
       scale_e1 <- 1 / x$model_output$`model summary`$BUGSoutput$sims.list$tau_e[, 1]
       scale_e2 <- 1 / x$model_output$`model summary`$BUGSoutput$sims.list$tau_e[, 2]
       for(i in 1 : nrow(e1_rep)){ e1_rep[i, ] <- rlogis(n1_e, location = mu_e1[i, ], scale = scale_e1 [i]) } 
       for(i in 1 : nrow(e2_rep)){ e2_rep[i, ] <- rlogis(n2_e, location = mu_e2[i, ], scale = scale_e2 [i]) }
      } else if(grepl("PATTERN", x$model_output$type) == TRUE) {
       scale_e1 <- 1 / x$model_output$`model summary`$BUGSoutput$sims.list$tau_e_p1
       scale_e2 <- 1 / x$model_output$`model summary`$BUGSoutput$sims.list$tau_e_p2
       for(i in 1 : nrow(e1_rep)){ e1_rep[i, ] <- rlogis(n1_e, location = mu_e1[i, ], scale = scale_e1 [i, ]) }
       for(i in 1 : nrow(e2_rep)){ e2_rep[i, ] <- rlogis(n2_e, location = mu_e2[i, ], scale = scale_e2 [i, ]) }
      } else if(grepl("HURDLE", x$model_output$type) == TRUE) {
       scale_e1 <- 1 / x$model_output$`model summary`$BUGSoutput$sims.list$tau_e1
       scale_e2 <- 1 / x$model_output$`model summary`$BUGSoutput$sims.list$tau_e2
       for(i in 1 : nrow(e1_rep)){ e1_rep[i, ] <- rlogis(n1_e, location = mu_e1[i, ], scale = scale_e1 [i, ]) }
       for(i in 1 : nrow(e2_rep)){ e2_rep[i, ] <- rlogis(n2_e, location = mu_e2[i, ], scale = scale_e2 [i, ]) }
      }
   } else if(x$model_output$dist_e == "nbinom") {
     mu_e1 <- x$model_output$`model summary`$BUGSoutput$sims.list$mu_e1[, -index_mis_e1]
     mu_e2 <- x$model_output$`model summary`$BUGSoutput$sims.list$mu_e2[, -index_mis_e2]
     if(is.null(dim(mu_e1)) == TRUE | is.null(dim(mu_e2)) == TRUE){
       stop("Not possible to plot posterior predictive checks with a single observed effect value in one arm.")
     }
     if(grepl("SELECTION", x$model_output$type) == TRUE){
       prob_e1 <- x$model_output$`model summary`$BUGSoutput$sims.list$tau_e[, 1] / (x$model_output$`model summary`$BUGSoutput$sims.list$tau_e[, 1] + x$model_output$`model summary`$BUGSoutput$sims.list$mu_e1[, -index_mis_e1])
       prob_e2 <- x$model_output$`model summary`$BUGSoutput$sims.list$tau_e[, 2] / (x$model_output$`model summary`$BUGSoutput$sims.list$tau_e[, 2] + x$model_output$`model summary`$BUGSoutput$sims.list$mu_e2[, -index_mis_e2])
       size_e1 <- x$model_output$`model summary`$BUGSoutput$sims.list$tau_e[, 1]
       size_e2 <- x$model_output$`model summary`$BUGSoutput$sims.list$tau_e[, 2]
       for(i in 1 : nrow(e1_rep)){ e1_rep[i, ] <- rnbinom(n1_e, prob = prob_e1[i, ], size = size_e1 [i]) } 
       for(i in 1 : nrow(e2_rep)){ e2_rep[i, ] <- rnbinom(n2_e, prob = prob_e2[i, ], size = size_e2 [i]) }
     } else if(grepl("PATTERN", x$model_output$type) == TRUE) {
       prob_e1 <- array(NA, dim = c(nrow(e1_rep), dim(mu_e1)[2], dim(x$model_output$`model summary`$BUGSoutput$sims.list$tau_e_p1)[2]))
        for(j in 1 : dim(x$model_output$`model summary`$BUGSoutput$sims.list$tau_e_p1)[2]) {
          prob_e1[, , j] <- x$model_output$`model summary`$BUGSoutput$sims.list$tau_e_p1[, j] / (x$model_output$`model summary`$BUGSoutput$sims.list$tau_e_p1[, j] + x$model_output$`model summary`$BUGSoutput$sims.list$mu_e1[, -index_mis_e1])
        }
       prob_e2 <- array(NA, dim = c(nrow(e2_rep), dim(mu_e2)[2], dim(x$model_output$`model summary`$BUGSoutput$sims.list$tau_e_p2)[2]))
       for(j in 1 : dim(x$model_output$`model summary`$BUGSoutput$sims.list$tau_e_p2)[2]) {
         prob_e2[, , j] <- x$model_output$`model summary`$BUGSoutput$sims.list$tau_e_p2[, j] / (x$model_output$`model summary`$BUGSoutput$sims.list$tau_e_p2[, j] + x$model_output$`model summary`$BUGSoutput$sims.list$mu_e2[, -index_mis_e2])
       }
       size_e1 <- x$model_output$`model summary`$BUGSoutput$sims.list$tau_e_p1
       size_e2 <- x$model_output$`model summary`$BUGSoutput$sims.list$tau_e_p2
       for(i in 1 : nrow(e1_rep)){ e1_rep[i, ] <- rnbinom(n1_e, prob = prob_e1[i, , ], size = size_e1 [i, ]) } 
       for(i in 1 : nrow(e2_rep)){ e2_rep[i, ] <- rnbinom(n2_e, prob = prob_e2[i, , ], size = size_e2 [i, ]) }
     } else if(grepl("HURDLE", x$model_output$type) == TRUE) {
       if(x$model_output$type %in% c("HURDLE_e", "HURDLE_ec")){
         prob_e1 <- array(NA, dim = c(nrow(e1_rep), dim(mu_e1)[2], dim(x$model_output$`model summary`$BUGSoutput$sims.list$tau_e1)[2]))
         for(j in 1 : dim(x$model_output$`model summary`$BUGSoutput$sims.list$tau_e1)[2]) {
           prob_e1[, , j] <- x$model_output$`model summary`$BUGSoutput$sims.list$tau_e1[, j] / (x$model_output$`model summary`$BUGSoutput$sims.list$tau_e1[, j] + x$model_output$`model summary`$BUGSoutput$sims.list$mu_e1[, -index_mis_e1])
         }
         prob_e2 <- array(NA, dim = c(nrow(e2_rep), dim(mu_e2)[2], dim(x$model_output$`model summary`$BUGSoutput$sims.list$tau_e2)[2]))
         for(j in 1 : dim(x$model_output$`model summary`$BUGSoutput$sims.list$tau_e2)[2]) {
           prob_e2[, , j] <- x$model_output$`model summary`$BUGSoutput$sims.list$tau_e2[, j] / (x$model_output$`model summary`$BUGSoutput$sims.list$tau_e2[, j] + x$model_output$`model summary`$BUGSoutput$sims.list$mu_e2[, -index_mis_e2])
         }
         size_e1 <- x$model_output$`model summary`$BUGSoutput$sims.list$tau_e1
         size_e2 <- x$model_output$`model summary`$BUGSoutput$sims.list$tau_e2
         for(i in 1 : nrow(e1_rep)){ e1_rep[i, ] <- rnbinom(n1_e, prob = prob_e1[i, , ], size = size_e1 [i, ]) } 
         for(i in 1 : nrow(e2_rep)){ e2_rep[i, ] <- rnbinom(n2_e, prob = prob_e2[i, , ], size = size_e2 [i, ]) }
       } else if(x$model_output$type %in% c("HURDLE_c")){
       prob_e1 <- as.vector(x$model_output$`model summary`$BUGSoutput$sims.list$tau_e1) / (as.vector(x$model_output$`model summary`$BUGSoutput$sims.list$tau_e1) + x$model_output$`model summary`$BUGSoutput$sims.list$mu_e1[, -index_mis_e1])
       prob_e2 <- as.vector(x$model_output$`model summary`$BUGSoutput$sims.list$tau_e2) / (as.vector(x$model_output$`model summary`$BUGSoutput$sims.list$tau_e2) + x$model_output$`model summary`$BUGSoutput$sims.list$mu_e2[, -index_mis_e2])
       size_e1 <- as.vector(x$model_output$`model summary`$BUGSoutput$sims.list$tau_e1)
       size_e2 <- as.vector(x$model_output$`model summary`$BUGSoutput$sims.list$tau_e2)
       for(i in 1 : nrow(e1_rep)){ e1_rep[i, ] <- rnbinom(n1_e, prob = prob_e1[i, ], size = size_e1 [i]) } 
       for(i in 1 : nrow(e2_rep)){ e2_rep[i, ] <- rnbinom(n2_e, prob = prob_e2[i, ], size = size_e2 [i]) }
       }
     }
   } else if(x$model_output$dist_e == "exp") {
     mu_e1 <- x$model_output$`model summary`$BUGSoutput$sims.list$mu_e1[, -index_mis_e1]
     mu_e2 <- x$model_output$`model summary`$BUGSoutput$sims.list$mu_e2[, -index_mis_e2]
     if(is.null(dim(mu_e1)) == TRUE | is.null(dim(mu_e2)) == TRUE){
       stop("Not possible to plot posterior predictive checks with a single observed effect value in one arm.")
     }
     rate_e1 <- 1 / mu_e1
     rate_e2 <- 1 / mu_e2
     for(i in 1 : nrow(e1_rep)){ e1_rep[i, ] <- rexp(n1_e, rate  = rate_e1[i, ]) }
     for(i in 1 : nrow(e2_rep)){ e2_rep[i, ] <- rexp(n2_e, rate  = rate_e2[i, ]) }
   } else if(x$model_output$dist_e == "bern" | x$model_output$dist_e == "pois") {
     mu_e1 <- x$model_output$`model summary`$BUGSoutput$sims.list$mu_e1[, -index_mis_e1]
     mu_e2 <- x$model_output$`model summary`$BUGSoutput$sims.list$mu_e2[, -index_mis_e2]
      if(is.null(dim(mu_e1)) == TRUE | is.null(dim(mu_e2)) == TRUE){
       stop("Not possible to plot posterior predictive checks with a single observed effect value in one arm.")
     }
      if(x$model_output$dist_e == "bern"){
       for(i in 1 : nrow(e1_rep)){ e1_rep[i, ] <- rbinom(n1_e, size = 1, prob = mu_e1[i, ]) }
       for(i in 1 : nrow(e2_rep)){ e2_rep[i, ] <- rbinom(n2_e, size = 1, prob = mu_e2[i, ]) }
     }
      if(x$model_output$dist_e == "pois"){
       for(i in 1 : nrow(e1_rep)){ e1_rep[i, ] <- rpois(n1_e, lambda = mu_e1[i, ]) }
       for(i in 1 : nrow(e2_rep)){ e2_rep[i, ] <- rpois(n2_e, lambda = mu_e2[i, ]) }
     }
   } else if(x$model_output$dist_e == "gamma") {
     shape_e1 <- x$model_output$`model summary`$BUGSoutput$sims.list$mu_e1[, -index_mis_e1] * x$model_output$`model summary`$BUGSoutput$sims.list$tau_e1[, -index_mis_e1]
     rate_e1 <-  x$model_output$`model summary`$BUGSoutput$sims.list$tau_e1[, -index_mis_e1]
     shape_e2 <- x$model_output$`model summary`$BUGSoutput$sims.list$mu_e2[, -index_mis_e2] * x$model_output$`model summary`$BUGSoutput$sims.list$tau_e2[, -index_mis_e2]
     rate_e2 <-  x$model_output$`model summary`$BUGSoutput$sims.list$tau_e2[, -index_mis_e2]
     if(is.null(dim(shape_e1)) == TRUE | is.null(dim(shape_e2)) == TRUE){
       stop("Not possible to plot posterior predictive checks with a single observed effect value in one arm.")
     }
     for(i in 1 : nrow(e1_rep)){ e1_rep[i, ] <- rgamma(n1_e, shape = shape_e1[i, ], rate = rate_e1 [i, ]) }
     for(i in 1 : nrow(e2_rep)){ e2_rep[i, ] <- rgamma(n2_e, shape = shape_e2[i, ], rate = rate_e2 [i, ]) }
   } else if(x$model_output$dist_e == "weibull") {
     shape_e1 <- x$model_output$`model summary`$BUGSoutput$sims.list$tau_e1[, -index_mis_e1]
     shape_e2 <- x$model_output$`model summary`$BUGSoutput$sims.list$tau_e2[, -index_mis_e2]
     scale_e1 <- x$model_output$`model summary`$BUGSoutput$sims.list$mu_e1[, -index_mis_e1] / exp(lgamma(1 + 1 / x$model_output$`model summary`$BUGSoutput$sims.list$tau_e1[, -index_mis_e1]))
     scale_e2 <- x$model_output$`model summary`$BUGSoutput$sims.list$mu_e2[, -index_mis_e2] / exp(lgamma(1 + 1 / x$model_output$`model summary`$BUGSoutput$sims.list$tau_e2[, -index_mis_e2]))
     if(is.null(dim(shape_e1)) == TRUE | is.null(dim(shape_e2)) == TRUE){
       stop("Not possible to plot posterior predictive checks with a single observed effect value in one arm.")
     }
     for(i in 1 : nrow(e1_rep)){ e1_rep[i, ] <- rweibull(n1_c, shape = shape_e1[i, ], scale = scale_e1 [i, ]) }
     for(i in 1 : nrow(e2_rep)){ e2_rep[i, ] <- rweibull(n2_c, shape = shape_e2[i, ], scale = scale_e2 [i, ]) }
   } else if(x$model_output$dist_e == "beta") {
    shape1_e1 <- x$model_output$`model summary`$BUGSoutput$sims.list$mu_e1[, -index_mis_e1] * x$model_output$`model summary`$BUGSoutput$sims.list$tau_e1[, -index_mis_e1]
    shape2_e1 <- (1 - x$model_output$`model summary`$BUGSoutput$sims.list$mu_e1[, -index_mis_e1]) * x$model_output$`model summary`$BUGSoutput$sims.list$tau_e1[, -index_mis_e1]
    shape1_e2 <- x$model_output$`model summary`$BUGSoutput$sims.list$mu_e2[, -index_mis_e2] * x$model_output$`model summary`$BUGSoutput$sims.list$tau_e2[, -index_mis_e2]
    shape2_e2 <- (1 - x$model_output$`model summary`$BUGSoutput$sims.list$mu_e2[, -index_mis_e2]) * x$model_output$`model summary`$BUGSoutput$sims.list$tau_e2[, -index_mis_e2]
    if(is.null(dim(shape1_e1)) == TRUE | is.null(dim(shape1_e2)) == TRUE){
      stop("Not possible to plot posterior predictive checks with a single observed effect value in one arm.")
    }
    for(i in 1 : nrow(e1_rep)){ e1_rep[i, ] <- rbeta(n1_e, shape1 = shape1_e1[i, ], shape2 = shape2_e1 [i, ]) }
    for(i in 1 : nrow(e2_rep)){ e2_rep[i, ] <- rbeta(n2_e, shape1 = shape1_e2[i, ], shape2 = shape2_e2 [i, ]) }
  }
  if(x$model_output$dist_c == "norm") {
    mu_c1 <- x$model_output$`model summary`$BUGSoutput$sims.list$mu_c1[, -index_mis_c1]
    mu_c2 <- x$model_output$`model summary`$BUGSoutput$sims.list$mu_c2[, -index_mis_c2]
    if(is.null(dim(mu_c1)) == TRUE | is.null(dim(mu_c2)) == TRUE){
      stop("Not possible to plot posterior predictive checks with a single observed cost value in one arm.")
    }
    if(grepl("SELECTION", x$model_output$type) == TRUE){
      sigma_c1 <- 1 / sqrt(x$model_output$`model summary`$BUGSoutput$sims.list$tau_c[, 1])
      sigma_c2 <- 1 / sqrt(x$model_output$`model summary`$BUGSoutput$sims.list$tau_c[, 2])
      for(i in 1 : nrow(c1_rep)){ c1_rep[i, ] <- rnorm(n1_c, mean = mu_c1[i, ], sd = sigma_c1 [i]) }
      for(i in 1 : nrow(c2_rep)){ c2_rep[i, ] <- rnorm(n2_c, mean = mu_c2[i, ], sd = sigma_c2 [i]) }
    } else if(grepl("PATTERN", x$model_output$type) == TRUE) {
      sigma_c1 <- 1 / sqrt(x$model_output$`model summary`$BUGSoutput$sims.list$tau_c_p1)
      sigma_c2 <- 1 / sqrt(x$model_output$`model summary`$BUGSoutput$sims.list$tau_c_p2)
      for(i in 1 : nrow(c1_rep)){ c1_rep[i, ] <- rnorm(n1_c, mean = mu_c1[i, ], sd = sigma_c1 [i, ]) }
      for(i in 1 : nrow(c2_rep)){ c2_rep[i, ] <- rnorm(n2_c, mean = mu_c2[i, ], sd = sigma_c2 [i, ]) }
    } else if(grepl("HURDLE", x$model_output$type) == TRUE) {
      sigma_c1 <- 1 / sqrt(x$model_output$`model summary`$BUGSoutput$sims.list$tau_c1)
      sigma_c2 <- 1 / sqrt(x$model_output$`model summary`$BUGSoutput$sims.list$tau_c2)
      for(i in 1 : nrow(c1_rep)){ c1_rep[i, ] <- rnorm(n1_c, mean = mu_c1[i, ], sd = sigma_c1 [i, ]) }
      for(i in 1 : nrow(c2_rep)){ c2_rep[i, ] <- rnorm(n2_c, mean = mu_c2[i, ], sd = sigma_c2 [i, ]) }
    }
  } else if(x$model_output$dist_c == "gamma") {
    shape_c1 <- x$model_output$`model summary`$BUGSoutput$sims.list$mu_c1[, -index_mis_c1] * x$model_output$`model summary`$BUGSoutput$sims.list$tau_c1[, -index_mis_c1]
    rate_c1 <-  x$model_output$`model summary`$BUGSoutput$sims.list$tau_c1[, -index_mis_c1]
    shape_c2 <- x$model_output$`model summary`$BUGSoutput$sims.list$mu_c2[, -index_mis_c2] * x$model_output$`model summary`$BUGSoutput$sims.list$tau_c2[, -index_mis_c2]
    rate_c2 <-  x$model_output$`model summary`$BUGSoutput$sims.list$tau_c2[, -index_mis_c2]
    if(is.null(dim(shape_c1)) == TRUE | is.null(dim(shape_c2)) == TRUE){
      stop("Not possible to plot posterior predictive checks with a single observed cost value in one arm.")
    }
    for(i in 1 : nrow(c1_rep)){ c1_rep[i, ] <- rgamma(n1_c, shape = shape_c1[i, ], rate = rate_c1 [i, ]) }
    for(i in 1 : nrow(c2_rep)){ c2_rep[i, ] <- rgamma(n2_c, shape = shape_c2[i, ], rate = rate_c2 [i, ]) }
  } else if(x$model_output$dist_c == "lnorm") {
    mu_c1 <- x$model_output$`model summary`$BUGSoutput$sims.list$lmu_c1[, -index_mis_c1]
    mu_c2 <- x$model_output$`model summary`$BUGSoutput$sims.list$lmu_c2[, -index_mis_c2]
    if(is.null(dim(mu_c1)) == TRUE | is.null(dim(mu_c2)) == TRUE){
      stop("Not possible to plot posterior predictive checks with a single observed cost value in one arm.")
    }
    if(grepl("SELECTION", x$model_output$type) == TRUE){
      sigma_c1 <- 1 / sqrt(x$model_output$`model summary`$BUGSoutput$sims.list$ltau_c[, 1])
      sigma_c2 <- 1 / sqrt(x$model_output$`model summary`$BUGSoutput$sims.list$ltau_c[, 2])
      for(i in 1 : nrow(c1_rep)){ c1_rep[i, ] <- rlnorm(n1_c, meanlog = mu_c1[i, ], sdlog = sigma_c1 [i]) }
      for(i in 1 : nrow(c2_rep)){ c2_rep[i, ] <- rlnorm(n2_c, meanlog = mu_c2[i, ], sdlog = sigma_c2 [i]) }
    } else if(grepl("PATTERN", x$model_output$type) == TRUE) {
      sigma_c1 <- 1 / sqrt(x$model_output$`model summary`$BUGSoutput$sims.list$ltau_c_p1)
      sigma_c2 <- 1 / sqrt(x$model_output$`model summary`$BUGSoutput$sims.list$ltau_c_p2)
      for(i in 1 : nrow(c1_rep)){ c1_rep[i, ] <- rlnorm(n1_c, meanlog = mu_c1[i, ], sdlog = sigma_c1 [i]) }
      for(i in 1 : nrow(c2_rep)){ c2_rep[i, ] <- rlnorm(n2_c, meanlog = mu_c2[i, ], sdlog = sigma_c2 [i]) }
    } else if(grepl("HURDLE", x$model_output$type) == TRUE) {
      sigma_c1 <- 1 / sqrt(x$model_output$`model summary`$BUGSoutput$sims.list$ltau_c1)
      sigma_c2 <- 1 / sqrt(x$model_output$`model summary`$BUGSoutput$sims.list$ltau_c2)
      for(i in 1 : nrow(c1_rep)){ c1_rep[i, ] <- rlnorm(n1_c, meanlog = mu_c1[i, ], sdlog = sigma_c1 [i]) }
      for(i in 1 : nrow(c2_rep)){ c2_rep[i, ] <- rlnorm(n2_c, meanlog = mu_c2[i, ], sdlog = sigma_c2 [i]) }
    }
  }
  bayesplot::color_scheme_set(scheme = scheme_set) 
  if(type == "histogram") {
    if(exists("binwidth_e", where = exArgs)) {binwidth_e = exArgs$binwidth_e} else {binwidth_e = 1 / 30 }
    if(exists("binwidth_c", where = exArgs)) {binwidth_c = exArgs$binwidth_c} else {binwidth_c = 30 }
    if(exists("breaks", where = exArgs)) {breaks = exArgs$breaks} else {breaks = NULL }
    if(exists("freq", where = exArgs)) {freq = exArgs$freq} else {freq = TRUE }
    ppc_plot_e1 <- bayesplot::ppc_hist(y = e1_obs, yrep = e1_rep[1 : ndisplay, ], binwidth = binwidth_e, breaks = breaks, freq = freq)
    ppc_plot_e2 <- bayesplot::ppc_hist(y = e2_obs, yrep = e2_rep[1 : ndisplay, ], binwidth = binwidth_e, breaks = breaks, freq = freq)
    ppc_plot_c1 <- bayesplot::ppc_hist(y = c1_obs, yrep = c1_rep[1 : ndisplay, ], binwidth = binwidth_c, breaks = breaks, freq = freq)
    ppc_plot_c2 <- bayesplot::ppc_hist(y = c2_obs, yrep = c2_rep[1 : ndisplay, ], binwidth = binwidth_c, breaks = breaks, freq = freq)
  } else if(type == "boxplot") {
    if(exists("notch", where = exArgs)) {notch = exArgs$notch} else {notch = FALSE }
    if(exists("size", where = exArgs)) {size = exArgs$size} else {size = 0.5 }
    if(exists("alpha", where = exArgs)) {alpha = exArgs$alpha} else {alpha = 1 }
    ppc_plot_e1 <- bayesplot::ppc_boxplot(y = e1_obs, yrep = e1_rep[1 : ndisplay, ], notch = notch, size = size, alpha = alpha)
    ppc_plot_e2 <- bayesplot::ppc_boxplot(y = e2_obs, yrep = e2_rep[1 : ndisplay, ], notch = notch, size = size, alpha = alpha)
    ppc_plot_c1 <- bayesplot::ppc_boxplot(y = c1_obs, yrep = c1_rep[1 : ndisplay, ], notch = notch, size = size, alpha = alpha)
    ppc_plot_c2 <- bayesplot::ppc_boxplot(y = c2_obs, yrep = c2_rep[1 : ndisplay, ], notch = notch, size = size, alpha = alpha)
  } else if(type == "freqpoly") {
    if(exists("binwidth_e", where = exArgs)) {binwidth_e = exArgs$binwidth_e} else {binwidth_e = 1 / 30 }
    if(exists("binwidth_c", where = exArgs)) {binwidth_c = exArgs$binwidth_c} else {binwidth_c = 30 }
    if(exists("freq", where = exArgs)) {freq = exArgs$freq} else {freq = TRUE }
    if(exists("size", where = exArgs)) {size = exArgs$size} else {size = 0.25 }
    if(exists("alpha", where = exArgs)) {alpha = exArgs$alpha} else {alpha = 1 }
    ppc_plot_e1 <- bayesplot::ppc_freqpoly(y = e1_obs, yrep = e1_rep[1 : ndisplay, ], binwidth = binwidth_e, freq = freq, size = size, alpha = alpha)
    ppc_plot_e2 <- bayesplot::ppc_freqpoly(y = e2_obs, yrep = e2_rep[1 : ndisplay, ], binwidth = binwidth_e, freq = freq, size = size, alpha = alpha)
    ppc_plot_c1 <- bayesplot::ppc_freqpoly(y = c1_obs, yrep = c1_rep[1 : ndisplay, ], binwidth = binwidth_c, freq = freq, size = size, alpha = alpha)
    ppc_plot_c2 <- bayesplot::ppc_freqpoly(y = c2_obs, yrep = c2_rep[1 : ndisplay, ], binwidth = binwidth_c, freq = freq, size = size, alpha = alpha)
  } else if(type == "dens") {
    if(exists("trim", where = exArgs)) {trim = exArgs$trim} else {trim = FALSE }
    if(exists("size", where = exArgs)) {size = exArgs$size} else {size = 0.5 }
    if(exists("alpha", where = exArgs)) {alpha = exArgs$alpha} else {alpha = 1 }
    ppc_plot_e1 <- bayesplot::ppc_dens(y = e1_obs, yrep = e1_rep[1 : ndisplay, ], size = size, alpha = alpha, trim = FALSE)
    ppc_plot_e2 <- bayesplot::ppc_dens(y = e2_obs, yrep = e2_rep[1 : ndisplay, ], size = size, alpha = alpha, trim = FALSE)
    ppc_plot_c1 <- bayesplot::ppc_dens(y = c1_obs, yrep = c1_rep[1 : ndisplay, ], size = size, alpha = alpha, trim = FALSE)
    ppc_plot_c2 <- bayesplot::ppc_dens(y = c2_obs, yrep = c2_rep[1 : ndisplay, ], size = size, alpha = alpha, trim = FALSE)
  } else if(type == "dens_overlay") {
    if(exists("adjust", where = exArgs)) {adjust = exArgs$adjust} else {adjust = 1 }
    if(exists("bw", where = exArgs)) {bw = exArgs$bw} else {bw = "nrd0" }
    if(exists("trim", where = exArgs)) {trim = exArgs$trim} else {trim = FALSE }
    if(exists("size", where = exArgs)) {size = exArgs$size} else {size = 0.25 }
    if(exists("alpha", where = exArgs)) {alpha = exArgs$alpha} else {alpha = 0.7 }
    if(exists("kernel", where = exArgs)) {kernel = exArgs$kernel} else {kernel = "gaussian" }
    if(exists("n_dens", where = exArgs)) {n_dens = exArgs$n_dens} else {n_dens = 1024 }
    ppc_plot_e1 <- bayesplot::ppc_dens_overlay(y = e1_obs, yrep = e1_rep[1 : ndisplay, ], size = size, alpha = alpha, trim = FALSE, bw = bw, adjust = adjust, kernel = kernel, n_dens = n_dens)
    ppc_plot_e2 <- bayesplot::ppc_dens_overlay(y = e2_obs, yrep = e2_rep[1 : ndisplay, ], size = size, alpha = alpha, trim = FALSE, bw = bw, adjust = adjust, kernel = kernel, n_dens = n_dens)
    ppc_plot_c1 <- bayesplot::ppc_dens_overlay(y = c1_obs, yrep = c1_rep[1 : ndisplay, ], size = size, alpha = alpha, trim = FALSE, bw = bw, adjust = adjust, kernel = kernel, n_dens = n_dens)
    ppc_plot_c2 <- bayesplot::ppc_dens_overlay(y = c2_obs, yrep = c2_rep[1 : ndisplay, ], size = size, alpha = alpha, trim = FALSE, bw = bw, adjust = adjust, kernel = kernel, n_dens = n_dens)
  } else if(type == "ecdf_overlay") {
    if(exists("size", where = exArgs)) {size = exArgs$size} else {size = 0.25 }
    if(exists("alpha", where = exArgs)) {alpha = exArgs$alpha} else {alpha = 0.7 }
    if(exists("pad", where = exArgs)) {pad = exArgs$pad} else {pad = TRUE }
    if(exists("discrete", where = exArgs)) {discrete = exArgs$discrete} else {discrete = FALSE }
    ppc_plot_e1 <- bayesplot::ppc_ecdf_overlay(y = e1_obs, yrep = e1_rep[1 : ndisplay, ], size = size, alpha = alpha, discrete = discrete, pad = pad)
    ppc_plot_e2 <- bayesplot::ppc_ecdf_overlay(y = e2_obs, yrep = e2_rep[1 : ndisplay, ], size = size, alpha = alpha, discrete = discrete, pad = pad)
    ppc_plot_c1 <- bayesplot::ppc_ecdf_overlay(y = c1_obs, yrep = c1_rep[1 : ndisplay, ], size = size, alpha = alpha, discrete = discrete, pad = pad)
    ppc_plot_c2 <- bayesplot::ppc_ecdf_overlay(y = c2_obs, yrep = c2_rep[1 : ndisplay, ], size = size, alpha = alpha, discrete = discrete, pad = pad)
  } else if(type == "stat" & is.null(corr) == TRUE) {
    if(bpp == FALSE){
    if(exists("binwidth_e", where = exArgs)) {binwidth_e = exArgs$binwidth_e} else {binwidth_e = 1 / 30 }
    if(exists("binwidth_c", where = exArgs)) {binwidth_c = exArgs$binwidth_c} else {binwidth_c = 30 }
    if(exists("freq", where = exArgs)) {freq = exArgs$freq} else {freq = TRUE }
    if(exists("breaks", where = exArgs)) {breaks = exArgs$breaks} else {breaks = NULL }
    ppc_plot_e1 <- bayesplot::ppc_stat(y = e1_obs, yrep = e1_rep[1 : ndisplay, ], stat = stat, binwidth = binwidth_e, breaks = breaks, freq = freq)
    ppc_plot_e2 <- bayesplot::ppc_stat(y = e2_obs, yrep = e2_rep[1 : ndisplay, ], stat = stat, binwidth = binwidth_e, breaks = breaks, freq = freq)
    ppc_plot_c1 <- bayesplot::ppc_stat(y = c1_obs, yrep = c1_rep[1 : ndisplay, ], stat = stat, binwidth = binwidth_c, breaks = breaks, freq = freq)
    ppc_plot_c2 <- bayesplot::ppc_stat(y = c2_obs, yrep = c2_rep[1 : ndisplay, ], stat = stat, binwidth = binwidth_c, breaks = breaks, freq = freq)
    }
    if(bpp == TRUE){
    if(exists("trim", where = exArgs)) {trim = exArgs$trim} else {trim = FALSE }
    if(exists("size", where = exArgs)) {size = exArgs$size} else {size = 0.5 }
    if(exists("alpha", where = exArgs)) {alpha = exArgs$alpha} else {alpha = 1 }    
      stat_e1 <- stat(e1_obs)
      stat_c1 <- stat(c1_obs)
      values_e1 <- apply(e1_rep[1 : ndisplay, ], 1, stat)
      values_c1 <- apply(c1_rep[1 : ndisplay, ], 1, stat)
      df1<-data.frame(values_e1, values_c1)
      stat_e2 <- stat(e2_obs)
      stat_c2 <- stat(c2_obs)
      values_e2 <- apply(e2_rep[1 : ndisplay, ], 1, stat)
      values_c2 <- apply(c2_rep[1 : ndisplay, ], 1, stat)
      df2 <- data.frame(values_e2, values_c2)
      pval_e1 <- sum(values_e1 > stat_e1) / length(values_e1)
      pval_c1 <- sum(values_c1 > stat_c1) / length(values_c1)
      pval_e2 <- sum(values_e2 > stat_e2) / length(values_e2)
      pval_c2 <- sum(values_c2 > stat_c2) / length(values_c2)
      color_sheme_rep <- bayesplot::color_scheme_get()[1]
      color_sheme_rep_border <- bayesplot::color_scheme_get()[2]
      color_sheme_obs <- bayesplot::color_scheme_get()[5]
      color_sheme_obs_border <- bayesplot::color_scheme_get()[6]
      set_theme <- bayesplot::bayesplot_theme_get()
      ggplot2::theme_set(set_theme)
      paste_p_e1 <- sprintf("Prob (T(y[rep]) > T(y)) == '%0.3f'", pval_e1)
      paste_p_c1 <- sprintf("Prob (T(y[rep]) > T(y)) == '%0.3f'", pval_c1)
      paste_p_e2 <- sprintf("Prob (T(y[rep]) > T(y)) == '%0.3f'", pval_e2)
      paste_p_c2 <- sprintf("Prob (T(y[rep]) > T(y)) == '%0.3f'", pval_c2)
      ppc_plot_e1 <- ggplot2::ggplot(df1, ggplot2::aes(x = values_e1)) + 
        ggplot2::geom_density(fill = color_sheme_rep$light, color = color_sheme_rep_border$light_highlight, 
                              trim = trim, size = size, alpha = alpha) +
        ggplot2::geom_vline(xintercept = stat_e1, color = color_sheme_obs$dark, size = 1, linetype = "dashed") +
        ggplot2::scale_y_continuous(expand = c(0, 0)) +
        ggplot2::theme(axis.text.y = ggplot2::element_blank(), axis.ticks.y = ggplot2::element_blank(), 
                       axis.title.y = ggplot2::element_blank(), axis.title.x = ggplot2::element_blank()) +
        ggplot2::annotate("text", x = -Inf, y = Inf, label = paste_p_e1, parse = TRUE, hjust = -0.1, vjust = 2, size = 4)
      ppc_plot_c1 <- ggplot2::ggplot(df1, ggplot2::aes(x = values_c1)) + 
        ggplot2::geom_density(fill = color_sheme_rep$light, color = color_sheme_rep_border$light_highlight, 
                              trim = trim, size = size, alpha = alpha) +
        ggplot2::geom_vline(xintercept = stat_c1, color = color_sheme_obs$dark, size = 1, linetype = "dashed") +
        ggplot2::scale_y_continuous(expand = c(0, 0)) +
        ggplot2::theme(axis.text.y = ggplot2::element_blank(), axis.ticks.y = ggplot2::element_blank(), 
                       axis.title.y = ggplot2::element_blank(), axis.title.x = ggplot2::element_blank()) +
        ggplot2::annotate("text", x = -Inf, y = Inf, label = paste_p_c1, parse = TRUE, hjust = -0.1, vjust = 2, size = 4)
      ppc_plot_e2 <- ggplot2::ggplot(df2, ggplot2::aes(x = values_e2)) + 
        ggplot2::geom_density(fill = color_sheme_rep$light, color = color_sheme_rep_border$light_highlight, 
                              trim = trim, size = size, alpha = alpha) +
        ggplot2::geom_vline(xintercept = stat_e2, color = color_sheme_obs$dark, size = 1, linetype = "dashed") +
        ggplot2::scale_y_continuous(expand = c(0, 0)) +
        ggplot2::theme(axis.text.y = ggplot2::element_blank(), axis.ticks.y = ggplot2::element_blank(), 
                       axis.title.y = ggplot2::element_blank(), axis.title.x = ggplot2::element_blank()) +
        ggplot2::annotate("text", x = -Inf, y = Inf, label = paste_p_e2, parse = TRUE, hjust = -0.1, vjust = 2, size = 4)
      ppc_plot_c2 <- ggplot2::ggplot(df2, ggplot2::aes(x = values_c2)) + 
        ggplot2::geom_density(fill = color_sheme_rep$light, color = color_sheme_rep_border$light_highlight, 
                              trim = trim, size = size, alpha = alpha) +
        ggplot2::geom_vline(xintercept = stat_c2, color = color_sheme_obs$dark, size = 1, linetype = "dashed") +
        ggplot2::scale_y_continuous(expand = c(0, 0)) +
        ggplot2::theme(axis.text.y = ggplot2::element_blank(), axis.ticks.y = ggplot2::element_blank(), 
                       axis.title.y = ggplot2::element_blank(), axis.title.x = ggplot2::element_blank()) +
        ggplot2::annotate("text", x = -Inf, y = Inf, label = paste_p_c2, parse = TRUE, hjust = -0.1, vjust = 2, size = 4)
     }
  } else if(type == "stat" & is.null(corr) == FALSE) {
      ec_arm1 <- cbind(x$data_set$effects$Control, x$data_set$costs$Control)[complete.cases(cbind(x$data_set$effects$Control, x$data_set$costs$Control)), ]
      stat_biv1 <- corr(ec_arm1[, 1], ec_arm1[, 2])
      array_1 <- array(data = NA, dim = c(nrow(e1_rep[1 : ndisplay, ]), nrow(ec_arm1), 2))
      array_1[, , 1] <- e1_rep[1 : ndisplay, 1:nrow(ec_arm1)]
      array_1[, , 2] <- c1_rep[1 : ndisplay, 1:nrow(ec_arm1)]
      values_biv1 <- apply(array_1, 1, corr)[2, ]
      df1 <- data.frame(values_biv1)
      ec_arm2 <- cbind(x$data_set$effects$Intervention, x$data_set$costs$Intervention)[complete.cases(cbind(x$data_set$effects$Intervention, x$data_set$costs$Intervention)), ]
      stat_biv2 <- corr(ec_arm2[, 1], ec_arm2[, 2])
      array_2 <- array(data = NA, dim = c(nrow(e2_rep[1 : ndisplay, ]), nrow(ec_arm2), 2))
      array_2[, , 1] <- e2_rep[1 : ndisplay, 1:nrow(ec_arm2)]
      array_2[, , 2] <- c2_rep[1 : ndisplay, 1:nrow(ec_arm2)]
      values_biv2 <- apply(array_2, 1, corr)[2, ]
      df2 <- data.frame(values_biv2)
      pval_1 <- sum(values_biv1 > stat_biv1) / length(values_biv1)
      pval_2 <- sum(values_biv2 > stat_biv2) / length(values_biv2)
      paste_p_biv1 <- sprintf("Prob (T(y[rep]) > T(y)) == '%0.3f'", pval_1)
      paste_p_biv2 <- sprintf("Prob (T(y[rep]) > T(y)) == '%0.3f'", pval_2)
      color_sheme_rep <- bayesplot::color_scheme_get()[1]
      color_sheme_rep_border <- bayesplot::color_scheme_get()[2]
      color_sheme_obs <- bayesplot::color_scheme_get()[5]
      color_sheme_obs_border <- bayesplot::color_scheme_get()[6]
      set_theme <- bayesplot::bayesplot_theme_get()
      if(bpp == FALSE) {
      if(exists("binwidth", where = exArgs)) {binwidth = exArgs$binwidth} else {binwidth = 1 / 30 }
      ppc_plot_e1 <- ggplot2::ggplot(df1, ggplot2::aes(x = values_biv1)) + 
        ggplot2::geom_histogram(fill = color_sheme_rep$light, 
                                color = color_sheme_rep_border$light_highlight, 
                                binwidth = binwidth) +
        ggplot2::geom_vline(xintercept = stat_biv1, color = color_sheme_obs$dark, size = 2) +
        ggplot2::scale_y_continuous(expand = c(0, 0)) +
        ggplot2::theme(axis.text.y = ggplot2::element_blank(), axis.ticks.y = ggplot2::element_blank(), 
                       axis.title.y = ggplot2::element_blank(), axis.title.x = ggplot2::element_blank()) 
      ppc_plot_c2 <- ggplot2::ggplot(df2, ggplot2::aes(x = values_biv2)) + 
        ggplot2::geom_histogram(fill = color_sheme_rep$light, 
                                color = color_sheme_rep_border$light_highlight, 
                                binwidth = binwidth) +
        ggplot2::geom_vline(xintercept = stat_biv2, color = color_sheme_obs$dark, size = 2) +
        ggplot2::scale_y_continuous(expand = c(0, 0)) +
        ggplot2::theme(axis.text.y = ggplot2::element_blank(), axis.ticks.y = ggplot2::element_blank(), 
                       axis.title.y = ggplot2::element_blank(), axis.title.x = ggplot2::element_blank())
      }
      if(bpp == TRUE) {
        if(exists("trim", where = exArgs)) {trim = exArgs$trim} else {trim = FALSE }
        if(exists("size", where = exArgs)) {size = exArgs$size} else {size = 0.5 }
        if(exists("alpha", where = exArgs)) {alpha = exArgs$alpha} else {alpha = 1 } 
        ppc_plot_e1 <- ggplot2::ggplot(df1, ggplot2::aes(x = values_biv1)) +
          ggplot2::geom_density(fill = color_sheme_rep$light, 
                                color = color_sheme_rep_border$light_highlight, trim = trim, size = size, alpha = alpha) +
          ggplot2::geom_vline(xintercept = stat_biv1, color = color_sheme_obs$dark, size = 1, linetype = "dashed") +
          ggplot2::scale_y_continuous(expand = c(0, 0)) +
          ggplot2::theme(axis.text.y = ggplot2::element_blank(), axis.ticks.y = ggplot2::element_blank(), 
                         axis.title.y = ggplot2::element_blank(), axis.title.x = ggplot2::element_blank()) +
          ggplot2::annotate("text", x = -Inf, y = Inf, 
                            label = paste_p_biv1, parse = TRUE, hjust = -0.1, vjust = 2, size=4)
        ppc_plot_c2 <- ggplot2::ggplot(df2, ggplot2::aes(x = values_biv2)) +
          ggplot2::geom_density(fill = color_sheme_rep$light, 
                                color = color_sheme_rep_border$light_highlight, trim = trim, size = size, alpha = alpha) +
          ggplot2::geom_vline(xintercept = stat_biv2, color = color_sheme_obs$dark, size = 1, linetype = "dashed") +
          ggplot2::scale_y_continuous(expand = c(0, 0)) +
          ggplot2::theme(axis.text.y = ggplot2::element_blank(), axis.ticks.y = ggplot2::element_blank(), 
                         axis.title.y = ggplot2::element_blank(), axis.title.x = ggplot2::element_blank()) +
          ggplot2::annotate("text", x = -Inf, y = Inf, 
                            label = paste_p_biv2, parse = TRUE, hjust = -0.1, vjust = 2, size=4)
      }
  } else if(type == "stat_2d") {
    if(exists("size", where = exArgs)) {size = exArgs$size} else {size = 2.5 }
    if(exists("alpha", where = exArgs)) {alpha = exArgs$alpha} else {alpha = 0.7 }
    if(exists("stat", where = exArgs)) {stat = exArgs$stat} else {stat = c("mean", "sd") }
    ppc_plot_e1 <- bayesplot::ppc_stat_2d(y = e1_obs, yrep = e1_rep[1 : ndisplay, ], stat = stat, size = size, alpha = alpha)
    ppc_plot_e2 <- bayesplot::ppc_stat_2d(y = e2_obs, yrep = e2_rep[1 : ndisplay, ], stat = stat, size = size, alpha = alpha)
    ppc_plot_c1 <- bayesplot::ppc_stat_2d(y = c1_obs, yrep = c1_rep[1 : ndisplay, ], stat = stat, size = size, alpha = alpha)
    ppc_plot_c2 <- bayesplot::ppc_stat_2d(y = c2_obs, yrep = c2_rep[1 : ndisplay, ], stat = stat, size = size, alpha = alpha)
  } else if(type == "error_hist") {
    if(exists("binwidth_e", where = exArgs)) {binwidth_e = exArgs$binwidth_e} else {binwidth_e = 1 / 30 }
    if(exists("binwidth_c", where = exArgs)) {binwidth_c = exArgs$binwidth_c} else {binwidth_c = 30 }
    if(exists("freq", where = exArgs)) {freq = exArgs$freq} else {freq = TRUE }
    if(exists("breaks", where = exArgs)) {breaks = exArgs$breaks} else {breaks = NULL }
    ppc_plot_e1 <- bayesplot::ppc_error_hist(y = e1_obs, yrep = e1_rep[1 : ndisplay, ], binwidth = binwidth_e, breaks = breaks, freq = freq)
    ppc_plot_e2 <- bayesplot::ppc_error_hist(y = e2_obs, yrep = e2_rep[1 : ndisplay, ], binwidth = binwidth_e, breaks = breaks, freq = freq)
    ppc_plot_c1 <- bayesplot::ppc_error_hist(y = c1_obs, yrep = c1_rep[1 : ndisplay, ], binwidth = binwidth_c, breaks = breaks, freq = freq)
    ppc_plot_c2 <- bayesplot::ppc_error_hist(y = c2_obs, yrep = c2_rep[1 : ndisplay, ], binwidth = binwidth_c, breaks = breaks, freq = freq)
  } else if(type == "error_scatter") {
    if(exists("size", where = exArgs)) {size = exArgs$size} else {size = 2.5 }
    if(exists("alpha", where = exArgs)) {alpha = exArgs$alpha} else {alpha = 0.8 }
    ppc_plot_e1<-bayesplot::ppc_error_scatter(y = e1_obs, yrep = e1_rep[1 : ndisplay, ], size = size, alpha = alpha)
    ppc_plot_e2<-bayesplot::ppc_error_scatter(y = e2_obs, yrep = e2_rep[1 : ndisplay, ], size = size, alpha = alpha)
    ppc_plot_c1<-bayesplot::ppc_error_scatter(y = c1_obs, yrep = c1_rep[1 : ndisplay, ], size = size, alpha = alpha)
    ppc_plot_c2<-bayesplot::ppc_error_scatter(y = c2_obs, yrep = c2_rep[1 : ndisplay, ], size = size, alpha = alpha)
  } else if(type == "error_scatter_avg") {
    if(exists("size", where = exArgs)) {size = exArgs$size} else {size = 2.5 }
    if(exists("alpha", where = exArgs)) {alpha = exArgs$alpha} else {alpha = 0.8 }
    ppc_plot_e1 <- bayesplot::ppc_error_scatter_avg(y = e1_obs, yrep = e1_rep[1 : ndisplay, ], size = size, alpha = alpha)
    ppc_plot_e2 <- bayesplot::ppc_error_scatter_avg(y = e2_obs, yrep = e2_rep[1 : ndisplay, ], size = size, alpha = alpha)
    ppc_plot_c1 <- bayesplot::ppc_error_scatter_avg(y = c1_obs, yrep = c1_rep[1 : ndisplay, ], size = size, alpha = alpha)
    ppc_plot_c2 <- bayesplot::ppc_error_scatter_avg(y = c2_obs, yrep = c2_rep[1 : ndisplay, ], size = size, alpha = alpha)
  } else if(type == "error_binned") {
    if(exists("bins", where = exArgs)) {bins = exArgs$bins} else {bins = NULL }
    if(exists("size", where = exArgs)) {size = exArgs$size} else {size = 1 }
    if(exists("alpha", where = exArgs)) {alpha = exArgs$alpha} else {alpha = 0.25 }
    ppc_plot_e1 <- bayesplot::ppc_error_binned(y = e1_obs, yrep = e1_rep[1 : ndisplay, ], bins = bins, size = size, alpha = alpha)
    ppc_plot_e2 <- bayesplot::ppc_error_binned(y = e2_obs, yrep = e2_rep[1 : ndisplay, ], bins = bins, size = size, alpha = alpha)
    ppc_plot_c1 <- bayesplot::ppc_error_binned(y = c1_obs, yrep = c1_rep[1 : ndisplay, ], bins = bins, size = size, alpha = alpha)
    ppc_plot_c2 <- bayesplot::ppc_error_binned(y = c2_obs, yrep = c2_rep[1 : ndisplay, ], bins = bins, size = size, alpha = alpha)
  } else if(type == "intervals") {
    if(exists("prob", where = exArgs)) {prob = exArgs$prob} else {prob = 0.5 }
    if(exists("prob_outer", where = exArgs)) {prob_outer = exArgs$prob_outer} else {prob_outer = 0.9 }
    if(exists("size", where = exArgs)) {size = exArgs$size} else {size = 1 }
    if(exists("fatten", where = exArgs)) {fatten = exArgs$fatten} else {fatten = 3 }
    ppc_plot_e1 <- bayesplot::ppc_intervals(y = e1_obs, yrep = e1_rep[1 : ndisplay, ], prob = prob, prob_outer = prob_outer, size = size, fatten = fatten)
    ppc_plot_e2 <- bayesplot::ppc_intervals(y = e2_obs, yrep = e2_rep[1 : ndisplay, ], prob = prob, prob_outer = prob_outer, size = size, fatten = fatten)
    ppc_plot_c1 <- bayesplot::ppc_intervals(y = c1_obs, yrep = c1_rep[1 : ndisplay, ], prob = prob, prob_outer = prob_outer, size = size, fatten = fatten)
    ppc_plot_c2 <- bayesplot::ppc_intervals(y = c2_obs, yrep = c2_rep[1 : ndisplay, ], prob = prob, prob_outer = prob_outer, size = size, fatten = fatten)
  } else if(type == "ribbon") {
    if(exists("prob", where = exArgs)) {prob = exArgs$prob} else {prob = 0.5 }
    if(exists("prob_outer", where = exArgs)) {prob_outer = exArgs$prob_outer} else {prob_outer = 0.9 }
    if(exists("size", where = exArgs)) {size = exArgs$size} else {size = 0.25 }
    if(exists("alpha", where = exArgs)) {alpha = exArgs$alpha} else {alpha = 0.33 }
    ppc_plot_e1 <- bayesplot::ppc_ribbon(y = e1_obs, yrep = e1_rep[1 : ndisplay, ], prob = prob, prob_outer = prob_outer, alpha = alpha, size = size)
    ppc_plot_e2 <- bayesplot::ppc_ribbon(y = e2_obs, yrep = e2_rep[1 : ndisplay, ], prob = prob, prob_outer = prob_outer, alpha = alpha, size = size)
    ppc_plot_c1 <- bayesplot::ppc_ribbon(y = c1_obs, yrep = c1_rep[1 : ndisplay, ], prob = prob, prob_outer = prob_outer, alpha = alpha, size = size)
    ppc_plot_c2 <- bayesplot::ppc_ribbon(y = c2_obs, yrep = c2_rep[1 : ndisplay, ], prob = prob, prob_outer = prob_outer, alpha = alpha, size = size)
  } else if(type == "scatter") {
    if(exists("size", where = exArgs)) {size = exArgs$size} else {size = 2.5 }
    if(exists("alpha", where = exArgs)) {alpha = exArgs$alpha} else {alpha = 0.8 }
    ppc_plot_e1 <- bayesplot::ppc_scatter(y = e1_obs, yrep = e1_rep[1 : ndisplay, ], alpha = alpha, size = size)
    ppc_plot_e2 <- bayesplot::ppc_scatter(y = e2_obs, yrep = e2_rep[1 : ndisplay, ], alpha = alpha, size = size)
    ppc_plot_c1 <- bayesplot::ppc_scatter(y = c1_obs, yrep = c1_rep[1 : ndisplay, ], alpha = alpha, size = size)
    ppc_plot_c2 <- bayesplot::ppc_scatter(y = c2_obs, yrep = c2_rep[1 : ndisplay, ], alpha = alpha, size = size)
  } else if(type == "scatter_avg") {
    if(exists("size", where = exArgs)) {size = exArgs$size} else {size = 2.5 }
    if(exists("alpha", where = exArgs)) {alpha = exArgs$alpha} else {alpha = 0.8 }
    ppc_plot_e1 <- bayesplot::ppc_scatter_avg(y = e1_obs, yrep = e1_rep[1 : ndisplay, ], alpha = alpha, size = size)
    ppc_plot_e2 <- bayesplot::ppc_scatter_avg(y = e2_obs, yrep = e2_rep[1 : ndisplay, ], alpha = alpha, size = size)
    ppc_plot_c1 <- bayesplot::ppc_scatter_avg(y = c1_obs, yrep = c1_rep[1 : ndisplay, ], alpha = alpha, size = size)
    ppc_plot_c2 <- bayesplot::ppc_scatter_avg(y = c2_obs, yrep = c2_rep[1 : ndisplay, ], alpha = alpha, size = size)
  }
  if(outcome == "effects_arm1") {
    ppc_plot_output <- ppc_plot_e1
  } else if(outcome == "effects_arm2"){
    ppc_plot_output <- ppc_plot_e2
  } else if(outcome == "costs_arm1"){
    ppc_plot_output <- ppc_plot_c1
  } else if(outcome == "costs_arm2"){
    ppc_plot_output <- ppc_plot_c2
  } else if(outcome == "effects") {
    ppc_plot_output <- ggpubr::ggarrange(ppc_plot_e1, ppc_plot_e2, labels = c("effects (t=1)", "effects (t=2)"), ncol = 2, nrow = 1, common.legend = TRUE, legend = legend)
  } else if(outcome == "costs") {
    ppc_plot_output <- ggpubr::ggarrange(ppc_plot_c1, ppc_plot_c2, labels = c("costs (t=1)", "costs (t=2)"), ncol = 2, nrow = 1, common.legend = TRUE, legend = legend)
  } else if(outcome == "all") {
    if(is.null(corr) == TRUE) {
    ppc_plot_output <- ggpubr::ggarrange(ppc_plot_e1, ppc_plot_e2, ppc_plot_c1, ppc_plot_c2, 
                                         labels = c("effects (t=1)", "effects (t=2)", "costs (t=1)", "costs (t=2)"), ncol = 2, nrow = 2, common.legend = TRUE, legend = legend)
    }
    if(is.null(corr) == FALSE) {
      ppc_plot_output <- ggpubr::ggarrange(ppc_plot_e1, ppc_plot_c2, labels = c("control", "intervention"), ncol = 2, nrow = 1, common.legend = TRUE, legend = legend)
    }
  }
  ppc_plot_output <- ppc_plot_output + ggplot2::theme(legend.position = legend)
  if(length(theme) != 0) {
    if(theme == "base") ppc_plot_output <- ppc_plot_output + ggthemes::theme_base()
    if(theme == "calc") ppc_plot_output <- ppc_plot_output + ggthemes::theme_calc()
    if(theme == "economist") ppc_plot_output <- ppc_plot_output + ggthemes::theme_economist()
    if(theme == "excel") ppc_plot_output <- ppc_plot_output + ggthemes::theme_excel()
    if(theme == "few") ppc_plot_output <- ppc_plot_output + ggthemes::theme_few()
    if(theme == "538") ppc_plot_output <- ppc_plot_output + ggthemes::theme_fivethirtyeight()
    if(theme == "gdocs") ppc_plot_output <- ppc_plot_output + ggthemes::theme_gdocs()
    if(theme == "hc") ppc_plot_output <- ppc_plot_output + ggthemes::theme_hc()
    if(theme == "par") ppc_plot_output <- ppc_plot_output + ggthemes::theme_par()
    if(theme == "solarized") ppc_plot_output <- ppc_plot_output + ggthemes::theme_solarized()
    if(theme == "pander") ppc_plot_output <- ppc_plot_output + ggthemes::theme_pander()
    if(theme == "stata") ppc_plot_output <- ppc_plot_output + ggthemes::theme_stata()
    if(theme == "tufte") ppc_plot_output <- ppc_plot_output + ggthemes::theme_tufte()
    if(theme == "wsj") ppc_plot_output <- ppc_plot_output + ggthemes::theme_wsj()
  }
  if(legend == FALSE) { 
    ppc_plot_output <- ppc_plot_output + legend_none() 
  }
  return(ppc_plot_output)
}

Try the missingHE package in your browser

Any scripts or data that you put into this service are public.

missingHE documentation built on July 1, 2020, 5:50 p.m.