Permalink
Switch branches/tags
Nothing to show
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
73 lines (61 sloc) 2.59 KB
# R function for summarizing MCMC output in a regression-style table
# Johannes Karreth
# I use this function mainly for teaching.
# The function produces a table with medians, SDs, credible intervals, and
# the % of posterior draws with the same sign as the median
# from MCMC output from
# R2jags, rjags, R2WinBUGS, R2OpenBUGS, MCMCpack, rstan, and rstanarm
# Arguments:
# sims: MCMC object generated by R2jags, rjags, R2WinBUGS, R2OpenBUGS, MCMCpack, rstan, rstanarm
# ci: desired level for credible intervals; defaults to 0.95
# pars: character vector of parameters to be printed; defaults to NULL (all parameters are printed)
# Pr: print % of posterior draws with same sign as median; defaults to TRUE
mcmctab <- function(sims, ci = 0.95, pars = NULL, Pr = FALSE)
{
if(class(sims)[1] == "jags" || class(sims)[1] == "rjags"){
sims <- as.matrix(as.mcmc(sims))
}
if(class(sims)[1] == "bugs"){
sims <- sims$sims.matrix
}
if(class(sims)[1] == "mcmc"){
sims <- as.matrix(sims)
}
if(class(sims)[1] == "mcmc.list"){
sims <- as.matrix(sims)
}
if(class(sims)[1] == "stanreg"){
sims <- as.matrix(sims)
}
if(class(sims)[1] == "stanfit"){
sims <- as.matrix(sims)
}
if(is.null(pars) == TRUE){
dat <- sims
}
if(is.null(pars) == FALSE & length(pars) == 1){
dat <- sims[, grepl(x = colnames(sims), pattern = pars)]
}
if(is.null(pars) == FALSE & length(pars) > 1){
dat <- sims[, pars]
}
dat_wide <- t(dat)
mcmctab <- apply(dat_wide, 1,
function(x) c(Median = round(median(x), digits = 3), # Posterior median
SD = round(sd(x), digits = 3), # Posterior SD
Lower = as.numeric(round(quantile(x, probs = c((1 - ci) / 2)), digits = 3)), # Lower CI of posterior
Upper = as.numeric(round(quantile(x, probs = c((1 + ci) / 2)), digits = 3)), # Upper CI of posterior
Pr = round(ifelse(median(x) > 0, length(x[x > 0]) / length(x), length(x[x < 0]) / length(x)), digits = 3) # % of posterior draws with same sign as median
))
if(Pr == FALSE){
mcmctab <- apply(dat_wide, 1,
function(x) c(Median = round(median(x), digits = 3), # Posterior median
SD = round(sd(x), digits = 3), # Posterior SD
Lower = as.numeric(round(quantile(x, probs = c((1 - ci) / 2)), digits = 3)), # Lower CI of posterior
Upper = as.numeric(round(quantile(x, probs = c((1 + ci) / 2)), digits = 3))))
}
# return(t(mcmctab))
out_dat <- data.frame("Variable" = colnames(mcmctab), t(mcmctab),
row.names = NULL)
return(out_dat)
}