R/mcmc.summary.R

In dosreislab/bppr: An R package for BPP

#' MCMC summary
#'
#' Create an MCMC summary from a BPP or MCMCTree analysis
#'
#' @param mcmc Data frame with the MCMC output of BPP or MCMCTree
#' @param prob Numeric, probability for credibility interval calculation
#'
#' @details \code{mcmc} should contain the output (say from file
#'   \code{mcmc.txt}) generated by a BPP A00 or MCMCtree analysis. The function
#'   will calculate the posterior (prior) means, the equal-tail credibility
#'   interval (CI), and the highest posterior (prior) density (HPD) CI. \code{prob} is
#'   used to calculate the CIs. For example, if \code{prob} = 95%, then you get the
#'   95% CIs. This function requires the CODA pacakge to calculate the HPD CIs.
#'
#' @return A list with elements \code{means}, \code{eq.ci}, and \code{hpd.ci}
#'   containing the posterior (or prior) means, equal tail CI and HPD CI.
#'
#' @examples
#' \dontrun{
#' mcmc.summary(hominids$mcmc[,-1])
#' }
#'
#' @author Mario dos Reis
#'
#' @export
mcmc.summary <- function(mcmc, prob=0.95) {
  prob <- 1 - prob
  p.means <- apply(mcmc, 2, mean) # means
  ci <- t( apply(mcmc, 2, quantile, probs=c(prob/2, 1 - prob/2)) )
  hpd <- coda::HPDinterval(coda::as.mcmc(mcmc))

  return(list(means=p.means, eq.ci=ci, hpd.ci=hpd))
}