Nothing
#' @title Create BayesTools ensemble summary tables
#'
#' @description Creates estimate summaries based on posterior
#' distributions created by [mix_posteriors], inference summaries
#' based on inference created by [ensemble_inference], or ensemble
#' summary/diagnostics based on a list of [models_inference] models
#' (or [marginal_inference] in case of [marginal_estimates_table]).
#'
#' @param samples posterior samples created by [mix_posteriors]
#' @param inference model inference created by [ensemble_inference]
#' @param parameters character vector of parameters (or a
#' named list with of character vectors for summary and
#' diagnostics tables) specifying the parameters
#' (and their grouping) for the summary table
#' @param models list of [models_inference] model objects,
#' each of which containing a list of \code{priors}
#' and \code{inference} object, The \code{inference} must be a
#' named list with information about the model: model number
#' \code{m_number}, marginal likelihood \code{marglik}, prior and
#' posterior probability \code{prior_prob} and \code{post_prob},
#' inclusion Bayes factor \code{inclusion_BF}, and fit summary
#' generated by [runjags_estimates_table] for the diagnostics
#' table
#' @param probs quantiles for parameter estimates
#' @param logBF whether the Bayes factor should be on log scale
#' @param BF01 whether the Bayes factor should be inverted
#' @param short_name whether the prior distribution names should be
#' shortened. Defaults to \code{FALSE}.
#' @param remove_spike_0 whether prior distributions equal to spike
#' at 0 should be removed from the \code{prior_list}
#' @param transform_factors whether factors with orthonormal/meandif
#' prior distribution should be transformed to differences from the
#' grand mean
#' @param transform_orthonormal (to be depreciated) whether factors
#' with orthonormal prior distributions should be transformed to
#' differences from the grand mean
#' @param title title to be added to the table
#' @param footnotes footnotes to be added to the table
#' @param warnings warnings to be added to the table
#' @param formula_prefix whether the parameter prefix from formula should
#' be printed. Defaults to \code{TRUE}.
#'
#'
#' @return \code{ensemble_estimates_table} returns a table with the
#' model-averaged estimates, \code{ensemble_inference_table} returns
#' a table with the prior and posterior probabilities and inclusion
#' Bayes factors, \code{ensemble_summary_table} returns a table with
#' overview of the models included in the ensemble, and
#' \code{ensemble_diagnostics_table} returns an overview of the MCMC
#' diagnostics for the models included in the ensemble. All of the
#' tables are objects of class 'BayesTools_table'.
#'
#' @export ensemble_estimates_table
#' @export ensemble_inference_table
#' @export ensemble_summary_table
#' @export ensemble_diagnostics_table
#' @export ensemble_estimates_empty_table
#' @export ensemble_inference_empty_table
#' @export ensemble_summary_empty_table
#' @export ensemble_diagnostics_empty_table
#' @export marginal_estimates_table
#' @name BayesTools_ensemble_tables
#'
#' @seealso [ensemble_inference] [mix_posteriors] [BayesTools_model_tables]
NULL
#' @rdname BayesTools_ensemble_tables
ensemble_estimates_table <- function(samples, parameters, probs = c(0.025, 0.95), title = NULL, footnotes = NULL, warnings = NULL, transform_factors = FALSE, transform_orthonormal = FALSE, formula_prefix = TRUE){
# check input
check_char(parameters, "parameters", check_length = 0)
check_list(samples, "samples", check_names = parameters, all_objects = TRUE, allow_other = TRUE)
check_real(probs, "probs", lower = 0, upper = 1, check_length = 0, allow_NULL = TRUE)
check_char(title, "title", allow_NULL = TRUE)
check_char(footnotes, "footnotes", check_length = 0, allow_NULL = TRUE)
check_char(warnings, "warnings", check_length = 0, allow_NULL = TRUE)
check_bool(transform_factors, "transform_factors")
check_bool(transform_orthonormal, "transform_orthonormal")
check_bool(formula_prefix, "formula_prefix")
# depreciate
transform_factors <- .depreciate.transform_orthonormal(transform_orthonormal, transform_factors)
# transform meandif/orthonormal posterior
if(transform_factors){
samples <- transform_factor_samples(samples)
}
# extract values
estimates_table <- NULL
for(parameter in parameters){
if(is.matrix(samples[[parameter]])){
par_summary <- cbind(
"Mean" = apply(samples[[parameter]], 2, mean),
"Median" = apply(samples[[parameter]], 2, stats::median)
)
for(i in seq_along(probs)){
par_summary <- cbind(par_summary, apply(samples[[parameter]], 2, stats::quantile, probs = probs[i]))
colnames(par_summary)[ncol(par_summary)] <- probs[i]
}
if(inherits(samples[[parameter]], "mixed_posteriors.formula")){
parameter_name <- format_parameter_names(colnames(samples[[parameter]]), formula_parameters = attr(samples[[parameter]], "formula_parameter"), formula_prefix = formula_prefix)
}else{
parameter_name <- colnames(samples[[parameter]])
}
rownames(par_summary) <- parameter_name
estimates_table <- rbind(estimates_table, par_summary)
}else if(is.numeric(samples[[parameter]])){
par_summary <- c(
"Mean" = mean(samples[[parameter]]),
"Median" = stats::median((samples[[parameter]]))
)
for(i in seq_along(probs)){
par_summary <- c(par_summary, stats::quantile(samples[[parameter]], probs = probs[i]))
names(par_summary)[length(par_summary)] <- probs[i]
}
estimates_table <- rbind(estimates_table, par_summary)
if(inherits(samples[[parameter]], "mixed_posteriors.formula")){
parameter_name <- gsub(
paste0(attr(samples[[parameter]], "formula_parameter"), "_"),
if(formula_prefix) paste0("(", attr(samples[[parameter]], "formula_parameter"), ") ") else "",
parameter)
parameter_name <- gsub("__xXx__", ":", parameter_name)
}else{
parameter_name <- parameter
}
rownames(estimates_table)[nrow(estimates_table)] <- parameter_name
}else{
stop("Uknown parameter type.")
}
}
# prepare output
estimates_table <- data.frame(estimates_table)
colnames(estimates_table) <- gsub("X", "", colnames(estimates_table))
class(estimates_table) <- c("BayesTools_table", "BayesTools_ensemble_summary", class(estimates_table))
attr(estimates_table, "type") <- rep("estimate", ncol(estimates_table))
attr(estimates_table, "rownames") <- TRUE
attr(estimates_table, "title") <- title
attr(estimates_table, "footnotes") <- footnotes
attr(estimates_table, "warnings") <- warnings
return(estimates_table)
}
#' @rdname BayesTools_ensemble_tables
ensemble_inference_table <- function(inference, parameters, logBF = FALSE, BF01 = FALSE, title = NULL, footnotes = NULL, warnings = NULL){
# check input
check_char(parameters, "parameters", check_length = 0)
check_list(inference, "inference", check_names = parameters, all_objects = TRUE, allow_other = TRUE)
check_bool(logBF, "logBF")
check_bool(BF01, "BF01")
check_char(title, "title", allow_NULL = TRUE)
check_char(footnotes, "footnotes", check_length = 0, allow_NULL = TRUE)
check_char(warnings, "warnings", check_length = 0, allow_NULL = TRUE)
if(attr(inference,"conditional"))
stop("The inference object cannot be 'conditional'.")
# extract values
inference_table <- NULL
n_models <- NULL
for(parameter in parameters){
inference_table <- rbind(inference_table, c(
"models" = sum(!attr(inference[[parameter]], "is_null")),
"prior_prob" = sum(inference[[parameter]][["prior_probs"]][!attr(inference[[parameter]], "is_null")]),
"post_prob" = sum(inference[[parameter]][["post_probs"]][!attr(inference[[parameter]], "is_null")] ),
"inclusion_BF" = inference[[parameter]][["BF"]]
))
rownames(inference_table)[nrow(inference_table)] <- attr(inference[[parameter]], "parameter_name")
n_models <- c(n_models, length(attr(inference[[parameter]], "is_null")))
}
inference_table <- data.frame(inference_table)
# format BF
inference_table[,"inclusion_BF"] <- format_BF(inference_table[,"inclusion_BF"], logBF = logBF, BF01 = BF01, inclusion = TRUE)
# prepare output
class(inference_table) <- c("BayesTools_table", "BayesTools_ensemble_summary", class(inference_table))
attr(inference_table, "type") <- c("n_models", "prior_prob", "post_prob", "inclusion_BF")
attr(inference_table, "rownames") <- TRUE
attr(inference_table, "n_models") <- n_models
attr(inference_table, "title") <- title
attr(inference_table, "footnotes") <- footnotes
attr(inference_table, "warnings") <- warnings
return(inference_table)
}
#' @rdname BayesTools_ensemble_tables
ensemble_summary_table <- function(models, parameters, logBF = FALSE, BF01 = FALSE, title = NULL, footnotes = NULL, warnings = NULL,
remove_spike_0 = TRUE, short_name = FALSE){
# check input
check_list(models, "models")
for(i in seq_along(models)){
model <- models[[i]]
check_list(model, "model", check_names = "inference", allow_other = TRUE, all_objects = TRUE)
prior_list <- attr(model[["fit"]], "prior_list")
check_list(prior_list, "model:priors")
if(!all(sapply(prior_list, is.prior)))
stop("'model:priors' must be a list of priors.")
model_inference <- model[["inference"]]
check_list(model_inference, "model:inference", check_names = c("m_number", "marglik", "prior_prob", "post_prob", "inclusion_BF"), allow_other = TRUE, all_objects = TRUE)
check_int(model_inference[["m_number"]], "model_inference:model_number")
check_real(model_inference[["marglik"]], "model_inference:marglik")
check_real(model_inference[["prior_prob"]], "model_inference:prior_prob", lower = 0, upper = 1)
check_real(model_inference[["post_prob"]], "model_inference:post_prob", lower = 0, upper = 1)
check_real(model_inference[["inclusion_BF"]], "model_inference:inclusion_BF", lower = 0)
}
if(is.list(parameters)){
check_list(parameters, "parameters")
check_char(names(parameters), "names(parameters)", check_length = FALSE, allow_NULL = length(parameters) == 0)
sapply(parameters, check_char, name = "parameters", check_length = FALSE)
}else{
check_char(parameters, "parameters", check_length = FALSE)
if(is.null(names(parameters))){
names(parameters) <- parameters
}
}
check_bool(logBF, "logBF")
check_bool(BF01, "BF01")
check_char(title, "title", allow_NULL = TRUE)
check_char(footnotes, "footnotes", check_length = 0, allow_NULL = TRUE)
check_char(warnings, "warnings", check_length = 0, allow_NULL = TRUE)
check_bool(short_name, "short_name")
# create the output
ensemble_table <- .ensemble_table_foundation(models, parameters, remove_spike_0, short_name)
ensemble_table <- cbind(
ensemble_table,
"prior_prob" = sapply(models, function(model)model[["inference"]][["prior_prob"]]),
"marglik" = sapply(models, function(model)model[["inference"]][["marglik"]]),
"post_prob" = sapply(models, function(model)model[["inference"]][["post_prob"]]),
"inclusion_BF" = sapply(models, function(model)model[["inference"]][["inclusion_BF"]])
)
ensemble_table[,"inclusion_BF"] <- format_BF(ensemble_table[,"inclusion_BF"], logBF = logBF, BF01 = BF01, inclusion = TRUE)
# prepare output
class(ensemble_table) <- c("BayesTools_table", "BayesTools_ensemble_summary", class(ensemble_table))
attr(ensemble_table, "type") <- c("integer", rep("prior", length(parameters)), "prior_prob", "marglik", "post_prob", "inclusion_BF")
attr(ensemble_table, "rownames") <- FALSE
attr(ensemble_table, "title") <- title
attr(ensemble_table, "footnotes") <- footnotes
attr(ensemble_table, "warnings") <- warnings
return(ensemble_table)
}
#' @rdname BayesTools_ensemble_tables
ensemble_diagnostics_table <- function(models, parameters, title = NULL, footnotes = NULL, warnings = NULL, remove_spike_0 = TRUE, short_name = FALSE){
# check input
check_list(models, "models")
for(i in seq_along(models)){
model <- models[[i]]
check_list(model, "model", check_names = c("fit_summary", "inference"), allow_other = TRUE, all_objects = TRUE)
prior_list <- attr(model[["fit"]], "prior_list")
check_list(prior_list, "model:priors")
if(!all(sapply(prior_list, is.prior)))
stop("'model:priors' must be a list of priors.")
if(!is.null(model[["fit_summary"]]) && !(inherits(model[["fit_summary"]], "BayesTools_runjags_summary") || inherits(model[["fit_summary"]], "BayesTools_stan_summary")))
stop("'fit_summary' must be a runjags or rstan summary generated by 'JAGS_estimates_table()' / 'stan_estimates_table()'.")
model_inference <- model[["inference"]]
check_list(model_inference, "model:inference", check_names = "m_number", allow_other = TRUE, all_objects = TRUE)
check_int(model_inference[["m_number"]], "model_inference:model_number")
}
if(is.list(parameters)){
check_list(parameters, "parameters")
check_char(names(parameters), "names(parameters)", check_length = FALSE, allow_NULL = length(parameters) == 0)
sapply(parameters, check_char, name = "parameters", check_length = FALSE)
}else{
check_char(parameters, "parameters", check_length = FALSE)
if(is.null(names(parameters))){
names(parameters) <- parameters
}
}
check_char(title, "title", allow_NULL = TRUE)
check_char(footnotes, "footnotes", check_length = 0, allow_NULL = TRUE)
check_char(warnings, "warnings", check_length = 0, allow_NULL = TRUE)
check_bool(short_name, "short_name")
# create the output
ensemble_table <- .ensemble_table_foundation(models, parameters, remove_spike_0, short_name)
ensemble_table <- cbind(
ensemble_table,
"max_MCMC_error" = sapply(models, function(model){
MCMC_error <- model[["fit_summary"]][,"MCMC_error"]
if(all(is.na(MCMC_error))){
return(NA)
}else{
max(MCMC_error, na.rm = TRUE)
}
}),
"max_MCMC_SD_error" = sapply(models, function(model){
MCMC_SD_error <- model[["fit_summary"]][,"MCMC_SD_error"]
if(all(is.na(MCMC_SD_error))){
return(NA)
}else{
max(MCMC_SD_error, na.rm = TRUE)
}
}),
"min_ESS" = sapply(models, function(model){
ESS <- model[["fit_summary"]][,"ESS"]
if(all(is.na(ESS))){
return(NA)
}else{
min(ESS, na.rm = TRUE)
}
}),
"max_R_hat" = sapply(models, function(model){
Rhat <- model[["fit_summary"]][,"R_hat"]
if(all(is.na(Rhat))){
return(NA)
}else{
max(Rhat, na.rm = TRUE)
}
})
)
# prepare output
class(ensemble_table) <- c("BayesTools_table", "BayesTools_ensemble_summary", class(ensemble_table))
attr(ensemble_table, "type") <- c("integer", rep("prior", length(parameters)), "max_MCMC_error", "max_MCMC_SD_error", "min_ESS", "max_R_hat")
attr(ensemble_table, "rownames") <- FALSE
attr(ensemble_table, "title") <- title
attr(ensemble_table, "footnotes") <- footnotes
attr(ensemble_table, "warnings") <- warnings
return(ensemble_table)
}
#' @rdname BayesTools_ensemble_tables
ensemble_estimates_empty_table <- function(probs = c(0.025, 0.95), title = NULL, footnotes = NULL, warnings = NULL){
empty_table <- data.frame(matrix(nrow = 0, ncol = 2 + length(probs)))
colnames(empty_table) <- c("Mean", "Median", probs)
class(empty_table) <- c("BayesTools_table", "BayesTools_ensemble_summary", class(empty_table))
attr(empty_table, "type") <- rep("estimate", ncol(empty_table))
attr(empty_table, "rownames") <- TRUE
attr(empty_table, "title") <- title
attr(empty_table, "footnotes") <- footnotes
attr(empty_table, "warnings") <- warnings
return(empty_table)
}
#' @rdname BayesTools_ensemble_tables
ensemble_inference_empty_table <- function(title = NULL, footnotes = NULL, warnings = NULL){
empty_table <- data.frame(matrix(nrow = 0, ncol = 4))
colnames(empty_table) <- c("models", "prior_prob", "post_prob", "inclusion_BF")
class(empty_table) <- c("BayesTools_table", "BayesTools_ensemble_summary", class(empty_table))
attr(empty_table, "type") <- c("n_models", "prior_prob", "post_prob", "inclusion_BF")
attr(empty_table, "rownames") <- TRUE
attr(empty_table, "title") <- title
attr(empty_table, "footnotes") <- footnotes
attr(empty_table, "warnings") <- warnings
return(empty_table)
}
#' @rdname BayesTools_ensemble_tables
ensemble_summary_empty_table <- function(title = NULL, footnotes = NULL, warnings = NULL){
empty_table <- data.frame(matrix(nrow = 0, ncol = 5))
colnames(empty_table) <- c("Model", "prior_prob", "marglik", "post_prob", "inclusion_BF")
# prepare output
class(empty_table) <- c("BayesTools_table", "BayesTools_ensemble_summary", class(empty_table))
attr(empty_table, "type") <- c("integer", "prior_prob", "marglik", "post_prob", "inclusion_BF")
attr(empty_table, "rownames") <- FALSE
attr(empty_table, "title") <- title
attr(empty_table, "footnotes") <- footnotes
attr(empty_table, "warnings") <- warnings
return(empty_table)
}
#' @rdname BayesTools_ensemble_tables
ensemble_diagnostics_empty_table <- function(title = NULL, footnotes = NULL, warnings = NULL){
empty_table <- data.frame(matrix(nrow = 0, ncol = 5))
colnames(empty_table) <- c("Model", "max_MCMC_error", "max_MCMC_SD_error", "min_ESS", "max_R_hat")
# prepare output
class(empty_table) <- c("BayesTools_table", "BayesTools_ensemble_summary", class(empty_table))
attr(empty_table, "type") <- c("integer", "max_MCMC_error", "max_MCMC_SD_error", "min_ESS", "max_R_hat")
attr(empty_table, "rownames") <- FALSE
attr(empty_table, "title") <- title
attr(empty_table, "footnotes") <- footnotes
attr(empty_table, "warnings") <- warnings
return(empty_table)
}
#' @rdname BayesTools_ensemble_tables
marginal_estimates_table <- function(samples, inference, parameters, probs = c(0.025, 0.95), logBF = FALSE, BF01 = FALSE, title = NULL, footnotes = NULL, warnings = NULL, formula_prefix = TRUE){
# check input
check_char(parameters, "parameters", check_length = 0)
check_list(samples, "samples", check_names = parameters, all_objects = TRUE, allow_other = TRUE)
check_list(inference, "inference", check_names = parameters, all_objects = TRUE, allow_other = TRUE)
check_real(probs, "probs", lower = 0, upper = 1, check_length = 0, allow_NULL = TRUE)
check_bool(logBF, "logBF")
check_bool(BF01, "BF01")
check_char(title, "title", allow_NULL = TRUE)
check_char(footnotes, "footnotes", check_length = 0, allow_NULL = TRUE)
check_char(warnings, "warnings", check_length = 0, allow_NULL = TRUE)
check_bool(formula_prefix, "formula_prefix")
# extract values
estimates_table <- NULL
for(parameter in parameters){
# extract the relevant information
if(is.list(samples[[parameter]]) && length(samples[[parameter]]) > 1){
temp_samples <- do.call(cbind, samples[[parameter]])
temp_BF <- do.call(c, inference[[parameter]])
temp_warnings <- do.call(c, lapply(names(inference[[parameter]]), function(lvl) {
if(is.null(attr(inference[[parameter]][[lvl]], "warnings"))){
return()
}else{
paste0("[", lvl, "]: ", attr(inference[[parameter]][[lvl]], "warnings"))
}
}))
}else{
temp_samples <- matrix(samples[[parameter]][[1]], ncol = 1)
temp_BF <- inference[[parameter]][[1]]
if(is.null(attr(inference[[parameter]][[1]], "warnings"))){
temp_warnings <- NULL
}else{
temp_warnings <- paste0(if(names(inference[[parameter]]) != "intercept") paste0("[", names(inference[[parameter]]), "]: ") else ": ", attr(inference[[parameter]][[1]], "warnings"))
}
}
# add estimates
par_summary <- cbind(
"Mean" = apply(temp_samples, 2, mean),
"Median" = apply(temp_samples, 2, stats::median)
)
for(i in seq_along(probs)){
par_summary <- cbind(par_summary, apply(temp_samples, 2, stats::quantile, probs = probs[i]))
colnames(par_summary)[ncol(par_summary)] <- probs[i]
}
# add BF
par_summary <- cbind(par_summary, "inclusion_BF" = temp_BF)
# format parameter names
if(inherits(samples[[parameter]], "marginal_posterior.formula")){
if(length(names(samples[[parameter]])) == 1 && names(samples[[parameter]]) == "intercept"){
parameter_name <- parameter
}else{
parameter_name <- paste0(parameter, "[", names(samples[[parameter]]), "]")
}
parameter_name <- format_parameter_names(parameter_name, formula_parameters = attr(samples[[parameter]], "formula_parameter"), formula_prefix = formula_prefix)
}else{
parameter_name <- paste0(parameter, "[", names(samples[[parameter]]), "]")
}
rownames(par_summary) <- parameter_name
estimates_table <- rbind(estimates_table, par_summary)
# add warnings
if(!is.null(temp_warnings)){
warnings <- c(warnings, paste0(parameter, temp_warnings))
}
}
estimates_table[,"inclusion_BF"] <- format_BF(estimates_table[,"inclusion_BF"], logBF = logBF, BF01 = BF01, inclusion = TRUE)
# prepare output
estimates_table <- data.frame(estimates_table)
colnames(estimates_table) <- gsub("X", "", colnames(estimates_table))
class(estimates_table) <- c("BayesTools_table", "BayesTools_marginal_estimates", class(estimates_table))
attr(estimates_table, "type") <- c(rep("estimate", ncol(estimates_table) - 1), "inclusion_BF")
attr(estimates_table, "rownames") <- TRUE
attr(estimates_table, "title") <- title
attr(estimates_table, "footnotes") <- footnotes
attr(estimates_table, "warnings") <- warnings
return(estimates_table)
}
.ensemble_table_foundation <- function(models, parameters, remove_spike_0, short_name){
model_rows <- list()
for(i in seq_along(models)){
model_row <- list()
model_row[["Model"]] <- models[[i]][["inference"]][["m_number"]]
temp_prior_list <- attr(models[[i]][["fit"]], "prior_list")
for(p in seq_along(parameters)){
if(is.list(parameters)){
if(sum(names(temp_prior_list) %in% parameters[[p]]) == 1){
temp_prior <- temp_prior_list[[parameters[[p]][parameters[[p]] %in% names(temp_prior_list)]]]
}else if(sum(names(temp_prior_list) %in% parameters[[p]]) == 0){
temp_prior <- prior_none()
}else{
stop("More than one prior matching the specified grouping.")
}
}else{
if(any(names(temp_prior_list) == parameters[p])){
temp_prior <- temp_prior_list[[parameters[p]]]
}else{
temp_prior <- prior_none()
}
}
if(remove_spike_0 && is.prior.point(temp_prior) && temp_prior[["parameters"]][["location"]] == 0){
model_row[[names(parameters)[p]]] <- ""
}else if(is.prior.none(temp_prior)){
model_row[[names(parameters)[p]]] <- ""
}else{
model_row[[names(parameters)[p]]] <- print(temp_prior, silent = TRUE, short_name = short_name)
}
}
model_rows[[i]] <- model_row
}
summary_table <- data.frame(do.call(rbind, model_rows))
for(i in 1:ncol(summary_table)){
summary_table[,i] <- unlist(summary_table[,i])
}
colnames(summary_table) <- c("Model", names(parameters))
parameters <- unlist(parameters) # deal with potential matching of multiple parameters withing a name
for(p in seq_along(parameters)){
parameter_name <- parameters[p]
formula_parameter <- unique(unlist(lapply(models, function(m) attr(attr(m[["fit"]], "prior_list")[[parameter_name]], "parameter"))))
if(!is.null(unlist(formula_parameter))){
parameter_name <- gsub(paste0(formula_parameter, "_"), paste0("(", formula_parameter, ") "), parameter_name)
parameter_name <- gsub("__xXx__", ":", parameter_name)
colnames(summary_table)[colnames(summary_table) == parameters[p]] <- parameter_name
}
}
return(summary_table)
}
#' @title Create BayesTools model tables
#'
#' @description Creates model summary based on a model objects or
#' provides estimates table for a runjags fit.
#'
#' @param model model object containing a list of \code{priors}
#' and \code{inference} object, The \code{inference} must be a
#' named list with information about the model: model number
#' \code{m_number}, marginal likelihood \code{marglik}, prior and
#' posterior probability \code{prior_prob} and \code{post_prob},
#' and model inclusion Bayes factor \code{inclusion_BF}
#' @param fit runjags model fit
#' @param conditional summarizes estimates conditional on being included
#' in the model for spike and slab priors. Defaults to \code{FALSE}.
#' @param transformations named list of transformations to be applied
#' to specific parameters
#' @param model_description named list with additional description
#' to be added to the table
#' @param remove_inclusion whether estimates of the inclusion probabilities
#' should be excluded from the summary table. Defaults to \code{FALSE}.
#' @param remove_parameters parameters to be removed from the summary. Defaults
#' to \code{NULL}, i.e., including all parameters.
#' @inheritParams BayesTools_ensemble_tables
#'
#'
#' @return \code{model_summary_table} returns a table with
#' overview of the fitted model, \code{runjags_estimates_table} returns
#' a table with MCMC estimates, and \code{runjags_estimates_empty_table}
#' returns an empty estimates table. All of the tables are objects of
#' class 'BayesTools_table'.
#'
#' @export JAGS_summary_table
#' @export JAGS_estimates_table
#' @export JAGS_inference_table
#' @export model_summary_table
#' @export runjags_estimates_table
#' @export runjags_inference_table
#' @export model_summary_empty_table
#' @export JAGS_inference_empty_table
#' @export JAGS_estimates_empty_table
#' @export runjags_estimates_empty_table
#' @export runjags_inference_empty_table
#' @export stan_estimates_table
#' @name BayesTools_model_tables
#'
#' @seealso [BayesTools_ensemble_tables]
NULL
#' @rdname BayesTools_model_tables
model_summary_table <- function(model, model_description = NULL, title = NULL, footnotes = NULL, warnings = NULL, remove_spike_0 = TRUE, short_name = FALSE, formula_prefix = TRUE, remove_parameters = NULL){
# check input
check_list(model, "model", check_names = "inference", allow_other = TRUE, all_objects = TRUE)
prior_list <- attr(model[["fit"]], "prior_list")
check_list(prior_list, "model:priors")
if(!all(sapply(prior_list, is.prior)))
stop("'model:priors' must be a list of priors.")
model_inference <- model[["inference"]]
check_list(model_inference, "model:inference", check_names = c("m_number", "marglik", "prior_prob", "post_prob", "inclusion_BF"), allow_other = TRUE, all_objects = TRUE)
check_int(model_inference[["m_number"]], "model_inference:model_number")
check_real(model_inference[["marglik"]], "model_inference:marglik")
check_real(model_inference[["prior_prob"]], "model_inference:prior_prob", lower = 0, upper = 1)
check_real(model_inference[["post_prob"]], "model_inference:post_prob", lower = 0, upper = 1)
check_real(model_inference[["inclusion_BF"]], "model_inference:inclusion_BF", lower = 0)
check_list(model_description, "model_description", allow_NULL = TRUE)
check_bool(short_name, "short_name")
check_char(title, "title", allow_NULL = TRUE)
check_char(footnotes, "footnotes", check_length = 0, allow_NULL = TRUE)
check_char(warnings, "warnings", check_length = 0, allow_NULL = TRUE)
check_bool(formula_prefix, "formula_prefix")
check_char(remove_parameters, "remove_parameters", allow_NULL = TRUE, check_length = 0)
# prepare the columns
summary_names <- c(
"Model",
if(!is.null(model_description)) names(model_description),
"Prior prob.",
"log(marglik)",
"Post. prob.",
"Inclusion BF")
summary_values <- c(
model_inference[["m_number"]],
if(!is.null(model_description)) unlist(model_description),
.format_column(model_inference[["prior_prob"]], "probability"),
.format_column(model_inference[["marglik"]], "marglik"),
.format_column(model_inference[["post_prob"]], "probability"),
.format_column(model_inference[["inclusion_BF"]], "BF"))
summary_priors <- "Parameter prior distributions"
for(i in seq_along(prior_list)){
# get the prior name
if(is.prior.none(prior_list[[i]])){
next
}else if(remove_spike_0 && is.prior.point(prior_list[[i]]) && prior_list[[i]][["parameters"]][["location"]] == 0 || (names(prior_list)[[i]] %in% remove_parameters)){
next
}else if(is.prior.weightfunction(prior_list[[i]]) | is.prior.PET(prior_list[[i]]) | is.prior.PEESE(prior_list[[i]])){
temp_prior <- print(prior_list[[i]], silent = TRUE, short_name = short_name)
}else if(is.prior.simple(prior_list[[i]]) | is.prior.vector(prior_list[[i]]) | is.prior.factor(prior_list[[i]]) | is.prior.spike_and_slab(prior_list[[i]])){
temp_prior <- paste0(names(prior_list)[i], " ~ " , print(prior_list[[i]], silent = TRUE, short_name = short_name))
}else if(is.prior.point(prior_list[[i]])){
temp_prior <- paste0(names(prior_list)[i], " = " , print(prior_list[[i]], silent = TRUE, short_name = short_name))
}
# change the formula formatting
if(!is.null(attr(prior_list[[i]], "parameter"))){
temp_prior <- gsub(
paste0(attr(prior_list[[i]], "parameter"), "_"),
if(formula_prefix) paste0("(", attr(prior_list[[i]], "parameter"), ") ") else "",
temp_prior)
temp_prior <- gsub("__xXx__", ":", temp_prior)
}
summary_priors <- c(summary_priors, temp_prior)
}
if(length(summary_names) > length(summary_priors)){
summary_priors <- c(summary_priors, rep("", length(summary_names) - length(summary_priors)))
}else if(length(summary_names) < length(summary_priors)){
summary_names <- c(summary_names, rep("", length(summary_priors) - length(summary_names)))
summary_values <- c(summary_values, rep("", length(summary_priors) - length(summary_values)))
}
summary_names <- paste0(summary_names, " ")
summary_table <- data.frame(cbind(
summary_names,
summary_values,
rep(" ", length(summary_names)),
summary_priors
))
names(summary_table) <- NULL
# prepare output
class(summary_table) <- c("BayesTools_table", class(summary_table))
attr(summary_table, "type") <- c("string_left", "string", "string", "prior")
attr(summary_table, "rownames") <- FALSE
attr(summary_table, "as.matrix") <- TRUE
attr(summary_table, "title") <- title
attr(summary_table, "footnotes") <- footnotes
attr(summary_table, "warnings") <- warnings
return(summary_table)
}
#' @rdname BayesTools_model_tables
runjags_estimates_table <- function(fit, transformations = NULL, title = NULL, footnotes = NULL, warnings = NULL, conditional = FALSE, remove_spike_0 = TRUE, transform_factors = FALSE, transform_orthonormal = FALSE, formula_prefix = TRUE, remove_inclusion = FALSE, remove_parameters = NULL){
.check_runjags()
# most of the code is shared with .diagnostics_plot_data function (keep them in sync on update)
# check fits
if(!inherits(fit, "runjags"))
stop("'fit' must be a runjags fit")
if(!inherits(fit, "BayesTools_fit"))
stop("'fit' must be a BayesTools fit")
prior_list <- attr(fit, "prior_list")
check_list(prior_list, "prior_list")
if(!all(sapply(prior_list, is.prior)))
stop("'prior_list' must be a list of priors.")
check_list(transformations, "transformations", allow_NULL = TRUE)
if(!is.null(transformations) && any(!sapply(transformations, function(trans)is.function(trans[["fun"]]))))
stop("'transformations' must be list of functions in the 'fun' element.")
check_char(title, "title", allow_NULL = TRUE)
check_char(footnotes, "footnotes", check_length = 0, allow_NULL = TRUE)
check_char(warnings, "warnings", check_length = 0, allow_NULL = TRUE)
check_bool(remove_spike_0, "remove_spike_0")
check_bool(conditional, "conditional")
check_bool(transform_factors, "transform_factors")
check_bool(transform_orthonormal, "transform_orthonormal")
check_bool(formula_prefix, "formula_prefix")
check_char(remove_parameters, "remove_parameters", allow_NULL = TRUE, check_length = 0)
# depreciate
transform_factors <- .depreciate.transform_orthonormal(transform_orthonormal, transform_factors)
# obtain model information
invisible(utils::capture.output(runjags_summary <- suppressWarnings(summary(fit, silent.jags = TRUE))))
runjags_summary <- data.frame(runjags_summary)
model_samples <- suppressWarnings(coda::as.mcmc(fit))
# change HPD to quantile intervals
for(par in rownames(runjags_summary)){
runjags_summary[par, "Lower95"] <- stats::quantile(model_samples[,par], .025, na.rm = TRUE)
runjags_summary[par, "Upper95"] <- stats::quantile(model_samples[,par], .975, na.rm = TRUE)
}
# deal with missing median in case of non-stochastic variables
if(!any(colnames(runjags_summary) == "Median")){
runjags_summary[,"Median"] <- NA
}
# remove un-wanted estimates (or support values) - spike and slab priors already dealt with later
# also remove the item from prior list
for(i in rev(seq_along(prior_list))){
if(is.prior.weightfunction(prior_list[[i]])){
# remove etas
if(prior_list[[i]][["distribution"]] %in% c("one.sided", "two.sided")){
runjags_summary <- runjags_summary[!grepl("eta", rownames(runjags_summary)),,drop=FALSE]
}
# remove wrong diagnostics for the constant
runjags_summary[max(grep("omega", rownames(runjags_summary))),c("MCerr", "MC.ofSD","SSeff","psrf")] <- NA
# reorder
runjags_summary[grep("omega", rownames(runjags_summary)),] <- runjags_summary[rev(grep("omega", rownames(runjags_summary))),]
# rename
omega_cuts <- weightfunctions_mapping(prior_list[i], cuts_only = TRUE)
omega_names <- sapply(1:(length(omega_cuts)-1), function(i)paste0("omega[",omega_cuts[i],",",omega_cuts[i+1],"]"))
rownames(runjags_summary)[grep("omega", rownames(runjags_summary))] <- omega_names
# remove if requested
if("omega" %in% remove_parameters){
prior_list[[i]] <- NULL
runjags_summary <- runjags_summary[,!rownames(runjags_summary) %in% omega_names]
}
}else if((remove_spike_0 && is.prior.point(prior_list[[i]]) && prior_list[[i]][["parameters"]][["location"]] == 0) || (names(prior_list)[[i]] %in% remove_parameters)){
if(is.prior.factor(prior_list[[i]])){
runjags_summary <- runjags_summary[!rownames(runjags_summary) %in% .JAGS_prior_factor_names(names(prior_list)[i], prior_list[[i]]),,drop=FALSE]
}else{
runjags_summary <- runjags_summary[rownames(runjags_summary) != names(prior_list)[i],,drop=FALSE]
}
if(prior_list[[i]][["distribution"]] == "invgamma"){
runjags_summary <- runjags_summary[rownames(runjags_summary) != paste0("inv_",names(prior_list)[i]),,drop=FALSE]
}
prior_list[i] <- NULL
}else if(is.prior.simple(prior_list[[i]]) && prior_list[[i]][["distribution"]] == "invgamma"){
runjags_summary <- runjags_summary[rownames(runjags_summary) != paste0("inv_",names(prior_list)[i]),,drop=FALSE]
}
}
# remove transformations for removed variables
if(!is.null(transformations)){
transformations <- transformations[names(transformations) %in% names(prior_list)]
}
# simplify spike and slab priors to simple priors -- the samples and summary can be dealt with as any other prior
for(par in names(prior_list)){
if(is.prior.spike_and_slab(prior_list[[par]])){
# prepare parameter names
if(is.prior.factor(prior_list[[par]][["variable"]])){
if(.get_prior_factor_levels(prior_list[[par]][["variable"]]) == 1){
par_names <- par
}else{
par_names <- paste0(par, "[", 1:.get_prior_factor_levels(prior_list[[par]][["variable"]]), "]")
}
}else{
par_names <- par
}
# change the samples between conditional/averaged based on the preferences
if(conditional){
# set the spike samples to NA
model_samples[
model_samples[,colnames(model_samples) == paste0(par, "_indicator")] == 0,
colnames(model_samples) %in% par_names] <- NA
# recompute the summaries
runjags_summary[par_names, "Mean"] <- mean(model_samples[,par_names], na.rm = TRUE)
runjags_summary[par_names, "Median"] <- stats::median(model_samples[,par_names], na.rm = TRUE)
runjags_summary[par_names, "SD"] <- sd(model_samples[,par_names], na.rm = TRUE)
runjags_summary[par_names, "Lower95"] <- stats::quantile(model_samples[,par_names], .025, na.rm = TRUE)
runjags_summary[par_names, "Upper95"] <- stats::quantile(model_samples[,par_names], .975, na.rm = TRUE)
}
# remove the indicator
runjags_summary <- runjags_summary[rownames(runjags_summary) != paste0(par, "_indicator"),,drop=FALSE]
model_samples <- model_samples[colnames(runjags_summary) != paste0(par, "_indicator"),,drop=FALSE]
# remove the latent variable
runjags_summary <- runjags_summary[!rownames(runjags_summary) %in% gsub(par, paste0(par, "_variable"), par_names),,drop=FALSE]
model_samples <- model_samples[!colnames(runjags_summary) %in% gsub(par, paste0(par, "_variable"), par_names),,drop=FALSE]
# remove/rename the inclusions probabilities
if(remove_inclusion){
runjags_summary <- runjags_summary[rownames(runjags_summary) != paste0(par, "_inclusion"),,drop=FALSE]
model_samples <- model_samples[colnames(runjags_summary) != paste0(par, "_inclusion"),,drop=FALSE]
}else{
rownames(runjags_summary)[rownames(runjags_summary) == paste0(par, "_inclusion")] <- paste0(par, " (inclusion)")
colnames(model_samples)[colnames(model_samples) == paste0(par, "_inclusion")] <- paste0(par, " (inclusion)")
}
# modify the parameter list
prior_list[[par]] <- prior_list[[par]]$variable
}
}
# apply transformations (not orthornormal if they are to be returned transformed to diffs)
if(!is.null(transformations)){
for(par in names(transformations)){
if(!is.prior.factor(prior_list[[par]])){
# non-factor priors
model_samples[,par] <- do.call(transformations[[par]][["fun"]], c(list(model_samples[,par]), transformations[[par]][["arg"]]))
runjags_summary[par, "Mean"] <- mean(model_samples[,par], na.rm = TRUE)
runjags_summary[par, "SD"] <- sd(model_samples[,par], na.rm = TRUE)
runjags_summary[par, "Lower95"] <- stats::quantile(model_samples[,par], .025, na.rm = TRUE)
runjags_summary[par, "Upper95"] <- stats::quantile(model_samples[,par], .975, na.rm = TRUE)
runjags_summary[par, "Median"] <- do.call(transformations[[par]][["fun"]], c(list(runjags_summary[par, "Median"]), transformations[[par]][["arg"]]))
runjags_summary[par, "MCerr"] <- do.call(transformations[[par]][["fun"]], c(list(runjags_summary[par, "MCerr"]), transformations[[par]][["arg"]]))
runjags_summary[par, "MC.ofSD"] <- 100 * runjags_summary[par, "MCerr"] / runjags_summary[par, "SD"]
}else if((!transform_factors && (is.prior.orthonormal(prior_list[[par]]) | is.prior.meandif(prior_list[[par]]))) || is.prior.treatment(prior_list[[par]])){
# treatment priors
par_names <- .JAGS_prior_factor_names(par, prior_list[[par]])
for(i in seq_along(par_names)){
model_samples[,par_names[i]] <- do.call(transformations[[par]][["fun"]], c(list(model_samples[,par_names[i]]), transformations[[par]][["arg"]]))
runjags_summary[par_names[i], "Mean"] <- mean(model_samples[,par_names[i]], na.rm = TRUE)
runjags_summary[par_names[i], "SD"] <- sd(model_samples[,par_names[i]], na.rm = TRUE)
runjags_summary[par_names[i], "Lower95"] <- stats::quantile(model_samples[,par_names[i]], .025, na.rm = TRUE)
runjags_summary[par_names[i], "Upper95"] <- stats::quantile(model_samples[,par_names[i]], .975, na.rm = TRUE)
runjags_summary[par_names[i], "Median"] <- do.call(transformations[[par]][["fun"]], c(list(runjags_summary[par_names[i], "Median"]), transformations[[par]][["arg"]]))
runjags_summary[par_names[i], "MCerr"] <- do.call(transformations[[par]][["fun"]], c(list(runjags_summary[par_names[i], "MCerr"]), transformations[[par]][["arg"]]))
runjags_summary[par_names[i], "MC.ofSD"] <- 100 * runjags_summary[par_names[i], "MCerr"] / runjags_summary[par_names[i], "SD"]
}
}
}
}
# transform orthonormal factors to differences from mean
if(transform_factors & any(sapply(prior_list, function(x) is.prior.orthonormal(x) | is.prior.meandif(x)))){
message("The transformation was applied to the differences from the mean. Note that non-linear transformations do not map from the orthonormal/meandif contrasts to the differences from the mean.")
for(par in names(prior_list)[sapply(prior_list, function(x) is.prior.orthonormal(x) | is.prior.meandif(x))]){
par_names <- .JAGS_prior_factor_names(par, prior_list[[par]])
original_samples <- model_samples[,par_names,drop = FALSE]
if(is.prior.orthonormal(prior_list[[par]])){
transformed_samples <- original_samples %*% t(contr.orthonormal(1:(.get_prior_factor_levels(prior_list[[par]])+1)))
}else if(is.prior.meandif(prior_list[[par]])){
transformed_samples <- original_samples %*% t(contr.meandif(1:(.get_prior_factor_levels(prior_list[[par]])+1)))
}
# apply transformation if specified
if(!is.null(transformations[par])){
for(i in 1:ncol(transformed_samples)){
transformed_samples[,i] <- do.call(transformations[[par]][["fun"]], c(list(transformed_samples[,i]), transformations[[par]][["arg"]]))
}
}
if(.is_prior_interaction(prior_list[[par]])){
if(length(.get_prior_factor_level_names(prior_list[[par]])) == 1){
transformed_names <- paste0(par, " [dif: ", .get_prior_factor_level_names(prior_list[[par]])[[1]],"]")
}else{
stop("orthonormal/meandif de-transformation for interaction of multiple factors is not implemented.")
}
}else{
transformed_names <- paste0(par, " [dif: ", .get_prior_factor_level_names(prior_list[[par]]),"]")
}
colnames(transformed_samples) <- transformed_names
# update samples
model_samples <- model_samples[,!colnames(model_samples) %in% par_names,drop=FALSE]
model_samples <- cbind(model_samples, transformed_samples)
# update summary
if(anyNA(transformed_samples)){
# remove NA's introduced by conditional models and spike & slab priors -- also removes the
transformed_chains <- lapply(split(data.frame(transformed_samples), sort(rep(1:length(fit[["mcmc"]]), fit[["sample"]]))), function(x) coda::mcmc(stats::na.omit(x)))
transformed_summary <- summary(runjags::combine.mcmc(transformed_chains, collapse.chains = FALSE))
transformed_summary <- cbind(
Lower95 = transformed_summary$quantiles[,"2.5%"],
Median = transformed_summary$quantiles[,"50%"],
Upper95 = transformed_summary$quantiles[,"97.5%"],
Mean = transformed_summary$statistics[,"Mean"],
SD = transformed_summary$statistics[,"SD"],
Mode = NA,
MCerr = NA,
MC.ofSD = NA,
SSeff = NA,
AC.10 = NA,
psrf = NA
)
}else{
transformed_chains <- lapply(split(data.frame(transformed_samples), sort(rep(1:length(fit[["mcmc"]]), fit[["sample"]]))), coda::mcmc)
transformed_summary <- summary(runjags::combine.mcmc(transformed_chains, collapse.chains = FALSE))
transformed_summary <- cbind(
Lower95 = transformed_summary$quantiles[,"2.5%"],
Median = transformed_summary$quantiles[,"50%"],
Upper95 = transformed_summary$quantiles[,"97.5%"],
Mean = transformed_summary$statistics[,"Mean"],
SD = transformed_summary$statistics[,"SD"],
Mode = NA,
MCerr = if(is.prior.point(prior_list[[par]])) NA else transformed_summary$statistics[,"Naive SE"],
MC.ofSD = if(is.prior.point(prior_list[[par]])) NA else 100 * transformed_summary$statistics[,"Naive SE"] / transformed_summary$statistics[,"SD"],
SSeff = if(is.prior.point(prior_list[[par]])) NA else unname(coda::effectiveSize(coda::as.mcmc(transformed_samples))),
AC.10 = if(is.prior.point(prior_list[[par]])) NA else coda::autocorr.diag(coda::as.mcmc(transformed_samples), lags = 10)[1,],
psrf = if(is.prior.point(prior_list[[par]])) NA else if(length(fit$mcmc) > 1) unname(coda::gelman.diag(transformed_chains, multivariate = FALSE)$psrf[,"Point est."]) else NA
)
}
rownames(transformed_summary) <- transformed_names
par_index <- which.max(rownames(runjags_summary) %in% par_names)
runjags_summary <- runjags_summary[!rownames(runjags_summary) %in% par_names,]
runjags_summary <- rbind(
if(par_index > 1) runjags_summary[1:(par_index-1),],
transformed_summary,
if(par_index <= nrow(runjags_summary)) runjags_summary[par_index:nrow(runjags_summary),]
)
}
}
# remove un-wanted columns
runjags_summary <- runjags_summary[,!colnames(runjags_summary) %in% c("Mode", "AC.10"),drop = FALSE]
# rename treatment factor levels
if(any(sapply(prior_list, is.prior.treatment))){
for(par in names(prior_list)[sapply(prior_list, is.prior.treatment)]){
if(!.is_prior_interaction(prior_list[[par]])){
if(.get_prior_factor_levels(prior_list[[par]]) == 1){
rownames(runjags_summary)[rownames(runjags_summary) == par] <-
paste0(par,"[",.get_prior_factor_level_names(prior_list[[par]])[-1], "]")
}else{
rownames(runjags_summary)[rownames(runjags_summary) %in% paste0(par,"[",1:.get_prior_factor_levels(prior_list[[par]]),"]")] <-
paste0(par,"[",.get_prior_factor_level_names(prior_list[[par]])[-1], "]")
}
}else if(length(attr(prior_list[[par]], "levels")) == 1){
rownames(runjags_summary)[rownames(runjags_summary) %in% paste0(par,"[",1:.get_prior_factor_levels(prior_list[[par]]),"]")] <-
paste0(par,"[",.get_prior_factor_level_names(prior_list[[par]])[[1]][-1], "]")
}
}
}
# rename independent factor levels
if(any(sapply(prior_list, is.prior.independent))){
for(par in names(prior_list)[sapply(prior_list, is.prior.independent)]){
if(!.is_prior_interaction(prior_list[[par]])){
if(.get_prior_factor_levels(prior_list[[par]]) == 1){
rownames(runjags_summary)[rownames(runjags_summary) == par] <-
paste0(par,"[",.get_prior_factor_level_names(prior_list[[par]]), "]")
}else{
rownames(runjags_summary)[rownames(runjags_summary) %in% paste0(par,"[",1:.get_prior_factor_levels(prior_list[[par]]),"]")] <-
paste0(par,"[",.get_prior_factor_level_names(prior_list[[par]]), "]")
}
}else if(length(attr(prior_list[[par]], "levels")) == 1){
rownames(runjags_summary)[rownames(runjags_summary) %in% paste0(par,"[",1:.get_prior_factor_levels(prior_list[[par]]),"]")] <-
paste0(par,"[",.get_prior_factor_level_names(prior_list[[par]])[[1]], "]")
}
}
}
# store parameter names before removing formula attachments
parameter_names <- rownames(runjags_summary)
# rename formula parameters
if(any(!sapply(lapply(prior_list, attr, which = "parameter"), is.null))){
rownames(runjags_summary) <- format_parameter_names(
parameters = rownames(runjags_summary),
formula_parameters = unique(unlist(lapply(prior_list, attr, which = "parameter"))),
formula_prefix = formula_prefix)
}
# rename the rest
colnames(runjags_summary)[colnames(runjags_summary) == "Lower95"] <- "lCI"
colnames(runjags_summary)[colnames(runjags_summary) == "Upper95"] <- "uCI"
colnames(runjags_summary)[colnames(runjags_summary) == "MCerr"] <- "MCMC_error"
colnames(runjags_summary)[colnames(runjags_summary) == "MC.ofSD"] <- "MCMC_SD_error"
colnames(runjags_summary)[colnames(runjags_summary) == "SSeff"] <- "ESS"
colnames(runjags_summary)[colnames(runjags_summary) == "psrf"] <- "R_hat"
# change the SD error to a fraction
runjags_summary[, "MCMC_SD_error"] <- runjags_summary[, "MCMC_SD_error"] / 100
# reorder the columns
runjags_summary <- runjags_summary[,c("Mean", "SD", "lCI", "Median", "uCI", "MCMC_error", "MCMC_SD_error", "ESS", "R_hat"), drop = FALSE]
runjags_summary <- data.frame(runjags_summary)
# prepare output
class(runjags_summary) <- c("BayesTools_table", "BayesTools_runjags_summary", class(runjags_summary))
attr(runjags_summary, "type") <- c(rep("estimate", 5), "MCMC_error", "MCMC_SD_error", "ESS", "R_hat")
attr(runjags_summary, "parameters") <- parameter_names
attr(runjags_summary, "rownames") <- TRUE
attr(runjags_summary, "title") <- title
attr(runjags_summary, "footnotes") <- footnotes
attr(runjags_summary, "warnings") <- warnings
return(runjags_summary)
}
#' @rdname BayesTools_model_tables
runjags_inference_table <- function(fit, title = NULL, footnotes = NULL, warnings = NULL, formula_prefix = TRUE){
# check fits
if(!inherits(fit, "runjags"))
stop("'fit' must be a runjags fit")
if(!inherits(fit, "BayesTools_fit"))
stop("'fit' must be a BayesTools fit")
prior_list <- attr(fit, "prior_list")
check_list(prior_list, "prior_list")
if(!all(sapply(prior_list, is.prior)))
stop("'prior_list' must be a list of priors.")
check_char(title, "title", allow_NULL = TRUE)
check_char(footnotes, "footnotes", check_length = 0, allow_NULL = TRUE)
check_char(warnings, "warnings", check_length = 0, allow_NULL = TRUE)
check_bool(formula_prefix, "formula_prefix")
# return empty table if none of the priors is spike and slab
if(!any(sapply(prior_list, is.prior.spike_and_slab))){
runjags_summary <- runjags_inference_empty_table(title = title, footnotes = footnotes, warnings = warnings)
return(runjags_summary)
}
# extract samples
model_samples <- suppressWarnings(coda::as.mcmc(fit))
runjags_summary <- data.frame(matrix(nrow = 0, ncol = 4))
for(par in names(prior_list)){
if(is.prior.spike_and_slab(prior_list[[par]])){
temp_prior_prob <- mean(prior_list[[par]][["inclusion"]])
temp_post_prob <- mean(model_samples[,paste0(par, "_indicator")])
runjags_summary <- rbind(runjags_summary, data.frame(
Parameter = par,
prior_prob = temp_prior_prob,
post_prob = temp_post_prob,
inclusion_BF = (temp_post_prob / (1-temp_post_prob)) / (temp_prior_prob / (1-temp_prior_prob))
))
}
}
runjags_summary$Parameter <- format_parameter_names(
parameters = runjags_summary$Parameter,
formula_parameters = unique(unlist(lapply(prior_list, attr, which = "parameter"))),
formula_prefix = formula_prefix)
class(runjags_summary) <- c("BayesTools_table", "BayesTools_runjags_summary", class(runjags_summary))
attr(runjags_summary, "type") <- c("string", "prior_prob", "post_prob", "inclusion_BF")
attr(runjags_summary, "rownames") <- FALSE
attr(runjags_summary, "title") <- title
attr(runjags_summary, "footnotes") <- footnotes
attr(runjags_summary, "warnings") <- warnings
return(runjags_summary)
}
#' @rdname BayesTools_model_tables
JAGS_estimates_table <- runjags_estimates_table
#' @rdname BayesTools_model_tables
JAGS_inference_table <- runjags_inference_table
#' @rdname BayesTools_model_tables
JAGS_summary_table <- model_summary_table
#' @rdname BayesTools_model_tables
model_summary_empty_table <- function(model_description = NULL, title = NULL, footnotes = NULL, warnings = NULL){
check_list(model_description, "model_description", allow_NULL = TRUE)
summary_names <- c(
"Model",
if(!is.null(model_description)) names(model_description),
"Prior prob.",
"log(marglik)",
"Post. prob.",
"Inclusion BF")
summary_names <- paste0(summary_names, " ")
empty_table <- data.frame(cbind(
summary_names,
rep("", length(summary_names)),
rep(" ", length(summary_names)),
c("Parameter prior distributions", rep("", length(summary_names) - 1))
))
names(empty_table) <- NULL
# prepare output
class(empty_table) <- c("BayesTools_table", class(empty_table))
attr(empty_table, "type") <- c("string_left", "string", "string", "prior")
attr(empty_table, "rownames") <- FALSE
attr(empty_table, "as.matrix") <- TRUE
attr(empty_table, "title") <- title
attr(empty_table, "footnotes") <- footnotes
attr(empty_table, "warnings") <- warnings
return(empty_table)
}
#' @rdname BayesTools_model_tables
runjags_estimates_empty_table <- function(title = NULL, footnotes = NULL, warnings = NULL){
empty_table <- data.frame(matrix(nrow = 0, ncol = 9))
colnames(empty_table) <- c("Mean", "SD", "lCI", "Median", "uCI", "MCMC_error", "MCMC_SD_error", "ESS", "R_hat")
class(empty_table) <- c("BayesTools_table", "BayesTools_runjags_summary", class(empty_table))
attr(empty_table, "type") <- c(rep("estimate", 5), "MCMC_error", "MCMC_SD_error", "ESS", "R_hat")
attr(empty_table, "rownames") <- FALSE
attr(empty_table, "title") <- title
attr(empty_table, "footnotes") <- footnotes
attr(empty_table, "warnings") <- warnings
return(empty_table)
}
#' @rdname BayesTools_model_tables
runjags_inference_empty_table <- function(title = NULL, footnotes = NULL, warnings = NULL){
empty_table <- data.frame(matrix(nrow = 0, ncol = 4))
colnames(empty_table) <- c("Parameter", "prior_prob", "post_prob", "inclusion_BF")
class(empty_table) <- c("BayesTools_table", "BayesTools_runjags_summary", class(empty_table))
attr(empty_table, "type") <- c("string", "prior_prob", "post_prob", "inclusion_BF")
attr(empty_table, "rownames") <- FALSE
attr(empty_table, "title") <- title
attr(empty_table, "footnotes") <- footnotes
attr(empty_table, "warnings") <- warnings
return(empty_table)
}
#' @rdname BayesTools_model_tables
JAGS_estimates_empty_table <- runjags_estimates_empty_table
#' @rdname BayesTools_model_tables
JAGS_inference_empty_table <- runjags_inference_empty_table
#' @rdname BayesTools_model_tables
stan_estimates_table <- function(fit, transformations = NULL, title = NULL, footnotes = NULL, warnings = NULL){
# this is a simplification of the runjags_estimates_table function for stan
.check_rstan()
# check fits
if(!inherits(fit, "stanfit"))
stop("'fit' must be a rstan fit")
prior_list <- attr(fit, "prior_list")
check_list(prior_list, "prior_list")
if(!all(sapply(prior_list, is.prior)))
stop("'prior_list' must be a list of priors.")
check_list(transformations, "transformations", allow_NULL = TRUE)
if(!is.null(transformations) && any(!sapply(transformations, function(trans)is.function(trans[["fun"]]))))
stop("'transformations' must be list of functions in the 'fun' element.")
check_char(title, "title", allow_NULL = TRUE)
check_char(footnotes, "footnotes", check_length = 0, allow_NULL = TRUE)
check_char(warnings, "warnings", check_length = 0, allow_NULL = TRUE)
# obtain model information
stan_summary <- data.frame(rstan::summary(fit)$summary)
model_samples <- .extract_stan(fit, drop = FALSE)
# remove un-wanted columns
stan_summary <- stan_summary[,!colnames(stan_summary) %in% c("X25.", "X75."), drop = FALSE]
# remove un-wanted rows
stan_summary <- stan_summary[-nrow(stan_summary),, drop = FALSE]
# rename columns to match runjags output
colnames(stan_summary) <- c("Mean", "MCerr", "SD", "Lower95", "Median", "Upper95", "SSeff", "psrf")
# add MC.ofSD estimates
stan_summary[, "MC.ofSD"] <- stan_summary[, "MCerr"] / stan_summary[, "SD"]
# apply transformations
if(!is.null(transformations)){
for(par in names(transformations)){
model_samples[,par] <- do.call(transformations[[par]][["fun"]], c(list(model_samples[,par]), transformations[[par]][["arg"]]))
stan_summary[par, "Mean"] <- mean(model_samples[,par], na.rm = TRUE)
stan_summary[par, "SD"] <- sd(model_samples[,par], na.rm = TRUE)
stan_summary[par, "Median"] <- do.call(transformations[[par]][["fun"]], c(list(stan_summary[par, "Median"]), transformations[[par]][["arg"]]))
stan_summary[par, "MCerr"] <- do.call(transformations[[par]][["fun"]], c(list(stan_summary[par, "MCerr"]), transformations[[par]][["arg"]]))
stan_summary[par, "MC.ofSD"] <- stan_summary[par, "MCerr"] / stan_summary[par, "SD"]
}
}
# rename the rest
colnames(stan_summary)[colnames(stan_summary) == "Lower95"] <- "lCI"
colnames(stan_summary)[colnames(stan_summary) == "Upper95"] <- "uCI"
colnames(stan_summary)[colnames(stan_summary) == "MCerr"] <- "MCMC_error"
colnames(stan_summary)[colnames(stan_summary) == "MC.ofSD"] <- "MCMC_SD_error"
colnames(stan_summary)[colnames(stan_summary) == "SSeff"] <- "ESS"
colnames(stan_summary)[colnames(stan_summary) == "psrf"] <- "R_hat"
# reorder the columns
stan_summary <- stan_summary[,c("Mean", "SD", "lCI", "Median", "uCI", "MCMC_error", "MCMC_SD_error", "ESS", "R_hat"), drop = FALSE]
# store parameter names
parameter_names <- rownames(stan_summary)
# prepare output
class(stan_summary) <- c("BayesTools_table", "BayesTools_stan_summary", class(stan_summary))
attr(stan_summary, "type") <- c(rep("estimate", 5), "MCMC_error", "MCMC_SD_error", "ESS", "R_hat")
attr(stan_summary, "parameters") <- parameter_names
attr(stan_summary, "rownames") <- TRUE
attr(stan_summary, "title") <- title
attr(stan_summary, "footnotes") <- footnotes
attr(stan_summary, "warnings") <- warnings
return(stan_summary)
}
#' @title Print a BayesTools table
#'
#' @param x a BayesTools_values_tables
#' @param ... additional arguments.
#'
#' @return \code{print.BayesTools_table} returns \code{NULL}.
#'
#' @exportS3Method
print.BayesTools_table <- function(x, ...){
# print formatting
for(i in seq_along(attr(x, "type"))){
colnames(x)[i] <- .format_column_names(colnames(x)[i], attr(x, "type")[i], x[,i])
x[,i] <- .format_column(x[,i], attr(x, "type")[i], attr(x, "n_models")[i])
}
# print title
if(!is.null(attr(x, "title"))){
cat(paste0(attr(x, "title"), "\n"))
}
# print the table
print(as.data.frame(x), quote = FALSE, right = TRUE, row.names = attr(x, "rownames"))
# print footnotes
for(i in seq_along(attr(x, "footnotes"))){
cat(paste0(attr(x, "footnotes")[i], "\n"))
}
# print warnings in red
for(i in seq_along(attr(x, "warnings"))){
cat(paste0("\033[0;31m", attr(x, "warnings")[i], "\033[0m\n"))
}
return(invisible())
}
#' @title Format Bayes factor
#'
#' @description Formats Bayes factor
#'
#' @param BF Bayes factor(s)
#' @param logBF log(BF)
#' @param BF01 1/BF
#' @param inclusion whether the Bayes factor is an inclusion BF (for naming purposes)
#'
#' @return \code{format_BF} returns a formatted Bayes factor.
#'
#' @export
format_BF <- function(BF, logBF = FALSE, BF01 = FALSE, inclusion = FALSE){
check_real(BF, "BF", lower = 0, check_length = FALSE)
check_bool(logBF, "logBF")
check_bool(BF01, "BF01")
if(BF01){
BF <- 1/BF
name <- ifelse(inclusion, "Exclusion BF", "1/BF")
}else{
name <- ifelse(inclusion, "Inclusion BF", "BF")
}
if(logBF){
BF <- log(BF)
name <- paste0("log(", name, ")")
}
attr(BF, "name") <- name
attr(BF, "logBF") <- logBF
attr(BF, "BF01") <- BF01
return(BF)
}
#' @title Adds column to BayesTools table
#'
#' @description Adds column to a BayesTools table while not
#' breaking formatting, attributes, etc...
#'
#' @param table BayesTools table
#' @param column_title title of the new column
#' @param column_values values of the new column
#' @param column_position position of the new column (defaults to \code{NULL} which
#' appends the column to the end)
#' @param column_type type of values of the new column table (important for formatting,
#' defaults to \code{NULL} = the function tries to guess numeric / character based on the
#' \code{column_values} but many more specific types are available)
#'
#' @return returns an object of 'BayesTools_table' class.
#'
#' @export
add_column <- function(table, column_title, column_values, column_position = NULL, column_type = NULL){
if(!inherits(table, "BayesTools_table"))
stop("The 'table' must be of class 'BayesTools_table'.")
check_char(column_title, "column_title")
if(!(is.vector(column_values) | is.numeric(column_values)) || length(column_values) != nrow(table))
stop("The 'column_values' must be a vector of the same length as has the table rows.")
check_int(column_position, "column_position", allow_NULL = TRUE, lower = 0, upper = ncol(table) + 1)
.check_table_types(column_type, "column_type", allow_NULL = TRUE)
# fill defaults
if(is.null(column_position)){
column_position <- ncol(table) + 1
}
if(is.null(column_type)){
if(all(.is.wholenumber(column_values))){
column_type <- "integer"
}else if(is.numeric(column_values)){
column_type <- "estimate"
}else if(is.character(column_values)){
column_type <- "string"
}else{
stop("The 'column_type' could not be guessed. Please, supply it manually.")
}
}
# add the column (must be constructed in this was as it can't contain NULL)
to_bind <- list(
if(column_position > 1) table[,1:(column_position - 1)],
column_values,
if(column_position <= ncol(table)) table[,column_position:ncol(table)]
)
new_table <- do.call(cbind, to_bind[!sapply(to_bind, is.null)])
# transfer column names
colnames(new_table) <- c(
if(column_position > 1) colnames(table)[1:(column_position - 1)],
column_title,
if(column_position <= ncol(table)) colnames(table)[column_position:ncol(table)]
)
# copy attributes
class(new_table) <- class(table)
attr(new_table, "type") <- c(
if(column_position > 1) attr(table, "type")[1:(column_position - 1)],
column_type,
if(column_position <= ncol(table)) attr(table, "type")[column_position:ncol(table)]
)
for(a in names(attributes(table))[!names(attributes(table)) %in% c("class", "type", "names")]){
attr(new_table, a) <- attr(table, a)
}
return(new_table)
}
#' @title Removes column to BayesTools table
#'
#' @description Removes column to a BayesTools table while not
#' breaking formatting, attributes, etc...
#'
#' @param table BayesTools table
#' @param column_position position of the to be removed column (defaults to \code{NULL} which
#' removes the last column)
#'
#' @return returns an object of 'BayesTools_table' class.
#'
#' @export
remove_column <- function(table, column_position = NULL){
if(!inherits(table, "BayesTools_table"))
stop("The 'table' must be of class 'BayesTools_table'.")
check_int(column_position, "column_position", allow_NULL = TRUE, lower = 1, upper = ncol(table))
# fill defaults
if(is.null(column_position)){
column_position <- ncol(table)
}
# add the column (must be constructed in this was as it can't contain NULL)
new_table <- table[,-column_position]
# transfer column names
colnames(new_table) <- colnames(table)[-column_position]
# copy attributes
class(new_table) <- class(table)
attr(new_table, "type") <- attr(table, "type")[-column_position]
for(a in names(attributes(table))[!names(attributes(table)) %in% c("class", "type", "names")]){
attr(new_table, a) <- attr(table, a)
}
return(new_table)
}
.format_column <- function(x, type, n_models){
if(is.null(x) || length(x) == 0){
return(x)
}else{
return(switch(
type,
"integer" = round(x),
"prior" = .center_priors(x),
"string_left" = .string_left(x),
"string" = x,
"estimate" = format(round(x, digits = 3), nsmall = 3),
"prior_prob" = format(round(x, digits = 3), nsmall = 3),
"post_prob" = format(round(x, digits = 3), nsmall = 3),
"probability" = format(round(x, digits = 3), nsmall = 3),
"marglik" = format(round(x, digits = 2), nsmall = 2),
"BF" = format(round(x, digits = 3), nsmall = 3),
"inclusion_BF" = format(round(x, digits = 3), nsmall = 3),
"n_models" = paste0(round(x), "/", n_models),
"ESS" = round(x),
"R_hat" = format(round(x, digits = 3), nsmall = 3),
"MCMC_error" = format(round(x, digits = 5), nsmall = 5),
"MCMC_SD_error" = format(round(x, digits = 3), nsmall = 3),
"min_ESS" = round(x),
"max_R_hat" = format(round(x, digits = 3), nsmall = 3),
"max_MCMC_error" = format(round(x, digits = 5), nsmall = 5),
"max_MCMC_SD_error" = format(round(x, digits = 3), nsmall = 3)
))
}
}
.format_column_names <- function(x, type, values){
if(is.null(x) || length(x) == 0){
return(x)
}else{
return(switch(
type,
"integer" = x,
"prior" = .string_center(paste0("Prior ", x), values),
"string_left" = .string_left(x, values),
"string" = x,
"estimate" = x,
"probability" = x,
"prior_prob" = "Prior prob.",
"post_prob" = "Post. prob.",
"marglik" = "log(marglik)",
"BF" = if(is.null(attr(values, "name"))) "BF" else attr(values, "name"),
"inclusion_BF" = if(is.null(attr(values, "name"))) "Inclusion BF" else attr(values, "name"),
"n_models" = "Models",
"ESS" = "ESS",
"R_hat" = "R-hat",
"MCMC_error" = "error(MCMC)",
"MCMC_SD_error" = "error(MCMC)/SD",
"min_ESS" = "min(ESS)",
"max_R_hat" = "max(R-hat)",
"max_MCMC_error" = "max[error(MCMC)]",
"max_MCMC_SD_error" = "max[error(MCMC)/SD]",
))
}
}
.center_priors <- function(x){
if(any(grepl("~", x) | grepl("=", x))){
position_tilda <- regexpr("~", x)
position_equal <- regexpr("=", x)
from_right <- sapply(seq_along(x), function(i){
if(position_tilda[i] != -1){
return(nchar(x[i]) - position_tilda[i])
}else if(position_equal[i] != -1){
return(nchar(x[i]) - position_equal[i])
}else{
return(0)
}
})
add_to_right <- ifelse(from_right == 0, 0, max(from_right) - from_right)
x <- paste0(x, sapply(seq_along(x), function(i)paste0(rep(" ", add_to_right[i]), collapse = "")))
}
return(x)
}
.string_left <- function(x, reference = x){
if(length(x) > 0){
add_to_right <- max(nchar(reference)) - nchar(x)
x <- paste0(x, sapply(seq_along(x), function(i)paste0(rep(" ", add_to_right[i]), collapse = "")))
}
return(x)
}
.string_center <- function(x, reference = x){
if(length(x) > 0){
add_to_sides <- max(nchar(reference)) - nchar(x)
add_to_sides <- ifelse(add_to_sides < 0, 0, add_to_sides)
x <- paste0( sapply(seq_along(x), function(i)paste0(rep(" ", round(add_to_sides[i]/2)), collapse = "")), x, sapply(seq_along(x), function(i)paste0(rep(" ", round(add_to_sides[i]/2)), collapse = "")))
}
return(x)
}
.check_table_types <- function(x, name, allow_NULL = FALSE){
check_char(x, name, allow_values = c(
"integer", "prior", "string_left", "string", "estimate",
"probability", "prior_prob", "post_prob",
"marglik", "BF", "inclusion_BF", "n_models",
"ESS", "R_hat", "MCMC_error", "MCMC_SD_error", "min_ESS", "max_R_hat", "max_MCMC_error", "max_MCMC_SD_error"),
allow_NULL = allow_NULL)
}
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.