# What is posterior predictive check, and how I can do that in R?

I am using Bayesian hierarchical modeling to predict an ordered categorical variable from a metric variable. For example, I want to regress Happiness (in 1-5 ratings) on Money (a metric variable):

Happiness∼log(Dollars)

After estimating posterior distribution using MCMC with RJags, I want to do a posterior predictive check, so I need to model a discrepency between posterior distribution and the data. But what should I choose as the discrepency?

This is the result of a sample from pp.check from jagsUI package which calculates Bayesian $p$-value for posterior predictive checking Bayesian analyses fit in JAGS.

• Have you looked into the Anderson-Darling goodness of fit test? Commented Oct 1, 2015 at 7:36
• @GuðmundurEinarsson Yes, but it is not Bayesian, and it uses $p$-value.
– Ho1
Commented Oct 1, 2015 at 7:52
• I suggest you look into this CV question. Commented Oct 1, 2015 at 7:58

The hierarchical model you describe is a generative model. The model you constructed can be used to generate "fake" data. This is a little different conceptually than using your model to make predictions.

The assumption underlying this concept is that a good model should generate fake data that is similar to the actual data set you used to make your model. A bad model will generate data that is in some way fundamentally or systematically different.

You can assess this visually or by using some metric, such as the pp.check method you tried in JAGS (I am not a JAGS user, so can't comment specifically on how this is implemented).

Procedurally how this works is:

1. You specify your model. In your case it looks like you want to do an ordinal regression. This looks like a similar example. Specifically I refer you to the chapter called "Ordinal Predicted Variable" in this book.

2. You sample and obtain posterior distributions for the parameters in your model. Looking at the figure in the linked example, these parameters are $\beta_0$, $\beta_1$ and $\sigma$.

3. Now draw posterior predictive samples. Over the range of your input (Dollars), draw many samples from the posteriors (or take the samples of your posteriors) of the parameters you estimated, then plug those samples into your model equation, the Happiness ~ log(Dollars) you wrote down.

You should end up with many samples of "Happiness" data at a given log(Dollars). From these samples you could, for instance, compute and plot 90% credible intervals across log(Dollar).

1. Plot actual data (on the y axis: Happiness, on the x axis: log(Dollars)), then overlay the draws and credible intervals of your posterior predictive samples.

Now check visually. Does your 90% credible interval contain 90% of the actual Happiness data points? Are there systematic departures of the true data from your model? Then resort to metrics such as pp.check.

This is one way of performing model validation, there are many others.

• Thanks for the answer, and yes, I have read the book throughly, including the chapter on "ordinal predicted variable". But I couldn't get my answer. First: page 38 of (mlg.eng.cam.ac.uk/zoubin/talks/lect1bayes.pdf) says: Myth: Bayesian Methods = Generative Models. So why do you say a Bayesian method would be only useful for generating data? I can use this model to do predtions. Both point forecast, and posterior distribution estimation. What's the problem here? A good forecasting model usually fits well to data, and this is not a problem.
– Ho1
Commented Oct 19, 2015 at 16:36
• Edited. yes not all Bayesian models are "generative", by the definition given in your link. I didn't say that Bayesian models are only useful for generating data, just that they could be used to do this. This underlies model validation based on the posterior predictive distribution. Commented Oct 19, 2015 at 16:50