R Academy · Lesson

Posterior Predictive Checks

Validate model fit by comparing simulated vs observed data distributions.

Posterior Predictive Checks is a free R Academy lesson on CoddyKit — lesson 4 of 4. You can read the complete lesson below for free — then practise it hands-on in the browser with a built-in code editor and a 24/7 AI tutor. It is part of the R Academy learning path, one of 4 lessons in the course, and your progress syncs across the web and the CoddyKit app.

What Are Posterior Predictive Checks?

After fitting a Bayesian model you need to ask: does this model generate data that looks like the data you observed? Posterior Predictive Checks (PPCs) answer this by simulating replicated datasets yrep from the posterior and comparing them graphically to the observed y.

Extracting Posterior Samples

extract(fit, pars='mu') returns a named list; $mu is a numeric vector of all post-warmup samples for that parameter. With 4 chains x 1000 post-warmup iterations you get 4000 samples.

# library(rstan)
# mu_samples    <- extract(fit, pars = 'mu')$mu
# sigma_samples <- extract(fit, pars = 'sigma')$sigma
#
# cat('Samples drawn:', length(mu_samples), '
')
# cat('Posterior mean of mu:', mean(mu_samples), '
')
# cat('90% CI:', quantile(mu_samples, c(0.05, 0.95)), '
')

Generating yrep from Posterior Samples

For each posterior draw (mu_s, sigma_s), simulate a replicated dataset of the same size as the original data. Store these in a matrix yrep where each row is one simulated dataset.

# y <- c(2.1, 1.8, 2.4, 1.9, 2.3, 2.0, 1.7, 2.2)
# n_obs <- length(y)
# S     <- length(mu_samples)  # 4000 posterior draws
#
# yrep <- matrix(NA, nrow = S, ncol = n_obs)
# for (s in seq_len(S)) {
#   yrep[s, ] <- rnorm(n_obs, mean = mu_samples[s], sd = sigma_samples[s])
# }
# dim(yrep)  # [4000, 8]

ppc_dens_overlay() — Density Comparison

bayesplot::ppc_dens_overlay(y, yrep[1:50,]) overlays the kernel density of the observed data (dark line) with densities from 50 randomly chosen simulated datasets (light lines). Good model fit means the dark line sits inside the cloud of light lines.

# library(bayesplot)
#
# ppc_dens_overlay(y, yrep[1:50, ])
#
# Interpretation:
# - Dark line (y_obs) surrounded by light lines (yrep): good fit
# - Dark line systematically outside the cloud: model misfit
# - Light lines much wider than dark: overdispersed model
# - Light lines much narrower than dark: underdispersed model

ppc_stat() — Test Statistic Check

ppc_stat(y, yrep, stat = 'mean') shows a histogram of the test statistic (e.g., mean) computed on each simulated dataset, with a vertical line at the observed statistic. If the observed value falls in the bulk of the histogram, the model captures that aspect of the data.

# library(bayesplot)
#
# ppc_stat(y, yrep, stat = 'mean')   # does model capture the mean?
# ppc_stat(y, yrep, stat = 'sd')     # does model capture spread?
# ppc_stat(y, yrep, stat = 'max')    # does model capture extremes?
#
# If observed stat is in the tail of the histogram,
# the model fails to reproduce that statistic.

Bayesian p-value

The Bayesian p-value (posterior predictive p-value) is the proportion of simulated datasets whose test statistic is more extreme than the observed value. Values near 0.5 indicate good calibration; values near 0 or 1 indicate model misfit for that statistic.

# Bayesian p-value for the mean:
# obs_mean <- mean(y)
# rep_means <- apply(yrep, 1, mean)
# pval <- mean(rep_means >= obs_mean)
# cat('Bayesian p-value (mean):', round(pval, 3), '
')
# # 0.5 is perfect; < 0.05 or > 0.95 suggests misfit

More bayesplot PPC Functions

bayesplot offers many PPC visualizations beyond density overlays:

ppc_hist(y, yrep[1:8,]) — histogram grid
ppc_scatter_avg(y, yrep) — observed vs mean of yrep scatter
ppc_intervals(y, yrep) — uncertainty intervals around each observation
ppc_rootogram(y, yrep) — for count data

# library(bayesplot)
#
# # Grid of 8 simulated histograms vs the observed
# ppc_hist(y, yrep[1:8, ])
#
# # Scatter: y_obs (x) vs mean of yrep (y) — should hug diagonal
# ppc_scatter_avg(y, yrep)
#
# # 50% and 90% posterior predictive intervals around each y_i
# ppc_intervals(y, yrep)

