Nothing
#' 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}}, \code{\link{hurdle}} or \code{\link{lmdm}}
#'
#' The focus is restricted to full Bayesian models in cost-effectiveness analyses based on the function \code{\link{selection}}, \code{\link{pattern}},
#' \code{\link{hurdle}} or \code{\link{lmdm}} 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}}, \code{\link{hurdle}} or \code{\link{lmdm}}.
#' @param type Type of posterior predictive check among some pre-defined types, mostly taken from the package \strong{bayesplot}. For a full list of available options see details.
#' @param outcome The outcome variables that should be displayed. Options are: 'both' (default) for both effects and costs; 'effects' or 'costs' for the effects
#' or costs separately.
#' @param ndisplay Number of posterior replications to be used to generate the plots.
#' @param trt treatment group for which plots should be displayed. Choices include: 'all' (default) for all groups; 'none' for results across all groups; any
#' character or numeric value denoting the treatment group name or index associated with the treatment variable in the original data set.
#' @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 ... 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{lmdm}} \code{\link{diagnostic}}
#' @details The function 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}}, \code{\link{hurdle}} or \code{\link{lmdm}}. 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 ggpubr
#' @importFrom stats rnorm rbeta rgamma rlnorm rweibull rnbinom rbinom rpois rlogis rexp complete.cases
#' @importFrom bayesplot color_scheme_set ppc_hist ppc_boxplot ppc_freqpoly ppc_dens ppc_dens_overlay ppc_ecdf_overlay ppc_pit_ecdf ppc_stat color_scheme_get bayesplot_theme_get ppc_stat_2d ppc_error_hist ppc_error_scatter ppc_error_scatter_avg ppc_error_binned ppc_intervals ppc_ribbon ppc_scatter ppc_scatter_avg legend_none ppc_loo_pit_ecdf ppc_loo_pit_overlay
#' @importFrom graphics text
#' @export
#' @examples
#' # For examples see the function \code{\link{selection}}, \code{\link{pattern}},
#' # \code{\link{hurdle}} or \code{\link{lmdm}}
#' #
ppc <- function(x, type = "histogram", outcome = "both", ndisplay = 15,
trt = "all", theme = NULL, scheme_set = NULL, ...) {
if(!inherits(x, "missingHE")) { stop("Only objects of class 'missingHE' can be used")}
if(!isTRUE(requireNamespace("ggplot2")) | !isTRUE(requireNamespace("bayesplot")) | !isTRUE(requireNamespace("ggpubr"))) {
stop("You need to install the R packages 'ggplot2', 'bayesplot' and 'ggpubr'.
Please run in your R terminal:\n install.packages('ggplot2', 'bayesplot', 'ggpubr')")
}
if(is.null(theme)) { theme = "classic"}
if(!theme %in% c("classic", "base", "calc", "economist", "excel", "few", "538", "gdocs",
"hc", "par", "pander", "solarized", "stata", "tufte", "wsj")) {
stop("Please provide a valid theme type among those available. Type ''help(ggthemes)'' for more details.")}
if(is.null(scheme_set)) { 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("Please provide a valid scheme set 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",
"pit_ecdf", "loo_pit_ecdf", "loo_pit_overlay")) {
stop("Please provide a plot type among those available.")}
if(!outcome %in% c("both", "effects", "costs")) {
stop("Please provide an outcome name among those available.")}
exArgs <- list(...)
if(exists("legend.position", where = exArgs)) { legend.pos = exArgs$legend.position } else { legend.pos = "top"}
if(exists("hjust", where = exArgs)) { hjust = exArgs$hjust } else { hjust = -0.5}
if(exists("vjust", where = exArgs)) { vjust = exArgs$vjust } else { vjust = 0.3}
if(exists("digits", where = exArgs)) { digits = exArgs$digits } else { digits = 7}
if(exists("corr", where = exArgs)) { corr = exArgs$corr} else { corr = NULL}
values <- NULL
if(!is.null(corr) & outcome != "both" | !is.null(corr) & !is.function(corr)) {
stop("The argument 'corr' must be a function and is available only when outcome is set to 'both'.")}
if(exists("stat", where = exArgs)) { stat = exArgs$stat} else { stat = NULL}
if(exists("bpp", where = exArgs)) { bpp = exArgs$bpp} else { bpp = FALSE}
if(bpp & type == "stat_2d") { stop("Bayesian p-values are not available for 'stat_2d' type.")}
if(type == "stat" & is.null(corr)) {
if(is.null(stat) | length(stat) !=1 | !is.function(stat)) {
stop("Please use the argument 'stat' to provide a single function for computing statistics.")}
}
if(type == "stat_2d" & is.null(stat) | type == "stat_2d" & length(stat) !=2 | type == "stat_2d" & !any(sapply(stat, is.function))) {
stop("Please use the argument 'stat' to provide a vector of two functions for computing statistics.")}
n_sim <- x$model_output$model$n.iter
n_trt <- length(x$data_set$data_raw$trt_index)
trt_lev <- names(x$data_set$data_raw$trt_index)
if(is.numeric(trt)) {
if(length(trt) == 1 & !any(trt %in% c(1:n_trt))) {
stop("Please provide a trt value among those available.")
} else if(length(trt) > 1 & !any(trt %in% c(2:n_trt))){
stop("Please provide a trt value among those available.")}
} else if(is.character(trt)) {
if(length(trt) == 1 & !any(trt %in% c("all", trt_lev, "none"))) {
stop("Please provide a trt name among those available.")
} else if(length(trt) > 1 & !length(trt) %in% c(2:n_trt)){
stop("Please provide a trt name among those available.")}
} else { stop("Please provide a trt name or value among those available.")}
if(!is.numeric(ndisplay) | ndisplay < 1 | ndisplay > n_sim) {
stop("Number of replications to display is not valid.")}
if(exists("time.plot", where = exArgs)) { time.plot = exArgs$time.plot } else { time.plot = 1}
if(x$data_format == "long") {
n_time <- max(x$data_set$data_raw$time_long, na.rm = TRUE)
} else { n_time <- 1}
if(is.numeric(time.plot)) {
if(length(time.plot) == 1 & !any(time.plot %in% c(1:n_time))) {
stop("Please provide a time.plot value among those available.")
} else if(length(time.plot) > 1 & !any(time.plot %in% c(2:n_time))){
stop("Please provide a time.plot value among those available.")}
} else if(is.character(time.plot)) {
if(length(time.plot) == 1 & !any(time.plot %in% c("all", 1:n_time))) {
stop("Please provide a time.plot name among those available.")
} else if(length(time.plot) > 1 & !length(time.plot) %in% c(2:n_time)){
stop("Please provide a time.plot name among those available.")}
} else { stop("Please provide a time.plot name or value among those available.")}
if(x$data_format == "wide") {
me_type <- ifelse(x$data_set$data_raw$me == 1, "missing", "observed")
mc_type <- ifelse(x$data_set$data_raw$mc == 1, "missing", "observed")
e_obs <- x$data_set$data_raw$e[me_type == "observed"]
c_obs <- x$data_set$data_raw$c[mc_type == "observed"]
trt_obs_e <- x$data_set$data_raw$data$trt[me_type == "observed"]
trt_obs_c <- x$data_set$data_raw$data$trt[mc_type == "observed"]
n_all <- sum(x$data_set$data_raw$n)
n_trt_obs_e <- x$data_set$data_raw$n_obs_e
n_trt_obs_c <- x$data_set$data_raw$n_obs_c
n_obs_e <- sum(n_trt_obs_e)
n_obs_c <- sum(n_trt_obs_c)
e_rep <- matrix(NA, nrow = n_sim, n_obs_e)
c_rep <- matrix(NA, nrow = n_sim, n_obs_c)
e_obs_trt <- lapply(x$data_set$data_raw$efft, function(x) x[!is.na(x)])
c_obs_trt <- lapply(x$data_set$data_raw$costt, function(x) x[!is.na(x)])
e_rep_trt <- lapply(x$data_set$data_raw$m_efft, function(x) matrix(NA, nrow = n_sim, ncol = length(x[x == 0])))
c_rep_trt <- lapply(x$data_set$data_raw$m_costt, function(x) matrix(NA, nrow = n_sim, ncol = length(x[x == 0])))
if(x$model_output$method == "SELECTION") {
cmu_e <- x$model_output$model$BUGSoutput$sims.list$cmu_e[, me_type == "observed"]
s_e <- x$model_output$model$BUGSoutput$sims.list$s_e
if(x$model_output$dist_c == "lnorm") {
clmu_c <- x$model_output$model$BUGSoutput$sims.list$clmu_c[, mc_type == "observed"]
ltau_c <- x$model_output$model$BUGSoutput$sims.list$ltau_c
} else {
cmu_c <- x$model_output$model$BUGSoutput$sims.list$cmu_c[, mc_type == "observed"]
s_c <- x$model_output$model$BUGSoutput$sims.list$s_c
}
if(x$model_output$dist_e == "norm") {
for(i in 1:n_sim){
e_rep[i, ] <- rnorm(n_obs_e, mean = cmu_e[i, ], sd = s_e[i])
for(trt_i in trt_lev){ e_rep_trt[[trt_i]][i, ] <- rnorm(n_trt_obs_e[trt_i], mean = cmu_e[i, trt_obs_e == trt_i], sd = s_e[i])}
}
}
if(x$model_output$dist_e == "logis") {
for(i in 1:n_sim){ e_rep[i, ] <- rlogis(n_obs_e, location = cmu_e[i, ], scale = s_e[i])
for(trt_i in trt_lev){ e_rep_trt[[trt_i]][i, ] <- rlogis(n_trt_obs_e[trt_i], location = cmu_e[i, trt_obs_e == trt_i], scale = s_e[i])}
}
}
if(x$model_output$dist_e == "negbin") {
tau_e <- x$model_output$model$BUGSoutput$sims.list$tau_e
for(i in 1:n_sim){ e_rep[i, ] <- rnbinom(n_obs_e, prob = tau_e[i] / (tau_e[i] + cmu_e[i, ]), size = tau_e[i])
for(trt_i in trt_lev){ e_rep_trt[[trt_i]][i, ] <- rnbinom(n_trt_obs_e[trt_i], prob = tau_e[i] / (tau_e[i] + cmu_e[i, trt_obs_e == trt_i]), size = tau_e[i])}
}
}
if(x$model_output$dist_e == "exp") {
for(i in 1:n_sim){ e_rep[i, ] <- rexp(n_obs_e, rate = 1 / cmu_e[i, ])
for(trt_i in trt_lev){ e_rep_trt[[trt_i]][i, ] <- rexp(n_trt_obs_e[trt_i], rate = 1 / cmu_e[i, trt_obs_e == trt_i])}
}
}
if(x$model_output$dist_e == "bern") {
for(i in 1:n_sim){ e_rep[i, ] <- rbinom(n_obs_e, size = 1, prob = cmu_e[i, ])
for(trt_i in trt_lev){ e_rep_trt[[trt_i]][i, ] <- rbinom(n_trt_obs_e[trt_i], size = 1, prob = cmu_e[i, trt_obs_e == trt_i])}
}
}
if(x$model_output$dist_e == "pois") {
for(i in 1:n_sim){ e_rep[i, ] <- rpois(n_obs_e, lambda = cmu_e[i, ])
for(trt_i in trt_lev){ e_rep_trt[[trt_i]][i, ] <- rpois(n_trt_obs_e[trt_i], lambda = cmu_e[i, trt_obs_e == trt_i])}
}
}
if(x$model_output$dist_e == "gamma") {
ctau_e <- x$model_output$model$BUGSoutput$sims.list$ctau_e[, me_type == "observed"]
for(i in 1:n_sim){ e_rep[i, ] <- rgamma(n_obs_e, shape = cmu_e[i, ] * ctau_e[i, ], rate = ctau_e[i, ])
for(trt_i in trt_lev){ e_rep_trt[[trt_i]][i, ] <- rgamma(n_trt_obs_e[trt_i], shape = cmu_e[i, trt_obs_e == trt_i] * ctau_e[i, trt_obs_e == trt_i], rate = ctau_e[i, trt_obs_e == trt_i])}
}
}
if(x$model_output$dist_e == "weib") {
ctau_e <- x$model_output$model$BUGSoutput$sims.list$ctau_e[, me_type == "observed"]
for(i in 1:n_sim){ e_rep[i, ] <- rweibull(n_obs_e, shape = ctau_e[i, ], scale = cmu_e[i, ] / exp(lgamma(1 + 1 / ctau_e[i, ])))
for(trt_i in trt_lev){ e_rep_trt[[trt_i]][i, ] <- rweibull(n_trt_obs_e[trt_i], shape = ctau_e[i, trt_obs_e == trt_i], scale = cmu_e[i, trt_obs_e == trt_i] / exp(lgamma(1 + 1 / ctau_e[i, trt_obs_e == trt_i])))}
}
}
if(x$model_output$dist_e == "beta") {
ctau_e <- x$model_output$model$BUGSoutput$sims.list$ctau_e[, me_type == "observed"]
for(i in 1:n_sim){ e_rep[i, ] <- rbeta(n_obs_e, shape1 = cmu_e[i, ] * ctau_e[i, ], shape2 = (1 - cmu_e[i, ]) * ctau_e[i, ])
for(trt_i in trt_lev){ e_rep_trt[[trt_i]][i, ] <- rbeta(n_trt_obs_e[trt_i], shape1 = cmu_e[i, trt_obs_e == trt_i] * ctau_e[i, trt_obs_e == trt_i], shape2 = (1 - cmu_e[i, trt_obs_e == trt_i]) * ctau_e[i, trt_obs_e == trt_i])}
}
}
if(x$model_output$dist_c == "norm") {
for(i in 1:n_sim){ c_rep[i, ] <- rnorm(n_obs_c, mean = cmu_c[i, ], sd = s_c[i])
for(trt_i in trt_lev){ c_rep_trt[[trt_i]][i, ] <- rnorm(n_trt_obs_c[trt_i], mean = cmu_c[i, trt_obs_c == trt_i], sd = s_c[i])}
}
}
if(x$model_output$dist_c == "gamma") {
ctau_c <- x$model_output$model$BUGSoutput$sims.list$ctau_c[, mc_type == "observed"]
for(i in 1:n_sim){ c_rep[i, ] <- rgamma(n_obs_c, shape = cmu_c[i, ] * ctau_c[i, ], rate = ctau_c[i, ])
for(trt_i in trt_lev){ c_rep_trt[[trt_i]][i, ] <- rgamma(n_trt_obs_c[trt_i], shape = cmu_c[i, trt_obs_c == trt_i] * ctau_c[i, trt_obs_c == trt_i], rate = ctau_c[i, trt_obs_c == trt_i])}
}
}
if(x$model_output$dist_c == "lnorm") {
for(i in 1:n_sim){ c_rep[i, ] <- rlnorm(n_obs_c, meanlog = clmu_c[i, ], sdlog = 1 / sqrt(ltau_c[i]))
for(trt_i in trt_lev){ c_rep_trt[[trt_i]][i, ] <- rlnorm(n_trt_obs_c[trt_i], meanlog = clmu_c[i, trt_obs_c == trt_i], sdlog = 1 / sqrt(ltau_c[i]))}
}
}
}
if(x$model_output$method == "HURDLE") {
cmu_e <- round(x$model_output$model$BUGSoutput$sims.list$cmu_e[, me_type == "observed"], digits = digits)
s_e <- x$model_output$model$BUGSoutput$sims.list$s_e
if(x$model_output$dist_c == "lnorm") {
clmu_c <- x$model_output$model$BUGSoutput$sims.list$clmu_c[, mc_type == "observed"]
ltau_c <- x$model_output$model$BUGSoutput$sims.list$ltau_c
} else {
cmu_c <- round(x$model_output$model$BUGSoutput$sims.list$cmu_c[, mc_type == "observed"], digits = digits)
s_c <- x$model_output$model$BUGSoutput$sims.list$s_c
}
if(!is.null(x$model_output$mean_nostr_effects)) {
str_e <- x$model_output$se
if(x$model_output$dist_e %in% c("gamma", "weib", "exp", "pois", "negbin")) {
str_e <- round(exp(str_e), digits = digits)
}
if(x$model_output$dist_e %in% c("beta")) {
str_e <- round(exp(str_e) / (1 + exp(str_e)), digits = digits)}
} else { str_e <- NULL}
if(!is.null(x$model_output$mean_nostr_costs)) {
str_c <- x$model_output$sc
if(x$model_output$dist_c %in% c("gamma", "lnorm")) {
str_c <- round(exp(str_c), digits = digits)}
} else { str_c <- NULL}
if(!is.null(str_e)) {
s_eff <- x$data_set$data_raw$se
s_eff[is.na(s_eff)] <- -999999
se_type <- ifelse(s_eff == 1, "structural", s_eff)
se_type <- ifelse(s_eff == 0, "non-structural", se_type)
e_obs_ns <- x$data_set$data_raw$e[se_type == "non-structural"]
e_obs_s <- x$data_set$data_raw$e[se_type == "structural"]
trt_obs_ns_e <- x$data_set$data_raw$data$trt[se_type == "non-structural"]
trt_obs_s_e <- x$data_set$data_raw$data$trt[se_type == "structural"]
n_trt_obs_ns_e <- x$data_set$data_raw$n_ns_eff
n_trt_obs_s_e <- x$data_set$data_raw$n_s_eff
n_obs_ns_e <- sum(n_trt_obs_ns_e)
n_obs_s_e <- sum(n_trt_obs_s_e)
e_ns_rep <- e_ns_rep <- matrix(NA, nrow = n_sim, n_obs_ns_e)
e_s_rep <- e_s_rep <- matrix(NA, nrow = n_sim, n_obs_s_e)
e_obs_ns_trt <- lapply(x$data_set$data_raw$s_efft, function(x) x[!is.na(x) & x == 0])
e_obs_s_trt <- lapply(x$data_set$data_raw$s_efft, function(x) x[!is.na(x) & x == 1])
e_rep_ns_trt <- lapply(x$data_set$data_raw$s_efft, function(x) matrix(NA, nrow = n_sim, ncol = length(x[ !is.na(x) & x == 0])))
e_rep_s_trt <- lapply(x$data_set$data_raw$s_efft, function(x) matrix(NA, nrow = n_sim, ncol = length(x[ !is.na(x) & x == 1])))
}
if(!is.null(str_c)) {
s_cost <- x$data_set$data_raw$sc
s_cost[is.na(s_cost)] <- -999999
sc_type <- ifelse(s_cost == 1, "structural", s_cost)
sc_type <- ifelse(s_cost == 0, "non-structural", sc_type)
c_obs_ns <- x$data_set$data_raw$c[sc_type == "non-structural"]
c_obs_s <- x$data_set$data_raw$c[sc_type == "structural"]
trt_obs_ns_c <- x$data_set$data_raw$data$trt[sc_type == "non-structural"]
trt_obs_s_c <- x$data_set$data_raw$data$trt[sc_type == "structural"]
n_trt_obs_ns_c <- x$data_set$data_raw$n_ns_cost
n_trt_obs_s_c <- x$data_set$data_raw$n_s_cost
n_obs_ns_c <- sum(n_trt_obs_ns_c)
n_obs_s_c <- sum(n_trt_obs_s_c)
c_ns_rep <- c_ns_rep <- matrix(NA, nrow = n_sim, n_obs_ns_c)
c_s_rep <- c_s_rep <- matrix(NA, nrow = n_sim, n_obs_s_c)
c_obs_ns_trt <- lapply(x$data_set$data_raw$s_costt, function(x) x[!is.na(x) & x == 0])
c_obs_s_trt <- lapply(x$data_set$data_raw$s_costt, function(x) x[!is.na(x) & x == 1])
c_rep_ns_trt <- lapply(x$data_set$data_raw$s_costt, function(x) matrix(NA, nrow = n_sim, ncol = length(x[ !is.na(x) & x == 0])))
c_rep_s_trt <- lapply(x$data_set$data_raw$s_costt, function(x) matrix(NA, nrow = n_sim, ncol = length(x[ !is.na(x) & x == 1])))
}
if(is.null(str_e)) {
if(x$model_output$dist_e == "norm") {
for(i in 1:n_sim){ e_rep[i, ] <- rnorm(n_obs_e, mean = cmu_e[i, ], sd = s_e[i])
for(trt_i in trt_lev){ e_rep_trt[[trt_i]][i, ] <- rnorm(n_trt_obs_e[trt_i], mean = cmu_e[i, trt_obs_e == trt_i], sd = s_e[i])}
}
}
if(x$model_output$dist_e == "logis") {
for(i in 1:n_sim){ e_rep[i, ] <- rlogis(n_obs_e, location = cmu_e[i, ], scale = s_e[i])
for(trt_i in trt_lev){ e_rep_trt[[trt_i]][i, ] <- rlogis(n_trt_obs_e[trt_i], location = cmu_e[i, trt_obs_e == trt_i], scale = s_e[i])}
}
}
if(x$model_output$dist_e == "negbin") {
tau_e <- x$model_output$model$BUGSoutput$sims.list$tau_e
for(i in 1:n_sim){ e_rep[i, ] <- rnbinom(n_obs_e, prob = tau_e[i] / (tau_e[i] + cmu_e[i, ]), size = tau_e[i])
for(trt_i in trt_lev){ e_rep_trt[[trt_i]][i, ] <- rnbinom(n_trt_obs_e[trt_i], prob = tau_e[i] / (tau_e[i] + cmu_e[i, trt_obs_e == trt_i]), size = tau_e[i])}
}
}
if(x$model_output$dist_e == "exp") {
for(i in 1:n_sim){ e_rep[i, ] <- rexp(n_obs_e, rate = 1 / cmu_e[i, ])
for(trt_i in trt_lev){ e_rep_trt[[trt_i]][i, ] <- rexp(n_trt_obs_e[trt_i], rate = 1 / cmu_e[i, trt_obs_e == trt_i])}
}
}
if(x$model_output$dist_e == "bern") {
for(i in 1:n_sim){ e_rep[i, ] <- rbinom(n_obs_e, size = 1, prob = cmu_e[i, ])
for(trt_i in trt_lev){ e_rep_trt[[trt_i]][i, ] <- rbinom(n_trt_obs_e[trt_i], size = 1, prob = cmu_e[i, trt_obs_e == trt_i])}
}
}
if(x$model_output$dist_e == "pois") {
for(i in 1:n_sim){ e_rep[i, ] <- rpois(n_obs_e, lambda = cmu_e[i, ])
for(trt_i in trt_lev){ e_rep_trt[[trt_i]][i, ] <- rpois(n_trt_obs_e[trt_i], lambda = cmu_e[i, trt_obs_e == trt_i])}
}
}
if(x$model_output$dist_e == "gamma") {
ctau_e <- x$model_output$model$BUGSoutput$sims.list$ctau_e[, me_type == "observed"]
for(i in 1:n_sim){ e_rep[i, ] <- rgamma(n_obs_e, shape = cmu_e[i, ] * ctau_e[i, ], rate = ctau_e[i, ])
for(trt_i in trt_lev){ e_rep_trt[[trt_i]][i, ] <- rgamma(n_trt_obs_e[trt_i], shape = cmu_e[i, trt_obs_e == trt_i] * ctau_e[i, trt_obs_e == trt_i], rate = ctau_e[i, trt_obs_e == trt_i])}
}
}
if(x$model_output$dist_e == "weib") {
ctau_e <- x$model_output$model$BUGSoutput$sims.list$ctau_e[, me_type == "observed"]
for(i in 1:n_sim){ e_rep[i, ] <- rweibull(n_obs_e, shape = ctau_e[i, ], scale = cmu_e[i, ] / exp(lgamma(1 + 1 / ctau_e[i, ])))
for(trt_i in trt_lev){ e_rep_trt[[trt_i]][i, ] <- rweibull(n_trt_obs_e[trt_i], shape = ctau_e[i, trt_obs_e == trt_i], scale = cmu_e[i, trt_obs_e == trt_i] / exp(lgamma(1 + 1 / ctau_e[i, trt_obs_e == trt_i])))}
}
}
if(x$model_output$dist_e == "beta") {
ctau_e <- x$model_output$model$BUGSoutput$sims.list$ctau_e[, me_type == "observed"]
for(i in 1:n_sim){ e_rep[i, ] <- rbeta(n_obs_e, shape1 = cmu_e[i, ] * ctau_e[i, ], shape2 = (1 - cmu_e[i, ]) * ctau_e[i, ])
for(trt_i in trt_lev){ e_rep_trt[[trt_i]][i, ] <- rbeta(n_trt_obs_e[trt_i], shape1 = cmu_e[i, trt_obs_e == trt_i] * ctau_e[i, trt_obs_e == trt_i], shape2 = (1 - cmu_e[i, trt_obs_e == trt_i]) * ctau_e[i, trt_obs_e == trt_i])}
}
}
}
if(is.null(str_c)) {
if(x$model_output$dist_c == "norm") {
for(i in 1:n_sim){ c_rep[i, ] <- rnorm(n_obs_c, mean = cmu_c[i, ], sd = s_c[i])
for(trt_i in trt_lev){ c_rep_trt[[trt_i]][i, ] <- rnorm(n_trt_obs_c[trt_i], mean = cmu_c[i, trt_obs_c == trt_i], sd = s_c[i])}
}
}
if(x$model_output$dist_c == "gamma") {
ctau_c <- x$model_output$model$BUGSoutput$sims.list$ctau_c[, mc_type == "observed"]
for(i in 1:n_sim){ c_rep[i, ] <- rgamma(n_obs_c, shape = cmu_c[i, ] * ctau_c[i, ], rate = ctau_c[i, ])
for(trt_i in trt_lev){ c_rep_trt[[trt_i]][i, ] <- rgamma(n_trt_obs_c[trt_i], shape = cmu_c[i, trt_obs_c == trt_i] * ctau_c[i, trt_obs_c == trt_i], rate = ctau_c[i, trt_obs_c == trt_i])}
}
}
if(x$model_output$dist_c == "lnorm") {
for(i in 1:n_sim){ c_rep[i, ] <- rlnorm(n_obs_c, meanlog = clmu_c[i, ], sdlog = 1 / sqrt(ltau_c[i]))
for(trt_i in trt_lev){ c_rep_trt[[trt_i]][i, ] <- rlnorm(n_trt_obs_c[trt_i], meanlog = clmu_c[i, trt_obs_c == trt_i], sdlog = 1 / sqrt(ltau_c[i]))}
}
}
}
if(!is.null(str_e)) {
if(x$model_output$dist_e == "norm") {
for(i in 1:n_sim){
e_ns_rep[i, ] <- rnorm(n_obs_ns_e, mean = cmu_e[i, ][cmu_e[i, ] != str_e], sd = s_e[i, 1])
e_s_rep[i, ] <- rnorm(n_obs_s_e, mean = cmu_e[i, ][cmu_e[i, ] == str_e], sd = s_e[i, 2])
for(trt_i in trt_lev){
e_rep_ns_trt[[trt_i]][i, ] <- rnorm(n_trt_obs_ns_e[trt_i], mean = cmu_e[i, trt_obs_ns_e == trt_i][cmu_e[i, trt_obs_ns_e == trt_i] != str_e], sd = s_e[i, 1])
e_rep_s_trt[[trt_i]][i, ] <- rnorm(n_trt_obs_s_e[trt_i], mean = cmu_e[i, trt_obs_s_e == trt_i][cmu_e[i, trt_obs_s_e == trt_i] != str_e], sd = s_e[i, 2])
}
}
}
if(x$model_output$dist_e == "logis") {
for(i in 1:n_sim){
e_ns_rep[i, ] <- rlogis(n_obs_ns_e, location = cmu_e[i, ][cmu_e[i, ] != str_e], scale = s_e[i, 1])
e_s_rep[i, ] <- rlogis(n_obs_s_e, location = cmu_e[i, ][cmu_e[i, ] == str_e], scale = s_e[i, 2])
for(trt_i in trt_lev){
e_rep_ns_trt[[trt_i]][i, ] <- rlogis(n_trt_obs_ns_e[trt_i], location = cmu_e[i, trt_obs_ns_e == trt_i][cmu_e[i, trt_obs_ns_e == trt_i] != str_e], scale = s_e[i, 1])
e_rep_s_trt[[trt_i]][i, ] <- rlogis(n_trt_obs_s_e[trt_i], location = cmu_e[i, trt_obs_s_e == trt_i][cmu_e[i, trt_obs_s_e == trt_i] == str_e], scale = s_e[i, 2])
}
}
}
if(x$model_output$dist_e == "negbin") {
tau_e <- x$model_output$model$BUGSoutput$sims.list$tau_e
for(i in 1:n_sim){
e_ns_rep[i, ] <- rnbinom(n_obs_ns_e, prob = tau_e[i, 1] / (tau_e[i, 1] + cmu_e[i, ][cmu_e[i, ] != str_e]), size = tau_e[i, 1])
e_s_rep[i, ] <- rnbinom(n_obs_s_e, prob = tau_e[i, 2] / (tau_e[i, 2] + cmu_e[i, ][cmu_e[i, ] == str_e]), size = tau_e[i, 2])
for(trt_i in trt_lev){
e_rep_ns_trt[[trt_i]][i, ] <- rnbinom(n_trt_obs_ns_e[trt_i], prob = tau_e[i, 1] / (tau_e[i, 1] + cmu_e[i, trt_obs_ns_e == trt_i][cmu_e[i, trt_obs_ns_e == trt_i] != str_e]), size = tau_e[i, 1])
e_rep_s_trt[[trt_i]][i, ] <- rnbinom(n_trt_obs_s_e[trt_i], prob = tau_e[i, 2] / (tau_e[i, 2] + cmu_e[i, trt_obs_s_e == trt_i][cmu_e[i, trt_obs_s_e == trt_i] == str_e]), size = tau_e[i, 2])
}
}
}
if(x$model_output$dist_e == "exp") {
for(i in 1:n_sim){
e_ns_rep[i, ] <- rexp(n_obs_ns_e, rate = 1 / cmu_e[i, ][cmu_e[i, ] != str_e])
e_s_rep[i, ] <- rexp(n_obs_s_e, rate = 1 / cmu_e[i, ][cmu_e[i, ] == str_e])
for(trt_i in trt_lev){
e_rep_ns_trt[[trt_i]][i, ] <- rexp(n_trt_obs_ns_e[trt_i], rate = 1 / cmu_e[i, trt_obs_ns_e == trt_i][cmu_e[i, trt_obs_ns_e == trt_i] != str_e])
e_rep_s_trt[[trt_i]][i, ] <- rexp(n_trt_obs_s_e[trt_i], rate = 1 / cmu_e[i, trt_obs_s_e == trt_i][cmu_e[i, trt_obs_s_e == trt_i] == str_e])
}
}
}
if(x$model_output$dist_e == "bern") {
for(i in 1:n_sim){
e_ns_rep[i, ] <- rbinom(n_obs_ns_e, size = 1, prob = cmu_e[i, ][cmu_e[i, ] != str_e])
e_s_rep[i, ] <- rbinom(n_obs_s_e, size = 1, prob = cmu_e[i, ][cmu_e[i, ] == str_e])
for(trt_i in trt_lev){
e_rep_ns_trt[[trt_i]][i, ] <- rbinom(n_trt_obs_ns_e[trt_i], size = 1, prob = cmu_e[i, trt_obs_ns_e == trt_i][cmu_e[i, trt_obs_ns_e == trt_i] != str_e])
e_rep_s_trt[[trt_i]][i, ] <- rbinom(n_trt_obs_s_e[trt_i], size = 1, prob = cmu_e[i, trt_obs_s_e == trt_i][cmu_e[i, trt_obs_s_e == trt_i] == str_e])
}
}
}
if(x$model_output$dist_e == "pois") {
for(i in 1:n_sim){
e_ns_rep[i, ] <- rpois(n_obs_ns_e, lambda = cmu_e[i, ][cmu_e[i, ] != str_e])
e_s_rep[i, ] <- rpois(n_obs_s_e, lambda = cmu_e[i, ][cmu_e[i, ] == str_e])
for(trt_i in trt_lev){
e_rep_ns_trt[[trt_i]][i, ] <- rpois(n_trt_obs_ns_e[trt_i], lambda = cmu_e[i, trt_obs_ns_e == trt_i][cmu_e[i, trt_obs_ns_e == trt_i] != str_e])
e_rep_s_trt[[trt_i]][i, ] <- rpois(n_trt_obs_s_e[trt_i], lambda = cmu_e[i, trt_obs_s_e == trt_i][cmu_e[i, trt_obs_s_e == trt_i] == str_e])
}
}
}
if(x$model_output$dist_e == "gamma") {
ctau_e <- x$model_output$model$BUGSoutput$sims.list$ctau_e[, me_type == "observed"]
for(i in 1:n_sim){
e_ns_rep[i, ] <- rgamma(n_obs_ns_e, shape = cmu_e[i, ][cmu_e[i, ] != str_e] * ctau_e[i, ][cmu_e[i, ] != str_e], rate = ctau_e[i, ][cmu_e[i, ] != str_e])
e_s_rep[i, ] <- rgamma(n_obs_s_e, shape = cmu_e[i, ][cmu_e[i, ] == str_e] * ctau_e[i, ][cmu_e[i, ] == str_e], rate = ctau_e[i, ][cmu_e[i, ] == str_e])
for(trt_i in trt_lev){
e_rep_ns_trt[[trt_i]][i, ] <- rgamma(n_trt_obs_ns_e[trt_i], shape = cmu_e[i, trt_obs_ns_e == trt_i][cmu_e[i, trt_obs_ns_e == trt_i] != str_e] * ctau_e[i, trt_obs_ns_e == trt_i][cmu_e[i, trt_obs_ns_e == trt_i] != str_e],
rate = ctau_e[i, trt_obs_ns_e == trt_i][cmu_e[i, trt_obs_ns_e == trt_i] != str_e])
e_rep_s_trt[[trt_i]][i, ] <- rgamma(n_trt_obs_s_e[trt_i], shape = cmu_e[i, trt_obs_s_e == trt_i][cmu_e[i, trt_obs_s_e == trt_i] == str_e] * ctau_e[i, trt_obs_s_e == trt_i][cmu_e[i, trt_obs_s_e == trt_i] == str_e],
rate = ctau_e[i, trt_obs_s_e == trt_i][cmu_e[i, trt_obs_s_e == trt_i] == str_e])
}
}
}
if(x$model_output$dist_e == "weib") {
ctau_e <- x$model_output$model$BUGSoutput$sims.list$ctau_e[, me_type == "observed"]
for(i in 1:n_sim){
e_ns_rep[i, ] <- rweibull(n_obs_ns_e, shape = ctau_e[i, ][cmu_e[i, ] != str_e], scale = cmu_e[i, ][cmu_e[i, ] != str_e] / exp(lgamma(1 + 1 / ctau_e[i, ][cmu_e[i, ] != str_e])))
e_s_rep[i, ] <- rweibull(n_obs_s_e, shape = ctau_e[i, ][cmu_e[i, ] == str_e], scale = cmu_e[i, ][cmu_e[i, ] == str_e] / exp(lgamma(1 + 1 / ctau_e[i, ][cmu_e[i, ] == str_e])))
for(trt_i in trt_lev){
e_rep_ns_trt[[trt_i]][i, ] <- rweibull(n_trt_obs_ns_e[trt_i], shape = ctau_e[i, trt_obs_ns_e == trt_i][cmu_e[i, trt_obs_ns_e == trt_i] != str_e],
scale = cmu_e[i, trt_obs_ns_e == trt_i][cmu_e[i, trt_obs_ns_e == trt_i] != str_e] / exp(lgamma(1 + 1 / ctau_e[i, trt_obs_ns_e == trt_i][cmu_e[i, trt_obs_ns_e == trt_i] != str_e])))
e_rep_s_trt[[trt_i]][i, ] <- rweibull(n_trt_obs_s_e[trt_i], shape = ctau_e[i, trt_obs_s_e == trt_i][cmu_e[i, trt_obs_s_e == trt_i] == str_e],
scale = cmu_e[i, trt_obs_s_e == trt_i][cmu_e[i, trt_obs_s_e == trt_i] == str_e] / exp(lgamma(1 + 1 / ctau_e[i, trt_obs_s_e == trt_i][cmu_e[i, trt_obs_s_e == trt_i] == str_e])))
}
}
}
if(x$model_output$dist_e == "beta") {
ctau_e <- x$model_output$model$BUGSoutput$sims.list$ctau_e[, me_type == "observed"]
for(i in 1:n_sim){
e_ns_rep[i, ] <- rbeta(n_obs_ns_e, shape1 = cmu_e[i, ][cmu_e[i, ] != str_e] * ctau_e[i, ][cmu_e[i, ] != str_e], shape2 = (1 - cmu_e[i, ][cmu_e[i, ] != str_e]) * ctau_e[i, ][cmu_e[i, ] != str_e])
e_s_rep[i, ] <- rbeta(n_obs_s_e, shape1 = cmu_e[i, ][cmu_e[i, ] == str_e] * ctau_e[i, ][cmu_e[i, ] == str_e], shape2 = (1 - cmu_e[i, ][cmu_e[i, ] == str_e]) * ctau_e[i, ][cmu_e[i, ] == str_e])
for(trt_i in trt_lev){
e_rep_ns_trt[[trt_i]][i, ] <- rbeta(n_trt_obs_ns_e[trt_i], shape1 = cmu_e[i, trt_obs_ns_e == trt_i][cmu_e[i, trt_obs_ns_e == trt_i] != str_e] * ctau_e[i, trt_obs_ns_e == trt_i][cmu_e[i, trt_obs_ns_e == trt_i] != str_e],
shape2 = (1 - cmu_e[i, trt_obs_ns_e == trt_i][cmu_e[i, trt_obs_ns_e == trt_i] != str_e]) * ctau_e[i, trt_obs_ns_e == trt_i][cmu_e[i, trt_obs_ns_e == trt_i] != str_e])
e_rep_s_trt[[trt_i]][i, ] <- rbeta(n_trt_obs_s_e[trt_i], shape1 = cmu_e[i, trt_obs_ns_e == trt_i][cmu_e[i, trt_obs_ns_e == trt_i] == str_e] * ctau_e[i, trt_obs_ns_e == trt_i][cmu_e[i, trt_obs_ns_e == trt_i] == str_e],
shape2 = (1 - cmu_e[i, trt_obs_ns_e == trt_i][cmu_e[i, trt_obs_ns_e == trt_i] == str_e]) * ctau_e[i, trt_obs_ns_e == trt_i][cmu_e[i, trt_obs_ns_e == trt_i] == str_e])
}
}
}
e_rep <- cbind(e_ns_rep, e_s_rep)
e_rep_trt <- mapply(cbind, e_rep_ns_trt, e_rep_s_trt)
}
if(!is.null(str_c)) {
if(x$model_output$dist_c == "norm") {
for(i in 1:n_sim){
c_ns_rep[i, ] <- rnorm(n_obs_ns_c, mean = cmu_c[i, ][cmu_c[i, ] != str_c], sd = s_c[i, 1])
c_s_rep[i, ] <- rnorm(n_obs_s_c, mean = cmu_c[i, ][cmu_c[i, ] == str_c], sd = s_c[i, 2])
for(trt_i in trt_lev){
c_rep_ns_trt[[trt_i]][i, ] <- rnorm(n_trt_obs_ns_c[trt_i], mean = cmu_c[i, trt_obs_ns_c == trt_i][cmu_c[i, trt_obs_ns_c == trt_i] != str_c], sd = s_c[i, 1])
c_rep_s_trt[[trt_i]][i, ] <- rnorm(n_trt_obs_s_c[trt_i], mean = cmu_c[i, trt_obs_s_c == trt_i][cmu_c[i, trt_obs_s_c == trt_i] != str_c], sd = s_c[i, 2])
}
}
}
if(x$model_output$dist_c == "gamma") {
ctau_c <- x$model_output$model$BUGSoutput$sims.list$ctau_c[, mc_type == "observed"]
for(i in 1:n_sim){
c_ns_rep[i, ] <- rgamma(n_obs_ns_c, shape = cmu_c[i, ][cmu_c[i, ] != str_c] * ctau_c[i, ][cmu_c[i, ] != str_c], rate = ctau_c[i, ][cmu_c[i, ] != str_c])
c_s_rep[i, ] <- rgamma(n_obs_s_c, shape = cmu_c[i, ][cmu_c[i, ] == str_c] * ctau_c[i, ][cmu_c[i, ] == str_c], rate = ctau_c[i, ][cmu_c[i, ] == str_c])
for(trt_i in trt_lev){
c_rep_ns_trt[[trt_i]][i, ] <- rgamma(n_trt_obs_ns_c[trt_i], shape = cmu_c[i, trt_obs_ns_c == trt_i][cmu_c[i, trt_obs_ns_c == trt_i] != str_c] * ctau_c[i, trt_obs_ns_c == trt_i][cmu_c[i, trt_obs_ns_c == trt_i] != str_c],
rate = ctau_c[i, trt_obs_ns_c == trt_i][cmu_c[i, trt_obs_ns_c == trt_i] != str_c])
c_rep_s_trt[[trt_i]][i, ] <- rgamma(n_trt_obs_s_c[trt_i], shape = cmu_c[i, trt_obs_s_c == trt_i][cmu_c[i, trt_obs_s_c == trt_i] == str_c] * ctau_c[i, trt_obs_s_c == trt_i][cmu_c[i, trt_obs_s_c == trt_i] == str_c],
rate = ctau_c[i, trt_obs_s_c == trt_i][cmu_c[i, trt_obs_s_c == trt_i] == str_c])
}
}
}
if(x$model_output$dist_c == "lnorm") {
for(i in 1:n_sim){
c_ns_rep[i, ] <- rlnorm(n_obs_ns_c, meanlog = clmu_c[i, ][clmu_c[i, ] != str_c], sdlog = 1 / sqrt(ltau_c[i, 1]))
c_s_rep[i, ] <- rlnorm(n_obs_s_c, meanlog = clmu_c[i, ][clmu_c[i, ] == str_c], sdlog = 1 / sqrt(ltau_c[i, 2]))
for(trt_i in trt_lev){
c_rep_ns_trt[[trt_i]][i, ] <- rlnorm(n_trt_obs_ns_c[trt_i], meanlog = clmu_c[i, trt_obs_ns_c == trt_i][clmu_c[i, trt_obs_ns_c == trt_i] != str_c], sdlog = 1 / sqrt(ltau_c[i, 1]))
c_rep_s_trt[[trt_i]][i, ] <- rlnorm(n_trt_obs_s_c[trt_i], meanlog = clmu_c[i, trt_obs_s_c == trt_i][clmu_c[i, trt_obs_s_c == trt_i] == str_c], sdlog = 1 / sqrt(ltau_c[i, 2]))
}
}
}
c_rep <- cbind(c_ns_rep, c_s_rep)
c_rep_trt <- mapply(cbind, c_rep_ns_trt, c_rep_s_trt)
}
}
if(x$model_output$method == "PATTERN") {
d_mod <- x$model_output$d_mod
d_mod_e_obs <- d_mod[me_type == "observed"]
d_mod_c_obs <- d_mod[mc_type == "observed"]
d_e_obs <- x$data_set$data_raw$d[me_type == "observed"]
d_c_obs <- x$data_set$data_raw$d[mc_type == "observed"]
pat_obs_e_val <- unique(d_mod_e_obs)
pat_obs_c_val <- unique(d_mod_c_obs)
n_pat_e_obs <- length(unique(d_mod_e_obs))
n_pat_c_obs <- length(unique(d_mod_c_obs))
if(1 %in% x$data_set$data_raw$d) {
e_rep_pat1_trt <- c_rep_pat1_trt <- lapply(x$data_set$data_raw$dt, function(x) matrix(NA, nrow = n_sim, ncol = length(x[x == 1])))
} else { e_rep_pat1_trt <- c_rep_pat1_trt <- NULL}
if(2 %in% x$data_set$data_raw$d) {
c_rep_pat2_trt <- lapply(x$data_set$data_raw$dt, function(x) matrix(NA, nrow = n_sim, ncol = length(x[x == 2])))
} else { c_rep_pat2_trt <- NULL}
if(3 %in% x$data_set$data_raw$d) {
e_rep_pat3_trt <- lapply(x$data_set$data_raw$dt, function(x) matrix(NA, nrow = n_sim, ncol = length(x[x == 3])))
} else { e_rep_pat3_trt <- NULL}
par_e_pat3_ind <- 3; par_c_pat2_ind <- 2
if(all(unique(x$data_set$data_raw$d) %in% c(3, 4)) & length(unique(x$data_set$data_raw$d)) == 2) { par_e_pat3_ind <- 1}
if(all(unique(x$data_set$data_raw$d) %in% c(1, 3)) & length(unique(x$data_set$data_raw$d)) == 2) { par_e_pat3_ind <- 2}
if(all(unique(x$data_set$data_raw$d) %in% c(2, 3)) & length(unique(x$data_set$data_raw$d)) == 2) { par_e_pat3_ind <- 2}
if(all(unique(x$data_set$data_raw$d) %in% c(1, 3, 4)) & length(unique(x$data_set$data_raw$d)) == 3) { par_e_pat3_ind <- 2}
if(all(unique(x$data_set$data_raw$d) %in% c(2, 4)) & length(unique(x$data_set$data_raw$d)) == 2) { par_e_pat2_ind <- 1}
if(all(unique(x$data_set$data_raw$d) %in% c(2, 3)) & length(unique(x$data_set$data_raw$d)) == 2) { par_e_pat2_ind <- 1}
if(all(unique(x$data_set$data_raw$d) %in% c(2, 3, 4)) & length(unique(x$data_set$data_raw$d)) == 3) { par_e_pat2_ind <- 1}
d_mod_e_obs_target <- d_mod_e_obs; d_mod_c_obs_target <- d_mod_c_obs
if(length(unique(d_mod_e_obs)) == 1) {
d_mod_e_obs_target[d_mod_e_obs == unique(d_mod_e_obs)[1]] <- 1}
if(length(unique(d_mod_e_obs)) == 2) {
d_mod_e_obs_target[d_mod_e_obs == unique(d_mod_e_obs)[1]] <- 1
d_mod_e_obs_target[d_mod_e_obs == unique(d_mod_e_obs)[2]] <- 2}
if(length(unique(d_mod_c_obs)) == 1) {
d_mod_c_obs_target[d_mod_c_obs == unique(d_mod_c_obs)[1]] <- 1}
if(length(unique(d_mod_e_obs)) == 2) {
d_mod_c_obs_target[d_mod_c_obs == unique(d_mod_c_obs)[1]] <- 1
d_mod_c_obs_target[d_mod_c_obs == unique(d_mod_c_obs)[2]] <- 2}
p_obs_e_index <- lapply(1:n_pat_e_obs, function(target) which(d_mod_e_obs_target == target))
p_obs_c_index <- lapply(1:n_pat_c_obs, function(target) which(d_mod_c_obs_target == target))
e_rep_pat <- lapply(p_obs_e_index, function(x) matrix(NA, nrow = n_sim, ncol = length(x)))
c_rep_pat <- lapply(p_obs_c_index, function(x) matrix(NA, nrow = n_sim, ncol = length(x)))
cmu_e <- x$model_output$model$BUGSoutput$sims.list$cmu_e[, me_type == "observed"]
s_e <- as.matrix(x$model_output$model$BUGSoutput$sims.list$s_e_p[, pat_obs_e_val])
if(x$model_output$dist_c == "lnorm") {
clmu_c <- x$model_output$model$BUGSoutput$sims.list$clmu_c[, mc_type == "observed"]
ltau_c <- as.matrix(x$model_output$model$BUGSoutput$sims.list$ltau_c_p[, pat_obs_c_val])
} else {
cmu_c <- x$model_output$model$BUGSoutput$sims.list$cmu_c[, mc_type == "observed"]
s_c <- as.matrix(x$model_output$model$BUGSoutput$sims.list$s_c_p[, pat_obs_c_val])
}
if(x$model_output$dist_e == "norm") {
for(i in 1:n_sim){
for(p in 1:n_pat_e_obs){ e_rep_pat[[p]][i, ] <- rnorm(dim(e_rep_pat[[p]])[2], mean = cmu_e[i, unlist(p_obs_e_index[p])], sd = s_e[i, p])}
for(trt_i in trt_lev){
if(!is.null(e_rep_pat1_trt)) { e_rep_pat1_trt[[trt_i]][i, ] <- rnorm(dim(e_rep_pat1_trt[[trt_i]])[2], mean = cmu_e[i, trt_obs_e == trt_i & d_mod_e_obs == 1], sd = s_e[i, 1])}
if(!is.null(e_rep_pat3_trt)) { e_rep_pat3_trt[[trt_i]][i, ] <- rnorm(dim(e_rep_pat3_trt[[trt_i]])[2], mean = cmu_e[i, trt_obs_e == trt_i & d_mod_e_obs == par_e_pat3_ind], sd = s_e[i, par_e_pat3_ind])}
}
}
}
if(x$model_output$dist_e == "logis") {
for(i in 1:n_sim){
for(p in 1:n_pat_e_obs){ e_rep_pat[[p]][i, ] <- rlogis(dim(e_rep_pat[[p]])[2], location = cmu_e[i, unlist(p_obs_e_index[p])], scale = s_e[i, p])}
for(trt_i in trt_lev){
if(!is.null(e_rep_pat1_trt)) { e_rep_pat1_trt[[trt_i]][i, ] <- rlogis(dim(e_rep_pat1_trt[[trt_i]])[2], location = cmu_e[i, trt_obs_e == trt_i & d_mod_e_obs == 1], scale = s_e[i, 1])}
if(!is.null(e_rep_pat3_trt)) { e_rep_pat3_trt[[trt_i]][i, ] <- rlogis(dim(e_rep_pat3_trt[[trt_i]])[2], location = cmu_e[i, trt_obs_e == trt_i & d_mod_e_obs == par_e_pat3_ind], scale = s_e[i, par_e_pat3_ind])}
}
}
}
if(x$model_output$dist_e == "negbin") {
tau_e <- as.matrix(x$model_output$model$BUGSoutput$sims.list$tau_e_p[, pat_obs_e_val])
for(i in 1:n_sim){
for(p in 1:n_pat_e_obs){ e_rep_pat[[p]][i, ] <- rnbinom(dim(e_rep_pat[[p]])[2], prob = tau_e[i, p] / (tau_e[i, p] + cmu_e[i, unlist(p_obs_e_index[p])]), size = tau_e[i, p])}
for(trt_i in trt_lev){
if(!is.null(e_rep_pat1_trt)) { e_rep_pat1_trt[[trt_i]][i, ] <- rnbinom(dim(e_rep_pat1_trt[[trt_i]])[2], prob = tau_e[i, 1] / (tau_e[i, 1] + cmu_e[i, trt_obs_e == trt_i & d_mod_e_obs == 1]), size = tau_e[i, 1])}
if(!is.null(e_rep_pat3_trt)) { e_rep_pat3_trt[[trt_i]][i, ] <- rnbinom(dim(e_rep_pat3_trt[[trt_i]])[2], prob = tau_e[i, par_e_pat3_ind] / (tau_e[i, par_e_pat3_ind] + cmu_e[i, trt_obs_e == trt_i & d_mod_e_obs == par_e_pat3_ind]), size = tau_e[i, par_e_pat3_ind])}
}
}
}
if(x$model_output$dist_e == "exp") {
for(i in 1:n_sim){
for(p in 1:n_pat_e_obs){ e_rep_pat[[p]][i, ] <- rexp(dim(e_rep_pat[[p]])[2], rate = 1 / cmu_e[i, unlist(p_obs_e_index[p])])}
for(trt_i in trt_lev){
if(!is.null(e_rep_pat1_trt)) { e_rep_pat1_trt[[trt_i]][i, ] <- rexp(dim(e_rep_pat1_trt[[trt_i]])[2], rate = 1 / cmu_e[i, trt_obs_e == trt_i & d_mod_e_obs == 1])}
if(!is.null(e_rep_pat3_trt)) { e_rep_pat3_trt[[trt_i]][i, ] <- rexp(dim(e_rep_pat3_trt[[trt_i]])[2], rate = 1 / cmu_e[i, trt_obs_e == trt_i & d_mod_e_obs == par_e_pat3_ind])}
}
}
}
if(x$model_output$dist_e == "bern") {
for(i in 1:n_sim){
for(p in 1:n_pat_e_obs){ e_rep_pat[[p]][i, ] <- rbinom(dim(e_rep_pat[[p]])[2], size = 1, prob = cmu_e[i, unlist(p_obs_e_index[p])])}
for(trt_i in trt_lev){
if(!is.null(e_rep_pat1_trt)) { e_rep_pat1_trt[[trt_i]][i, ] <- rbinom(dim(e_rep_pat1_trt[[trt_i]])[2], size = 1, prob = cmu_e[i, trt_obs_e == trt_i & d_mod_e_obs == 1])}
if(!is.null(e_rep_pat3_trt)) { e_rep_pat3_trt[[trt_i]][i, ] <- rbinom(dim(e_rep_pat3_trt[[trt_i]])[2], size = 1, prob = cmu_e[i, trt_obs_e == trt_i & d_mod_e_obs == par_e_pat3_ind])}
}
}
}
if(x$model_output$dist_e == "pois") {
for(i in 1:n_sim){
for(p in 1:n_pat_e_obs){ e_rep_pat[[p]][i, ] <- rpois(dim(e_rep_pat[[p]])[2], lambda = cmu_e[i, unlist(p_obs_e_index[p])])}
for(trt_i in trt_lev){
if(!is.null(e_rep_pat1_trt)) { e_rep_pat1_trt[[trt_i]][i, ] <- rpois(dim(e_rep_pat1_trt[[trt_i]])[2], lambda = cmu_e[i, trt_obs_e == trt_i & d_mod_e_obs == 1])}
if(!is.null(e_rep_pat3_trt)) { e_rep_pat3_trt[[trt_i]][i, ] <- rpois(dim(e_rep_pat3_trt[[trt_i]])[2], lambda = cmu_e[i, trt_obs_e == trt_i & d_mod_e_obs == par_e_pat3_ind])}
}
}
}
if(x$model_output$dist_e == "gamma") {
ctau_e <- x$model_output$model$BUGSoutput$sims.list$ctau_e[, me_type == "observed"]
for(i in 1:n_sim){
for(p in 1:n_pat_e_obs){ e_rep_pat[[p]][i, ] <- rgamma(dim(e_rep_pat[[p]])[2], shape = cmu_e[i, unlist(p_obs_e_index[p])] * ctau_e[i, unlist(p_obs_e_index[p])], rate = ctau_e[i, unlist(p_obs_e_index[p])])}
for(trt_i in trt_lev){
if(!is.null(e_rep_pat1_trt)) { e_rep_pat1_trt[[trt_i]][i, ] <- rgamma(dim(e_rep_pat1_trt[[trt_i]])[2], shape = cmu_e[i, trt_obs_e == trt_i & d_mod_e_obs == 1] * ctau_e[i, trt_obs_e == trt_i & d_mod_e_obs == 1], rate = ctau_e[i, trt_obs_e == trt_i & d_mod_e_obs == 1])}
if(!is.null(e_rep_pat3_trt)) { e_rep_pat3_trt[[trt_i]][i, ] <- rgamma(dim(e_rep_pat3_trt[[trt_i]])[2], shape = cmu_e[i, trt_obs_e == trt_i & d_mod_e_obs == par_e_pat3_ind] * ctau_e[i, trt_obs_e == trt_i & d_mod_e_obs == par_e_pat3_ind], rate = ctau_e[i, trt_obs_e == trt_i & d_mod_e_obs == par_e_pat3_ind])}
}
}
}
if(x$model_output$dist_e == "weib") {
ctau_e <- x$model_output$model$BUGSoutput$sims.list$ctau_e[, me_type == "observed"]
for(i in 1:n_sim){
for(p in 1:n_pat_e_obs){ e_rep_pat[[p]][i, ] <- rweibull(dim(e_rep_pat[[p]])[2], shape = ctau_e[i, unlist(p_obs_e_index[p])], scale = cmu_e[i, unlist(p_obs_e_index[p])] / exp(lgamma(1 + 1 / ctau_e[i, unlist(p_obs_e_index[p])])))}
for(trt_i in trt_lev){
if(!is.null(e_rep_pat1_trt)) { e_rep_pat1_trt[[trt_i]][i, ] <- rweibull(dim(e_rep_pat1_trt[[trt_i]])[2], shape = ctau_e[i, trt_obs_e == trt_i & d_mod_e_obs == 1], scale = cmu_e[i, trt_obs_e == trt_i & d_mod_e_obs == 1] / exp(lgamma(1 + 1 / ctau_e[i, trt_obs_e == trt_i & d_mod_e_obs == par_e_pat3_ind])))}
if(!is.null(e_rep_pat3_trt)) { e_rep_pat3_trt[[trt_i]][i, ] <- rweibull(dim(e_rep_pat3_trt[[trt_i]])[2], shape = ctau_e[i, trt_obs_e == trt_i & d_mod_e_obs == par_e_pat3_ind], scale = cmu_e[i, trt_obs_e == trt_i & d_mod_e_obs == par_e_pat3_ind] / exp(lgamma(1 + 1 / ctau_e[i, trt_obs_e == trt_i & d_mod_e_obs == par_e_pat3_ind])))}
}
}
}
if(x$model_output$dist_e == "beta") {
ctau_e <- x$model_output$model$BUGSoutput$sims.list$ctau_e[, me_type == "observed"]
for(i in 1:n_sim){
for(p in 1:n_pat_e_obs){ e_rep_pat[[p]][i, ] <- rbeta(dim(e_rep_pat[[p]])[2], shape1 = cmu_e[i, unlist(p_obs_e_index[p])] * ctau_e[i, unlist(p_obs_e_index[p])], shape2 = (1 - cmu_e[i, unlist(p_obs_e_index[p])]) * ctau_e[i, unlist(p_obs_e_index[p])])}
for(trt_i in trt_lev){
if(!is.null(e_rep_pat1_trt)) { e_rep_pat1_trt[[trt_i]][i, ] <- rbeta(dim(e_rep_pat1_trt[[trt_i]])[2], shape1 = cmu_e[i, trt_obs_e == trt_i & d_mod_e_obs == 1] * ctau_e[i, trt_obs_e == trt_i & d_mod_e_obs == 1], shape2 = (1 - cmu_e[i, trt_obs_e == trt_i & d_mod_e_obs == 1]) * ctau_e[i, trt_obs_e == trt_i & d_mod_e_obs == 1])}
if(!is.null(e_rep_pat3_trt)) { e_rep_pat3_trt[[trt_i]][i, ] <- rbeta(dim(e_rep_pat3_trt[[trt_i]])[2], shape1 = cmu_e[i, trt_obs_e == trt_i & d_mod_e_obs == par_e_pat3_ind] * ctau_e[i, trt_obs_e == trt_i & d_mod_e_obs == par_e_pat3_ind], shape2 = (1 - cmu_e[i, trt_obs_e == trt_i & d_mod_e_obs == par_e_pat3_ind]) * ctau_e[i, trt_obs_e == trt_i & d_mod_e_obs == par_e_pat3_ind])}
}
}
}
if(x$model_output$dist_c == "norm") {
for(i in 1:n_sim){
for(p in 1:n_pat_c_obs){ c_rep_pat[[p]][i, ] <- rnorm(dim(c_rep_pat[[p]])[2], mean = cmu_c[i, unlist(p_obs_c_index[p])], sd = s_c[i, p])}
for(trt_i in trt_lev){
if(!is.null(c_rep_pat1_trt)) { c_rep_pat1_trt[[trt_i]][i, ] <- rnorm(dim(c_rep_pat1_trt[[trt_i]])[2], mean = cmu_c[i, trt_obs_c == trt_i & d_mod_c_obs == 1], sd = s_c[i, 1])}
if(!is.null(c_rep_pat2_trt)) { c_rep_pat2_trt[[trt_i]][i, ] <- rnorm(dim(c_rep_pat2_trt[[trt_i]])[2], mean = cmu_c[i, trt_obs_c == trt_i & d_mod_c_obs == par_c_pat2_ind], sd = s_c[i, par_c_pat2_ind])}
}
}
}
if(x$model_output$dist_c == "gamma") {
ctau_c <- x$model_output$model$BUGSoutput$sims.list$ctau_c[, mc_type == "observed"]
for(i in 1:n_sim){
for(p in 1:n_pat_c_obs){ c_rep_pat[[p]][i, ] <- rgamma(dim(c_rep_pat[[p]])[2], shape = cmu_c[i, unlist(p_obs_c_index[p])] * ctau_c[i, unlist(p_obs_c_index[p])], rate = ctau_c[i, unlist(p_obs_c_index[p])])}
for(trt_i in trt_lev){
if(!is.null(c_rep_pat1_trt)) { c_rep_pat1_trt[[trt_i]][i, ] <- rgamma(dim(c_rep_pat1_trt[[trt_i]])[2], shape = cmu_c[i, trt_obs_c == trt_i & d_mod_c_obs == 1] * ctau_c[i, trt_obs_c == trt_i & d_mod_c_obs == 1], rate = ctau_c[i, trt_obs_c == trt_i & d_mod_c_obs == 1])}
if(!is.null(c_rep_pat2_trt)) { c_rep_pat2_trt[[trt_i]][i, ] <- rgamma(dim(c_rep_pat2_trt[[trt_i]])[2], shape = cmu_c[i, trt_obs_c == trt_i & d_mod_c_obs == par_c_pat2_ind] * ctau_c[i, trt_obs_c == trt_i & d_mod_c_obs == par_c_pat2_ind], rate = ctau_c[i, trt_obs_c == trt_i & d_mod_c_obs == par_c_pat2_ind])}
}
}
}
if(x$model_output$dist_c == "lnorm") {
for(i in 1:n_sim){
for(p in 1:n_pat_c_obs){ c_rep_pat[[p]][i, ] <- rlnorm(dim(c_rep_pat[[p]])[2], meanlog = clmu_c[i, unlist(p_obs_c_index[p])], sdlog = 1 / sqrt(ltau_c[i, p]))}
for(trt_i in trt_lev){
if(!is.null(c_rep_pat1_trt)) { c_rep_pat1_trt[[trt_i]][i, ] <- rlnorm(dim(c_rep_pat1_trt[[trt_i]])[2], meanlog = clmu_c[i, trt_obs_c == trt_i & d_mod_c_obs == 1], sdlog = 1 / sqrt(ltau_c[i, 1]))}
if(!is.null(c_rep_pat2_trt)) { c_rep_pat2_trt[[trt_i]][i, ] <- rlnorm(dim(c_rep_pat2_trt[[trt_i]])[2], meanlog = clmu_c[i, trt_obs_c == trt_i & d_mod_c_obs == par_c_pat2_ind], sdlog = 1 / sqrt(ltau_c[i, par_c_pat2_ind]))}
}
}
}
e_rep <- do.call(cbind, e_rep_pat)
c_rep <- do.call(cbind, c_rep_pat)
if(!is.null(e_rep_pat1_trt) & !is.null(e_rep_pat3_trt)) { e_rep_trt <- mapply(cbind, e_rep_pat1_trt, e_rep_pat3_trt)}
if(!is.null(e_rep_pat1_trt) & is.null(e_rep_pat3_trt)) { e_rep_trt <- e_rep_pat1_trt}
if(is.null(e_rep_pat1_trt) & !is.null(e_rep_pat3_trt)) { e_rep_trt <- e_rep_pat3_trt}
if(!is.null(c_rep_pat1_trt) & !is.null(c_rep_pat2_trt)) { c_rep_trt <- mapply(cbind, c_rep_pat1_trt, c_rep_pat2_trt)}
if(!is.null(c_rep_pat1_trt) & is.null(c_rep_pat2_trt)) { c_rep_trt <- c_rep_pat1_trt}
if(is.null(c_rep_pat1_trt) & !is.null(c_rep_pat2_trt)) { c_rep_trt <- c_rep_pat2_trt}
}
}
if(x$data_format == "long") {
me_type <- ifelse(x$data_set$data_raw$me_long == 1, "missing", "observed")
mc_type <- ifelse(x$data_set$data_raw$mc_long == 1, "missing", "observed")
e_obs <- x$data_set$data_raw$e_long[me_type == "observed"]
c_obs <- x$data_set$data_raw$c_long[mc_type == "observed"]
trt_obs_e <- x$data_set$data_raw$data_long$trt[me_type == "observed"]
trt_obs_c <- x$data_set$data_raw$data_long$trt[mc_type == "observed"]
time_obs_e <- x$data_set$data_raw$time_long[me_type == "observed"]
time_obs_c <- x$data_set$data_raw$time_long[mc_type == "observed"]
time_trt_obs_e <- x$data_set$data_raw$data_long$time_trt[me_type == "observed"]
time_trt_obs_c <- x$data_set$data_raw$data_long$time_trt[mc_type == "observed"]
n_all <- sum(x$data_set$data_raw$n_long)
time_trt_lev <- unique(x$data_set$data_raw$data_long$time_trt)
time_trt_index_time <- rep(1:n_time, each = n_trt)
names(time_trt_index_time) <- time_trt_lev
n_trt_obs_e <- x$data_set$data_raw$n_obs_e
n_trt_obs_c <- x$data_set$data_raw$n_obs_c
n_time_obs_e <- x$data_set$data_raw$n_obs_time_e
n_time_obs_c <- x$data_set$data_raw$n_obs_time_c
n_time_trt_obs_e <- x$data_set$data_raw$n_obs_time_trt_e
n_time_trt_obs_c <- x$data_set$data_raw$n_obs_time_trt_c
n_obs_e <- sum(n_trt_obs_e)
n_obs_c <- sum(n_trt_obs_c)
n_obs_e_time <- sum(n_time_obs_e)
n_obs_c_time <- sum(n_time_obs_c)
n_obs_e_time_trt <- sum(n_time_trt_obs_e)
n_obs_c_time_trt <- sum(n_time_trt_obs_c)
e_obs_trt <- lapply(x$data_set$data_raw$efft, function(x) x[!is.na(x)])
c_obs_trt <- lapply(x$data_set$data_raw$costt, function(x) x[!is.na(x)])
e_obs_time <- lapply(x$data_set$data_raw$eff_time_list, function(x) x[!is.na(x)])
c_obs_time <- lapply(x$data_set$data_raw$cost_time_list, function(x) x[!is.na(x)])
e_obs_time_trt <- lapply(x$data_set$data_raw$eff_time_trt_list, function(x) x[!is.na(x)])
c_obs_time_trt <- lapply(x$data_set$data_raw$cost_time_trt_list, function(x) x[!is.na(x)])
e_rep_trt <- lapply(x$data_set$data_raw$m_efft, function(x) matrix(NA, nrow = n_sim, ncol = length(x[x == 0])))
c_rep_trt <- lapply(x$data_set$data_raw$m_costt, function(x) matrix(NA, nrow = n_sim, ncol = length(x[x == 0])))
e_rep_time <- lapply(x$data_set$data_raw$eff_time_list, function(x) matrix(NA, nrow = n_sim, ncol = length(x[!is.na(x)])))
c_rep_time <- lapply(x$data_set$data_raw$cost_time_list, function(x) matrix(NA, nrow = n_sim, ncol = length(x[!is.na(x)])))
e_rep_time_trt <- lapply(x$data_set$data_raw$eff_time_trt_list, function(x) matrix(NA, nrow = n_sim, ncol = length(x[!is.na(x)])))
c_rep_time_trt <- lapply(x$data_set$data_raw$cost_time_trt_list, function(x) matrix(NA, nrow = n_sim, ncol = length(x[!is.na(x)])))
if(x$model_output$method == "LMDM") {
cmu_e <- as.matrix(as.data.frame(x$model_output$model$BUGSoutput$sims.list$cmu_e))[, me_type == "observed"]
s_e <- x$model_output$model$BUGSoutput$sims.list$s_e
if(x$model_output$dist_c == "lnorm") {
clmu_c <- as.matrix(as.data.frame(x$model_output$model$BUGSoutput$sims.list$clmu_c))[, mc_type == "observed"]
ltau_c <- x$model_output$model$BUGSoutput$sims.list$ltau_c
} else {
cmu_c <- as.matrix(as.data.frame(x$model_output$model$BUGSoutput$sims.list$cmu_c))[, mc_type == "observed"]
s_c <- x$model_output$model$BUGSoutput$sims.list$s_c
}
if(x$model_output$dist_e == "norm") {
for(i in 1:n_sim){
for(time_i in 1:n_time){ e_rep_time[[time_i]][i, ] <- rnorm(n_time_obs_e[time_i], mean = cmu_e[i, time_obs_e == time_i], sd = s_e[i, time_i])}
for(time_trt_i in time_trt_lev){ e_rep_time_trt[[time_trt_i]][i, ] <- rnorm(n_time_trt_obs_e[time_trt_i], mean = cmu_e[i, time_trt_obs_e == time_trt_i], sd = s_e[i, time_trt_index_time[time_trt_i]])}
}
}
if(x$model_output$dist_e == "logis") {
for(i in 1:n_sim){
for(time_i in 1:n_time){ e_rep_time[[time_i]][i, ] <- rlogis(n_time_obs_e[time_i], location = cmu_e[i, time_obs_e == time_i], scale = s_e[i, time_i])}
for(time_trt_i in time_trt_lev){ e_rep_time_trt[[time_trt_i]][i, ] <- rlogis(n_time_trt_obs_e[time_trt_i], location = cmu_e[i, time_trt_obs_e == time_trt_i], scale = s_e[i, time_trt_index_time[time_trt_i]])}
}
}
if(x$model_output$dist_e == "negbin") {
tau_e <- x$model_output$model$BUGSoutput$sims.list$tau_e
for(i in 1:n_sim){
for(time_i in 1:n_time){ e_rep_time[[time_i]][i, ] <- rnbinom(n_time_obs_e[time_i], prob = tau_e[i, time_i] / (tau_e[i, time_i] + cmu_e[i, time_obs_e == time_i]), size = tau_e[i, time_i])}
for(time_trt_i in time_trt_lev){ e_rep_time_trt[[time_trt_i]][i, ] <- rnbinom(n_time_trt_obs_e[time_trt_i], prob = tau_e[i, time_trt_index_time[time_trt_i]] / (tau_e[i, time_trt_index_time[time_trt_i]] + cmu_e[i, time_trt_obs_e == time_trt_i]), size = tau_e[i, time_trt_index_time[time_trt_i]])}
}
}
if(x$model_output$dist_e == "exp") {
for(i in 1:n_sim){
for(time_i in 1:n_time){ e_rep_time[[time_i]][i, ] <- rexp(n_time_obs_e[time_i], rate = 1 / cmu_e[i, time_obs_e == time_i])}
for(time_trt_i in time_trt_lev){ e_rep_time_trt[[time_trt_i]][i, ] <- rexp(n_time_trt_obs_e[time_trt_i], rate = 1 / cmu_e[i, time_trt_obs_e == time_trt_i])}
}
}
if(x$model_output$dist_e == "bern") {
for(i in 1:n_sim){
for(time_i in 1:n_time){ e_rep_time[[time_i]][i, ] <- rbinom(n_time_obs_e[time_i], size = 1, prob = cmu_e[i, time_obs_e == time_i])}
for(time_trt_i in time_trt_lev){ e_rep_time_trt[[time_trt_i]][i, ] <- rbinom(n_time_trt_obs_e[time_trt_i], size = 1, prob = cmu_e[i, time_trt_obs_e == time_trt_i])}
}
}
if(x$model_output$dist_e == "pois") {
for(i in 1:n_sim){
for(time_i in 1:n_time){ e_rep_time[[time_i]][i, ] <- rpois(n_time_obs_e[time_i], lambda = cmu_e[i, time_obs_e == time_i])}
for(time_trt_i in time_trt_lev){ e_rep_time_trt[[time_trt_i]][i, ] <- rpois(n_time_trt_obs_e[time_trt_i], lambda = cmu_e[i, time_trt_obs_e == time_trt_i])}
}
}
if(x$model_output$dist_e == "gamma") {
ctau_e <- as.matrix(as.data.frame(x$model_output$model$BUGSoutput$sims.list$ctau_e))[, me_type == "observed"]
for(i in 1:n_sim){
for(time_i in 1:n_time){ e_rep_time[[time_i]][i, ] <- rgamma(n_time_obs_e[time_i], shape = cmu_e[i, time_obs_e == time_i], rate = ctau_e[i, time_obs_e == time_i])}
for(time_trt_i in time_trt_lev){ e_rep_time_trt[[time_trt_i]][i, ] <- rgamma(n_time_trt_obs_e[time_trt_i], shape = cmu_e[i, time_trt_obs_e == time_trt_i] * ctau_e[i, time_trt_obs_e == time_trt_i], rate = ctau_e[i, time_trt_obs_e == time_trt_i])}
}
}
if(x$model_output$dist_e == "weib") {
ctau_e <- as.matrix(as.data.frame(x$model_output$model$BUGSoutput$sims.list$ctau_e))[, me_type == "observed"]
for(i in 1:n_sim){
for(time_i in 1:n_time){ e_rep_time[[time_i]][i, ] <- rweibull(n_time_obs_e[time_i], shape = ctau_e[i, time_obs_e == time_i], scale = cmu_e[i, time_obs_e == time_i] / exp(lgamma(1 + 1 / ctau_e[i, time_obs_e == time_i])))}
for(time_trt_i in time_trt_lev){ e_rep_time_trt[[time_trt_i]][i, ] <- rweibull(n_time_trt_obs_e[time_trt_i], shape = ctau_e[i, time_trt_obs_e == time_trt_i], scale = cmu_e[i, time_trt_obs_e == time_trt_i] / exp(lgamma(1 + 1 / ctau_e[i, time_trt_obs_e == time_trt_i])))}
}
}
if(x$model_output$dist_e == "beta") {
ctau_e <- as.matrix(as.data.frame(x$model_output$model$BUGSoutput$sims.list$ctau_e))[, me_type == "observed"]
for(i in 1:n_sim){
for(time_i in 1:n_time){ e_rep_time[[time_i]][i, ] <- rbeta(n_time_obs_e[time_i], shape1 = cmu_e[i, time_obs_e == time_i] * ctau_e[i, time_obs_e == time_i], shape2 = (1 - cmu_e[i, time_obs_e == time_i]) * ctau_e[i, time_obs_e == time_i])}
for(time_trt_i in time_trt_lev){ e_rep_time_trt[[time_trt_i]][i, ] <- rbeta(n_time_trt_obs_e[time_trt_i], shape1 = cmu_e[i, time_trt_obs_e == time_trt_i] * ctau_e[i, time_trt_obs_e == time_trt_i], shape2 = (1 - cmu_e[i, time_trt_obs_e == time_trt_i]) * ctau_e[i, time_trt_obs_e == time_trt_i])}
}
}
if(x$model_output$dist_c == "norm") {
for(i in 1:n_sim){
for(time_i in 1:n_time){ c_rep_time[[time_i]][i, ] <- rnorm(n_time_obs_c[time_i], mean = cmu_c[i, time_obs_c == time_i], sd = s_c[i, time_i])}
for(time_trt_i in time_trt_lev){ c_rep_time_trt[[time_trt_i]][i, ] <- rnorm(n_time_trt_obs_c[time_trt_i], mean = cmu_c[i, time_trt_obs_c == time_trt_i], sd = s_c[i, time_trt_index_time[time_trt_i]])}
}
}
if(x$model_output$dist_c == "gamma") {
ctau_c <- as.matrix(as.data.frame(x$model_output$model$BUGSoutput$sims.list$ctau_c))[, mc_type == "observed"]
for(i in 1:n_sim){
for(time_i in 1:n_time){ c_rep_time[[time_i]][i, ] <- rgamma(n_time_obs_c[time_i], shape = cmu_c[i, time_obs_c == time_i], rate = ctau_c[i, time_obs_c == time_i])}
for(time_trt_i in time_trt_lev){ c_rep_time_trt[[time_trt_i]][i, ] <- rgamma(n_time_trt_obs_c[time_trt_i], shape = cmu_c[i, time_trt_obs_c == time_trt_i] * ctau_c[i, time_trt_obs_c == time_trt_i], rate = ctau_c[i, time_trt_obs_c == time_trt_i])}
}
}
if(x$model_output$dist_c == "lnorm") {
for(i in 1:n_sim){
for(time_i in 1:n_time){ c_rep_time[[time_i]][i, ] <- rlnorm(n_time_obs_c[time_i], meanlog = clmu_c[i, time_obs_c == time_i], sdlog = 1 / sqrt(ltau_c[i, time_i]))}
for(time_trt_i in time_trt_lev){ c_rep_time_trt[[time_trt_i]][i, ] <- rlnorm(n_time_trt_obs_c[time_trt_i], meanlog = clmu_c[i, time_trt_obs_c == time_trt_i], sdlog = 1 / sqrt(ltau_c[i, time_trt_index_time[time_trt_i]]))}
}
}
}
}
if(length(trt) == 1) {
if(trt == "none") {
list.trt_name <- "none"}
if(trt == "all") {
list.trt_name <- trt_lev}
if(!trt %in% c("all", "none") & is.character(trt)) {
list.trt_name <- trt}
}
if(length(trt) > 1 & is.character(trt)) {
list.trt_name <- trt}
if(is.numeric(trt)) {
list.trt_name <- trt_lev[trt]
}
if(length(time.plot) == 1) {
if(time.plot == "all") {
list.time_name <- as.character(1:n_time)}
if(!time.plot %in% c("all") & is.character(time.plot)) {
list.time_name <- time.plot}
}
if(length(time.plot) > 1 & is.character(time.plot)) {
list.time_name <- time.plot}
if(is.numeric(time.plot)) {
list.time_name <- as.character(1:n_time)[time.plot]}
color_scheme_set(scheme = scheme_set)
gplot_e <- gplot_c <- gplot_e_trt <- gplot_c_trt <- list()
gplot_e_time <- gplot_c_time <- gplot_e_time_trt <- gplot_c_time_trt <- list()
if(type == "histogram") {
if(exists("bins", where = exArgs)) { bins = exArgs$bins} else { bins = 30}
if(exists("breaks", where = exArgs)) { breaks = exArgs$breaks} else { breaks = NULL}
if(exists("freq", where = exArgs)) { freq = exArgs$freq} else { freq = TRUE}
if(x$data_format == "wide") {
gplot_e <- ppc_hist(y = e_obs, yrep = e_rep[1 : ndisplay, ], bins = bins, breaks = breaks, freq = freq)
gplot_c <- ppc_hist(y = c_obs, yrep = c_rep[1 : ndisplay, ], bins = bins, breaks = breaks, freq = freq)
for(trt_i in trt_lev) {
gplot_e_trt[[trt_i]] <- ppc_hist(y = e_obs_trt[[trt_i]], yrep = e_rep_trt[[trt_i]][1 : ndisplay, ], bins = bins, breaks = breaks, freq = freq)
gplot_c_trt[[trt_i]] <- ppc_hist(y = c_obs_trt[[trt_i]], yrep = c_rep_trt[[trt_i]][1 : ndisplay, ], bins = bins, breaks = breaks, freq = freq)
}
}
if(x$data_format == "long") {
for(time_i in 1:n_time) {
gplot_e_time[[time_i]] <- ppc_hist(y = e_obs_time[[time_i]], yrep = e_rep_time[[time_i]][1 : ndisplay, ], bins = bins, breaks = breaks, freq = freq)
gplot_c_time[[time_i]] <- ppc_hist(y = c_obs_time[[time_i]], yrep = c_rep_time[[time_i]][1 : ndisplay, ], bins = bins, breaks = breaks, freq = freq)
}
for(time_trt_i in time_trt_lev) {
gplot_e_time_trt[[time_trt_i]] <- ppc_hist(y = e_obs_time_trt[[time_trt_i]], yrep = e_rep_time_trt[[time_trt_i]][1 : ndisplay, ], bins = bins, breaks = breaks, freq = freq)
gplot_c_time_trt[[time_trt_i]] <- ppc_hist(y = c_obs_time_trt[[time_trt_i]], yrep = c_rep_time_trt[[time_trt_i]][1 : ndisplay, ], bins = bins, breaks = breaks, freq = freq)
}
}
}
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}
if(x$data_format == "wide") {
gplot_e <- ppc_boxplot(y = e_obs, yrep = e_rep[1 : ndisplay, ], notch = notch, size = size, alpha = alpha)
gplot_c <- ppc_boxplot(y = c_obs, yrep = c_rep[1 : ndisplay, ], notch = notch, size = size, alpha = alpha)
for(trt_i in trt_lev) {
gplot_e_trt[[trt_i]] <- ppc_boxplot(y = e_obs_trt[[trt_i]], yrep = e_rep_trt[[trt_i]][1 : ndisplay, ], notch = notch, size = size, alpha = alpha)
gplot_c_trt[[trt_i]] <- ppc_boxplot(y = c_obs_trt[[trt_i]], yrep = c_rep_trt[[trt_i]][1 : ndisplay, ], notch = notch, size = size, alpha = alpha)
}
}
if(x$data_format == "long") {
for(time_i in 1:n_time) {
gplot_e_time[[time_i]] <- ppc_boxplot(y = e_obs_time[[time_i]], yrep = e_rep_time[[time_i]][1 : ndisplay, ], notch = notch, size = size, alpha = alpha)
gplot_c_time[[time_i]] <- ppc_boxplot(y = c_obs_time[[time_i]], yrep = c_rep_time[[time_i]][1 : ndisplay, ], notch = notch, size = size, alpha = alpha)
}
for(time_trt_i in time_trt_lev) {
gplot_e_time_trt[[time_trt_i]] <- ppc_boxplot(y = e_obs_time_trt[[time_trt_i]], yrep = e_rep_time_trt[[time_trt_i]][1 : ndisplay, ], notch = notch, size = size, alpha = alpha)
gplot_e_time_trt[[time_trt_i]] <- ppc_boxplot(y = c_obs_time_trt[[time_trt_i]], yrep = c_rep_time_trt[[time_trt_i]][1 : ndisplay, ], notch = notch, size = size, alpha = alpha)
}
}
}
if(type == "freqpoly") {
if(exists("bins", where = exArgs)) { bins = exArgs$bins} else { bins = 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}
if(x$data_format == "wide") {
gplot_e <- ppc_freqpoly(y = e_obs, yrep = e_rep[1 : ndisplay, ], bins = bins, freq = freq, size = size, alpha = alpha)
gplot_c <- ppc_freqpoly(y = c_obs, yrep = c_rep[1 : ndisplay, ], bins = bins, freq = freq, size = size, alpha = alpha)
for(trt_i in trt_lev) {
gplot_e_trt[[trt_i]] <- ppc_freqpoly(y = e_obs_trt[[trt_i]], yrep = e_rep_trt[[trt_i]][1 : ndisplay, ], bins = bins, freq = freq, size = size, alpha = alpha)
gplot_c_trt[[trt_i]] <- ppc_freqpoly(y = c_obs_trt[[trt_i]], yrep = c_rep_trt[[trt_i]][1 : ndisplay, ], bins = bins, freq = freq, size = size, alpha = alpha)
}
}
if(x$data_format == "long") {
for(time_i in 1:n_time) {
gplot_e_time[[time_i]] <- ppc_freqpoly(y = e_obs_time[[time_i]], yrep = e_rep_time[[time_i]][1 : ndisplay, ], bins = bins, freq = freq, size = size, alpha = alpha)
gplot_c_time[[time_i]] <- ppc_freqpoly(y = c_obs_time[[time_i]], yrep = c_rep_time[[time_i]][1 : ndisplay, ], bins = bins, freq = freq, size = size, alpha = alpha)
}
for(time_trt_i in time_trt_lev) {
gplot_e_time_trt[[time_trt_i]] <- ppc_freqpoly(y = e_obs_time_trt[[time_trt_i]], yrep = e_rep_time_trt[[time_trt_i]][1 : ndisplay, ], bins = bins, freq = freq, size = size, alpha = alpha)
gplot_e_time_trt[[time_trt_i]] <- ppc_freqpoly(y = c_obs_time_trt[[time_trt_i]], yrep = c_rep_time_trt[[time_trt_i]][1 : ndisplay, ], bins = bins, freq = freq, size = size, alpha = alpha)
}
}
}
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}
if(x$data_format == "wide") {
gplot_e <- ppc_dens(y = e_obs, yrep = e_rep[1 : ndisplay, ], size = size, alpha = alpha, trim = trim)
gplot_c <- ppc_dens(y = c_obs, yrep = c_rep[1 : ndisplay, ], size = size, alpha = alpha, trim = trim)
for(trt_i in trt_lev) {
gplot_e_trt[[trt_i]] <- ppc_dens(y = e_obs_trt[[trt_i]], yrep = e_rep_trt[[trt_i]][1 : ndisplay, ], size = size, alpha = alpha, trim = trim)
gplot_c_trt[[trt_i]] <- ppc_dens(y = c_obs_trt[[trt_i]], yrep = c_rep_trt[[trt_i]][1 : ndisplay, ], size = size, alpha = alpha, trim = trim)
}
}
if(x$data_format == "long") {
for(time_i in 1:n_time) {
gplot_e_time[[time_i]] <- ppc_dens(y = e_obs_time[[time_i]], yrep = e_rep_time[[time_i]][1 : ndisplay, ], size = size, alpha = alpha, trim = trim)
gplot_c_time[[time_i]] <- ppc_dens(y = c_obs_time[[time_i]], yrep = c_rep_time[[time_i]][1 : ndisplay, ], size = size, alpha = alpha, trim = trim)
}
for(time_trt_i in time_trt_lev) {
gplot_e_time_trt[[time_trt_i]] <- ppc_dens(y = e_obs_time_trt[[time_trt_i]], yrep = e_rep_time_trt[[time_trt_i]][1 : ndisplay, ], size = size, alpha = alpha, trim = trim)
gplot_c_time_trt[[time_trt_i]] <- ppc_dens(y = c_obs_time_trt[[time_trt_i]], yrep = c_rep_time_trt[[time_trt_i]][1 : ndisplay, ], size = size, alpha = alpha, trim = trim)
}
}
}
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}
if(x$data_format == "wide") {
gplot_e <- ppc_dens_overlay(y = e_obs, yrep = e_rep[1 : ndisplay, ], size = size, alpha = alpha, trim = trim, bw = bw, adjust = adjust, kernel = kernel, n_dens = n_dens)
gplot_c <- ppc_dens_overlay(y = c_obs, yrep = c_rep[1 : ndisplay, ], size = size, alpha = alpha, trim = trim, bw = bw, adjust = adjust, kernel = kernel, n_dens = n_dens)
for(trt_i in trt_lev) {
gplot_e_trt[[trt_i]] <- ppc_dens_overlay(y = e_obs_trt[[trt_i]], yrep = e_rep_trt[[trt_i]][1 : ndisplay, ], size = size, alpha = alpha, trim = trim, bw = bw, adjust = adjust, kernel = kernel, n_dens = n_dens)
gplot_c_trt[[trt_i]] <- ppc_dens_overlay(y = c_obs_trt[[trt_i]], yrep = c_rep_trt[[trt_i]][1 : ndisplay, ], size = size, alpha = alpha, trim = trim, bw = bw, adjust = adjust, kernel = kernel, n_dens = n_dens)
}
}
if(x$data_format == "long") {
for(time_i in 1:n_time) {
gplot_e_time[[time_i]] <- ppc_dens_overlay(y = e_obs_time[[time_i]], yrep = e_rep_time[[time_i]][1 : ndisplay, ], size = size, alpha = alpha, trim = trim, bw = bw, adjust = adjust, kernel = kernel, n_dens = n_dens)
gplot_c_time[[time_i]] <- ppc_dens_overlay(y = c_obs_time[[time_i]], yrep = c_rep_time[[time_i]][1 : ndisplay, ], size = size, alpha = alpha, trim = trim, bw = bw, adjust = adjust, kernel = kernel, n_dens = n_dens)
}
for(time_trt_i in time_trt_lev) {
gplot_e_time_trt[[time_trt_i]] <- ppc_dens_overlay(y = e_obs_time_trt[[time_trt_i]], yrep = e_rep_time_trt[[time_trt_i]][1 : ndisplay, ], size = size, alpha = alpha, trim = trim, bw = bw, adjust = adjust, kernel = kernel, n_dens = n_dens)
gplot_c_time_trt[[time_trt_i]] <- ppc_dens_overlay(y = c_obs_time_trt[[time_trt_i]], yrep = c_rep_time_trt[[time_trt_i]][1 : ndisplay, ], size = size, alpha = alpha, trim = trim, bw = bw, adjust = adjust, kernel = kernel, n_dens = n_dens)
}
}
}
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}
if(x$data_format == "wide") {
gplot_e <- ppc_ecdf_overlay(y = e_obs, yrep = e_rep[1 : ndisplay, ], size = size, alpha = alpha, discrete = discrete, pad = pad)
gplot_c <- ppc_ecdf_overlay(y = c_obs, yrep = c_rep[1 : ndisplay, ], size = size, alpha = alpha, discrete = discrete, pad = pad)
for(trt_i in trt_lev) {
gplot_e_trt[[trt_i]] <- ppc_ecdf_overlay(y = e_obs_trt[[trt_i]], yrep = e_rep_trt[[trt_i]][1 : ndisplay, ], size = size, alpha = alpha, discrete = discrete, pad = pad)
gplot_c_trt[[trt_i]] <- ppc_ecdf_overlay(y = c_obs_trt[[trt_i]], yrep = c_rep_trt[[trt_i]][1 : ndisplay, ], size = size, alpha = alpha, discrete = discrete, pad = pad)
}
}
if(x$data_format == "long") {
for(time_i in 1:n_time) {
gplot_e_time[[time_i]] <- ppc_ecdf_overlay(y = e_obs_time[[time_i]], yrep = e_rep_time[[time_i]][1 : ndisplay, ], size = size, alpha = alpha, discrete = discrete, pad = pad)
gplot_c_time[[time_i]] <- ppc_ecdf_overlay(y = c_obs_time[[time_i]], yrep = c_rep_time[[time_i]][1 : ndisplay, ], size = size, alpha = alpha, discrete = discrete, pad = pad)
}
for(time_trt_i in time_trt_lev) {
gplot_e_time_trt[[time_trt_i]] <- ppc_ecdf_overlay(y = e_obs_time_trt[[time_trt_i]], yrep = e_rep_time_trt[[time_trt_i]][1 : ndisplay, ], size = size, alpha = alpha, discrete = discrete, pad = pad)
gplot_c_time_trt[[time_trt_i]] <- ppc_ecdf_overlay(y = c_obs_time_trt[[time_trt_i]], yrep = c_rep_time_trt[[time_trt_i]][1 : ndisplay, ], size = size, alpha = alpha, discrete = discrete, pad = pad)
}
}
}
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")}
if(x$data_format == "wide") {
gplot_e <- ppc_stat_2d(y = e_obs, yrep = e_rep[1 : ndisplay, ], stat = stat, size = size, alpha = alpha)
gplot_c <- ppc_stat_2d(y = c_obs, yrep = c_rep[1 : ndisplay, ], stat = stat, size = size, alpha = alpha)
for(trt_i in trt_lev) {
gplot_e_trt[[trt_i]] <- ppc_stat_2d(y = e_obs_trt[[trt_i]], yrep = e_rep_trt[[trt_i]][1 : ndisplay, ], stat = stat, size = size, alpha = alpha)
gplot_c_trt[[trt_i]] <- ppc_stat_2d(y = c_obs_trt[[trt_i]], yrep = c_rep_trt[[trt_i]][1 : ndisplay, ], stat = stat, size = size, alpha = alpha)
}
}
if(x$data_format == "long") {
for(time_i in 1:n_time) {
gplot_e_time[[time_i]] <- ppc_stat_2d(y = e_obs_time[[time_i]], yrep = e_rep_time[[time_i]][1 : ndisplay, ], stat = stat, size = size, alpha = alpha)
gplot_c_time[[time_i]] <- ppc_stat_2d(y = c_obs_time[[time_i]], yrep = c_rep_time[[time_i]][1 : ndisplay, ], stat = stat, size = size, alpha = alpha)
}
for(time_trt_i in time_trt_lev) {
gplot_e_time_trt[[time_trt_i]] <- ppc_stat_2d(y = e_obs_time_trt[[time_trt_i]], yrep = e_rep_time_trt[[time_trt_i]][1 : ndisplay, ], stat = stat, size = size, alpha = alpha)
gplot_c_time_trt[[time_trt_i]] <- ppc_stat_2d(y = c_obs_time_trt[[time_trt_i]], yrep = c_rep_time_trt[[time_trt_i]][1 : ndisplay, ], stat = stat, size = size, alpha = alpha)
}
}
}
if(type == "error_hist") {
if(exists("bins", where = exArgs)) { bins = exArgs$bins} else { bins = 30}
if(exists("freq", where = exArgs)) { freq = exArgs$freq} else { freq = TRUE}
if(exists("breaks", where = exArgs)) { breaks = exArgs$breaks} else { breaks = NULL}
if(x$data_format == "wide") {
gplot_e <- ppc_error_hist(y = e_obs, yrep = e_rep[1 : ndisplay, ], bins = bins, breaks = breaks, freq = freq)
gplot_c <- ppc_error_hist(y = c_obs, yrep = c_rep[1 : ndisplay, ], bins = bins, breaks = breaks, freq = freq)
for(trt_i in trt_lev) {
gplot_e_trt[[trt_i]] <- ppc_error_hist(y = e_obs_trt[[trt_i]], yrep = e_rep_trt[[trt_i]][1 : ndisplay, ], bins = bins, breaks = breaks, freq = freq)
gplot_c_trt[[trt_i]] <- ppc_error_hist(y = c_obs_trt[[trt_i]], yrep = c_rep_trt[[trt_i]][1 : ndisplay, ], bins = bins, breaks = breaks, freq = freq)
}
}
if(x$data_format == "long") {
for(time_i in 1:n_time) {
gplot_e_time[[time_i]] <- ppc_error_hist(y = e_obs_time[[time_i]], yrep = e_rep_time[[time_i]][1 : ndisplay, ], bins = bins, breaks = breaks, freq = freq)
gplot_c_time[[time_i]] <- ppc_error_hist(y = c_obs_time[[time_i]], yrep = c_rep_time[[time_i]][1 : ndisplay, ], bins = bins, breaks = breaks, freq = freq)
}
for(time_trt_i in time_trt_lev) {
gplot_e_time_trt[[time_trt_i]] <- ppc_error_hist(y = e_obs_time_trt[[time_trt_i]], yrep = e_rep_time_trt[[time_trt_i]][1 : ndisplay, ], bins = bins, breaks = breaks, freq = freq)
gplot_c_time_trt[[time_trt_i]] <- ppc_error_hist(y = c_obs_time_trt[[time_trt_i]], yrep = c_rep_time_trt[[time_trt_i]][1 : ndisplay, ], bins = bins, breaks = breaks, freq = freq)
}
}
}
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}
if(x$data_format == "wide") {
gplot_e <- ppc_error_scatter(y = e_obs, yrep = e_rep[1 : ndisplay, ], size = size, alpha = alpha)
gplot_c <- ppc_error_scatter(y = c_obs, yrep = c_rep[1 : ndisplay, ], size = size, alpha = alpha)
for(trt_i in trt_lev) {
gplot_e_trt[[trt_i]] <- ppc_error_scatter(y = e_obs_trt[[trt_i]], yrep = e_rep_trt[[trt_i]][1 : ndisplay, ], size = size, alpha = alpha)
gplot_c_trt[[trt_i]] <- ppc_error_scatter(y = c_obs_trt[[trt_i]], yrep = c_rep_trt[[trt_i]][1 : ndisplay, ], size = size, alpha = alpha)
}
}
if(x$data_format == "long") {
for(time_i in 1:n_time) {
gplot_e_time[[time_i]] <- ppc_error_scatter(y = e_obs_time[[time_i]], yrep = e_rep_time[[time_i]][1 : ndisplay, ], size = size, alpha = alpha)
gplot_c_time[[time_i]] <- ppc_error_scatter(y = c_obs_time[[time_i]], yrep = c_rep_time[[time_i]][1 : ndisplay, ], size = size, alpha = alpha)
}
for(time_trt_i in time_trt_lev) {
gplot_e_time_trt[[time_trt_i]] <- ppc_error_scatter(y = e_obs_time_trt[[time_trt_i]], yrep = e_rep_time_trt[[time_trt_i]][1 : ndisplay, ], size = size, alpha = alpha)
gplot_c_time_trt[[time_trt_i]] <- ppc_error_scatter(y = c_obs_time_trt[[time_trt_i]], yrep = c_rep_time_trt[[time_trt_i]][1 : ndisplay, ], size = size, alpha = alpha)
}
}
}
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}
if(x$data_format == "wide") {
gplot_e <- ppc_error_scatter_avg(y = e_obs, yrep = e_rep[1 : ndisplay, ], size = size, alpha = alpha)
gplot_c <- ppc_error_scatter_avg(y = c_obs, yrep = c_rep[1 : ndisplay, ], size = size, alpha = alpha)
for(trt_i in trt_lev) {
gplot_e_trt[[trt_i]] <- ppc_error_scatter_avg(y = e_obs_trt[[trt_i]], yrep = e_rep_trt[[trt_i]][1 : ndisplay, ], size = size, alpha = alpha)
gplot_c_trt[[trt_i]] <- ppc_error_scatter_avg(y = c_obs_trt[[trt_i]], yrep = c_rep_trt[[trt_i]][1 : ndisplay, ], size = size, alpha = alpha)
}
}
if(x$data_format == "long") {
for(time_i in 1:n_time) {
gplot_e_time[[time_i]] <- ppc_error_scatter_avg(y = e_obs_time[[time_i]], yrep = e_rep_time[[time_i]][1 : ndisplay, ], size = size, alpha = alpha)
gplot_c_time[[time_i]] <- ppc_error_scatter_avg(y = c_obs_time[[time_i]], yrep = c_rep_time[[time_i]][1 : ndisplay, ], size = size, alpha = alpha)
}
for(time_trt_i in time_trt_lev) {
gplot_e_time_trt[[time_trt_i]] <- ppc_error_scatter_avg(y = e_obs_time_trt[[time_trt_i]], yrep = e_rep_time_trt[[time_trt_i]][1 : ndisplay, ], size = size, alpha = alpha)
gplot_c_time_trt[[time_trt_i]] <- ppc_error_scatter_avg(y = c_obs_time_trt[[time_trt_i]], yrep = c_rep_time_trt[[time_trt_i]][1 : ndisplay, ], size = size, alpha = alpha)
}
}
}
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}
if(x$data_format == "wide") {
gplot_e <- ppc_error_binned(y = e_obs, yrep = e_rep[1 : ndisplay, ], bins = bins, size = size, alpha = alpha)
gplot_c <- ppc_error_binned(y = c_obs, yrep = c_rep[1 : ndisplay, ], bins = bins, size = size, alpha = alpha)
for(trt_i in trt_lev) {
gplot_e_trt[[trt_i]] <- ppc_error_binned(y = e_obs_trt[[trt_i]], yrep = e_rep_trt[[trt_i]][1 : ndisplay, ], bins = bins, size = size, alpha = alpha)
gplot_c_trt[[trt_i]] <- ppc_error_binned(y = c_obs_trt[[trt_i]], yrep = c_rep_trt[[trt_i]][1 : ndisplay, ], bins = bins, size = size, alpha = alpha)
}
}
if(x$data_format == "long") {
for(time_i in 1:n_time) {
gplot_e_time[[time_i]] <- ppc_error_binned(y = e_obs_time[[time_i]], yrep = e_rep_time[[time_i]][1 : ndisplay, ], bins = bins, size = size, alpha = alpha)
gplot_c_time[[time_i]] <- ppc_error_binned(y = c_obs_time[[time_i]], yrep = c_rep_time[[time_i]][1 : ndisplay, ], bins = bins, size = size, alpha = alpha)
}
for(time_trt_i in time_trt_lev) {
gplot_e_time_trt[[time_trt_i]] <- ppc_error_binned(y = e_obs_time_trt[[time_trt_i]], yrep = e_rep_time_trt[[time_trt_i]][1 : ndisplay, ], bins = bins, size = size, alpha = alpha)
gplot_c_time_trt[[time_trt_i]] <- ppc_error_binned(y = c_obs_time_trt[[time_trt_i]], yrep = c_rep_time_trt[[time_trt_i]][1 : ndisplay, ], bins = bins, size = size, alpha = alpha)
}
}
}
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}
if(x$data_format == "wide") {
gplot_e <- ppc_intervals(y = e_obs, yrep = e_rep[1 : ndisplay, ], prob = prob, prob_outer = prob_outer, size = size, fatten = fatten)
gplot_c <- ppc_intervals(y = c_obs, yrep = c_rep[1 : ndisplay, ], prob = prob, prob_outer = prob_outer, size = size, fatten = fatten)
for(trt_i in trt_lev) {
gplot_e_trt[[trt_i]] <- ppc_intervals(y = e_obs_trt[[trt_i]], yrep = e_rep_trt[[trt_i]][1 : ndisplay, ], prob = prob, prob_outer = prob_outer, size = size, fatten = fatten)
gplot_c_trt[[trt_i]] <- ppc_intervals(y = c_obs_trt[[trt_i]], yrep = c_rep_trt[[trt_i]][1 : ndisplay, ], prob = prob, prob_outer = prob_outer, size = size, fatten = fatten)
}
}
if(x$data_format == "long") {
for(time_i in 1:n_time) {
gplot_e_time[[time_i]] <- ppc_intervals(y = e_obs_time[[time_i]], yrep = e_rep_time[[time_i]][1 : ndisplay, ], prob = prob, prob_outer = prob_outer, size = size, fatten = fatten)
gplot_c_time[[time_i]] <- ppc_intervals(y = c_obs_time[[time_i]], yrep = c_rep_time[[time_i]][1 : ndisplay, ], prob = prob, prob_outer = prob_outer, size = size, fatten = fatten)
}
for(time_trt_i in time_trt_lev) {
gplot_e_time_trt[[time_trt_i]] <- ppc_intervals(y = e_obs_time_trt[[time_trt_i]], yrep = e_rep_time_trt[[time_trt_i]][1 : ndisplay, ], prob = prob, prob_outer = prob_outer, size = size, fatten = fatten)
gplot_c_time_trt[[time_trt_i]] <- ppc_intervals(y = c_obs_time_trt[[time_trt_i]], yrep = c_rep_time_trt[[time_trt_i]][1 : ndisplay, ], prob = prob, prob_outer = prob_outer, size = size, fatten = fatten)
}
}
}
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}
if(x$data_format == "wide") {
gplot_e <- ppc_ribbon(y = e_obs, yrep = e_rep[1 : ndisplay, ], prob = prob, prob_outer = prob_outer, alpha = alpha, size = size)
gplot_c <- ppc_ribbon(y = c_obs, yrep = c_rep[1 : ndisplay, ], prob = prob, prob_outer = prob_outer, alpha = alpha, size = size)
for(trt_i in trt_lev) {
gplot_e_trt[[trt_i]] <- ppc_ribbon(y = e_obs_trt[[trt_i]], yrep = e_rep_trt[[trt_i]][1 : ndisplay, ], prob = prob, prob_outer = prob_outer, alpha = alpha, size = size)
gplot_c_trt[[trt_i]] <- ppc_ribbon(y = c_obs_trt[[trt_i]], yrep = c_rep_trt[[trt_i]][1 : ndisplay, ], prob = prob, prob_outer = prob_outer, alpha = alpha, size = size)
}
}
if(x$data_format == "long") {
for(time_i in 1:n_time) {
gplot_e_time[[time_i]] <- ppc_ribbon(y = e_obs_time[[time_i]], yrep = e_rep_time[[time_i]][1 : ndisplay, ], prob = prob, prob_outer = prob_outer, alpha = alpha, size = size)
gplot_c_time[[time_i]] <- ppc_ribbon(y = c_obs_time[[time_i]], yrep = c_rep_time[[time_i]][1 : ndisplay, ], prob = prob, prob_outer = prob_outer, alpha = alpha, size = size)
}
for(time_trt_i in time_trt_lev) {
gplot_e_time_trt[[time_trt_i]] <- ppc_ribbon(y = e_obs_time_trt[[time_trt_i]], yrep = e_rep_time_trt[[time_trt_i]][1 : ndisplay, ], prob = prob, prob_outer = prob_outer, alpha = alpha, size = size)
gplot_c_time_trt[[time_trt_i]] <- ppc_ribbon(y = c_obs_time_trt[[time_trt_i]], yrep = c_rep_time_trt[[time_trt_i]][1 : ndisplay, ], prob = prob, prob_outer = prob_outer, alpha = alpha, size = size)
}
}
}
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}
if(x$data_format == "wide") {
gplot_e <- ppc_scatter(y = e_obs, yrep = e_rep[1 : ndisplay, ], alpha = alpha, size = size)
gplot_c <- ppc_scatter(y = c_obs, yrep = c_rep[1 : ndisplay, ], alpha = alpha, size = size)
for(trt_i in trt_lev) {
gplot_e_trt[[trt_i]] <- ppc_scatter(y = e_obs_trt[[trt_i]], yrep = e_rep_trt[[trt_i]][1 : ndisplay, ], alpha = alpha, size = size)
gplot_c_trt[[trt_i]] <- ppc_scatter(y = c_obs_trt[[trt_i]], yrep = c_rep_trt[[trt_i]][1 : ndisplay, ], alpha = alpha, size = size)
}
}
if(x$data_format == "long") {
for(time_i in 1:n_time) {
gplot_e_time[[time_i]] <- ppc_scatter(y = e_obs_time[[time_i]], yrep = e_rep_time[[time_i]][1 : ndisplay, ], alpha = alpha, size = size)
gplot_c_time[[time_i]] <- ppc_scatter(y = c_obs_time[[time_i]], yrep = c_rep_time[[time_i]][1 : ndisplay, ], alpha = alpha, size = size)
}
for(time_trt_i in time_trt_lev) {
gplot_e_time_trt[[time_trt_i]] <- ppc_scatter(y = e_obs_time_trt[[time_trt_i]], yrep = e_rep_time_trt[[time_trt_i]][1 : ndisplay, ], alpha = alpha, size = size)
gplot_c_time_trt[[time_trt_i]] <- ppc_scatter(y = c_obs_time_trt[[time_trt_i]], yrep = c_rep_time_trt[[time_trt_i]][1 : ndisplay, ], alpha = alpha, size = size)
}
}
}
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}
if(x$data_format == "wide") {
gplot_e <- ppc_scatter_avg(y = e_obs, yrep = e_rep[1 : ndisplay, ], alpha = alpha, size = size)
gplot_c <- ppc_scatter_avg(y = c_obs, yrep = c_rep[1 : ndisplay, ], alpha = alpha, size = size)
for(trt_i in trt_lev) {
gplot_e_trt[[trt_i]] <- ppc_scatter_avg(y = e_obs_trt[[trt_i]], yrep = e_rep_trt[[trt_i]][1 : ndisplay, ], alpha = alpha, size = size)
gplot_c_trt[[trt_i]] <- ppc_scatter_avg(y = c_obs_trt[[trt_i]], yrep = c_rep_trt[[trt_i]][1 : ndisplay, ], alpha = alpha, size = size)
}
}
if(x$data_format == "long") {
for(time_i in 1:n_time) {
gplot_e_time[[time_i]] <- ppc_scatter_avg(y = e_obs_time[[time_i]], yrep = e_rep_time[[time_i]][1 : ndisplay, ], alpha = alpha, size = size)
gplot_c_time[[time_i]] <- ppc_scatter_avg(y = c_obs_time[[time_i]], yrep = c_rep_time[[time_i]][1 : ndisplay, ], alpha = alpha, size = size)
}
for(time_trt_i in time_trt_lev) {
gplot_e_time_trt[[time_trt_i]] <- ppc_scatter_avg(y = e_obs_time_trt[[time_trt_i]], yrep = e_rep_time_trt[[time_trt_i]][1 : ndisplay, ], alpha = alpha, size = size)
gplot_c_time_trt[[time_trt_i]] <- ppc_scatter_avg(y = c_obs_time_trt[[time_trt_i]], yrep = c_rep_time_trt[[time_trt_i]][1 : ndisplay, ], alpha = alpha, size = size)
}
}
}
if(type == "pit_ecdf") {
if(exists("pit", where = exArgs)) { pit = exArgs$pit} else { pit = NULL}
if(exists("K", where = exArgs)) { K = exArgs$K} else { K = NULL}
if(exists("prob", where = exArgs)) { prob = exArgs$prob} else { prob = 0.99}
if(exists("plot_diff", where = exArgs)) { plot_diff = exArgs$plot_diff} else { plot_diff = FALSE}
if(exists("interpolate_adj", where = exArgs)) { interpolate_adj = exArgs$interpolate_adj} else { interpolate_adj = NULL}
if(x$data_format == "wide") {
gplot_e <- ppc_pit_ecdf(y = e_obs, yrep = e_rep[1 : ndisplay, ], pit = pit, K = K, prob = prob, plot_diff = plot_diff, interpolate_adj = interpolate_adj)
gplot_c <- ppc_pit_ecdf(y = c_obs, yrep = c_rep[1 : ndisplay, ], pit = pit, K = K, prob = prob, plot_diff = plot_diff, interpolate_adj = interpolate_adj)
for(trt_i in trt_lev) {
gplot_e_trt[[trt_i]] <- ppc_pit_ecdf(y = e_obs_trt[[trt_i]], yrep = e_rep_trt[[trt_i]][1 : ndisplay, ], pit = pit, K = K, prob = prob, plot_diff = plot_diff, interpolate_adj = interpolate_adj)
gplot_c_trt[[trt_i]] <- ppc_pit_ecdf(y = c_obs_trt[[trt_i]], yrep = c_rep_trt[[trt_i]][1 : ndisplay, ], pit = pit, K = K, prob = prob, plot_diff = plot_diff, interpolate_adj = interpolate_adj)
}
}
if(x$data_format == "long") {
for(time_i in 1:n_time) {
gplot_e_time[[time_i]] <- ppc_pit_ecdf(y = e_obs_time[[time_i]], yrep = e_rep_time[[time_i]][1 : ndisplay, ], pit = pit, K = K, prob = prob, plot_diff = plot_diff, interpolate_adj = interpolate_adj)
gplot_c_time[[time_i]] <- ppc_pit_ecdf(y = c_obs_time[[time_i]], yrep = c_rep_time[[time_i]][1 : ndisplay, ], pit = pit, K = K, prob = prob, plot_diff = plot_diff, interpolate_adj = interpolate_adj)
}
for(time_trt_i in time_trt_lev) {
gplot_e_time_trt[[time_trt_i]] <- ppc_pit_ecdf(y = e_obs_time_trt[[time_trt_i]], yrep = e_rep_time_trt[[time_trt_i]][1 : ndisplay, ], pit = pit, K = K, prob = prob, plot_diff = plot_diff, interpolate_adj = interpolate_adj)
gplot_c_time_trt[[time_trt_i]] <- ppc_pit_ecdf(y = c_obs_time_trt[[time_trt_i]], yrep = c_rep_time_trt[[time_trt_i]][1 : ndisplay, ], pit = pit, K = K, prob = prob, plot_diff = plot_diff, interpolate_adj = interpolate_adj)
}
}
}
if(type == "loo_pit_overlay") {
if(exists("lw", where = exArgs)) { lw = exArgs$lw} else { lw = NULL}
if(exists("psis_object", where = exArgs)) { psis_object = exArgs$psis_object} else { psis_object = NULL}
if(exists("pit ", where = exArgs)) { pit = exArgs$pit } else { pit = NULL}
if(exists("samples", where = exArgs)) { samples = exArgs$samples} else { samples = 100}
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("boundary_correction", where = exArgs)) { boundary_correction = exArgs$boundary_correction} else { boundary_correction = TRUE}
if(exists("grid_len", where = exArgs)) { grid_len = exArgs$grid_len} else { grid_len = 512}
if(exists("bw", where = exArgs)) { bw = exArgs$bw} else { bw = "nrd0"}
if(exists("trim", where = exArgs)) { trim = exArgs$trim} else { trim = FALSE}
if(exists("adjust", where = exArgs)) { adjust = exArgs$adjust} else { adjust = 1}
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}
if(is.null(lw) & is.null(psis_object)) { stop("Please provide either a matrix of (smoothed) log weights or a 'psis_object' when computing LOO predictive checks.")}
if(!is.null(lw)) { lw = lw[1 : ndisplay, ]}
if(!is.null(psis_object)) { psis_object$log_weights = psis_object$log_weights[1 : ndisplay, ]}
if(x$data_format == "wide") {
if(outcome == "effects") {
gplot_e <- gplot_c <- ppc_loo_pit_overlay(y = e_obs, yrep = e_rep[1 : ndisplay, ], lw = lw, psis_object = psis_object, pit = pit, samples = samples, size = size, alpha = alpha, boundary_correction = boundary_correction,
grid_len = grid_len, bw = bw, trim = trim, adjust = adjust, kernel = kernel, n_dens = n_dens)
for(trt_i in trt_lev) {
if(!is.null(lw)) { lwt <- lw[, trt_obs_e == trt_i]} else {lwt <- lw}
if(!is.null(psis_object)) {
psis_objectt <- psis_object
psis_objectt$log_weights <- psis_objectt$log_weights[, trt_obs_e == trt_i]
psis_objectt$diagnostics$pareto_k <- psis_object$diagnostics$pareto_k[trt_obs_e == trt_i]
psis_objectt$diagnostics$n_eff <- psis_object$diagnostics$n_eff[trt_obs_e == trt_i]
psis_objectt$diagnostics$r_eff <- psis_object$diagnostics$r_eff[trt_obs_e == trt_i]} else { psis_objectt <- psis_object}
gplot_e_trt[[trt_i]] <- gplot_c_trt[[trt_i]] <- ppc_loo_pit_overlay(y = e_obs_trt[[trt_i]], yrep = e_rep_trt[[trt_i]][1 : ndisplay, ], lw = lwt, psis_object = psis_objectt, pit = pit, samples = samples, size = size, alpha = alpha, boundary_correction = boundary_correction,
grid_len = grid_len, bw = bw, trim = trim, adjust = adjust, kernel = kernel, n_dens = n_dens)
}
}
if(outcome == "costs") {
gplot_c <- gplot_e <- ppc_loo_pit_overlay(y = c_obs, yrep = c_rep[1 : ndisplay, ], lw = lw, psis_object = psis_object, pit = pit, samples = samples, size = size, alpha = alpha, boundary_correction = boundary_correction,
grid_len = grid_len, bw = bw, trim = trim, adjust = adjust, kernel = kernel, n_dens = n_dens)
for(trt_i in trt_lev) {
if(!is.null(lw)) { lwt <- lw[, trt_obs_c == trt_i]} else {lwt <- lw}
if(!is.null(psis_object)) {
psis_objectt <- psis_object
psis_objectt$log_weights <- psis_objectt$log_weights[, trt_obs_c == trt_i]
psis_objectt$diagnostics$pareto_k <- psis_object$diagnostics$pareto_k[trt_obs_c == trt_i]
psis_objectt$diagnostics$n_eff <- psis_object$diagnostics$n_eff[trt_obs_c == trt_i]
psis_objectt$diagnostics$r_eff <- psis_object$diagnostics$r_eff[trt_obs_c == trt_i]} else { psis_objectt <- psis_object}
gplot_c_trt[[trt_i]] <- gplot_e_trt[[trt_i]] <- ppc_loo_pit_overlay(y = c_obs_trt[[trt_i]], yrep = c_rep_trt[[trt_i]][1 : ndisplay, ], lw = lwt, psis_object = psis_objectt, pit = pit, samples = samples, size = size, alpha = alpha, boundary_correction = boundary_correction,
grid_len = grid_len, bw = bw, trim = trim, adjust = adjust, kernel = kernel, n_dens = n_dens)
}
}
if(outcome == "both") {
gplot_e <- ppc_loo_pit_overlay(y = e_obs, yrep = e_rep[1 : ndisplay, ], lw = lw, psis_object = psis_object, pit = pit, samples = samples, size = size, alpha = alpha, boundary_correction = boundary_correction,
grid_len = grid_len, bw = bw, trim = trim, adjust = adjust, kernel = kernel, n_dens = n_dens)
gplot_c <- ppc_loo_pit_overlay(y = c_obs, yrep = c_rep[1 : ndisplay, ], lw = lw, psis_object = psis_object, pit = pit, samples = samples, size = size, alpha = alpha, boundary_correction = boundary_correction,
grid_len = grid_len, bw = bw, trim = trim, adjust = adjust, kernel = kernel, n_dens = n_dens)
for(trt_i in trt_lev) {
if(!is.null(lw)) { lwte <- lw[, trt_obs_e == trt_i]; lwtc <- lw[, trt_obs_c == trt_i]} else {lwte <- lw; lwtc <- lw}
if(!is.null(psis_object)) {
psis_objectt <- psis_object
psis_objectt$log_weights <- psis_objectt$log_weights[, trt_obs_e == trt_i & trt_obs_c == trt_i]
psis_objectt$diagnostics$pareto_k <- psis_object$diagnostics$pareto_k[trt_obs_e == trt_i & trt_obs_c == trt_i]
psis_objectt$diagnostics$n_eff <- psis_object$diagnostics$n_eff[trt_obs_e == trt_i & trt_obs_c == trt_i]
psis_objectt$diagnostics$r_eff <- psis_object$diagnostics$r_eff[trt_obs_e == trt_i & trt_obs_c == trt_i]} else { psis_objectt <- psis_object}
gplot_e_trt[[trt_i]] <- ppc_loo_pit_overlay(y = e_obs_trt[[trt_i]], yrep = e_rep_trt[[trt_i]][1 : ndisplay, ], lw = lwte, psis_object = psis_objectt, pit = pit, samples = samples, size = size, alpha = alpha, boundary_correction = boundary_correction,
grid_len = grid_len, bw = bw, trim = trim, adjust = adjust, kernel = kernel, n_dens = n_dens)
gplot_c_trt[[trt_i]] <- ppc_loo_pit_overlay(y = c_obs_trt[[trt_i]], yrep = c_rep_trt[[trt_i]][1 : ndisplay, ], lw = lwtc, psis_object = psis_objectt, pit = pit, samples = samples, size = size, alpha = alpha, boundary_correction = boundary_correction,
grid_len = grid_len, bw = bw, trim = trim, adjust = adjust, kernel = kernel, n_dens = n_dens)
}
}
}
if(x$data_format == "long") {
if(outcome == "effects") {
for(time_i in 1:n_time) {
if(!is.null(lw)) { lwt <- lw[, time_obs_e == time_i]} else {lwt <- lw}
if(!is.null(psis_object)) {
psis_objectt <- psis_object
psis_objectt$log_weights <- psis_objectt$log_weights[, time_obs_e == time_i]
psis_objectt$diagnostics$pareto_k <- psis_object$diagnostics$pareto_k[time_obs_e == time_i]
psis_objectt$diagnostics$n_eff <- psis_object$diagnostics$n_eff[time_obs_e == time_i]
psis_objectt$diagnostics$r_eff <- psis_object$diagnostics$r_eff[time_obs_e == time_i]} else { psis_objectt <- psis_object}
gplot_e_time[[time_i]] <- gplot_c_time[[time_i]] <- ppc_loo_pit_overlay(y = e_obs_time[[time_i]], yrep = e_rep_time[[time_i]][1 : ndisplay, ], lw = lwt, psis_object = psis_objectt, pit = pit, samples = samples, size = size, alpha = alpha, boundary_correction = boundary_correction,
grid_len = grid_len, bw = bw, trim = trim, adjust = adjust, kernel = kernel, n_dens = n_dens)
}
for(time_trt_i in time_trt_lev) {
if(!is.null(lw)) { lwt <- lw[, time_trt_obs_e == time_trt_i]} else {lwt <- lw}
if(!is.null(psis_object)) {
psis_objectt <- psis_object
psis_objectt$log_weights <- psis_objectt$log_weights[, time_trt_obs_e == time_trt_i]
psis_objectt$diagnostics$pareto_k <- psis_object$diagnostics$pareto_k[time_trt_obs_e == time_trt_i]
psis_objectt$diagnostics$n_eff <- psis_object$diagnostics$n_eff[time_trt_obs_e == time_trt_i]
psis_objectt$diagnostics$r_eff <- psis_object$diagnostics$r_eff[time_trt_obs_e == time_trt_i]} else { psis_objectt <- psis_object}
gplot_e_time_trt[[time_trt_i]] <- gplot_c_time_trt[[time_trt_i]] <- ppc_loo_pit_overlay(y = e_obs_time_trt[[time_trt_i]], yrep = e_rep_time_trt[[time_trt_i]][1 : ndisplay, ], lw = lwt, psis_object = psis_objectt, pit = pit, samples = samples, size = size, alpha = alpha, boundary_correction = boundary_correction,
grid_len = grid_len, bw = bw, trim = trim, adjust = adjust, kernel = kernel, n_dens = n_dens)
}
}
if(outcome == "costs") {
for(time_i in 1:n_time) {
if(!is.null(lw)) { lwt <- lw[, time_obs_c == time_i]} else {lwt <- lw}
if(!is.null(psis_object)) {
psis_objectt <- psis_object
psis_objectt$log_weights <- psis_objectt$log_weights[, time_obs_c == time_i]
psis_objectt$diagnostics$pareto_k <- psis_object$diagnostics$pareto_k[time_obs_c == time_i]
psis_objectt$diagnostics$n_eff <- psis_object$diagnostics$n_eff[time_obs_c == time_i]
psis_objectt$diagnostics$r_eff <- psis_object$diagnostics$r_eff[time_obs_c == time_i]} else { psis_objectt <- psis_object}
gplot_c_time[[time_i]] <- gplot_e_time[[time_i]] <- ppc_loo_pit_overlay(y = c_obs_time[[time_i]], yrep = c_rep_time[[time_i]][1 : ndisplay, ], lw = lwt, psis_object = psis_objectt, pit = pit, samples = samples, size = size, alpha = alpha, boundary_correction = boundary_correction,
grid_len = grid_len, bw = bw, trim = trim, adjust = adjust, kernel = kernel, n_dens = n_dens)
}
for(time_trt_i in time_trt_lev) {
if(!is.null(lw)) { lwt <- lw[, time_trt_obs_c == time_trt_i]} else {lwt <- lw}
if(!is.null(psis_object)) {
psis_objectt <- psis_object
psis_objectt$log_weights <- psis_objectt$log_weights[, time_trt_obs_c == time_trt_i]
psis_objectt$diagnostics$pareto_k <- psis_object$diagnostics$pareto_k[time_trt_obs_c == time_trt_i]
psis_objectt$diagnostics$n_eff <- psis_object$diagnostics$n_eff[time_trt_obs_c == time_trt_i]
psis_objectt$diagnostics$r_eff <- psis_object$diagnostics$r_eff[time_trt_obs_c == time_trt_i]} else { psis_objectt <- psis_object}
gplot_c_time_trt[[time_trt_i]] <- gplot_e_time_trt[[time_trt_i]] <- ppc_loo_pit_overlay(y = c_obs_time_trt[[time_trt_i]], yrep = c_rep_time_trt[[time_trt_i]][1 : ndisplay, ], lw = lwt, psis_object = psis_objectt, pit = pit, samples = samples, size = size, alpha = alpha, boundary_correction = boundary_correction,
grid_len = grid_len, bw = bw, trim = trim, adjust = adjust, kernel = kernel, n_dens = n_dens)
}
}
if(outcome == "both") {
for(time_i in 1:n_time) {
if(!is.null(lw)) { lwte <- lw[, time_obs_e == time_i]; lwtc <- lw[, time_obs_c == time_i]} else {lwte <- lw; lwtc <- lw}
if(!is.null(psis_object)) {
psis_objectt <- psis_object
psis_objectt$log_weights <- psis_objectt$log_weights[, time_obs_e == time_i & time_obs_c == time_i]
psis_objectt$diagnostics$pareto_k <- psis_object$diagnostics$pareto_k[time_obs_e == time_i & time_obs_c == time_i]
psis_objectt$diagnostics$n_eff <- psis_object$diagnostics$n_eff[time_obs_e == time_i & time_obs_c == time_i]
psis_objectt$diagnostics$r_eff <- psis_object$diagnostics$r_eff[time_obs_e == time_i & time_obs_c == time_i]} else { psis_objectt <- psis_object}
gplot_e_time[[time_i]] <- ppc_loo_pit_overlay(y = e_obs_time[[time_i]], yrep = e_rep_time[[time_i]][1 : ndisplay, ], lw = lwte, psis_object = psis_objectt, pit = pit, samples = samples, size = size, alpha = alpha, boundary_correction = boundary_correction,
grid_len = grid_len, bw = bw, trim = trim, adjust = adjust, kernel = kernel, n_dens = n_dens)
gplot_c_time[[time_i]] <- ppc_loo_pit_overlay(y = c_obs_time[[time_i]], yrep = c_rep_time[[time_i]][1 : ndisplay, ], lw = lwtc, psis_object = psis_objectt, pit = pit, samples = samples, size = size, alpha = alpha, boundary_correction = boundary_correction,
grid_len = grid_len, bw = bw, trim = trim, adjust = adjust, kernel = kernel, n_dens = n_dens)
}
for(time_trt_i in time_trt_lev) {
if(!is.null(lw)) { lwte <- lw[, time_trt_obs_e == time_trt_i]; lwtc <- lw[, time_trt_obs_c == time_trt_i]} else {lwte <- lw; lwtc <- lw}
if(!is.null(psis_object)) {
psis_objectt <- psis_object
psis_objectt$log_weights <- psis_objectt$log_weights[, time_trt_obs_e == time_trt_i & time_trt_obs_c == time_trt_i]
psis_objectt$diagnostics$pareto_k <- psis_object$diagnostics$pareto_k[time_trt_obs_e == time_trt_i & time_trt_obs_c == time_trt_i]
psis_objectt$diagnostics$n_eff <- psis_object$diagnostics$n_eff[time_trt_obs_e == time_trt_i & time_trt_obs_c == time_trt_i]
psis_objectt$diagnostics$r_eff <- psis_object$diagnostics$r_eff[time_trt_obs_e == time_trt_i & time_trt_obs_c == time_trt_i]} else { psis_objectt <- psis_object}
gplot_e_time_trt[[time_trt_i]] <- ppc_loo_pit_overlay(y = e_obs_time_trt[[time_trt_i]], yrep = e_rep_time_trt[[time_trt_i]][1 : ndisplay, ], lw = lwte, psis_object = psis_objectt, pit = pit, samples = samples, size = size, alpha = alpha, boundary_correction = boundary_correction,
grid_len = grid_len, bw = bw, trim = trim, adjust = adjust, kernel = kernel, n_dens = n_dens)
gplot_c_time_trt[[time_trt_i]] <- ppc_loo_pit_overlay(y = c_obs_time_trt[[time_trt_i]], yrep = c_rep_time_trt[[time_trt_i]][1 : ndisplay, ], lw = lwtc, psis_object = psis_objectt, pit = pit, samples = samples, size = size, alpha = alpha, boundary_correction = boundary_correction,
grid_len = grid_len, bw = bw, trim = trim, adjust = adjust, kernel = kernel, n_dens = n_dens)
}
}
}
}
if(type == "loo_pit_ecdf") {
if(exists("lw", where = exArgs)) { lw = exArgs$lw} else { lw = NULL}
if(exists("psis_object", where = exArgs)) { psis_object = exArgs$psis_object} else { psis_object = NULL}
if(exists("pit", where = exArgs)) { pit = exArgs$pit} else { pit = NULL}
if(exists("K", where = exArgs)) { K = exArgs$K} else { K = NULL}
if(exists("prob", where = exArgs)) { prob = exArgs$prob} else { prob = 0.99}
if(exists("plot_diff", where = exArgs)) { plot_diff = exArgs$plot_diff} else { plot_diff = FALSE}
if(exists("interpolate_adj", where = exArgs)) { interpolate_adj = exArgs$interpolate_adj} else { interpolate_adj = NULL}
if(is.null(lw) & is.null(psis_object)) { stop("Please provide either a matrix of (smoothed) log weights or a 'psis_object' when computing LOO predictive checks.")}
if(!is.null(lw)) { lw = lw[1 : ndisplay, ]}
if(!is.null(psis_object)) { psis_object$log_weights = psis_object$log_weights[1 : ndisplay, ]}
if(x$data_format == "wide") {
if(outcome == "effects") {
gplot_e <- gplot_c <- ppc_loo_pit_ecdf(y = e_obs, yrep = e_rep[1 : ndisplay, ], lw = lw, psis_object = psis_object, pit = pit, K = K, prob = prob, plot_diff = plot_diff, interpolate_adj = interpolate_adj)
for(trt_i in trt_lev) {
if(!is.null(lw)) { lwt <- lw[, trt_obs_e == trt_i]} else {lwt <- lw}
if(!is.null(psis_object)) {
psis_objectt <- psis_object
psis_objectt$log_weights <- psis_objectt$log_weights[, trt_obs_e == trt_i]
psis_objectt$diagnostics$pareto_k <- psis_object$diagnostics$pareto_k[trt_obs_e == trt_i]
psis_objectt$diagnostics$n_eff <- psis_object$diagnostics$n_eff[trt_obs_e == trt_i]
psis_objectt$diagnostics$r_eff <- psis_object$diagnostics$r_eff[trt_obs_e == trt_i]} else { psis_objectt <- psis_object}
gplot_e_trt[[trt_i]] <- gplot_c_trt[[trt_i]] <- ppc_loo_pit_ecdf(y = e_obs_trt[[trt_i]], yrep = e_rep_trt[[trt_i]][1 : ndisplay, ], lw = lwt, psis_object = psis_objectt, pit = pit, K = K, prob = prob, plot_diff = plot_diff, interpolate_adj = interpolate_adj)
}
}
if(outcome == "costs") {
gplot_c <- gplot_e <- ppc_loo_pit_ecdf(y = c_obs, yrep = c_rep[1 : ndisplay, ], lw = lw, psis_object = psis_object, pit = pit, K = K, prob = prob, plot_diff = plot_diff, interpolate_adj = interpolate_adj)
for(trt_i in trt_lev) {
if(!is.null(lw)) { lwt <- lw[, trt_obs_c == trt_i]} else {lwt <- lw}
if(!is.null(psis_object)) {
psis_objectt <- psis_object
psis_objectt$log_weights <- psis_objectt$log_weights[, trt_obs_c == trt_i]
psis_objectt$diagnostics$pareto_k <- psis_object$diagnostics$pareto_k[trt_obs_c == trt_i]
psis_objectt$diagnostics$n_eff <- psis_object$diagnostics$n_eff[trt_obs_c == trt_i]
psis_objectt$diagnostics$r_eff <- psis_object$diagnostics$r_eff[trt_obs_c == trt_i]} else { psis_objectt <- psis_object}
gplot_c_trt[[trt_i]] <- gplot_e_trt[[trt_i]] <- ppc_loo_pit_ecdf(y = c_obs_trt[[trt_i]], yrep = c_rep_trt[[trt_i]][1 : ndisplay, ], lw = lwt, psis_object = psis_objectt, pit = pit, K = K, prob = prob, plot_diff = plot_diff, interpolate_adj = interpolate_adj)
}
}
if(outcome == "both") {
gplot_e <- ppc_loo_pit_ecdf(y = e_obs, yrep = e_rep[1 : ndisplay, ], lw = lw, psis_object = psis_object, pit = pit, K = K, prob = prob, plot_diff = plot_diff, interpolate_adj = interpolate_adj)
gplot_c <- ppc_loo_pit_ecdf(y = c_obs, yrep = c_rep[1 : ndisplay, ], lw = lw, psis_object = psis_object, pit = pit, K = K, prob = prob, plot_diff = plot_diff, interpolate_adj = interpolate_adj)
for(trt_i in trt_lev) {
if(!is.null(lw)) { lwte <- lw[, trt_obs_e == trt_i]; lwtc <- lw[, trt_obs_c == trt_i]} else {lwte <- lw; lwtc <- lw}
if(!is.null(psis_object)) {
psis_objectt <- psis_object
psis_objectt$log_weights <- psis_objectt$log_weights[, trt_obs_e == trt_i & trt_obs_c == trt_i]
psis_objectt$diagnostics$pareto_k <- psis_object$diagnostics$pareto_k[trt_obs_e == trt_i & trt_obs_c == trt_i]
psis_objectt$diagnostics$n_eff <- psis_object$diagnostics$n_eff[trt_obs_e == trt_i & trt_obs_c == trt_i]
psis_objectt$diagnostics$r_eff <- psis_object$diagnostics$r_eff[trt_obs_e == trt_i & trt_obs_c == trt_i]} else { psis_objectt <- psis_object}
gplot_e_trt[[trt_i]] <- ppc_loo_pit_ecdf(y = e_obs_trt[[trt_i]], yrep = e_rep_trt[[trt_i]][1 : ndisplay, ], lw = lwte, psis_object = psis_objectt, pit = pit, K = K, prob = prob, plot_diff = plot_diff, interpolate_adj = interpolate_adj)
gplot_c_trt[[trt_i]] <- ppc_loo_pit_ecdf(y = c_obs_trt[[trt_i]], yrep = c_rep_trt[[trt_i]][1 : ndisplay, ], lw = lwtc, psis_object = psis_objectt, pit = pit, K = K, prob = prob, plot_diff = plot_diff, interpolate_adj = interpolate_adj)
}
}
}
if(x$data_format == "long") {
if(outcome == "effects") {
for(time_i in 1:n_time) {
if(!is.null(lw)) { lwt <- lw[, time_obs_e == time_i]} else {lwt <- lw}
if(!is.null(psis_object)) {
psis_objectt <- psis_object
psis_objectt$log_weights <- psis_objectt$log_weights[, time_obs_e == time_i]
psis_objectt$diagnostics$pareto_k <- psis_object$diagnostics$pareto_k[time_obs_e == time_i]
psis_objectt$diagnostics$n_eff <- psis_object$diagnostics$n_eff[time_obs_e == time_i]
psis_objectt$diagnostics$r_eff <- psis_object$diagnostics$r_eff[time_obs_e == time_i]} else { psis_objectt <- psis_object}
gplot_e_time[[time_i]] <- gplot_c_time[[time_i]] <- ppc_loo_pit_ecdf(y = e_obs_time[[time_i]], yrep = e_rep_time[[time_i]][1 : ndisplay, ], lw = lwt, psis_object = psis_objectt, pit = pit, K = K, prob = prob, plot_diff = plot_diff, interpolate_adj = interpolate_adj)
}
for(time_trt_i in time_trt_lev) {
if(!is.null(lw)) { lwt <- lw[, time_trt_obs_e == time_trt_i]} else {lwt <- lw}
if(!is.null(psis_object)) {
psis_objectt <- psis_object
psis_objectt$log_weights <- psis_objectt$log_weights[, time_trt_obs_e == time_trt_i]
psis_objectt$diagnostics$pareto_k <- psis_object$diagnostics$pareto_k[time_trt_obs_e == time_trt_i]
psis_objectt$diagnostics$n_eff <- psis_object$diagnostics$n_eff[time_trt_obs_e == time_trt_i]
psis_objectt$diagnostics$r_eff <- psis_object$diagnostics$r_eff[time_trt_obs_e == time_trt_i]} else { psis_objectt <- psis_object}
gplot_e_time_trt[[time_trt_i]] <- gplot_c_time_trt[[time_trt_i]] <- ppc_loo_pit_ecdf(y = e_obs_time_trt[[time_trt_i]], yrep = e_rep_time_trt[[time_trt_i]][1 : ndisplay, ], lw = lwt, psis_object = psis_objectt, pit = pit, K = K, prob = prob, plot_diff = plot_diff, interpolate_adj = interpolate_adj)
}
}
if(outcome == "costs") {
for(time_i in 1:n_time) {
if(!is.null(lw)) { lwt <- lw[, time_obs_c == time_i]} else {lwt <- lw}
if(!is.null(psis_object)) {
psis_objectt <- psis_object
psis_objectt$log_weights <- psis_objectt$log_weights[, time_obs_c == time_i]
psis_objectt$diagnostics$pareto_k <- psis_object$diagnostics$pareto_k[time_obs_c == time_i]
psis_objectt$diagnostics$n_eff <- psis_object$diagnostics$n_eff[time_obs_c == time_i]
psis_objectt$diagnostics$r_eff <- psis_object$diagnostics$r_eff[time_obs_c == time_i]} else { psis_objectt <- psis_object}
gplot_c_time[[time_i]] <- gplot_e_time[[time_i]] <- ppc_loo_pit_ecdf(y = c_obs_time[[time_i]], yrep = c_rep_time[[time_i]][1 : ndisplay, ], lw = lwt, psis_object = psis_objectt, pit = pit, K = K, prob = prob, plot_diff = plot_diff, interpolate_adj = interpolate_adj)
}
for(time_trt_i in time_trt_lev) {
if(!is.null(lw)) { lwt <- lw[, time_trt_obs_c == time_trt_i]} else {lwt <- lw}
if(!is.null(psis_object)) {
psis_objectt <- psis_object
psis_objectt$log_weights <- psis_objectt$log_weights[, time_trt_obs_c == time_trt_i]
psis_objectt$diagnostics$pareto_k <- psis_object$diagnostics$pareto_k[time_trt_obs_c == time_trt_i]
psis_objectt$diagnostics$n_eff <- psis_object$diagnostics$n_eff[time_trt_obs_c == time_trt_i]
psis_objectt$diagnostics$r_eff <- psis_object$diagnostics$r_eff[time_trt_obs_c == time_trt_i]} else { psis_objectt <- psis_object}
gplot_c_time_trt[[time_trt_i]] <- gplot_e_time_trt[[time_trt_i]] <- ppc_loo_pit_ecdf(y = c_obs_time_trt[[time_trt_i]], yrep = c_rep_time_trt[[time_trt_i]][1 : ndisplay, ], lw = lwt, psis_object = psis_objectt, pit = pit, K = K, prob = prob, plot_diff = plot_diff, interpolate_adj = interpolate_adj)
}
}
if(outcome == "both") {
for(time_i in 1:n_time) {
if(!is.null(lw)) { lwte <- lw[, time_obs_e == time_i]; lwtc <- lw[, time_obs_c == time_i]} else {lwte <- lw; lwtc <- lw}
if(!is.null(psis_object)) {
psis_objectt <- psis_object
psis_objectt$log_weights <- psis_objectt$log_weights[, time_obs_e == time_i & trt_obs_c == time_i]
psis_objectt$diagnostics$pareto_k <- psis_object$diagnostics$pareto_k[time_obs_e == time_i & trt_obs_c == time_i]
psis_objectt$diagnostics$n_eff <- psis_object$diagnostics$n_eff[time_obs_e == time_i & trt_obs_c == time_i]
psis_objectt$diagnostics$r_eff <- psis_object$diagnostics$r_eff[time_obs_e == time_i & trt_obs_c == time_i]} else { psis_objectt <- psis_object}
gplot_e_time[[time_i]] <- ppc_loo_pit_ecdf(y = e_obs_time[[time_i]], yrep = e_rep_time[[time_i]][1 : ndisplay, ], lw = lwte, psis_object = psis_objectt, pit = pit, K = K, prob = prob, plot_diff = plot_diff, interpolate_adj = interpolate_adj)
gplot_c_time[[time_i]] <- ppc_loo_pit_ecdf(y = c_obs_time[[time_i]], yrep = c_rep_time[[time_i]][1 : ndisplay, ], lw = lwtc, psis_object = psis_objectt, pit = pit, K = K, prob = prob, plot_diff = plot_diff, interpolate_adj = interpolate_adj)
}
for(time_trt_i in time_trt_lev) {
if(!is.null(lw)) { lwte <- lw[, time_trt_obs_e == time_trt_i]; lwtc <- lw[, time_trt_obs_c == time_trt_i]} else {lwte <- lw; lwtc <- lw}
if(!is.null(psis_object)) {
psis_objectt <- psis_object
psis_objectt$log_weights <- psis_objectt$log_weights[, time_trt_obs_e == time_trt_i & time_trt_obs_c == time_trt_i]
psis_objectt$diagnostics$pareto_k <- psis_object$diagnostics$pareto_k[time_trt_obs_e == time_trt_i & time_trt_obs_c == time_trt_i]
psis_objectt$diagnostics$n_eff <- psis_object$diagnostics$n_eff[time_trt_obs_e == time_trt_i & time_trt_obs_c == time_trt_i]
psis_objectt$diagnostics$r_eff <- psis_object$diagnostics$r_eff[time_trt_obs_e == time_trt_i & time_trt_obs_c == time_trt_i]} else { psis_objectt <- psis_object}
gplot_e_time_trt[[time_trt_i]] <- ppc_loo_pit_ecdf(y = e_obs_time_trt[[time_trt_i]], yrep = e_rep_time_trt[[time_trt_i]][1 : ndisplay, ], lw = lwte, psis_object = psis_objectt, pit = pit, K = K, prob = prob, plot_diff = plot_diff, interpolate_adj = interpolate_adj)
gplot_c_time_trt[[time_trt_i]] <- ppc_loo_pit_ecdf(y = c_obs_time_trt[[time_trt_i]], yrep = c_rep_time_trt[[time_trt_i]][1 : ndisplay, ], lw = lwtc, psis_object = psis_objectt, pit = pit, K = K, prob = prob, plot_diff = plot_diff, interpolate_adj = interpolate_adj)
}
}
}
}
if(type == "stat" & is.null(corr)) {
if(!bpp){
if(exists("bins", where = exArgs)) { bins = exArgs$bins} else { bins = 30}
if(exists("freq", where = exArgs)) { freq = exArgs$freq} else { freq = TRUE}
if(exists("breaks", where = exArgs)) { breaks = exArgs$breaks} else { breaks = NULL}
if(x$data_format == "wide") {
gplot_e <- ppc_stat(y = e_obs, yrep = e_rep[1 : ndisplay, ], stat = stat, bins = bins, breaks = breaks, freq = freq)
gplot_c <- ppc_stat(y = c_obs, yrep = c_rep[1 : ndisplay, ], stat = stat, bins = bins, breaks = breaks, freq = freq)
for(trt_i in trt_lev) {
gplot_e_trt[[trt_i]] <- ppc_stat(y = e_obs_trt[[trt_i]], yrep = e_rep_trt[[trt_i]][1 : ndisplay, ], stat = stat, bins = bins, breaks = breaks, freq = freq)
gplot_c_trt[[trt_i]] <- ppc_stat(y = c_obs_trt[[trt_i]], yrep = c_rep_trt[[trt_i]][1 : ndisplay, ], stat = stat, bins = bins, breaks = breaks, freq = freq)
}
}
if(x$data_format == "long") {
for(time_i in 1:n_time) {
gplot_e_time[[time_i]] <- ppc_stat(y = e_obs_time[[time_i]], yrep = e_rep_time[[time_i]][1 : ndisplay, ], stat = stat, bins = bins, breaks = breaks, freq = freq)
gplot_c_time[[time_i]] <- ppc_stat(y = c_obs_time[[time_i]], yrep = c_rep_time[[time_i]][1 : ndisplay, ], stat = stat, bins = bins, breaks = breaks, freq = freq)
}
for(time_trt_i in time_trt_lev) {
gplot_e_time_trt[[time_trt_i]] <- ppc_stat(y = e_obs_time_trt[[time_trt_i]], yrep = e_rep_time_trt[[time_trt_i]][1 : ndisplay, ], stat = stat, bins = bins, breaks = breaks, freq = freq)
gplot_c_time_trt[[time_trt_i]] <- ppc_stat(y = c_obs_time_trt[[time_trt_i]], yrep = c_rep_time_trt[[time_trt_i]][1 : ndisplay, ], stat = stat, bins = bins, breaks = breaks, freq = freq)
}
}
}
if(bpp){
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("size.annotate", where = exArgs)) { size.annotate = exArgs$size.annotate} else { size.annotate = 3.5}
if(exists("alpha", where = exArgs)) { alpha = exArgs$alpha} else { alpha = 1}
if(exists("vjust", where = exArgs)) { vjust = exArgs$vjust } else { vjust = 1}
if(x$data_format == "wide") {
stat_e <- stat(e_obs); stat_c <- stat(c_obs)
stat_trt_e <- lapply(e_obs_trt, stat); stat_trt_c <- lapply(c_obs_trt, stat)
values_e <- apply(e_rep[1 : ndisplay, ], 1, stat)
values_c <- apply(c_rep[1 : ndisplay, ], 1, stat)
pval_e <- sum(values_e > stat_e) / length(values_e)
pval_c <- sum(values_c > stat_c) / length(values_c)
paste_p_e <- sprintf("Prob (T(y[rep]) > T(y)) = %0.3f", pval_e)
paste_p_c <- sprintf("Prob (T(y[rep]) > T(y)) = %0.3f", pval_c)
values_trt_e <- values_trt_c <- pval_trt_e <- pval_trt_c <- paste_p_trt_e <- paste_p_trt_c <- list()
for(trt_i in trt_lev) {
values_trt_e[[trt_i]] <- apply(e_rep_trt[[trt_i]][1 : ndisplay, ], 1, stat)
values_trt_c[[trt_i]] <- apply(c_rep_trt[[trt_i]][1 : ndisplay, ], 1, stat)
pval_trt_e[[trt_i]] <- sum(values_trt_e[[trt_i]] > stat_trt_e[[trt_i]]) / length(values_trt_e[[trt_i]])
pval_trt_c[[trt_i]] <- sum(values_trt_c[[trt_i]] > stat_trt_c[[trt_i]]) / length(values_trt_c[[trt_i]])
paste_p_trt_e[[trt_i]] <- sprintf("Prob (T(y[rep]) > T(y)) = %0.3f", pval_trt_e[[trt_i]])
paste_p_trt_c[[trt_i]] <- sprintf("Prob (T(y[rep]) > T(y)) = %0.3f", pval_trt_c[[trt_i]])
}
df_stat0 <- data.frame(values = c(stat_e, stat_c),
outcome = factor(c("effects", "costs"), levels = c("effects", "costs")),
text = c(paste_p_e, paste_p_c))
df_trt_stat0 <- data.frame(values = c(unlist(stat_trt_e), unlist(stat_trt_c)),
outcome = factor(rep(c("effects", "costs"), each = n_trt), levels = c("effects", "costs")),
text = c(unlist(paste_p_trt_e), unlist(paste_p_trt_c)),
trt = factor(rep(trt_lev, 2), levels = trt_lev))
df_stat <- data.frame(values = c(values_e, values_c),
outcome = factor(c(rep("effects", ndisplay), rep("costs", ndisplay)),
levels = c("effects", "costs")))
df_stat_trt <- data.frame(values = c(unlist(values_trt_e), unlist(values_trt_c)),
outcome = factor(c(rep("effects", ndisplay * n_trt), rep("costs", ndisplay * n_trt)),
levels = c("effects", "costs")),
trt = factor(c(rep(trt_lev, each = ndisplay), rep(trt_lev, each = ndisplay)),
levels = trt_lev))
if(outcome == "both") {
df_stat <- df_stat; df_stat_trt <- df_stat_trt
df_stat0 <- df_stat0; df_trt_stat0 <- df_trt_stat0}
if(outcome == "effects") {
df_stat <- df_stat[df_stat$outcome == "effects", ]
df_stat_trt <- df_stat_trt[df_stat_trt$outcome == "effects", ]
df_stat0 <- df_stat0[df_stat0$outcome == "effects"]
df_trt_stat0 <- df_trt_stat0[df_trt_stat0$outcome == "effects", ]}
if(outcome == "costs") {
df_stat <- df_stat[df_stat$outcome == "costs", ]
df_stat_trt <- df_stat_trt[df_stat_trt$outcome == "costs", ]
df_stat0 <- df_stat0[df_stat0$outcome == "costs"]
df_trt_stat0 <- df_trt_stat0[df_trt_stat0$outcome == "costs", ]}
if(length(trt) == 1) {
if(trt == "none") {
df_stat_gplot <- df_stat
df_stat0_gplot <- df_stat0
df.trt_name <- "none"}
if(trt == "all") {
df_stat_gplot <- df_stat_trt
df_stat0_gplot <- df_trt_stat0
df.trt_name <- trt_lev}
if(!trt %in% c("all", "none") & is.character(trt)) {
df_stat_gplot <- df_stat_trt[df_stat_trt$trt %in% trt, ]
df_stat0_gplot <- df_trt_stat0[df_trt_stat0$trt %in% trt, ]
df.trt_name <- trt
df_stat_gplot$trt <- factor(df_stat_gplot$trt, levels = df.trt_name)
df_stat0_gplot$trt <- factor(df_stat0_gplot$trt, levels = df.trt_name)}
}
if(length(trt) > 1 & is.character(trt)) {
df_stat_gplot <- df_stat_trt[df_stat_trt$trt %in% trt, ]
df_stat0_gplot <- df_trt_stat0[df_trt_stat0$trt %in% trt, ]
df.trt_name <- trt
df_stat_gplot$trt <- factor(df_stat_gplot$trt, levels = df.trt_name)
df_stat0_gplot$trt <- factor(df_stat0_gplot$trt, levels = df.trt_name)}
if(is.numeric(trt)) {
df.trt_name <- trt_lev[trt]
df_stat_gplot <- df_stat_trt[df_stat_trt$trt %in% df.trt_name, ]
df_stat0_gplot <- df_trt_stat0[df_trt_stat0$trt %in% df.trt_name, ]
df_stat_gplot$trt <- factor(df_stat_gplot$trt, levels = df.trt_name)
df_stat0_gplot$trt <- factor(df_stat0_gplot$trt, levels = df.trt_name)}
}
if(x$data_format == "long") {
stat_time_e <- lapply(e_obs_time, stat); stat_time_c <- lapply(c_obs_time, stat)
stat_time_trt_e <- lapply(e_obs_time_trt, stat); stat_time_trt_c <- lapply(c_obs_time_trt, stat)
values_time_e <- values_time_c <- pval_time_e <- pval_time_c <- paste_p_time_e <- paste_p_time_c <- list()
values_time_trt_e <- values_time_trt_c <- pval_time_trt_e <- pval_time_trt_c <- paste_p_time_trt_e <- paste_p_time_trt_c <- list()
for(time_i in 1:n_time) {
values_time_e[[time_i]] <- apply(e_rep_time[[time_i]][1 : ndisplay, ], 1, stat)
values_time_c[[time_i]] <- apply(c_rep_time[[time_i]][1 : ndisplay, ], 1, stat)
pval_time_e[[time_i]] <- sum(values_time_e[[time_i]] > stat_time_e[[time_i]]) / length(values_time_e[[time_i]])
pval_time_c[[time_i]] <- sum(values_time_c[[time_i]] > stat_time_c[[time_i]]) / length(values_time_c[[time_i]])
paste_p_time_e[[time_i]] <- sprintf("Prob (T(y[rep]) > T(y)) = %0.3f", pval_time_e[[time_i]])
paste_p_time_c[[time_i]] <- sprintf("Prob (T(y[rep]) > T(y)) = %0.3f", pval_time_c[[time_i]])
}
for(time_trt_i in time_trt_lev) {
values_time_trt_e[[time_trt_i]] <- apply(e_rep_time_trt[[time_trt_i]][1 : ndisplay, ], 1, stat)
values_time_trt_c[[time_trt_i]] <- apply(c_rep_time_trt[[time_trt_i]][1 : ndisplay, ], 1, stat)
pval_time_trt_e[[time_trt_i]] <- sum(values_time_trt_e[[time_trt_i]] > stat_time_trt_e[[time_trt_i]]) / length(values_time_trt_e[[time_trt_i]])
pval_time_trt_c[[time_trt_i]] <- sum(values_time_trt_c[[time_trt_i]] > stat_time_trt_c[[time_trt_i]]) / length(values_time_trt_c[[time_trt_i]])
paste_p_time_trt_e[[time_trt_i]] <- sprintf("Prob (T(y[rep]) > T(y)) = %0.3f", pval_time_trt_e[[time_trt_i]])
paste_p_time_trt_c[[time_trt_i]] <- sprintf("Prob (T(y[rep]) > T(y)) = %0.3f", pval_time_trt_c[[time_trt_i]])
}
df_time_stat0 <- data.frame(values = c(unlist(stat_time_e), unlist(stat_time_c)),
outcome = factor(rep(c("effects", "costs"), each = n_time), levels = c("effects", "costs")),
text = c(unlist(paste_p_time_e), unlist(paste_p_time_c)),
time = factor(rep(as.character(1:n_time), 2), levels = as.character(1:n_time)))
df_time_trt_stat0 <- data.frame(values = c(unlist(stat_time_trt_e), unlist(stat_time_trt_c)),
outcome = factor(rep(c("effects", "costs"), each = length(time_trt_lev)), levels = c("effects", "costs")),
text = c(unlist(paste_p_time_trt_e), unlist(paste_p_time_trt_c)),
time_trt = factor(rep(time_trt_lev, 2), levels = time_trt_lev))
df_stat_time <- data.frame(values = c(unlist(values_time_e), unlist(values_time_c)),
outcome = factor(c(rep("effects", ndisplay * n_time), rep("costs", ndisplay * n_time)),
levels = c("effects", "costs")),
time = factor(c(rep(as.character(1:n_time), each = ndisplay), rep(as.character(1:n_time), each = ndisplay)),
levels = as.character(1:n_time)))
df_stat_time_trt <- data.frame(values = c(unlist(values_time_trt_e), unlist(values_time_trt_c)),
outcome = factor(c(rep("effects", ndisplay * length(time_trt_lev)), rep("costs", ndisplay * length(time_trt_lev))),
levels = c("effects", "costs")),
time_trt = factor(c(rep(time_trt_lev, each = ndisplay), rep(time_trt_lev, each = ndisplay)),
levels = time_trt_lev))
if(outcome == "both") {
df_stat_time <- df_stat_time; df_stat_time_trt <- df_stat_time_trt
df_time_stat0 <- df_time_stat0; df_time_trt_stat0 <- df_time_trt_stat0}
if(outcome == "effects") {
df_stat_time <- df_stat_time[df_stat_time$outcome == "effects", ]
df_stat_time_trt <- df_stat_time_trt[df_stat_time_trt$outcome == "effects", ]
df_time_stat0 <- df_time_stat0[df_time_stat0$outcome == "effects", ]
df_time_trt_stat0 <- df_time_trt_stat0[df_time_trt_stat0$outcome == "effects", ]}
if(outcome == "costs") {
df_stat_time <- df_stat_time[df_stat_time$outcome == "costs", ]
df_stat_time_trt <- df_stat_time_trt[df_stat_time_trt$outcome == "costs", ]
df_time_stat0 <- df_time_stat0[df_time_stat0$outcome == "costs", ]
df_time_trt_stat0 <- df_time_trt_stat0[df_time_trt_stat0$outcome == "costs", ]}
if(all(length(list.trt_name) == 1 & length(list.time_name) == 1)) {
if(list.time_name %in% as.character(1:n_time)) {
if(list.trt_name == "none"){
df_stat_gplot <- df_stat_time[df_stat_time$time %in% as.character(list.time_name), ]
df_stat0_gplot <- df_time_stat0[df_time_stat0$time %in% as.character(list.time_name), ]
df_stat_gplot$time <- factor(df_stat_gplot$time, levels = list.time_name)
df_stat0_gplot$time <- factor(df_stat0_gplot$time, levels = list.time_name)
}
if(list.trt_name %in% trt_lev){
df_stat_gplot <- df_stat_time_trt[df_stat_time_trt$time_trt %in% c(sapply(as.character(list.time_name), paste, as.character(list.trt_name), sep = ".")), ]
df_stat0_gplot <- df_time_trt_stat0[df_time_trt_stat0$time_trt %in% c(sapply(as.character(list.time_name), paste, as.character(list.trt_name), sep = ".")), ]
df_stat_gplot$time_trt <- factor(df_stat_gplot$time_trt, levels = c(sapply(as.character(list.time_name), paste, as.character(list.trt_name), sep = ".")))
df_stat0_gplot$time_trt <- factor(df_stat0_gplot$time_trt, levels = c(sapply(as.character(list.time_name), paste, as.character(list.trt_name), sep = ".")))
}
}
}
if(all(length(list.trt_name) == 1 & length(list.time_name) > 1)) {
if(all(list.time_name %in% as.character(1:n_time))){
if(list.trt_name == "none"){
df_stat_gplot <- df_stat_time[df_stat_time$time %in% as.character(list.time_name), ]
df_stat0_gplot <- df_time_stat0[df_time_stat0$time %in% as.character(list.time_name), ]
df_stat_gplot$time <- factor(df_stat_gplot$time, levels = list.time_name)
df_stat0_gplot$time <- factor(df_stat0_gplot$time, levels = list.time_name)
}
if(list.trt_name %in% trt_lev){
df_stat_gplot <- df_stat_time_trt[df_stat_time_trt$time_trt %in% c(sapply(as.character(list.time_name), paste, as.character(list.trt_name), sep = ".")), ]
df_stat0_gplot <- df_time_trt_stat0[df_time_trt_stat0$time_trt %in% c(sapply(as.character(list.time_name), paste, as.character(list.trt_name), sep = ".")), ]
df_stat_gplot$time_trt <- factor(df_stat_gplot$time_trt, levels = c(sapply(as.character(list.time_name), paste, as.character(list.trt_name), sep = ".")))
df_stat0_gplot$time_trt <- factor(df_stat0_gplot$time_trt, levels = c(sapply(as.character(list.time_name), paste, as.character(list.trt_name), sep = ".")))
}
}
if(!any(list.time_name %in% as.character(1:n_time))){
stop("Please provide only valid time.plot values.")
}
}
if(all(length(list.trt_name) > 1 & length(list.time_name) == 1)) {
if(list.time_name %in% as.character(1:n_time)){
if(all(list.trt_name %in% trt_lev)){
df_stat_gplot <- df_stat_time_trt[df_stat_time_trt$time_trt %in% c(sapply(as.character(list.time_name), paste, as.character(list.trt_name), sep = ".")), ]
df_stat0_gplot <- df_time_trt_stat0[df_time_trt_stat0$time_trt %in% c(sapply(as.character(list.time_name), paste, as.character(list.trt_name), sep = ".")), ]
df_stat_gplot$time_trt <- factor(df_stat_gplot$time_trt, levels = c(sapply(as.character(list.time_name), paste, as.character(list.trt_name), sep = ".")))
df_stat0_gplot$time_trt <- factor(df_stat0_gplot$time_trt, levels = c(sapply(as.character(list.time_name), paste, as.character(list.trt_name), sep = ".")))
}
if(!any(as.character(list.trt_name) %in% trt_lev)){
stop("Please provide only valid treatment names or values.")
}
}
}
if(all(length(list.trt_name) > 1 & length(list.time_name) > 1)) {
if(all(list.time_name %in% as.character(1:n_time))){
if(all(list.trt_name %in% trt_lev)){
df_stat_gplot <- df_stat_time_trt[df_stat_time_trt$time_trt %in% c(sapply(as.character(list.time_name), paste, as.character(list.trt_name), sep = ".")), ]
df_stat0_gplot <- df_time_trt_stat0[df_time_trt_stat0$time_trt %in% c(sapply(as.character(list.time_name), paste, as.character(list.trt_name), sep = ".")), ]
df_stat_gplot$time_trt <- factor(df_stat_gplot$time_trt, levels = c(sapply(as.character(list.time_name), paste, as.character(list.trt_name), sep = ".")))
df_stat0_gplot$time_trt <- factor(df_stat0_gplot$time_trt, levels = c(sapply(as.character(list.time_name), paste, as.character(list.trt_name), sep = ".")))
}
if(!any(list.trt_name %in% trt_lev)){
stop("Please provide only valid treatment names or values.")
}
}
if(!any(list.time_name %in% as.character(1:n_time))){
stop("Please provide only valid time.plot values.")
}
}
}
color_sheme_rep <- color_scheme_get()[1]
color_sheme_rep_border <- color_scheme_get()[2]
color_sheme_obs <- color_scheme_get()[5]
color_sheme_obs_border <- color_scheme_get()[6]
set_theme <- bayesplot_theme_get()
ggplot2::theme_set(set_theme)
gplot_ppc <- ggplot2::ggplot(df_stat_gplot, ggplot2::aes(x = values)) +
ggplot2::geom_density(fill = color_sheme_rep$light, color = color_sheme_rep_border$light_highlight,
trim = trim, linewidth = size, alpha = alpha) +
ggplot2::geom_vline(data = df_stat0_gplot, ggplot2::aes(xintercept = values), color = color_sheme_obs$dark,
linewidth = size * 1.5, linetype = "dashed") +
ggplot2::geom_text(data = df_stat0_gplot, mapping = ggplot2::aes(x = -Inf, y = Inf, label = text),
hjust = hjust, vjust = vjust, size = size.annotate)
if(x$data_format == "wide") {
if(length(df.trt_name) == 1) {
if(df.trt_name == "none"){
gplot_ppc <- gplot_ppc + ggplot2::facet_wrap( ~ outcome, scales = "free")}
if(df.trt_name %in% trt_lev){
gplot_ppc <- gplot_ppc + ggplot2::facet_wrap( ~ outcome, scales = "free")}
}
if(length(df.trt_name) > 1) {
if(all(df.trt_name %in% trt_lev)){
gplot_ppc <- gplot_ppc + ggplot2::facet_wrap(outcome ~ trt, scales = "free")
}
if(!any(df.trt_name %in% trt_lev)){
stop("Please provide only valid treatment names or values.")}
}
}
if(x$data_format == "long") {
if(all(length(list.trt_name) == 1 & length(list.time_name) == 1)) {
if(all(list.trt_name == "none" & list.time_name %in% as.character(1:n_time))){
gplot_ppc <- gplot_ppc + ggplot2::facet_wrap( ~ outcome, scales = "free")}
if(all(list.trt_name %in% trt_lev & list.time_name %in% as.character(1:n_time))){
gplot_ppc <- gplot_ppc + ggplot2::facet_wrap(outcome ~ time_trt, scales = "free")}
}
if(all(length(list.trt_name) == 1 & length(list.time_name) > 1)) {
if(all(list.trt_name == "none" & list.time_name %in% as.character(1:n_time))){
gplot_ppc <- gplot_ppc + ggplot2::facet_wrap(outcome ~ time, scales = "free")}
if(all(list.trt_name %in% trt_lev & all(list.time_name %in% as.character(1:n_time)))){
gplot_ppc <- gplot_ppc + ggplot2::facet_wrap(outcome ~ time_trt, scales = "free")}
}
if(all(length(list.trt_name) > 1 & length(list.time_name) == 1)) {
if(all(all(list.trt_name %in% trt_lev) & list.time_name %in% as.character(1:n_time))){
gplot_ppc <- gplot_ppc + ggplot2::facet_wrap(outcome ~ time_trt, scales = "free")}
}
if(all(length(list.trt_name) > 1 & length(list.time_name) > 1)) {
if(all(all(list.trt_name %in% trt_lev) & all(list.time_name %in% as.character(1:n_time)))){
gplot_ppc <- gplot_ppc + ggplot2::facet_wrap(outcome ~ time_trt, scales = "free")}
}
}
}
}
if(type == "stat" & !is.null(corr)) {
if(x$data_format == "wide") {
df_ec <- data.frame(effects = x$data_set$data_raw$e, costs = x$data_set$data_raw$c,
trt = x$data_set$data_raw$data$trt)
df_ec_cc <- df_ec[complete.cases(df_ec), ]
stat_biv <- corr(df_ec_cc$effects, df_ec_cc$costs)
array_ec <- array(NA, dim = c(ndisplay, nrow(df_ec_cc), 2),
dimnames = list(NULL, NULL, c("effects", "costs")))
array_ec[, , "effects"] <- e_rep[1 : ndisplay, 1:nrow(df_ec_cc)]
array_ec[, , "costs"] <- c_rep[1 : ndisplay, 1:nrow(df_ec_cc)]
df.values_biv <- data.frame(values = apply(array_ec, 1, corr)[2, ])
p_biv <- sum(df.values_biv$values > stat_biv) / length(df.values_biv$values)
paste_p_biv <- sprintf("Prob (T(y[rep]) > T(y)) = %0.3f", p_biv)
df_trt_ec <- df_trt_ec_cc <- stat_trt_biv <- array_trt_ec <- df.values_trt_biv <- p_trt_biv <- paste_p_trt_biv <- list()
for(trt_i in trt_lev) {
df_trt_ec[[trt_i]] <- df_ec[df_ec$trt == trt_i, ]
df_trt_ec_cc[[trt_i]] <- df_trt_ec[[trt_i]][complete.cases(df_trt_ec[[trt_i]]), ]
stat_trt_biv[[trt_i]] <- corr(df_trt_ec_cc[[trt_i]]$effects, df_trt_ec_cc[[trt_i]]$costs)
array_trt_ec[[trt_i]] <- array(NA, dim = c(ndisplay, nrow(df_trt_ec_cc[[trt_i]]), 2),
dimnames = list(NULL, NULL, c("effects", "costs")))
array_trt_ec[[trt_i]][, , "effects"] <- e_rep_trt[[trt_i]][1 : ndisplay, 1:nrow(df_trt_ec_cc[[trt_i]])]
array_trt_ec[[trt_i]][, , "costs"] <- c_rep_trt[[trt_i]][1 : ndisplay, 1:nrow(df_trt_ec_cc[[trt_i]])]
df.values_trt_biv[[trt_i]] <- data.frame(values = apply(array_trt_ec[[trt_i]], 1, corr)[2, ])
p_trt_biv[[trt_i]] <- sum(df.values_trt_biv[[trt_i]]$values > stat_trt_biv[[trt_i]]) / length(df.values_trt_biv[[trt_i]]$values)
paste_p_trt_biv[[trt_i]] <- sprintf("Prob (T(y[rep]) > T(y)) = %0.3f", p_trt_biv[[trt_i]])
}
df_biv0 <- data.frame(values = stat_biv,
text = paste_p_biv)
df_trt_biv0 <- data.frame(values = unlist(stat_trt_biv),
text = c(unlist(paste_p_trt_biv)),
trt = factor(trt_lev, levels = trt_lev),
row.names = NULL)
df_biv <- df.values_biv
df_biv_trt <- data.frame(values = unlist(df.values_trt_biv),
trt = factor(rep(trt_lev, each = ndisplay), levels = trt_lev),
row.names = NULL)
if(length(trt) == 1) {
if(trt == "none") {
df_biv_gplot <- df_biv
df_biv0_gplot <- df_biv0
df.trt_name <- "none"}
if(trt == "all") {
df_biv_gplot <- df_biv_trt
df_biv0_gplot <- df_trt_biv0
df.trt_name <- trt_lev}
if(!trt %in% c("all", "none") & is.character(trt)) {
df_biv_gplot <- df_biv_trt[df_biv_trt$trt %in% trt, ]
df_biv0_gplot <- df_trt_biv0[df_trt_biv0$trt %in% trt, ]
df.trt_name <- trt
df_biv_gplot$trt <- factor(df_biv_gplot$trt, levels = df.trt_name)
df_biv0_gplot$trt <- factor(df_biv0_gplot$trt, levels = df.trt_name)}
}
if(length(trt) > 1 & is.character(trt)) {
df_biv_gplot <- df_biv_trt[df_biv_trt$trt %in% trt, ]
df_biv0_gplot <- df_trt_biv0[df_trt_biv0$trt %in% trt, ]
df.trt_name <- trt
df_biv_gplot$trt <- factor(df_biv_gplot$trt, levels = df.trt_name)
df_biv0_gplot$trt <- factor(df_biv0_gplot$trt, levels = df.trt_name)}
if(is.numeric(trt)) {
df.trt_name <- trt_lev[trt]
df_biv_gplot <- df_biv_trt[df_biv_trt$trt %in% df.trt_name, ]
df_biv0_gplot <- df_trt_biv0[df_trt_biv0$trt %in% df.trt_name, ]
df_biv_gplot$trt <- factor(df_biv_gplot$trt, levels = df.trt_name)
df_biv0_gplot$trt <- factor(df_biv0_gplot$trt, levels = df.trt_name)}
}
if(x$data_format == "long") {
df_ec <- data.frame(effects = x$data_set$data_raw$data_long$e, costs = x$data_set$data_raw$data_long$c,
time_trt = x$data_set$data_raw$data_long$time_trt, time = x$data_set$data_raw$data_long$time)
df_time_ec <- df_time_ec_cc <- stat_time_biv <- array_time_ec <- df.values_time_biv <- p_time_biv <- paste_p_time_biv <- list()
for(time_i in as.character(1:n_time)) {
df_time_ec[[time_i]] <- df_ec[df_ec$time == time_i, ]
df_time_ec_cc[[time_i]] <- df_time_ec[[time_i]][complete.cases(df_time_ec[[time_i]]), ]
stat_time_biv[[time_i]] <- corr(df_time_ec_cc[[time_i]]$effects, df_time_ec_cc[[time_i]]$costs)
array_time_ec[[time_i]] <- array(NA, dim = c(ndisplay, nrow(df_time_ec_cc[[time_i]]), 2),
dimnames = list(NULL, NULL, c("effects", "costs")))
array_time_ec[[time_i]][, , "effects"] <- e_rep_time[[time_i]][1 : ndisplay, 1:nrow(df_time_ec_cc[[time_i]])]
array_time_ec[[time_i]][, , "costs"] <- c_rep_time[[time_i]][1 : ndisplay, 1:nrow(df_time_ec_cc[[time_i]])]
df.values_time_biv[[time_i]] <- data.frame(values = apply(array_time_ec[[time_i]], 1, corr)[2, ])
p_time_biv[[time_i]] <- sum(df.values_time_biv[[time_i]]$values > stat_time_biv[[time_i]]) / length(df.values_time_biv[[time_i]]$values)
paste_p_time_biv[[time_i]] <- sprintf("Prob (T(y[rep]) > T(y)) = %0.3f", p_time_biv[[time_i]])
}
df_time_trt_ec <- df_time_trt_ec_cc <- stat_time_trt_biv <- array_time_trt_ec <- df.values_time_trt_biv <- p_time_trt_biv <- paste_p_time_trt_biv <- list()
for(time_trt_i in time_trt_lev) {
df_time_trt_ec[[time_trt_i]] <- df_ec[df_ec$time_trt == time_trt_i, ]
df_time_trt_ec_cc[[time_trt_i]] <- df_time_trt_ec[[time_trt_i]][complete.cases(df_time_trt_ec[[time_trt_i]]), ]
stat_time_trt_biv[[time_trt_i]] <- corr(df_time_trt_ec_cc[[time_trt_i]]$effects, df_time_trt_ec_cc[[time_trt_i]]$costs)
array_time_trt_ec[[time_trt_i]] <- array(NA, dim = c(ndisplay, nrow(df_time_trt_ec_cc[[time_trt_i]]), 2),
dimnames = list(NULL, NULL, c("effects", "costs")))
array_time_trt_ec[[time_trt_i]][, , "effects"] <- e_rep_time_trt[[time_trt_i]][1 : ndisplay, 1:nrow(df_time_trt_ec_cc[[time_trt_i]])]
array_time_trt_ec[[time_trt_i]][, , "costs"] <- c_rep_time_trt[[time_trt_i]][1 : ndisplay, 1:nrow(df_time_trt_ec_cc[[time_trt_i]])]
df.values_time_trt_biv[[time_trt_i]] <- data.frame(values = apply(array_time_trt_ec[[time_trt_i]], 1, corr)[2, ])
p_time_trt_biv[[time_trt_i]] <- sum(df.values_time_trt_biv[[time_trt_i]]$values > stat_time_trt_biv[[time_trt_i]]) / length(df.values_time_trt_biv[[time_trt_i]]$values)
paste_p_time_trt_biv[[time_trt_i]] <- sprintf("Prob (T(y[rep]) > T(y)) = %0.3f", p_time_trt_biv[[time_trt_i]])
}
df_time_biv0 <- data.frame(values = unlist(stat_time_biv),
text = c(unlist(paste_p_time_biv)),
time = factor(as.character(1:n_time), levels = as.character(1:n_time)),
row.names = NULL)
df_biv_time <- data.frame(values = unlist(df.values_time_biv),
time = factor(rep(as.character(1:n_time), each = ndisplay), levels = as.character(1:n_time)),
row.names = NULL)
df_time_trt_biv0 <- data.frame(values = unlist(stat_time_trt_biv),
text = c(unlist(paste_p_time_trt_biv)),
time_trt = factor(time_trt_lev, levels = time_trt_lev),
row.names = NULL)
df_biv_time_trt <- data.frame(values = unlist(df.values_time_trt_biv),
time_trt = factor(rep(time_trt_lev, each = ndisplay), levels = time_trt_lev),
row.names = NULL)
if(all(length(list.trt_name) == 1 & length(list.time_name) == 1)) {
if(list.time_name %in% as.character(1:n_time)) {
if(list.trt_name == "none"){
df_biv_gplot <- df_biv_time[df_biv_time$time %in% as.character(list.time_name), ]
df_biv0_gplot <- df_time_biv0[df_time_biv0$time %in% as.character(list.time_name), ]
df_biv_gplot$time <- factor(df_biv_gplot$time, levels = list.time_name)
df_biv0_gplot$time <- factor(df_biv0_gplot$time, levels = list.time_name)
}
if(list.trt_name %in% trt_lev){
df_biv_gplot <- df_biv_time_trt[df_biv_time_trt$time_trt %in% c(sapply(as.character(list.time_name), paste, as.character(list.trt_name), sep = ".")), ]
df_biv0_gplot <- df_time_trt_biv0[df_time_trt_biv0$time_trt %in% c(sapply(as.character(list.time_name), paste, as.character(list.trt_name), sep = ".")), ]
df_biv_gplot$time_trt <- factor(df_biv_gplot$time_trt, levels = c(sapply(as.character(list.time_name), paste, as.character(list.trt_name), sep = ".")))
df_biv0_gplot$time_trt <- factor(df_biv0_gplot$time_trt, levels = c(sapply(as.character(list.time_name), paste, as.character(list.trt_name), sep = ".")))
}
}
}
if(all(length(list.trt_name) == 1 & length(list.time_name) > 1)) {
if(all(list.trt_name == "none" & list.time_name %in% as.character(1:n_time))){
df_biv_gplot <- df_biv_time[df_biv_time$time %in% as.character(list.time_name), ]
df_biv0_gplot <- df_time_biv0[df_time_biv0$time %in% as.character(list.time_name), ]
df_biv_gplot$time <- factor(df_biv_gplot$time, levels = list.time_name)
df_biv0_gplot$time <- factor(df_biv0_gplot$time, levels = list.time_name)}
if(all(list.trt_name %in% trt_lev & all(list.time_name %in% as.character(1:n_time)))){
df_biv_gplot <- df_biv_time_trt[df_biv_time_trt$time_trt %in% c(sapply(as.character(list.time_name), paste, as.character(list.trt_name), sep = ".")), ]
df_biv0_gplot <- df_time_trt_biv0[df_time_trt_biv0$time_trt %in% c(sapply(as.character(list.time_name), paste, as.character(list.trt_name), sep = ".")), ]
df_biv_gplot$time_trt <- factor(df_biv_gplot$time_trt, levels = c(sapply(as.character(list.time_name), paste, as.character(list.trt_name), sep = ".")))
df_biv0_gplot$time_trt <- factor(df_biv0_gplot$time_trt, levels = c(sapply(as.character(list.time_name), paste, as.character(list.trt_name), sep = ".")))}
}
if(all(length(list.trt_name) > 1 & length(list.time_name) == 1)) {
if(all(all(list.trt_name %in% trt_lev) & list.time_name %in% as.character(1:n_time))){
df_biv_gplot <- df_biv_time_trt[df_biv_time_trt$time_trt %in% c(sapply(as.character(list.time_name), paste, as.character(list.trt_name), sep = ".")), ]
df_biv0_gplot <- df_time_trt_biv0[df_time_trt_biv0$time_trt %in% c(sapply(as.character(list.time_name), paste, as.character(list.trt_name), sep = ".")), ]
df_biv_gplot$time_trt <- factor(df_biv_gplot$time_trt, levels = c(sapply(as.character(list.time_name), paste, as.character(list.trt_name), sep = ".")))
df_biv0_gplot$time_trt <- factor(df_biv0_gplot$time_trt, levels = c(sapply(as.character(list.time_name), paste, as.character(list.trt_name), sep = ".")))}
}
if(all(length(list.trt_name) > 1 & length(list.time_name) > 1)) {
if(all(all(list.trt_name %in% trt_lev) & all(list.time_name %in% as.character(1:n_time)))){
df_biv_gplot <- df_biv_time_trt[df_biv_time_trt$time_trt %in% c(sapply(as.character(list.time_name), paste, as.character(list.trt_name), sep = ".")), ]
df_biv0_gplot <- df_time_trt_biv0[df_time_trt_biv0$time_trt %in% c(sapply(as.character(list.time_name), paste, as.character(list.trt_name), sep = ".")), ]
df_biv_gplot$time_trt <- factor(df_biv_gplot$time_trt, levels = c(sapply(as.character(list.time_name), paste, as.character(list.trt_name), sep = ".")))
df_biv0_gplot$time_trt <- factor(df_biv0_gplot$time_trt, levels = c(sapply(as.character(list.time_name), paste, as.character(list.trt_name), sep = ".")))}
}
}
color_sheme_rep <- color_scheme_get()[1]
color_sheme_rep_border <- color_scheme_get()[2]
color_sheme_obs <- color_scheme_get()[5]
color_sheme_obs_border <- color_scheme_get()[6]
set_theme <- bayesplot_theme_get()
if(!bpp) {
if(exists("bins", where = exArgs)) {bins = exArgs$bins} else { bins = 30}
if(exists("size", where = exArgs)) { size = exArgs$size} else { size = 0.5}
gplot_ppc <- ggplot2::ggplot(df_biv_gplot, ggplot2::aes(x = values)) +
ggplot2::geom_histogram(fill = color_sheme_rep$light,
color = color_sheme_rep_border$light_highlight, bins = bins) +
ggplot2::geom_vline(data = df_biv0_gplot, ggplot2::aes(xintercept = values), color = color_sheme_obs$dark,
linewidth = size * 1.5, linetype = "dashed")
}
if(bpp) {
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("size.annotate", where = exArgs)) { size.annotate = exArgs$size.annotate} else { size.annotate = 3.5}
if(exists("alpha", where = exArgs)) { alpha = exArgs$alpha} else { alpha = 1}
if(exists("vjust", where = exArgs)) { vjust = exArgs$vjust } else { vjust = 1}
gplot_ppc <- ggplot2::ggplot(df_biv_gplot, ggplot2::aes(x = values)) +
ggplot2::geom_density(fill = color_sheme_rep$light, color = color_sheme_rep_border$light_highlight,
trim = trim, linewidth = size, alpha = alpha) +
ggplot2::geom_vline(data = df_biv0_gplot, ggplot2::aes(xintercept = values), color = color_sheme_obs$dark,
linewidth = size * 1.5, linetype = "dashed") +
ggplot2::geom_text(data = df_biv0_gplot, mapping = ggplot2::aes(x = -Inf, y = Inf, label = text),
hjust = hjust, vjust = vjust, size = size.annotate)
}
}
gg_leg <- ggplot2::theme(legend.position = legend.pos)
if(theme == "classic") {
gg_par <- ggplot2::theme(panel.grid.major = ggplot2::element_blank(),
panel.grid.minor = ggplot2::element_blank(),
panel.background = ggplot2::element_blank(),
axis.line = ggplot2::element_line(colour = "black"),
strip.background = ggplot2::element_rect(fill="white"),
strip.text.x = ggplot2::element_text(size = 12, color = "black", face = "bold"),
strip.text.y = ggplot2::element_text(size = 12, color = "black", face = "bold"))
}
if(theme == "base") gg_par <- ggthemes::theme_base()
if(theme == "calc") gg_par <- ggthemes::theme_calc()
if(theme == "economist") gg_par <- ggthemes::theme_economist()
if(theme == "excel") gg_par <- ggthemes::theme_excel()
if(theme == "few") gg_par <- ggthemes::theme_few()
if(theme == "538") gg_par <- ggthemes::theme_fivethirtyeight()
if(theme == "gdocs") gg_par <- ggthemes::theme_gdocs()
if(theme == "hc") gg_par <- ggthemes::theme_hc()
if(theme == "par") gg_par <- ggthemes::theme_par()
if(theme == "solarized") gg_par <- ggthemes::theme_solarized()
if(theme == "pander") gg_par <- ggthemes::theme_pander()
if(theme == "stata") gg_par <- ggthemes::theme_stata()
if(theme == "tufte") gg_par <- ggthemes::theme_tufte()
if(theme == "wsj") gg_par <- ggthemes::theme_wsj()
if(x$data_format == "wide") {
if(length(gplot_e) != 0) {
gplot_e <- gplot_e + gg_leg + gg_par
gplot_c <- gplot_c + gg_leg + gg_par
gplot_e_trt <- lapply(gplot_e_trt, function(x) x + gg_leg + gg_par)
gplot_c_trt <- lapply(gplot_c_trt, function(x) x + gg_leg + gg_par)
if(length(trt) == 1) {
if(trt == "none") {
if(outcome == "both") {
gplot_ppc <- ggpubr::ggarrange(gplot_e, gplot_c, labels = c("effects", "costs"),
ncol = 2, nrow = 1, common.legend = TRUE,
legend = legend.pos, hjust = hjust, vjust = vjust)
}
if(outcome == "effects") {
gplot_ppc <- ggpubr::ggarrange(gplot_e, labels = "effects",
ncol = 1, nrow = 1, common.legend = TRUE,
legend = legend.pos, hjust = hjust, vjust = vjust)
}
if(outcome == "costs") {
gplot_ppc <- ggpubr::ggarrange(gplot_c, labels = "costs",
ncol = 1, nrow = 1, common.legend = TRUE,
legend = legend.pos, hjust = hjust, vjust = vjust)
}
} else if(trt == "all") {
if(outcome == "both") {
gplot_ppc <- ggpubr::ggarrange(plotlist = c(gplot_e_trt, gplot_c_trt),
labels = c(paste("effects", names(gplot_e_trt), sep = " "),
paste("costs", names(gplot_e_trt), sep = " ")),
ncol = length(names(gplot_e_trt)), nrow = 2, common.legend = TRUE,
legend = legend.pos, hjust = hjust, vjust = vjust)
}
if(outcome == "effects") {
gplot_ppc <- ggpubr::ggarrange(plotlist = gplot_e_trt,
labels = paste("effects", names(gplot_e_trt), sep = " "),
ncol = length(names(gplot_e_trt)), nrow = 1, common.legend = TRUE,
legend = legend.pos, hjust = hjust, vjust = vjust)
}
if(outcome == "costs") {
gplot_ppc <- ggpubr::ggarrange(plotlist = gplot_c_trt,
labels = paste("costs", names(gplot_c_trt), sep = " "),
ncol = length(names(gplot_c_trt)), nrow = 1, common.legend = TRUE,
legend = legend.pos, hjust = hjust, vjust = vjust)
}
} else {
gplot_e_trt <- gplot_e_trt[list.trt_name]; gplot_c_trt <- gplot_c_trt[list.trt_name]
if(outcome == "both") {
gplot_ppc <- ggpubr::ggarrange(plotlist = c(gplot_e_trt, gplot_c_trt),
labels = c(paste("effects", names(gplot_e_trt), sep = " "),
paste("costs", names(gplot_e_trt), sep = " ")),
ncol = length(names(gplot_e_trt)), nrow = 2, common.legend = TRUE,
legend = legend.pos, hjust = hjust, vjust = vjust)
}
if(outcome == "effects") {
gplot_ppc <- ggpubr::ggarrange(plotlist = gplot_e_trt,
labels = paste("effects", names(gplot_e_trt), sep = " "),
ncol = length(names(gplot_e_trt)), nrow = 1, common.legend = TRUE,
legend = legend.pos, hjust = hjust, vjust = vjust)
}
if(outcome == "costs") {
gplot_ppc <- ggpubr::ggarrange(plotlist = gplot_c_trt,
labels = paste("costs", names(gplot_c_trt), sep = " "),
ncol = length(names(gplot_c_trt)), nrow = 1, common.legend = TRUE,
legend = legend.pos, hjust = hjust, vjust = vjust)
}
}
}
if(length(trt) > 1) {
gplot_e_trt <- gplot_e_trt[list.trt_name]; gplot_c_trt <- gplot_c_trt[list.trt_name]
if(outcome == "both") {
gplot_ppc <- ggpubr::ggarrange(plotlist = c(gplot_e_trt, gplot_c_trt),
labels = c(paste("effects", names(gplot_e_trt), sep = " "),
paste("costs", names(gplot_e_trt), sep = " ")),
ncol = length(names(gplot_e_trt)), nrow = 2, common.legend = TRUE,
legend = legend.pos, hjust = hjust, vjust = vjust)
}
if(outcome == "effects") {
gplot_ppc <- ggpubr::ggarrange(plotlist = gplot_e_trt,
labels = paste("effects", names(gplot_e_trt), sep = " "),
ncol = length(names(gplot_e_trt)), nrow = 1, common.legend = TRUE,
legend = legend.pos, hjust = hjust, vjust = vjust)
}
if(outcome == "costs") {
gplot_ppc <- ggpubr::ggarrange(plotlist = gplot_c_trt,
labels = paste("costs", names(gplot_c_trt), sep = " "),
ncol = length(names(gplot_c_trt)), nrow = 1, common.legend = TRUE,
legend = legend.pos, hjust = hjust, vjust = vjust)
}
}
}
if(length(gplot_e) == 0) {
if(length(list.trt_name) == 1) {
if(list.trt_name == "none"){
if(is.null(corr)) {
gplot_ppc <- gplot_ppc + ggplot2::facet_wrap( ~ outcome, scales = "free")}
}
if(list.trt_name %in% trt_lev){
if(is.null(corr)) {
gplot_ppc <- gplot_ppc + ggplot2::facet_wrap( ~ outcome, scales = "free")}
}
}
if(length(list.trt_name) > 1) {
if(all(list.trt_name %in% trt_lev)){
if(is.null(corr)) {
gplot_ppc <- gplot_ppc + ggplot2::facet_wrap(outcome ~ trt, scales = "free")
} else {
gplot_ppc <- gplot_ppc + ggplot2::facet_wrap( ~ trt, scales = "free")
}
}
if(!any(list.trt_name %in% trt_lev)){
stop("Please provide only valid treatment names or values.")
}
}
}
}
if(x$data_format == "long") {
if(length(gplot_e_time) != 0) {
names(gplot_e_time) <- names(gplot_c_time) <- as.character(1:n_time)
gplot_e_time <- lapply(gplot_e_time, function(x) x + gg_leg + gg_par)
gplot_c_time <- lapply(gplot_c_time, function(x) x + gg_leg + gg_par)
gplot_e_time_trt <- lapply(gplot_e_time_trt, function(x) x + gg_leg + gg_par)
gplot_c_time_trt <- lapply(gplot_c_time_trt, function(x) x + gg_leg + gg_par)
if(all(length(trt) == 1 & length(list.time_name) == 1) | all(length(trt) == 1 & length(list.time_name) > 1)) {
if(trt == "none"){
gplot_e_time <- gplot_e_time[list.time_name]; gplot_c_time <- gplot_c_time[list.time_name]
if(outcome == "both") {
gplot_ppc <- ggpubr::ggarrange(plotlist = c(gplot_e_time, gplot_c_time),
labels = c(paste("effects", names(gplot_e_time), sep = " "),
paste("costs", names(gplot_e_time), sep = " ")),
ncol = length(names(gplot_e_time)), nrow = 2, common.legend = TRUE,
legend = legend.pos, hjust = hjust, vjust = vjust)
}
if(outcome == "effects") {
gplot_ppc <- ggpubr::ggarrange(plotlist = gplot_e_time,
labels = paste("effects", names(gplot_e_time), sep = " "),
ncol = length(names(gplot_e_time)), nrow = 1, common.legend = TRUE,
legend = legend.pos, hjust = hjust, vjust = vjust)
}
if(outcome == "costs") {
gplot_ppc <- ggpubr::ggarrange(plotlist = gplot_c_time,
labels = paste("costs", names(gplot_c_time), sep = " "),
ncol = length(names(gplot_c_time)), nrow = 1, common.legend = TRUE,
legend = legend.pos, hjust = hjust, vjust = vjust)
}
} else if(trt == "all") {
list.trt_name_select <- c(sapply(as.character(list.time_name), paste, as.character(list.trt_name), sep = "."))
gplot_e_time_trt <- gplot_e_time_trt[list.trt_name_select]; gplot_c_time_trt <- gplot_c_time_trt[list.trt_name_select]
if(outcome == "both") {
gplot_ppc <- ggpubr::ggarrange(plotlist = c(gplot_e_time_trt, gplot_c_time_trt),
labels = c(paste("effects", names(gplot_e_time_trt), sep = " "),
paste("costs", names(gplot_e_time_trt), sep = " ")),
ncol = length(names(gplot_e_time_trt)), nrow = 2, common.legend = TRUE,
legend = legend.pos, hjust = hjust, vjust = vjust)
}
if(outcome == "effects") {
gplot_ppc <- ggpubr::ggarrange(plotlist = gplot_e_time_trt,
labels = paste("effects", names(gplot_e_time_trt), sep = " "),
ncol = length(names(gplot_e_time_trt)), nrow = 1, common.legend = TRUE,
legend = legend.pos, hjust = hjust, vjust = vjust)
}
if(outcome == "costs") {
gplot_ppc <- ggpubr::ggarrange(plotlist = gplot_c_time_trt,
labels = paste("costs", names(gplot_c_time_trt), sep = " "),
ncol = length(names(gplot_c_time_trt)), nrow = 1, common.legend = TRUE,
legend = legend.pos, hjust = hjust, vjust = vjust)
}
} else {
list.trt_name_select <- c(sapply(as.character(list.time_name), paste, as.character(list.trt_name), sep = "."))
gplot_e_time_trt <- gplot_e_time_trt[list.trt_name_select]; gplot_c_time_trt <- gplot_c_time_trt[list.trt_name_select]
if(outcome == "both") {
gplot_ppc <- ggpubr::ggarrange(plotlist = c(gplot_e_time_trt, gplot_c_time_trt),
labels = c(paste("effects", names(gplot_e_time_trt), sep = " "),
paste("costs", names(gplot_e_time_trt), sep = " ")),
ncol = length(names(gplot_e_time_trt)), nrow = 2, common.legend = TRUE,
legend = legend.pos, hjust = hjust, vjust = vjust)
}
if(outcome == "effects") {
gplot_ppc <- ggpubr::ggarrange(plotlist = gplot_e_time_trt,
labels = paste("effects", names(gplot_e_time_trt), sep = " "),
ncol = length(names(gplot_e_time_trt)), nrow = 1, common.legend = TRUE,
legend = legend.pos, hjust = hjust, vjust = vjust)
}
if(outcome == "costs") {
gplot_ppc <- ggpubr::ggarrange(plotlist = gplot_c_time_trt,
labels = paste("costs", names(gplot_c_time_trt), sep = " "),
ncol = length(names(gplot_c_time_trt)), nrow = 1, common.legend = TRUE,
legend = legend.pos, hjust = hjust, vjust = vjust)
}
}
}
if(all(length(trt) > 1 & length(list.time_name) == 1) | all(length(trt) > 1 & length(list.time_name) > 1)) {
list.trt_name_select <- c(sapply(as.character(list.time_name), paste, as.character(list.trt_name), sep = "."))
gplot_e_time_trt <- gplot_e_time_trt[list.trt_name_select]; gplot_c_time_trt <- gplot_c_time_trt[list.trt_name_select]
if(outcome == "both") {
gplot_ppc <- ggpubr::ggarrange(plotlist = c(gplot_e_time_trt, gplot_c_time_trt),
labels = c(paste("effects", names(gplot_e_time_trt), sep = " "),
paste("costs", names(gplot_e_time_trt), sep = " ")),
ncol = length(names(gplot_e_time_trt)), nrow = 2, common.legend = TRUE,
legend = legend.pos, hjust = hjust, vjust = vjust)
}
if(outcome == "effects") {
gplot_ppc <- ggpubr::ggarrange(plotlist = gplot_e_time_trt,
labels = paste("effects", names(gplot_e_time_trt), sep = " "),
ncol = length(names(gplot_e_time_trt)), nrow = 1, common.legend = TRUE,
legend = legend.pos, hjust = hjust, vjust = vjust)
}
if(outcome == "costs") {
gplot_ppc <- ggpubr::ggarrange(plotlist = gplot_c_time_trt,
labels = paste("costs", names(gplot_c_time_trt), sep = " "),
ncol = length(names(gplot_c_time_trt)), nrow = 1, common.legend = TRUE,
legend = legend.pos, hjust = hjust, vjust = vjust)
}
}
}
if(length(gplot_e_time) == 0) {
if(all(length(list.trt_name) == 1 & length(list.time_name) == 1)) {
if(list.time_name %in% as.character(1:n_time)) {
if(list.trt_name == "none"){
if(is.null(corr)) {
gplot_ppc <- gplot_ppc + ggplot2::facet_wrap( ~ outcome, scales = "free")
}
}
if(list.trt_name %in% trt_lev){
if(is.null(corr)) {
gplot_ppc <- gplot_ppc + ggplot2::facet_wrap(outcome ~ time_trt, scales = "free")
} else {
gplot_ppc <- gplot_ppc + ggplot2::facet_wrap( ~ time_trt, scales = "free")
}
}
}
}
if(all(length(list.trt_name) == 1 & length(list.time_name) > 1)) {
if(all(list.trt_name == "none" & list.time_name %in% as.character(1:n_time))){
if(is.null(corr)) {
gplot_ppc <- gplot_ppc + ggplot2::facet_wrap(outcome ~ time, scales = "free")
} else {
gplot_ppc <- gplot_ppc + ggplot2::facet_wrap( ~ time, scales = "free")
}
}
if(all(list.trt_name %in% trt_lev & all(list.time_name %in% as.character(1:n_time)))){
if(is.null(corr)) {
gplot_ppc <- gplot_ppc + ggplot2::facet_wrap(outcome ~ time_trt, scales = "free")
} else {
gplot_ppc <- gplot_ppc + ggplot2::facet_wrap( ~ time_trt, scales = "free")
}
}
}
if(all(length(list.trt_name) > 1 & length(list.time_name) == 1)) {
if(all(all(list.trt_name %in% trt_lev) & list.time_name %in% as.character(1:n_time))){
if(is.null(corr)) {
gplot_ppc <- gplot_ppc + ggplot2::facet_wrap(outcome ~ time_trt, scales = "free")
} else {
gplot_ppc <- gplot_ppc + ggplot2::facet_wrap( ~ time_trt, scales = "free")
}
}
}
if(all(length(list.trt_name) > 1 & length(list.time_name) > 1)) {
if(all(all(list.trt_name %in% trt_lev) & all(list.time_name %in% as.character(1:n_time)))){
if(is.null(corr)) {
gplot_ppc <- gplot_ppc + ggplot2::facet_wrap(outcome ~ time_trt, scales = "free")
} else {
gplot_ppc <- gplot_ppc + ggplot2::facet_wrap( ~ time_trt, scales = "free")
}
}
}
}
}
gplot_ppc <- gplot_ppc + gg_leg + gg_par
return(gplot_ppc)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.