The purpose of this tutorial for our open lab is to show you some options to work with and efficiently present output from Bayesian models in article manuscripts: regression tables, regression plots, predicted probabilities, marginal effects, and others. This also reinforces the mechanics of working with MCMC output after fitting a Bayesian model in R (see the previous lab). While there are a number of dedicated R packages to generate similar output,1 this tutorial introduces the BayesPostEst package (Scogin et al. 2019) and aims to give you the tools to produce this output on your own and customize it, with the hope that you can use this in your own applied work.
For consistency with course content, we begin with a model fit using rstanarm, but the tutorial below also works with models fit with rstan, JAGS, or BUGS.
This tutorial should be read in conjunction with Lecture #10 and the two articles on postestimation listed on the syllabus:
Next, we extract the MCMC output from the stanreg object created by stan_glm(). Here, we want to keep only the columns with draws from posterior distributions of our parameters of interest (regression coefficients). Hence, we subset the matrix of MCMC output to contain only the columns with the same names as the variables in the model matrix, i.e. the intercept and the predictors: (Intercept), agegroup(30,45], agegroup(45,65], agegroup(65,96], female, education, income, demsatis, information.
Bayesian variants of the regression model don’t differ much in their interpretation from frequentist models—aside from how you interpret uncertainty around your estimate. Hence, if you want to present a regression table in your paper, a “Bayesian” regression table will look very similar to the “frequentist” regression tables you are all familiar with. For Bayesian models, you want to display the standard deviation of the posterior estimates in lieu of coefficient standard errors Alternatively, you may choose to display the 95% credible interval besides (or below) the posterior mean.
As before, you can use the mcmcTab() function to create a table:
sims: the posterior output from your model. mcmcTab() automatically recognizes posterior distributions that were produced by R2jags, rjags, R2WinBUGS, R2OpenBUGS, MCMCpack, rstan, and rstanarm.
ci: desired level for credible intervals; defaults to c(0.025, 0.975), i.e. a 95% credible interval
pars: character vector of parameters to be printed; defaults to NULL (all parameters are printed)
Pr: print the percent of posterior draws that have the same sign as the median of the posterior distribution; defaults to FALSE
ROPE defaults to NULL. If not NULL, a vector of two values defining the region of practical equivalence (“ROPE”); returns % of posterior draws to the left/right of ROPE. For this quantity to be meaningful, all parameters must be on the same scale (e.g. standardized coefficients or first differences). See Kruschke (2013) and Kruschke (2018) for more on the ROPE.
You can also export a conventional regression table to LaTeX or Word using mcmcReg. See the package website or Lab 2 for more details.
Regression dot plots
Coefficient dot plots (or caterpillar plots) are a more intuitive way to present regression results. Some readers can evaluate the size and significance of relationships more easily.
You have a variety of options to generate these plots. Here are a two, one from the BayesPostEst package, and one “by hand”:
As an alternative, you can show the full posterior distribution of each coefficient estimate. To do so, we first reorganize the MCMC output. Note that I’m using the dplyr::select() function below to remove columns I don’t want to retain for this plot:
For nonlinear models, predicted probabilities are usually an intuitive option to present interpretable results to readers. Two popular approaches plot the predicted probability of the outcome across the simulated range of a predictor while holding all other covariates constant at typical values (the average-case approach) or across all observations (the observed-values approach). To do this using the output of a Bayesian model for binary outcomes, follow the steps below, which are also documented in Lecture 10.2
2 The marginaleffects package also offers a very accessible way to compute most of the quantities and graphs in this section - check it out in action below!
Generate predicted probabilities for a simulated “average” case
The following approach (which you also saw in Lecture 10) follows the practice in King, Tomz, and Wittenberg (2000) and calculates the predicted probability for each value of a focal predictor with all other predictors set to one “typical” value (usually a mean or a median).
The relevant function in the BayesPostEst package is mcmcAveProb; you can find more information about it here.
Code
library("BayesPostEst")
Note: as with every R function, you can type the name of the function into the console to see what the function actually does.
Some of the objects to feed to these arguments are best created before using the function:
Generate average predicted probabilities for all observed cases
The following approach (which you also saw in Lecture 10) follows the points raised by Hanmer and Kalkan (2013) that the previous approach is somewhat limited in exploring the full range of cases. In the words of the authors, “this approach creates a weaker connection between the results and the larger goals of the research enterprise … rather than seeking to understand the effect for the average case, the goal is to obtain an estimate of the average effect in the population.”
BayesPostEst also provides a function for this approach: mcmcObsProb.
Because we already created the input objects above, we can move directly to calling the function. Note that this function can be quite slow as it calculates probabilities for \(n \times k \times z\) cases, where \(n\) is the number of observations in the data, \(k\) is the number of posterior draws, and \(z\) is the number of distinct values of the simulated variable.
Code
dissatis_pp_obs <-mcmcObsProb(modelmatrix = swiss_mm, mcmcout = swiss_mcmc, xcol =8, xrange = dissatis_sim, link ="logit", ci =c(0.025, 0.975))
Comparing the two side by side is useful to adjudicate the sensitivity of the results to either approach:
First differences
You can also easily use the posterior distribution to calculate first differences, i.e. the difference in predicted probabilities for two types of cases - one with low and one with high values in a predictor of interest.
First differences “by hand”
First, here is some step-by-step code; this echoes the walkthrough for predicted probabilities in Lecture 10. I calculate here the first difference in Pr(Voted) between men and women.
You can now plot the posterior distribution of the difference in voting propensity between men and women:
Code
ggplot(data = swiss_fd, aes(x = fd)) +geom_density(fill ="gray", color =NA) +xlab("First difference in voting propensity: Women - Men\n(negative numbers indicate a lower voting propensity for women)") +ylab("") +scale_y_continuous(breaks =FALSE) +geom_vline(xintercept =0, linetype ="dashed") +theme_bw()
Using a function to calculate first differences
To facilitate the calculation and plotting of first differences for all explanatory variables, BayesPostEst also contains a function: mcmcFD. It requires the following arguments, most of which are familiar from above:
modelmatrix: model matrix, including intercept. Create with model.matrix(rstanarm object)
mcmcout: posterior distributions of all coefficients in matrix form - can easily be created from rstan, MCMCpack, R2jags, R2OpenBUGS, etc.
link: type of generalized linear model; a character vector set to “logit” (default) or “probit”
ci: the desired credible interval for the first differences; defaults to c(0.025, 0.975)
percentiles: the percentiles of each explanatory variable to compare; defaults to c(0.25, 0.75) for continuous variables. Binary variables are automatically set to 0 and 1. Note: for explanatory variables that are categorical with few values (e.g., 1, 2, 3, 4) and do not vary much across categories, the default c(0.25, 0.75) may return a first difference estimate of 0. This is because in that case, the 25th and 75th percentile of the variable are identical. You need to expand the range of percentiles in that case, e.g. to c(0.1, 0.9). I will eventually update the function to handle this case better. (And of course, you’re always welcome to fork the repository on Github and rework the function yourself!)
fullsims: should the function return the full distribution of first differences or only the median and credible interval? Defaults to FALSE.
The resulting first differences can then be plotted:
Code
ggplot(data = swiss_fd, aes(x = median_fd, y = VarName)) +geom_segment(aes(x = lower_fd, xend = upper_fd, yend = VarName)) +geom_point() +theme_bw() +geom_vline(xintercept =0, linetype ="dashed") +xlab("Difference in Pr(Voted) between low and high values") +ylab("")
Alternatively, you can plot the full density of the distribution of each first difference. Simply using plot() does this for you once you created a full posterior distribution of first differences with mcmcFD():
Note that this plotting function sets the x-axis to the percentage point change in \(\text{Pr}(y = 1)\), e.g. those in the age group of 45-64 years voted at a 20% higher rate than those in the baseline (18-29).
This function also offers a visualization of the “region of practical equivalence” (ROPE) to visualize what a meaningful first difference would look like. See the package vignette for details.
Postestimation with convenient functions from the modelsummary and marginaleffects packages
You can also speed up many of the operations before by using the marginaleffects and modelsummary packages as long as your Bayesian model was estimated with rstanarm or brms. As a general guide, please read through the excellent documentation of both packages and try out some of the code below!3
3 Please note that not all of the functions below reproduce the exact output generated by the BayesPostEst functions above. That’s because some of the defaults and estimation methods in marginaleffects differ. See the documentation of each package for details.
Regression tables
Use the modelsummary function from the modelsummary package. The documentation is rich and should help you customize this table in almost any way possible!
Turning to the marginaleffects package, please review the documentation carefully – the package is under active development and the functions shown here might change over time.
To plot predicted probabilities across the range of a variable, use plot_predictions. Please refer to the documentation for more details on how thse probabilities are calculated. One of the main advantages of plot_predictions (and its cousin, predictions) is that it allows for specifying values for all other covariates or combinations of covariates, even in the case of multilevel models.
First, I plot a similar quantity to “predictions for the average case”; marginaleffects refers to these as “adjusted predictions at representative values”:
To plot first differences (as defined above), use the comparisons function. Note that I am spelling out the differences of the explanatory variables (unlike mcmcFD above, which automatically chooses the 25th and 75th percentile for numeric variables, or the percentiles supplied by the user, as I did above). For other reasons for slight differences in the numbers below compared to mcmcFD above, please refer to the documentation of the comparisons function.
ggplot(data = swiss_fd2, aes(x = estimate, y =paste(term, contrast, sep =" "))) +geom_segment(aes(x = conf.low, xend = conf.high, yend =paste(term, contrast, sep =" "))) +geom_point() +theme_bw() +geom_vline(xintercept =0, linetype ="dashed") +xlab("Difference in Pr(Voted) between low and high values") +ylab("")
Predicted probabilities after ordered or multinomial logit/probit models
To generate similar figures, please refer to the code in the slides for Lectures 11-13. Here are some examples of useful presentations of predicted probabilities across the range of predictors.
Probability of each outcome in an ordered logistic regression.
Probability of each outcome in an ordered logistic regression.
I can replicate a similar version of the plot above in just one line of marginaleffects code:
Code
plot_predictions(vote_mod, condition =c("econselfworse")) +facet_wrap(~ group) +ylab("Probability of voting") +xlab("Personal economic situation is bad") +scale_x_discrete(labels =c("No", "Yes"))
To add the party names after estimating the model, I use this trick:
Differences in the probability of each choice in a multinomial logistic regression.
In addition to the more laborious code in my slides, the marginaleffects package offers some functions to create similar plots with ease. Please refer to the expanding documentation for examples.
Model fit: ROC curves and other assessments
Regression models for binary and ordinal/categorical outcomes offer a variety of tools for assessing model fit. Many of these measures evaluate how well the model classifies cases by the outcome variable. The BayesPostEst package offers a function to plot precision-recall curves and receiver operating characteristic curves: mcmcRocPrcGen. See the package vignette for details.
The authors of rstan and rstanarm also implemented a popular tool: approximate leave-one-out cross-validation (LOO, LOOIC). See the vignette for the rstanarm package.
References
Hanmer, Michael J., and Kerem Ozan Kalkan. 2013. “Behind the Curve: Clarifying the Best Approach to Calculating Predicted Probabilities and Marginal Effects from Limited Dependent Variable Models.”American Journal of Political Science 57 (1): 263–77. https://doi.org/10.1111/j.1540-5907.2012.00602.x.
King, Gary, Michael Tomz, and Jason Wittenberg. 2000. “Making the Most of Statistical Analyses: Improving Interpretation and Presentation.”American Journal of Political Science 44 (2): 347–61. http://www.jstor.org/stable/2669316.
Kruschke, John K. 2013. “Bayesian Estimation Supersedes the t-Test.”Journal of Experimental Psychology: General 142 (2): 573–603. https://doi.org/10.1037/a0029146.
———. 2018. “Rejecting or Accepting Parameter Values in Bayesian Estimation.”Advances in Methods and Practices in Psychological Science 1 (2): 270–80. https://doi.org/10.1177/2515245918771304.
Mize, Trenton D. 2019. “Best Practices for Estimating, Interpreting, and Presenting Nonlinear Interaction Effects.”Sociological Science 6 (4): 81–117. https://doi.org/10.15195/v6.a4.
Mood, Carina. 2010. “Logistic Regression: Why We Cannot Do What We Think We Can Do, and What We Can Do about It.”European Sociological Review 26 (1): 67–82. https://doi.org/10.1093/esr/jcp006.
Muller, Clemma J, and Richard F MacLehose. 2014. “Estimating predicted probabilities from logistic regression: different methods correspond to different target populations.”International Journal of Epidemiology 43 (3): 962–70. https://doi.org/10.1093/ije/dyu029.
Scogin, Shana, Johannes Karreth, Andreas Beger, and Rob Williams. 2019. “BayesPostEst: An r Package to Generate Postestimation Quantities for Bayesian MCMC Estimation.”Journal of Open Source Software 4 (42): 1722. https://doi.org/10.21105/joss.01722.
I can replicate a similar version of the plot above in just one line of 
In addition to the more laborious code in my slides, the