0
$\begingroup$

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"

Model deviance vs simulated deviance?

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.

$\endgroup$

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.