Check your prior posterior overlap (PPO) – MCMC wrangling in R made easy with `MCMCvis`

When fitting a Bayesian model using MCMC (often via JAGS/BUGS/Stan), a number of checks are typically performed to make sure your model is worth interpreting without further manipulation (remember: all models are wrong, some are useful!):

  • R-hat (AKA Gelman-Rubin statistic) – used to assess convergence of chains in the model
  • Visual assessment of chains – used to assess whether posterior chains mixed well (convergence)
  • Visual assessment of posterior distribution shape – used to determine if the posterior distribution is constrained
  • Posterior predictive check (predicting data using estimated parameters) – used to make sure that the model can generate the data used in the model

PPO

One check, however, is often missing: a robust assessment of the degree to which the prior is informing the posterior distribution. Substantial influence of the prior on the posterior may not be apparent through the use of R-hat and visual checks alone. Version 0.9.2 of MCMCvis (now available on CRAN), makes quantifying and plotting the prior posterior overlap (PPO) simple.

MCMCvis is an R package designed to streamline analysis of Bayesian model results derived from MCMC samplers (e.g., JAGS, BUGS, Stan). It can be used to easily visualize, manipulate, and summarize MCMC output. The newest version is full of new features – a full tutorial can be found here.

An example

To check PPO for a model, we will use the function MCMCtrace. As the function is used to generate trace and density plots, checking for PPO is barely more work than just doing the routine checks that one would ordinarily perform. The function plots trace plots on the left and density plots for both the posterior (black) and prior (red) distributions on the right. The function calculates the percent overlap between the prior and posterior and prints this value on the plot. See ?MCMCtrace in R for details regarding the syntax.

#install package
install.packages('MCMCvis', repos = "http://cran.case.edu")

#load package
require(MCMCvis)

#load example data
data(MCMC_data)

#simulate data from the prior used in your model
#number of iterations should equal the number of draws times the number of chains (although the function will adjust if the correct number of iterations is not specified)
#in JAGS: parameter ~ dnorm(0, 0.001)
PR <- rnorm(15000, 0, 32)

#run the function for just beta parameters
MCMCtrace(MCMC_data, params = 'beta', priors = PR, pdf = FALSE)

unnamed-chunk-1-1

Why check?

Checking the PPO has particular utility when trying to determine if the parameters in your model are identifiable. If substantial PPO exists, the prior may simply be dictating the posterior distribution – the data may have little influence on the results. If a small degree of PPO exists, the data was informative enough to overcome the influence of the prior. In the field of ecology, nonidentifiability is a particular concern in some types of mark-recapture models. Gimenez (2009) developed quantitative guidelines to determine when parameters are robustly identifiable using PPO.

While a large degree of PPO is not always a bad thing (e.g., substantial prior knowledge about the system may result in very informative priors used in the model), it is important to know where data was and was not informative for parameter estimation. The degree of PPO that is acceptable for a particular model will depend on a great number of factors, and may be somewhat subjective (but see Gimenez [2009] for a less subjective case). Like other checks, PPO is just one of many tools to be used for model assessment. Finding substantial PPO when unexpected may suggest that further model manipulation is needed. Happy model building!

Other MCMCvis improvements

Check out the rest of the new package freatures, including the option to calculate the number of effective samples for each parameter, ability to take arguments in the form of a ‘regular expression’ for the params argument, ability to retain the structure of all parameters in model output (e.g., parameters specified as matrices in the model are summarized as matrices).

Follow Casey Youngflesh on Twitter @caseyyoungflesh. The MCMCvis source code can be found on GitHub.

Visualizing and wrangling MCMC output in R with `MCMCvis`

Model results can be thought of as a reward for the many hours of model design, troubleshooting, re-design, etc. that analyses often require. Following the potentially exhausting mental exercise to acquire these results, I think we’d all like the interpretation to be as straightforward as possible. Analyzing MCMC output from Bayesian analyses, which may include hundreds of parameters and/or derived quantities, however can often require a fair amount of code and (more importantly) time.

The MCMCvis package was designed to alleviate this problem, and streamline the analysis of Bayesian model results. The latest version (0.7.1) is now available on CRAN with bug fixes, speed improvements, and added functionality.

 

Why MCMCvis?

Using MCMCvis provides three principal benefits:

1) MCMC output fit with Stan, JAGS, or other MCMC samplers can be fed into all package functions as an argument with no further manipulation needed. No need to specify the type of object or how it was fit; the package does all of that behind the scenes.

