BayesPostEst contains functions to generate postestimation quantities after estimating Bayesian regression models. The package was inspired by a set of functions written originally for Johannes Karreth’s workshop on Bayesian modeling at the ICPSR Summer program. It has grown to include new functions (see mcmcReg
) and will continue to grow to support Bayesian postestimation. For now, the package focuses mostly on generalized linear regression models for binary outcomes (logistic and probit regression). More details on the package philosophy, its functions, and related packages can be found in Scogin et al. (2019).
To install the latest release on CRAN:
install.packages("BayesPostEst")
The latest development version on GitHub can be installed with:
library("devtools")
install_github("ShanaScogin/BayesPostEst")
Once you have installed the package, you can access it by calling:
library("BayesPostEst")
After the package is loaded, check out the ?BayesPostEst
to see a help file.
Most functions in this package work with posterior distributions of parameters. These distributions need to be converted into a matrix. All functions in the package do this automatically for posterior draws generated by JAGS, BUGS, MCMCpack, rstan, and rstanarm. For posterior draws generated by other tools, users must convert these objects into a matrix, where rows represent iterations and columns represent parameters.
This vignette uses the Cowles
dataset (Cowles and Davis 1987) from the carData package (Fox, Weisberg, and Price 2018).
<- carData::Cowles df
This data frame contains information on 1421 individuals in the following variables:
Before proceeding, we convert the two factor variables sex
and volunteer
into numeric variables. We also means-center and standardize the two continuous variables by dividing each by two standard deviations (Gelman and Hill 2007).
$female <- (as.numeric(df$sex) - 2) * (-1)
df$volunteer <- as.numeric(df$volunteer) - 1
df$extraversion <- (df$extraversion - mean(df$extraversion)) / (2 * sd(df$extraversion))
df$neuroticism <- (df$neuroticism - mean(df$neuroticism)) / (2 * sd(df$neuroticism)) df
We estimate a Bayesian generalized linear model with the inverse logit link function, where
\[ Pr(\text{Volunteering}_i) = \text{logit}^{-1}(\beta_1 + \beta_2 \text{Female}_i + \beta_3 \text{Neuroticism}_i + \beta_4 \text{Extraversion}_i) \]
BayesPostEst
functions accommodate GLM estimates for both logit and probit link functions. The examples proceed with the logit link function. If we had estimated a probit regression, the corresponding argument link
in relevant function calls would need to be set to link = "probit"
. Otherwise, it is set to link = "logit"
by default.
To use BayesPostEst
, we first estimate a Bayesian regression model. This vignette demonstrates five tools for doing so: JAGS (via the R2jags and rjags packages), MCMCpack, and the two Stan interfaces rstan and rstanarm.
First, we prepare the data for JAGS (Plummer 2017). Users need to combine all variables into a list and specify any other elements, like in this case N, the number of observations.
<- as.list(df[, c("volunteer", "female", "neuroticism", "extraversion")])
dl $N <- nrow(df) dl
We then write the JAGS model into the working directory.
<- paste("
mod.jags model {
for (i in 1:N){
volunteer[i] ~ dbern(p[i])
logit(p[i]) <- mu[i]
mu[i] <- b[1] + b[2] * female[i] + b[3] * neuroticism[i] + b[4] * extraversion[i]
}
for(j in 1:4){
b[j] ~ dnorm(0, 0.1)
}
}
")
writeLines(mod.jags, "mod.jags")
We then define the parameters for which we wish to retain posterior distributions and provide starting values.
<- c("b")
params.jags <- list("b" = rep(0, 4))
inits1.jags <- list(inits1.jags, inits1.jags, inits1.jags, inits1.jags) inits.jags
Now, fit the model using the R2jags package.
library("R2jags")
set.seed(123)
<- jags(data = dl, inits = inits.jags,
fit.jags parameters.to.save = params.jags, n.chains = 4, n.iter = 2000,
n.burnin = 1000, model.file = "mod.jags")
#> Compiling model graph
#> Resolving undeclared variables
#> Allocating nodes
#> Graph information:
#> Observed stochastic nodes: 1421
#> Unobserved stochastic nodes: 4
#> Total graph size: 6864
#>
#> Initializing model
The same data and model can be used to fit the model using the rjags package:
library("rjags")
<- jags.model(file = "mod.jags", data = dl, inits = inits.jags,
mod.rjags n.chains = 4, n.adapt = 1000)
#> Compiling model graph
#> Resolving undeclared variables
#> Allocating nodes
#> Graph information:
#> Observed stochastic nodes: 1421
#> Unobserved stochastic nodes: 4
#> Total graph size: 6864
#>
#> Initializing model
<- coda.samples(model = mod.rjags,
fit.rjags variable.names = params.jags,
n.iter = 2000)
We estimate the same model using MCMCpack (Martin, Quinn, and Park 2011).
library("MCMCpack")
<- MCMClogit(volunteer ~ female + neuroticism + extraversion,
fit.MCMCpack data = df, burning = 1000, mcmc = 2000, seed = 123,
b0 = 0, B0 = 0.1)
We write the same model in Stan language.
<- paste("
mod.stan data {
int<lower=0> N;
int<lower=0,upper=1> volunteer[N];
vector[N] female;
vector[N] neuroticism;
vector[N] extraversion;
}
parameters {
vector[4] b;
}
model {
volunteer ~ bernoulli_logit(b[1] + b[2] * female + b[3] * neuroticism + b[4] * extraversion);
for(i in 1:4){
b[i] ~ normal(0, 3);
}
}
")
writeLines(mod.stan, "mod.stan")
We then load rstan (Stan Development Team 2019)…
library("rstan")
rstan_options(auto_write = TRUE)
options(mc.cores = 2)
… and estimate the model, re-using the data in list format created for JAGS earlier.
<- stan(file = "mod.stan",
fit.stan data = dl,
pars = c("b"),
chains = 4,
iter = 2000,
seed = 123)
Lastly, we use the rstanarm interface (Goodrich et al. 2019) to estimate the same model again.
library("rstanarm")
<- stan_glm(volunteer ~ female + neuroticism + extraversion,
fit.rstanarm data = df, family = binomial(link = "logit"),
prior = normal(0, 3),
prior_intercept = normal(0, 3),
chains = 4,
iter = 2000,
seed = 123)
#>
#> SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 1).
#> Chain 1:
#> Chain 1: Gradient evaluation took 0.000113 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 1.13 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1:
#> Chain 1:
#> Chain 1: Iteration: 1 / 2000 [ 0%] (Warmup)
#> Chain 1: Iteration: 200 / 2000 [ 10%] (Warmup)
#> Chain 1: Iteration: 400 / 2000 [ 20%] (Warmup)
#> Chain 1: Iteration: 600 / 2000 [ 30%] (Warmup)
#> Chain 1: Iteration: 800 / 2000 [ 40%] (Warmup)
#> Chain 1: Iteration: 1000 / 2000 [ 50%] (Warmup)
#> Chain 1: Iteration: 1001 / 2000 [ 50%] (Sampling)
#> Chain 1: Iteration: 1200 / 2000 [ 60%] (Sampling)
#> Chain 1: Iteration: 1400 / 2000 [ 70%] (Sampling)
#> Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
#> Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
#> Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
#> Chain 1:
#> Chain 1: Elapsed Time: 0.392322 seconds (Warm-up)
#> Chain 1: 0.419418 seconds (Sampling)
#> Chain 1: 0.81174 seconds (Total)
#> Chain 1:
#>
#> SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 2).
#> Chain 2:
#> Chain 2: Gradient evaluation took 5.8e-05 seconds
#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.58 seconds.
#> Chain 2: Adjust your expectations accordingly!
#> Chain 2:
#> Chain 2:
#> Chain 2: Iteration: 1 / 2000 [ 0%] (Warmup)
#> Chain 2: Iteration: 200 / 2000 [ 10%] (Warmup)
#> Chain 2: Iteration: 400 / 2000 [ 20%] (Warmup)
#> Chain 2: Iteration: 600 / 2000 [ 30%] (Warmup)
#> Chain 2: Iteration: 800 / 2000 [ 40%] (Warmup)
#> Chain 2: Iteration: 1000 / 2000 [ 50%] (Warmup)
#> Chain 2: Iteration: 1001 / 2000 [ 50%] (Sampling)
#> Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
#> Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
#> Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
#> Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
#> Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
#> Chain 2:
#> Chain 2: Elapsed Time: 0.388864 seconds (Warm-up)
#> Chain 2: 0.409005 seconds (Sampling)
#> Chain 2: 0.797869 seconds (Total)
#> Chain 2:
#>
#> SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 3).
#> Chain 3:
#> Chain 3: Gradient evaluation took 5.8e-05 seconds
#> Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.58 seconds.
#> Chain 3: Adjust your expectations accordingly!
#> Chain 3:
#> Chain 3:
#> Chain 3: Iteration: 1 / 2000 [ 0%] (Warmup)
#> Chain 3: Iteration: 200 / 2000 [ 10%] (Warmup)
#> Chain 3: Iteration: 400 / 2000 [ 20%] (Warmup)
#> Chain 3: Iteration: 600 / 2000 [ 30%] (Warmup)
#> Chain 3: Iteration: 800 / 2000 [ 40%] (Warmup)
#> Chain 3: Iteration: 1000 / 2000 [ 50%] (Warmup)
#> Chain 3: Iteration: 1001 / 2000 [ 50%] (Sampling)
#> Chain 3: Iteration: 1200 / 2000 [ 60%] (Sampling)
#> Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
#> Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
#> Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
#> Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
#> Chain 3:
#> Chain 3: Elapsed Time: 0.39014 seconds (Warm-up)
#> Chain 3: 0.412336 seconds (Sampling)
#> Chain 3: 0.802476 seconds (Total)
#> Chain 3:
#>
#> SAMPLING FOR MODEL 'bernoulli' NOW (CHAIN 4).
#> Chain 4:
#> Chain 4: Gradient evaluation took 5.7e-05 seconds
#> Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.57 seconds.
#> Chain 4: Adjust your expectations accordingly!
#> Chain 4:
#> Chain 4:
#> Chain 4: Iteration: 1 / 2000 [ 0%] (Warmup)
#> Chain 4: Iteration: 200 / 2000 [ 10%] (Warmup)
#> Chain 4: Iteration: 400 / 2000 [ 20%] (Warmup)
#> Chain 4: Iteration: 600 / 2000 [ 30%] (Warmup)
#> Chain 4: Iteration: 800 / 2000 [ 40%] (Warmup)
#> Chain 4: Iteration: 1000 / 2000 [ 50%] (Warmup)
#> Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
#> Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
#> Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
#> Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
#> Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
#> Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
#> Chain 4:
#> Chain 4: Elapsed Time: 0.390436 seconds (Warm-up)
#> Chain 4: 0.406683 seconds (Sampling)
#> Chain 4: 0.797119 seconds (Total)
#> Chain 4:
BayesPostEst contains functions to generate regression tables from objects created by the following packages: R2jags, runjags, rjags, R2WinBUGS, MCMCpack, rstan, rstanarm, and brms. This includes the following object classes: jags
, rjags
, bugs
, mcmc
, mcmc.list
, stanreg
, stanfit
, brmsfit
. The package contains two different functions to produce regression tables:
mcmcTab
mcmcReg
Each has its own advantages which we discuss in depth below.
mcmcTab
generates a table summarizing the posterior distributions of all parameters contained in the model object. This table can then be used to summarize parameter quantities. By default, mcmcTab
generates a dataframe with one row per parameter and columns containing the median, standard deviation, and 95% credible interval of each parameter’s posterior distribution.
mcmcTab(fit.jags)
#> Variable Median SD Lower Upper
#> 1 b[1] -0.462 0.082 -0.622 -0.298
#> 2 b[2] 0.238 0.112 0.017 0.455
#> 3 b[3] 0.063 0.113 -0.158 0.283
#> 4 b[4] 0.515 0.112 0.291 0.729
#> 5 deviance 1909.455 2.875 1906.574 1917.541
mcmcTab(fit.rjags)
#> Variable Median SD Lower Upper
#> 1 b[1] -0.458 0.082 -0.623 -0.298
#> 2 b[2] 0.234 0.110 0.022 0.454
#> 3 b[3] 0.062 0.111 -0.154 0.287
#> 4 b[4] 0.519 0.111 0.306 0.739
mcmcTab(fit.MCMCpack)
#> Variable Median SD Lower Upper
#> 1 (Intercept) -0.463 0.083 -0.612 -0.304
#> 2 female 0.231 0.110 0.016 0.435
#> 3 neuroticism 0.058 0.112 -0.147 0.321
#> 4 extraversion 0.509 0.102 0.320 0.718
mcmcTab(fit.stan)
mcmcTab(fit.rstanarm)
#> Variable Median SD Lower Upper
#> 1 (Intercept) -0.458 0.081 -0.623 -0.306
#> 2 female 0.237 0.111 0.023 0.462
#> 3 neuroticism 0.055 0.108 -0.144 0.276
#> 4 extraversion 0.518 0.112 0.296 0.733
Users can add a column to the table that calculates the percent of posterior draws that have the same sign as the median of the posterior distribution.
mcmcTab(fit.jags, Pr = TRUE)
#> Variable Median SD Lower Upper Pr
#> 1 b[1] -0.462 0.082 -0.622 -0.298 1.000
#> 2 b[2] 0.238 0.112 0.017 0.455 0.984
#> 3 b[3] 0.063 0.113 -0.158 0.283 0.712
#> 4 b[4] 0.515 0.112 0.291 0.729 1.000
#> 5 deviance 1909.455 2.875 1906.574 1917.541 1.000
Users can also define a “region of practical equivalence” (ROPE; J. K. Kruschke (2013); J. K. Kruschke (2018)). This region is a band of values around 0 that are “practically equivalent” to 0 or no effect. For this to be useful, all parameters (e.g. regression coefficients) must be on the same scale because mcmcTab accepts only one definition of ROPE for all parameters. Users can standardize regression coefficients to achieve this. Because we standardized variables earlier, the coefficients (except the intercept) are on a similar scale and we define the ROPE to be between -0.1 and 0.1.
mcmcTab(fit.jags, pars = c("b[2]", "b[3]", "b[4]"), ROPE = c(-0.1, 0.1))
#> Variable Median SD Lower Upper PrOutROPE
#> 1 b[2] 0.238 0.112 0.017 0.455 0.885
#> 2 b[3] 0.063 0.113 -0.158 0.283 0.376
#> 3 b[4] 0.515 0.112 0.291 0.729 1.000
The mcmcReg
function serves as an interface to texreg
and produces more polished and publication-ready tables than mcmcTab
. mcmcReg
writes tables in HTML or LaTeX format. mcmcReg
can produce tables with multiple models with each model in a column and supports flexible renaming of parameters. However, these tables are more similar to standard frequentist regression tables, so they do not have a way to incorporate the percent of posterior draws that have the same sign as the median of the posterior distribution or a ROPE like mcmcTab
is able to. Uncertainty intervals can be either standard credible intervals or highest posterior density intervals (J. Kruschke 2015) using the hpdi
argument, and their level can be set with the ci
argument (default 95%). Separately calculated goodness of fit statistics can be included with the gof
argument.
mcmcReg(fit.jags, format = 'html', doctype = F)
Model 1 | |
---|---|
b[1] | -0.46^{*} |
[ -0.62; -0.30] | |
b[2] | 0.24^{*} |
[ 0.02; 0.46] | |
b[3] | 0.06 |
[ -0.16; 0.28] | |
b[4] | 0.51^{*} |
[ 0.29; 0.73] | |
deviance | 1910.11^{*} |
[1906.57; 1917.54] | |
^{*} 0 outside 95% credible interval. |
mcmcReg
supports limiting the parameters included in the table via the pars
argument. By default, all parameters saved in the model object will be included. In the case of fit.jags
, this include the deviance estimate. If we wish to exclude it, we can specify pars = 'b'
which will capture b[1]
-b[4]
using regular expression matching.
mcmcReg(fit.jags, pars = 'b', format = 'html', regex = T, doctype = F)
Model 1 | |
---|---|
b[1] | -0.46^{*} |
[-0.62; -0.30] | |
b[2] | 0.24^{*} |
[ 0.02; 0.46] | |
b[3] | 0.06 |
[-0.16; 0.28] | |
b[4] | 0.51^{*} |
[ 0.29; 0.73] | |
^{*} 0 outside 95% credible interval. |
If we only wish to exclude the intercept, we can do this by explicitly specifying the parameters we wish to include as a vector. Note that in this example we have to escape the []
s in pars
because they are a reserved character in regular expressions.
mcmcReg(fit.jags, pars = c('b\\[1\\]', 'b\\[3\\]', 'b\\[4\\]'),
format = 'html', regex = T, doctype = F)
Model 1 | |
---|---|
b[1] | -0.46^{*} |
[-0.62; -0.30] | |
b[3] | 0.06 |
[-0.16; 0.28] | |
b[4] | 0.51^{*} |
[ 0.29; 0.73] | |
^{*} 0 outside 95% credible interval. |
mcmcReg
also supports partial regular expression matching of multiple parameter family names as demonstrated below.
mcmcReg(fit.jags, pars = c('b', 'dev'),
format = 'html', regex = T, doctype = F)
Model 1 | |
---|---|
b[1] | -0.46^{*} |
[ -0.62; -0.30] | |
b[2] | 0.24^{*} |
[ 0.02; 0.46] | |
b[3] | 0.06 |
[ -0.16; 0.28] | |
b[4] | 0.51^{*} |
[ 0.29; 0.73] | |
deviance | 1910.11^{*} |
[1906.57; 1917.54] | |
^{*} 0 outside 95% credible interval. |
mcmcReg
supports custom coefficient names to support publication-ready tables. The simplest option is via the coefnames
argument. Note that the number of parameters and the number of custom coefficient names must match, so it is a good idea to use pars
in tandem with coefnames
.
mcmcReg(fit.jags, pars = 'b',
coefnames = c('(Constant)', 'Female', 'Neuroticism', 'Extraversion'),
format = 'html', regex = T, doctype = F)
Model 1 | |
---|---|
(Constant) | -0.46^{*} |
[-0.62; -0.30] | |
Female | 0.24^{*} |
[ 0.02; 0.46] | |
Neuroticism | 0.06 |
[-0.16; 0.28] | |
Extraversion | 0.51^{*} |
[ 0.29; 0.73] | |
^{*} 0 outside 95% credible interval. |
A more flexible way to include custom coefficient names is via the custom.coef.map
argument, which accepts a named list, with names as parameter names in the model and values as the custom coefficient names.
mcmcReg(fit.jags, pars = 'b',
custom.coef.map = list('b[1]' = '(Constant)',
'b[2]' = 'Female',
'b[3]' = 'Nueroticism',
'b[4]' = 'Extraversion'),
format = 'html', regex = T, doctype = F)
Model 1 | |
---|---|
(Constant) | -0.46^{*} |
[-0.62; -0.30] | |
Female | 0.24^{*} |
[ 0.02; 0.46] | |
Nueroticism | 0.06 |
[-0.16; 0.28] | |
Extraversion | 0.51^{*} |
[ 0.29; 0.73] | |
^{*} 0 outside 95% credible interval. |
The advantage of custom.coef.map
is that it can flexibly reorder and omit coefficients from the table based on their positions within the list. Notice in the code below that deviance does not have to be included in pars
because its absence from custom.coef.map
omits it from the resulting table.
mcmcReg(fit.jags,
custom.coef.map = list('b[2]' = 'Female',
'b[4]' = 'Extraversion',
'b[1]' = '(Constant)'),
format = 'html', doctype = F)
Model 1 | |
---|---|
Female | 0.24^{*} |
[ 0.02; 0.46] | |
Extraversion | 0.51^{*} |
[ 0.29; 0.73] | |
(Constant) | -0.46^{*} |
[-0.62; -0.30] | |
^{*} 0 outside 95% credible interval. |
However, it is important to remember that mcmcReg
will look for the parameter names in the model object, so be sure to inspect it for the correct parameter names. This is important because stan_glm
will produce a model object with variable names instead of indexed parameter names.
mcmcReg
accepts multiple model objects and will produce a table with one model per column. To produce a table from multiple models, pass a list of models as the mod
argument to mcmcReg
.
mcmcReg(list(fit.stan, fit.stan), format = 'html', doctype = F)
Note, however, that all model objects must be of the same class, so it is not possible to generate a table from a jags
object and a stanfit
object.
mcmcReg(list(fit.jags, fit.stan), format = 'html', doctype = F)
When including multiple models, supplying scalars or vectors to arguments will result in them being applied to each model equally. Treating models differentially is possible by supplying a list of scalars or vectors instead.
mcmcReg(list(fit.rstanarm, fit.rstanarm),
pars = list(c('female', 'extraversion'), 'neuroticism'),
format = 'html', doctype = F)
Model 1 | Model 2 | |
---|---|---|
female | 0.24^{*} | |
[0.02; 0.46] | ||
extraversion | 0.52^{*} | |
[0.30; 0.73] | ||
neuroticism | 0.06 | |
[-0.14; 0.28] | ||
^{*} Null hypothesis value outside 95% credible interval. |
texreg
argumentsAlthough custom.coef.map
is not an argument to mcmcReg
, it works because mcmcReg
supports all standard texreg
arguments (a few have been overridden, but they are explicit arguments to mcmcReg
). This introduces a high level of control over the output of mcmcReg
, as e.g. models can be renamed.
mcmcReg(fit.rstanarm, custom.model.names = 'Binary Outcome',
format = 'html', doctype = F)
Binary Outcome | |
---|---|
(Intercept) | -0.46^{*} |
[-0.62; -0.31] | |
female | 0.24^{*} |
[ 0.02; 0.46] | |
neuroticism | 0.06 |
[-0.14; 0.28] | |
extraversion | 0.52^{*} |
[ 0.30; 0.73] | |
^{*} 0 outside 95% credible interval. |
mcmcAveProb
To evaluate the relationship between covariates and a binary outcome, this function calculates the predicted probability (\(Pr(y = 1)\)) at pre-defined values of one covariate of interest (\(x\)), while all other covariates are held at a “typical” value. This follows suggestions outlined in King, Tomz, and Wittenberg (2000) and elsewhere, which are commonly adopted by users of GLMs. The mcmcAveProb
function by default calculates the median value of all covariates other than \(x\) as “typical” values.
Before moving on, we show how create a matrix of posterior draws of coefficients to pass onto these functions. Eventually, each function will contain code similar to the first section of mcmcTab
to do this as part of the function.
<- as.matrix(coda::as.mcmc(fit.jags))
mcmcmat.jags <- as.matrix(fit.MCMCpack)
mcmcmat.MCMCpack <- as.matrix(fit.rstanarm) mcmcmat.rstanarm
<- as.matrix(fit.stan) mcmcmat.stan
Next, we generate the model matrix to pass on to the function. A model matrix contains as many columns as estimated regression coefficients. The first column is a vector of 1s (corresponding to the intercept); the remaining columns are the observed values of covariates in the model. Note: the order of columns in the model matrix must correspond to the order of columns in the matrix of posterior draws.
<- model.matrix(volunteer ~ female + neuroticism + extraversion,
mm data = df)
We can now generate predicted probabilities for different values of a covariate of interest.
First, we generate full posterior distributions of the predicted probability of volunteering for a typical female and a typical male. In this function and mcmcObsProb
, users specify the range of \(x\) (here 0 and 1) as well as the number of the column of \(x\) in the matrix of posterior draws as well as the model matrix.
<- mcmcAveProb(modelmatrix = mm,
aveprob.female.jags mcmcout = mcmcmat.jags[, 1:ncol(mm)],
xcol = 2,
xrange = c(0, 1),
link = "logit",
ci = c(0.025, 0.975),
fullsims = TRUE)
Users can then visualize this posterior distribution using the ggplot2 and ggridges
packages.
library("ggplot2")
library("ggridges")
ggplot(data = aveprob.female.jags,
aes(y = factor(x), x = pp)) +
stat_density_ridges(quantile_lines = TRUE,
quantiles = c(0.025, 0.5, 0.975), vline_color = "white") +
scale_y_discrete(labels = c("Male", "Female")) +
ylab("") +
xlab("Estimated probability of volunteering") +
labs(title = "Probability based on average-case approach") +
theme_minimal()
For continuous variables of interest, users may want to set fullsims = FALSE
to obtain the median predicted probability along the range of \(x\) as well as a lower and upper bound of choice (here, the 95% credible interval).
<- mcmcAveProb(modelmatrix = mm,
aveprob.extra.jags mcmcout = mcmcmat.jags[, 1:ncol(mm)],
xcol = 4,
xrange = seq(min(df$extraversion), max(df$extraversion), length.out = 20),
link = "logit",
ci = c(0.025, 0.975),
fullsims = FALSE)
Users can then plot the resulting probabilities using any plotting functions, such as ggplot2.
ggplot(data = aveprob.extra.jags,
aes(x = x, y = median_pp)) +
geom_ribbon(aes(ymin = lower_pp, ymax = upper_pp), fill = "gray") +
geom_line() +
xlab("Extraversion") +
ylab("Estimated probability of volunteering") +
ylim(0, 1) +
labs(title = "Probability based on average-case approach") +
theme_minimal()
mcmcObsProb
As an alternative to probabilities for “typical” cases, Hanmer and Kalkan (2013) suggest to calculate predicted probabilities for all observed cases and then derive an “average effect.” In their words, the goal of this postestimation “is to obtain an estimate of the average effect in the population … rather than seeking to understand the effect for the average case.”
We first calculate the average “effect” of sex on volunteering, again generating a full posterior distribution. Again, xcol
represents the position of the covariate of interest, and xrange
specifies the values for which \(Pr(y = 1)\) is to be calculated.
<- mcmcObsProb(modelmatrix = mm,
obsprob.female.jags mcmcout = mcmcmat.jags[, 1:ncol(mm)],
xcol = 2,
xrange = c(0, 1),
link = "logit",
ci = c(0.025, 0.975),
fullsims = TRUE)
Users can again plot the resulting densities.
ggplot(data = obsprob.female.jags,
aes(y = factor(x), x = pp)) +
stat_density_ridges(quantile_lines = TRUE,
quantiles = c(0.025, 0.5, 0.975), vline_color = "white") +
scale_y_discrete(labels = c("Male", "Female")) +
ylab("") +
xlab("Estimated probability of volunteering") +
labs(title = "Probability based on observed-case approach") +
theme_minimal()
For this continuous predictor, we use fullsims = FALSE
.
<- mcmcObsProb(modelmatrix = mm,
obsprob.extra.jags mcmcout = mcmcmat.jags[, 1:ncol(mm)],
xcol = 4,
xrange = seq(min(df$extraversion), max(df$extraversion), length.out = 20),
link = "logit",
ci = c(0.025, 0.975),
fullsims = FALSE)
We then plot the resulting probabilities across observed cases.
ggplot(data = obsprob.extra.jags,
aes(x = x, y = median_pp)) +
geom_ribbon(aes(ymin = lower_pp, ymax = upper_pp), fill = "gray") +
geom_line() +
xlab("Extraversion") +
ylab("Estimated probability of volunteering") +
ylim(0, 1) +
labs(title = "Probability based on observed-case approach") +
theme_minimal()
mcmcFD
To summarize typical effects across covariates, we generate “first differences” (Long (1997), King, Tomz, and Wittenberg (2000)). This quantity represents, for each covariate, the difference in predicted probabilities for cases with low and high values of the respective covariate. For each of these differences, all other variables are held constant at their median.
<- mcmcFD(modelmatrix = mm,
fdfull.jags mcmcout = mcmcmat.jags[, 1:ncol(mm)],
link = "logit",
ci = c(0.025, 0.975),
fullsims = TRUE)
summary(fdfull.jags)
#> female neuroticism extraversion
#> Min. :-0.04432 Min. :-0.055522 Min. :0.01957
#> 1st Qu.: 0.03972 1st Qu.:-0.002247 1st Qu.:0.06957
#> Median : 0.05795 Median : 0.011224 Median :0.08132
#> Mean : 0.05779 Mean : 0.011263 Mean :0.08124
#> 3rd Qu.: 0.07619 3rd Qu.: 0.024511 3rd Qu.:0.09313
#> Max. : 0.15245 Max. : 0.080707 Max. :0.14269
The posterior distribution can be summarized as above, or users can directly obtain a summary when setting fullsims
to FALSE.
<- mcmcFD(modelmatrix = mm,
fdsum.jags mcmcout = mcmcmat.jags[, 1:ncol(mm)],
link = "logit",
ci = c(0.025, 0.975),
fullsims = FALSE)
fdsum.jags#> median_fd lower_fd upper_fd VarName VarID
#> female 0.05794760 0.004299115 0.11073472 female 1
#> neuroticism 0.01122442 -0.028097810 0.05008092 neuroticism 2
#> extraversion 0.08132383 0.045716391 0.11482003 extraversion 3
Users can plot the median and credible intervals of the summary of the first differences.
ggplot(data = fdsum.jags,
aes(x = median_fd, y = VarName)) +
geom_point() +
geom_segment(aes(x = lower_fd, xend = upper_fd, yend = VarName)) +
geom_vline(xintercept = 0) +
xlab("Change in Pr(Volunteering)") +
ylab("") +
theme_minimal()
mcmcFD
objectsTo make use of the full posterior distribution of first differences, we provide a dedicated plotting method, plot.mcmcFD
, which returns a ggplot2 object that can be further customized. The function is modeled after Figure 1 in Karreth (2018). Users can specify a region of practical equivalence and print the percent of posterior draws to the right or left of the ROPE. If ROPE is not specified, the figure automatically prints the percent of posterior draws to the left or right of 0.
plot(fdfull.jags, ROPE = c(-0.01, 0.01))
The user can further customize the plot.
<- plot(fdfull.jags, ROPE = c(-0.01, 0.01))
p + labs(title = "First differences") +
p ::theme_ridges() ggridges
mcmcRocPrc
One way to assess model fit is to calculate the area under the Receiver Operating Characteristic (ROC) and Precision-Recall curves. A short description of these curves and their utility for model assessment is provided in Beger (2016). The mcmcRocPrc
function produces an object with four elements: the area under the ROC curve, the area under the PR curve, and two dataframes to plot each curve. When fullsims
is set to FALSE
, the elements represent the median of the posterior distribution of each quantity.
mcmcRocPrc
currently requires an “rjags” object (a model fitted in R2jags) as input. Future package versions will generalize this input to allow for model objects fit with any of the other packages used in BayesPostEst.
<- mcmcRocPrc(object = fit.jags,
fitstats yname = "volunteer",
xnames = c("female", "neuroticism", "extraversion"),
curves = TRUE,
fullsims = FALSE)
Users can then print the area under the each curve:
$area_under_roc
fitstats#> V1
#> 0.5840611
$area_under_prc
fitstats#> V1
#> 0.4867844
Users can also plot the ROC curve…
ggplot(data = as.data.frame(fitstats, what = "roc"), aes(x = x, y = y)) +
geom_line() +
geom_abline(intercept = 0, slope = 1, color = "gray") +
labs(title = "ROC curve") +
xlab("1 - Specificity") +
ylab("Sensitivity") +
theme_minimal()
… as well as the precision-recall curve.
ggplot(data = as.data.frame(fitstats, what = "prc"), aes(x = x, y = y)) +
geom_line() +
labs(title = "Precision-Recall curve") +
xlab("Recall") +
ylab("Precision") +
theme_minimal()
To plot the posterior distribution of the area under the curves, users set the fullsims
argument to TRUE
. Unless a user wishes to plot credible intervals around the ROC and PR curves themselves, we recommend keeping curves
at FALSE
to avoid long computation time.
<- mcmcRocPrc(object = fit.jags,
fitstats.fullsims yname = "volunteer",
xnames = c("female", "neuroticism", "extraversion"),
curves = FALSE,
fullsims = TRUE)
We can then plot the posterior density of the area under each curve.
ggplot(as.data.frame(fitstats.fullsims),
aes(x = area_under_roc)) +
geom_density() +
labs(title = "Area under the ROC curve") +
theme_minimal()
ggplot(as.data.frame(fitstats.fullsims),
aes(x = area_under_prc)) +
geom_density() +
labs(title = "Area under the Precision-Recall curve") +
theme_minimal()