Posterior Predictive Checks

Lesson 5 ended with partial pooling cutting the league table's error in half, and with one loose thread deliberately left hanging: every diagnostic you own so far checks the sampling, not the model. R-hat, effective sample size, divergences, the whole Lesson 4 dashboard can glow green while the model itself is quietly wrong about the world. A dashboard tells you the engine ran; it never tells you the map is right.

This lesson builds the tool that checks the map: make the fitted model simulate the data it claims to explain, then compare its fakes to the real thing. A model that cannot re-create its own data has no business forecasting anyone's future.

By the end of this lesson you will be able to:

• Simulate replicated datasets from a fitted model in base R, carrying both the posterior's doubt about the parameter and the day-to-day noise of the data
• Choose test statistics with teeth: features the model was not tuned to reproduce, picked from the decision at stake
• Compute and read a posterior predictive p-value, including why values near 0 and near 1 are both alarms
• Place the check inside the Bayesian workflow loop and state its limits: the data is used twice, so a pass is never proof

Prerequisites: Lessons 1 to 5 of this course (the prior-times-likelihood update, the Beta-Binomial and Normal-Normal conjugate pairs, credible intervals from draws, the Lesson 4 trust dashboard, and the hierarchical model), plus base R functions, replicate() and apply().

Below is the whole lesson in one picture. A model was fitted to sixty days of sales at Asha's store (you will meet the data on the next step), then asked to invent sixty-day runs of its own, two thousand times. The bars show how many zero-sale days each invented run contained; the red line marks the 15 zero-sale days that actually happened. Toggle which model does the inventing, and watch the red line walk out of the crowd, or into it.

Sixty days of terrarium kits

Asha added a terrarium kit to the store at the start of spring: a glass bowl, live moss, three tiny ferns, $39. For sixty days she has logged how many kits sold each day. Here is the log, built right here so every line on this page runs in interactive R. And, exactly as in Lesson 5, simulating the data hands us an answer key real life never shows: the counts truly come from a Poisson process selling 1.5 kits a day on average. Put that face down. The analyst fitting the model never sees it.

One hundred kits in sixty days, about 1.7 a day, and one striking feature: on 15 of the 60 days, not a single kit sold. Those blank days are not trivia. Every kit is assembled with live moss, and a kit that sits unsold needs re-misting every few days, so the pace of zero-sale days decides real work. Asha's planning question is literally a question about zeros: how often does a whole day pass with nothing sold?

Her analyst reaches for the default that has served this course since Lesson 2: a Normal model. Daily counts scatter around some unknown mean; put a prior on the mean and update. The prior is Asha's honest guess from similar launches, about 2 kits a day, give or take 1. For the day-to-day spread, hand the model the sample value, 1.43, as if known, the same simplification Lesson 2 made. The Lesson 2 update then gives the posterior in closed form:

Mean demand 1.68 kits a day, give or take 0.18. The posterior is tight, the interval is sensible, and nothing anywhere will complain. There are no chains to diagnose, the update was closed-form, so Lesson 4's dashboard has nothing to inspect and nothing to flag. Yet this model is about to tell Asha something absurd about her blank days, and no number printed so far can reveal it.