2) Specific parameters or derived quantities of interest can be specified within each function, to avoid additional steps of data processing. This works using a grep like call for optimal efficiency.

3) The package creates ‘publication-ready’ posterior estimate visualizations (below). Parameters can now be plotted vertically or horizontally.

rplot

The package has four functions for basic MCMC output tasks:

MCMCsummary – summarize MCMC output for particular parameters of interest

MCMCtrace – create trace and density plots of MCMC chains for particular parameters of interest

MCMCchains – easily extract posterior chains from MCMC output for particular parameters of interest

MCMCplot – create caterpillar plots from MCMC output for particular parameters of interest

The vignette can be found here.

 

An example workflow may go as follows:

– summarize posterior estimates for just beta parameters

#install package
install.packages('MCMCvis')

#load package
require(MCMCvis)

#load example data
data(MCMC_data)

#run summary function
MCMCsummary(MCMC_data, 
             params = 'beta')
##            mean   2.5%    50% 97.5% Rhat
## beta[1]    0.16   0.06   0.15  0.25    1
## beta[2]   -7.77 -25.82  -7.68  9.78    1
## beta[3]   -5.64 -28.53  -5.76 17.23    1
## beta[4]  -10.39 -25.98 -10.63  5.27    1
## beta[5]    7.52   6.03   7.52  9.05    1
## beta[6]   10.89  10.10  10.89 11.68    1
## beta[7]   -1.91  -4.83  -1.92  1.08    1
## beta[8]    5.38  -6.86   5.45 17.67    1
## beta[9]   13.39   3.28  13.38 23.60    1
## beta[10]  17.63  14.41  17.63 20.86    1

– check posteriors for convergence

MCMCtrace(MCMC_data, 
           params = c('beta[1]', 'beta[2]', 'beta[3]'), 
           ind = TRUE)

rplot01

– extract chains for beta parameters so that they can be manipulated directly

ex <- MCMCchains(MCMC_data, params = 'beta')
#find 22nd quantile for all beta parameters
apply(ex, 2, function(x){round(quantile(x, probs = 0.22), digits = 2)})
##  beta[1]  beta[2]  beta[3]  beta[4]  beta[5]  beta[6]  beta[7]  beta[8]  beta[9] beta[10] 
##    0.12   -14.86   -14.80   -16.48     6.91    10.58    -3.09     0.68     9.29    16.36

– create caterpillar plots for posterior estimates. Shading represents whether 50% CI (gray with open circle), 95% CI (gray with closed circle), or neither (black) overlap 0. This option can be turned off (as below). A variety of options exist, including the ability to plot posteriors vertically rather than horizontally

MCMCplot(MCMC_data, 
          params = 'beta', 
          horiz = FALSE,
          rank = TRUE,
          ref_ovl = FALSE,
          xlab = 'My x-axis label', 
          main = 'MCMCvis plot', 
          labels = c('First param', 'Second param', 'Third param', 
          'Fourth param', 'Fifth param', 'Sixth param', 'Seventh param', 
          'Eighth param', 'Nineth param', 'Tenth param'), 
          labels_sz = 1.5, med_sz = 2, thick_sz = 7, thin_sz = 3, 
          ax_sz = 4, main_text_sz = 2)

 rplot03

Follow Casey Youngflesh on Twitter @caseyyoungflesh. The MCMCvis source code can be found on GitHub.

R package `MCMCvis` now on CRAN

MCMCvis, an R package for visualizing, manipulating, and summarizing MCMC output is now available on CRAN!

It can be installed with: install.packages('MCMCvis')

The vignette (tutorial) can be run with: vignette('MCMCvis')

MCMCvis was designed to perform key functions for MCMC analysis using minimal code, in order to free up time/brainpower for interpretation of analysis results. Functions support simple and straightforward subsetting of model parameters within the calls, and produce presentable and ‘publication-ready’ output. Model output can be from JAGS, Stan, or other MCMC samplers. The package deals with different data types behind the scenes so you don’t have to think about it. You can specify which parameters you want to visualize or extract within the functions.

For example, this plot can be made using just one line of code (after loading package/data, of course)!

library('MCMCvis') #load package

data(MCMC_data) #load example data

MCMCplot(MCMC_data, params = 'beta') #create plot

 

rplot02

More information can be found within the package vignette. The package was created and authored by Lynch Lab PhD candidate, Casey Youngflesh. He can be found here on Github and here on Twitter.