R/summary-tables.R

In BayesTools: Tools for Bayesian Analyses

#' @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)
}