#' @title Create a summary table of an mcmc.list object
#'
#' @description
#' \code{coda_table} takes as input an mcmc.list object from the coda
#' package and outputs a dataframe giving the mean, standard
#' deviation, and median for each parameter estimate along with quantiles
#' for 95\% and 80\% credibility intervals.
#'
#' @param coda.object
#' an mcmc.list object
#'
#' @param parameters
#' character vector of parameter names to include in table. If none are
#' supplied all monitored parameters are included.
#'
#' @param gelman.diag
#' logical, should Rhat be computed for each parameter, if TRUE input
#' `coda.object` should have more than one chain.
#'
#' @param quantiles
#' vector of probability quantiles to compute
#'
#' @return
#' A dataframe
#'
#' @seealso \code{\link{coda_df}}
#'
#' @export
#'
#' @author Michael Malick
#'
#' @examples
#' library(coda)
#' data(line)
#'
#' coda_table(line)
#'
#' coda_table(line, quantiles = c(0.05, 0.95))
#'
#' coda_table(line, parameters = "alpha")
#'
#' coda_table(line, parameters = c("alpha", "beta"))
#'
#' coda_table(line, parameters = grep("sig", coda::varnames(line), value = TRUE))
#'
#' coda_table(line, parameters = grep("a", coda::varnames(line), value = TRUE))
#'
#' coda_table(line, parameters = c("alpha", grep("sig", coda::varnames(line),
#' value = TRUE)))
#'
coda_table <- function(coda.object,
parameters = NULL,
gelman.diag = TRUE,
quantiles = c(0.025, 0.1, 0.9, 0.975)) {
sim.dat <- coda_df(coda.object, parameters = parameters)
sim.dat$chain <- NULL
sim.dat$iter <- NULL
parms <- names(sim.dat)
n.parms <- length(parms)
ll <- list(mean = apply(as.matrix(sim.dat), 2, mean),
sd = apply(as.matrix(sim.dat), 2, stats::sd),
median = apply(as.matrix(sim.dat), 2, stats::median),
p.conf = apply(as.matrix(sim.dat), 2, stats::quantile,
probs = quantiles),
n.iter = rep(coda::nchain(coda.object) * coda::niter(coda.object),
n.parms),
n.eff = coda::effectiveSize(coda.object)[parms])
if(gelman.diag) {
ll$Rhat <- as.vector(coda::gelman.diag(coda.object)$psrf[,1][parms])
}
out <- do.call("rbind", ll)
out <- as.data.frame(t(round(out, 3)))
return(out)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.