Interpreting PPC Plots — Model Misfit

Common misfit patterns and their causes:

yrep too wide — prior is too diffuse or model is overdispersed
yrep shifted — wrong likelihood family (e.g., normal for skewed data)
yrep misses multimodality — mixture model needed
yrep fails for extreme values — heavy-tailed distribution needed

# Example: if data has a long right tail but yrep does not,
# consider switching:
# y ~ normal(mu, sigma)  =>  y ~ student_t(nu, mu, sigma)
#
# Or for count data:
# y ~ poisson(lambda)  =>  y ~ neg_binomial_2(mu, phi)  (overdispersion)
cat('PPCs guide model improvement by revealing specific failure modes
')

PPCs in the Stan Model Block

You can generate yrep directly in Stan using the generated quantities block. This avoids re-extracting parameters in R and is computationally equivalent.

# Stan model with generated quantities:
# '
# generated quantities {
#   array[N] real y_rep;
#   for (n in 1:N) {
#     y_rep[n] = normal_rng(mu, sigma);
#   }
# }
# '
# Then extract in R:
# yrep <- extract(fit, pars = 'y_rep')$y_rep  # [S, N] matrix

Leave-One-Out Cross-Validation

Beyond PPCs, loo::loo(fit) computes leave-one-out cross-validation to compare competing models. The model with higher ELPD (expected log predictive density) is preferred. Use loo::loo_compare(loo1, loo2) to rank models.

# library(loo)
# loo1 <- loo(fit1)  # normal model
# loo2 <- loo(fit2)  # student-t model
#
# comparison <- loo_compare(loo1, loo2)
# print(comparison)
#
# Model with elpd_diff > 0 is preferred
# se_diff > |elpd_diff| means difference is not reliable

PPC Best Practices

Follow these practices for rigorous posterior predictive checking:

Always start with ppc_dens_overlay() as a global sanity check
Follow up with domain-specific statistics (ppc_stat()) relevant to your analysis goals
Use at least 50 yrep draws for visual checks; all 4000 for Bayesian p-values
Failed PPCs guide model improvement — they are diagnostic, not failure

# Workflow:
# 1. Fit model -> extract() -> generate yrep matrix
# 2. ppc_dens_overlay(y, yrep[1:50,])  -- visual global check
# 3. ppc_stat(y, yrep, stat='mean')    -- check mean
# 4. ppc_stat(y, yrep, stat='sd')      -- check spread
# 5. ppc_stat(y, yrep, stat='max')     -- check tails
# 6. If misfit found -> revise model -> refit -> re-check

Quick Check: Bayesian p-value Interpretation

A Bayesian p-value of 0.03 for the maximum statistic means what about the model?

Posterior Predictive Checks Recap

PPC workflow in RStan and bayesplot:

Extract samples: extract(fit, pars='mu')$mu
Generate yrep: loop over posterior draws calling rnorm(n, mu_s, sigma_s)
Global check: ppc_dens_overlay(y, yrep[1:50,])
Statistic checks: ppc_stat(y, yrep, stat='mean')
Bayesian p-value: mean(apply(yrep,1,stat) >= stat(y)) — near 0.5 is good
For model comparison use loo::loo_compare()

Frequently asked questions

Is the “Posterior Predictive Checks” lesson free?

Yes — the full text of “Posterior Predictive Checks” is free to read here on the web, and the R Academy course includes 4 lessons in total. To practise it interactively (a built-in code editor and a 24/7 AI tutor) and unlock the rest of the R Academy course, upgrade to CoddyKit PRO.

What will I learn in “Posterior Predictive Checks”?

Validate model fit by comparing simulated vs observed data distributions. You practise R Academy with hands-on code you run directly in the browser, and a 24/7 AI tutor answers your questions as you work through the lesson.

Do I need any experience to start R Academy?

No prior experience is required. R Academy on CoddyKit is structured for beginners through advanced learners; this is — lesson 4 of 4, so you can start here or from the beginning and move at your own pace.

How long does the “Posterior Predictive Checks” lesson take?

Most CoddyKit lessons take about 5–10 minutes. Each one is bite-sized and interactive, so you make steady progress and pick up exactly where you left off across the web and the app.

Can I write and run code in this R Academy lesson?

Yes. Every R Academy lesson includes a built-in code editor, so you write and run real code right in your browser and get instant AI feedback — no local setup required.

All lessons in this course

← Back to R Academy