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.
x
andy
. JAGS will automatically sample from the PP distribution. Set the monitoring onx[11]
andy[11]
(for a sample size of 10) to get the PP distribution forx
andy
. Then compare those with the actual values ofx
andy
. That would be a posterior predictive check in this case, if I'm not mistaken. $\endgroup$