mcmcRocPrc.R is quite long
andybega
| 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 .