1
#' Summarize Bayesian MCMC Output
2
#'  
3
#' R function for summarizing MCMC output in a regression-style table.
4
#'  
5
#' @param sims Bayesian model object generated by R2jags, rjags, R2WinBUGS, 
6
#'   R2OpenBUGS, MCMCpack, rstan, and rstanarm.
7
#' @param ci desired level for credible intervals; defaults to c(0.025, 0.975).
8
#' @param pars character vector of parameters to be printed; defaults to \code{NULL} 
9
#'   (all parameters are printed). If not \code{NULL}, the user can either specify the exact names of 
10
#'   parameters to be printed (e.g. \code{c("alpha", "beta1", "beta2")}) or part of a name 
11
#'   so that all parameters containing that name will be printed (e.g. \code{"beta"} will print \code{beta1}, \code{beta2}, etc.).
12
#' @param Pr print percent of posterior draws with same sign as median; defaults to \code{FALSE}.
13
#' @param ROPE defaults to \code{NULL}. If not \code{NULL}, a vector of two values defining the region of 
14
#'   practical equivalence ("ROPE"); returns \% of posterior draws to the left/right of ROPE. For this quantity 
15
#'   to be meaningful, all parameters must be on the same scale (e.g. standardized coefficients 
16
#'   or first differences). See Kruschke (2013, Journal of Experimental 
17
#'   Psychology 143(2): 573-603) for more on the ROPE.
18
#' @param regex use regular expression matching with \code{pars}?
19
#'
20
#' @references Kruschke, John K. 2013. “Bayesian Estimation Supersedes the T-Test.” Journal of 
21
#'   Experimental Psychology: General 142 (2): 573–603. https://doi.org/10.1037/a0029146.
22
#'
23
#' @return a data frame containing MCMC summary statistics.
24
#'
25
#' @examples
26
#' \dontshow{.old_wd <- setwd(tempdir())}
27
#' \donttest{
28
#' data("jags_logit")
29
#' 
30
#' ## printing out table
31
#' object <- mcmcTab(jags_logit, 
32
#'           ci = c(0.025, 0.975), 
33
#'           pars = NULL, 
34
#'           Pr = FALSE,
35
#'           ROPE = NULL)
36
#' object
37
#' }
38
#' 
39
#' \dontshow{setwd(.old_wd)}
40
#' @export
41

42
mcmcTab <- function(sims, 
43
                    ci = c(0.025, 0.975), 
44
                    pars = NULL, 
45
                    Pr = FALSE,
46
                    ROPE = NULL,
47
                    regex = FALSE) {
48
  
49 1
  if (inherits(sims, what = c("jags", "rjags"))) {
50 1
    sims <- as.matrix(coda::as.mcmc(sims))
51
  }
52 1
  if (inherits(sims, what = "bugs")) {
53 1
    sims <- sims$sims.matrix
54
  }
55 1
  if (inherits(sims, what = c("mcmc", "mcmc.list", "stanfit", "stanreg",
56 1
                              "brmsfit"))) {
57 1
    sims <- as.matrix(sims)
58
  }
59
  
60 1
  ROPE <- check_ROPE_argument(ROPE)
61
  
62 1
  if (is.null(pars)) {
63 1
    dat <- sims
64 1
  } else if (regex) {
65 1
    dat <- sims[, grepl(x = colnames(sims), pattern = paste(pars, collapse = "|"))]
66
  } else {
67 1
    dat <- matrix(sims[, pars], nrow = nrow(sims), byrow = FALSE,
68 1
                  dimnames = list(NULL, pars))
69
  }
70
  
71 1
  dat_wide <- t(dat)
72
  
73 1
  mcmctab <- apply(dat_wide, 1, 
74 1
                   function(x) c(Median = round(median(x), digits = 3), # Posterior median
75 1
                                 SD = round(sd(x), digits = 3), # Posterior SD
76 1
                                 Lower = as.numeric(round(quantile(x, probs = ci[1]), digits = 3)), # Lower CI of posterior
77 1
                                 Upper = as.numeric(round(quantile(x, probs = ci[2]), digits = 3)), # Upper CI of posterior
78 1
                                 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
79
                   ))
80
  
81 1
  if(Pr == FALSE){
82 1
    mcmctab <- apply(dat_wide, 1, 
83 1
                     function(x) c(Median = round(median(x), digits = 3), # Posterior median
84 1
                                   SD = round(sd(x), digits = 3), # Posterior SD
85 1
                                   Lower = as.numeric(round(quantile(x, probs = ci[1]), digits = 3)), # Lower CI of posterior
86 1
                                   Upper = as.numeric(round(quantile(x, probs = ci[2]), digits = 3))))
87
  }
88
  
89 1
  if(!is.null(ROPE)){
90 1
    message("This table contains an estimate for parameter values outside of the region of 
91 1
          practical equivalence (ROPE). For this quantity to be meaningful, all parameters 
92 1
          must be on the same scale (e.g. standardized coefficients or first differences).")
93
    
94 1
    mcmctab <- apply(dat_wide, 1, 
95 1
                     function(x) c(Median = round(median(x), digits = 3), # Posterior median
96 1
                                   SD = round(sd(x), digits = 3), # Posterior SD
97 1
                                   Lower = as.numeric(round(quantile(x, probs = ci[1]), digits = 3)), # Lower CI of posterior
98 1
                                   Upper = as.numeric(round(quantile(x, probs = ci[2]), digits = 3)),
99 1
                                   PrOutROPE = round(ifelse(median(x) > 0, length(x[x > ROPE[2]]) / length(x), length(x[x < ROPE[1]]) / length(x)), digits = 3)))
100
  }
101
  
102
  # return(t(mcmctab))
103 1
  out_dat <- data.frame("Variable" = colnames(mcmctab), 
104 1
                        t(mcmctab),
105 1
                        row.names = NULL,
106 1
                        stringsAsFactors = TRUE) # check this, new with R 4.0.0
107
                                                 # recommended if sort order used 
108
                                                 # in the string to factor conversion
109
                                                 # does not matter
110
  
111 1
  return(out_dat)
112
  
113
} 

Read our documentation on viewing source code .

Loading