# What does it mean to compare deviance of a model based on observed values to the deviance for predicted values (deviance posterior predictive check)?

For an intro to Bayesian statistics course we were asked to compute a "Bayesian P-value" (what I understand is a posterior predictive check) by comparing the model deviance based on observed values to the deviance of model based on posterior predictions. The former I believe is "regular deviance." The latter though has me confused, because deviance already seems to me a measure of how well our model fits the data.

Posterior predictive checks of data summaries such as the grand mean or standard deviation make sense to me, but I wasn't able to find any resources about how this specific posterior predictive check works, or what it means.

Here is a reproducible example in Jags using the jagsUI package for how we are calculating our posterior predictive checks. As mentioned, I'm confused about the one calculated for deviance.

library(jagsUI)
set.seed(5454)

x1 <- rgamma(100, shape = 2, rate = 2)
x2 <- runif(100, min = 2, max = 8)
mu <- -3 + -2 * x1 + 0.75 * x2
Y <- rnorm(100, mean = mu, sd = 1)

jags_mod <- tempfile()

writeLines("
model {

# Priors
b0 ~ dnorm(0, .0001)
b1 ~ dnorm(0, .0001)
b2 ~ dnorm(0, .0001)

SD ~ dunif(0, 100)
ep.sig <- 1 / (SD^2)

# Likelihood:
for (i in 1:n){
Y[i] ~ dnorm(mu[i], ep.sig)
mu[i] <- b0 + b1 * x1[i] + b2 * x2[i]

# For GOF/ Bayesian p-value
Y_sim[i] ~ dnorm(mu[i], ep.sig)

# Here is where we calculate deviance for observed data then
# also for simulated (predicted) data
dev_unit[i] <- -2 * logdensity.norm(Y[i], mu[i], ep.sig)
dev_sim_unit[i] <- -2 * logdensity.norm(Y_sim[i], mu[i], ep.sig)
}

dev <- sum(dev_unit)
dev_sim <- sum(dev_sim_unit)
sim_mean <- mean(Y_sim)
sim_sd <- sd(Y_sim)

}", jags_mod)

d_jags <- list(
x1 = x1,
x2 = x2,
Y = Y,
n = 100
)

params <- c("dev", "dev_sim", "sim_mean", "sim_sd")

fit <- jags(data = d_jags,
parameters.to.save = params,
model.file = jags_mod,
n.chains = 3,
n.iter = 20000,
n.thin = 5,
n.burnin = 10000,
parallel = TRUE,
verbose = FALSE
)

s <- fit\$sims.list
d_mean <- mean(mu)
d_sd <- sd(mu)

plot(s$$dev, s$$dev_sim,
col = scales::alpha("forest green", .4),
main = "GOF",
ylab = "Simulated",
xlab = "Observed")
abline(0, 1, lwd = 2)

# Here I calculate check values, closer to 0.5 means better fit.
dev_p <- mean(s$$dev_sim > s$$dev)
sd_p <- mean(s$$sim_sd > d_sd) mean_p <- mean(s$$sim_mean > d_mean)

print(paste("Posterior predictive deviance:", round(dev_p, 2)))
print(paste("Posterior predictive sd:", round(sd_p, 2)))
print(paste("Posterior predictive mean:", round(mean_p, 2)))



Output

[1] "Posterior predictive deviance: 0.52"
[1] "Posterior predictive sd: 0.99"
[1] "Posterior predictive mean: 0.38"


I guess this means that our model is fitting the mean ok, but sd of the simulated data is often larger. I've noticed that the "deviance" metric almost always comes out the same: around 0.5.