# How do I perform an actual "posterior predictive check"?

This question is the follow-up of this previous question: Bayesian inference and testable implications.

For concreteness, consider the following bayesian model. This model is not to be taken literally, it is simply suppose to stand for a model that cannot capture the DGP but we do not know that a priori. However, I very much would like an answer that takes this concrete model to perform an actual posterior predictive check, so we avoid generic answers.

$$\text{Likelihood:}\\ \\ y \sim \mathcal{N}(\mu_1, \sigma_1)\\ x \sim \mathcal{N}(\mu_2, \sigma_2)\\ \text{Prior:}\\ \\ \mu_1 \sim \mathcal{N}(0, 1000)\\ a \sim \mathcal{U(0,2)}\\ \mu_2 \leftarrow \mu_1 + a\\ \sigma_1 \sim \mathcal{U}(0, 100)\\ \sigma_2 \sim \mathcal{U}(0, 100)$$

Where $$\mathcal{N}()$$ denotes a gaussian and $$\mathcal{U}()$$ denotes a uniform distribution. Here is the implementation in rjags:

library(rjags)
model <- "
model {
for (i in 1:length(x)){
x[i] ~ dnorm(mu1, tau1)
}

for (i in 1:length(y)){
y[i] ~ dnorm(mu2, tau2)
}

mu1 ~ dnorm(0, .00001)
a ~ dunif(0, 2)
mu2 <- mu1 + a

sigma1 ~ dunif(0,100)
tau1 <- pow(sigma1, -2)

sigma2 ~ dunif(0,100)
tau2 <- pow(sigma2, -2)
}
"


And here is the model fitted to some simulated data that does not conform to the model's assumptions.

n <- 10
dat <- list(x = rnorm(n, mean = 2, sd = 2),
y = rnorm(n, mean = 10, sd = 10))

jags.model   <- jags.model(textConnection(model), data =dat)
#> Compiling model graph
#>    Resolving undeclared variables
#>    Allocating nodes
#> Graph information:
#>    Observed stochastic nodes: 20
#>    Unobserved stochastic nodes: 4
#>    Total graph size: 32
#>
#> Initializing model
samp <- coda.samples(jags.model, n.iter = 1e4,
variable.names = c("mu1", "mu2", "sigma1", "sigma2"))
post  <- as.data.frame(samp[[1]])
summary(post$$mu1) #> Min. 1st Qu. Median Mean 3rd Qu. Max. #> -1.732 1.456 1.977 2.004 2.526 6.897 summary(post$$mu2)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
#> -0.9573  2.4011  3.0740  3.0808  3.7376  8.2234


Now, how do I formally perform a "posterior predictive check" in this model with this data? And how do I formally decide, using the posterior predictive check, that the model misfit is "bad enough" so that I "reject" this model? What "test statistic" would you use? Which "threshold" for decision would you use? And so on. If there are missing details that are required for solving this problem (like, say, you need a cost or loss function) please feel free to add those details in your answer as needed; these details are part of a good answer, since they clarify what we need to know to actually perform the check.

Finally, please try to provide an actual solution to this toy problem. It doesn't need to be code, if you can derive the numerical results by hand that works as well. But the main idea is to have this toy problem actually solved.

• A simple way to sample from the posterior predictive is to include a missing value in x and y. JAGS will automatically sample from the PP distribution. Set the monitoring on x[11] and y[11] (for a sample size of 10) to get the PP distribution for x and y. Then compare those with the actual values of x and y. That would be a posterior predictive check in this case, if I'm not mistaken. May 24 '20 at 7:45
• I think this is a reasonable question and don't quite understand the downvotes. But the request for an implementation is off-topic here, and I'd recommend you remove it.
– mkt
May 26 '20 at 14:51
• @mkt-ReinstateMonica think of it as just a small cost to avoid those people who are tempted to give generic answers like "there are several ways to do it, you could do it like this, or like that." If the person knows how to do posterior predictive checks, it should be trivial to do it in this example. No?
– user272422
May 26 '20 at 14:53
• @mkt-ReinstateMonica thanks I just reworded the question, hope it is a bit better.
– user272422
May 26 '20 at 15:03
• Too lazy to construct an actual answer, but have you consulted Gelman's Bayesian Data Analysis? The PDF is available for free online, and chapter 6.3 is devoted to posterior predictive checks with some examples partially worked out May 26 '20 at 15:43