###########################
# Bayesian p-value example (see Example 2.17, Carlin and Louis, 2008)
# data is from the "stacks" example (see WinBUGS Examples Vol I)
###########################
# The trick: we add a new replicate data vector "Ystar", which is not observed.
# WinBUGS will estimate Ystar during the run, and this is used as the predicted
# Y* vector in the Bayesian p-value calculation.
# We use three different "discrepancy measures":
# Dp is the "omnibus goodness-of-fit" measure of Gelman et al. (1995, p.172);
# Sp sets D(y, theta)= min(Y), to check fit in the lower tail of the distribution
# Bp sets D(y, theta)= max(Y), to check fit in the upper tail of the distribution
model
{
# standardise x's and coefficients
for (j in 1 : p) {
b[j] <- beta[j] / sd(x[ , j ])
for (i in 1 : N) {z[i,j] <- (x[i,j] - mean(x[,j])) / sd(x[,j])}
}
b0 <- beta0 - b[1]*mean(x[,1]) - b[2]*mean(x[,2]) - b[3]*mean(x[,3])
# statistical model
for (i in 1 : N) {
Y[i] ~ dnorm(mu[i], tau)
mu[i] <- beta0 + beta[1]*z[i,1] + beta[2]*z[i,2] + beta[3]*z[i,3]
Ystar[i] ~ dnorm(mu[i], tau)
D[i]<- (Y[i]-mu[i])* (Y[i]-mu[i]) *tau
Dstar[i] <- (Ystar[i]-mu[i])* (Ystar[i]-mu[i])*tau
p_val[i] <- step(Dstar[i] - D[i])
}
sumDstar<- sum(Dstar[])
sumD<- sum(D[])
Dp<- step(sumDstar-sumD)
Sp<-step(ranked(Y[], 1) - ranked(Ystar[], 1)) #min(Y) =7
Bp<-step(ranked(Ystar[], N) - ranked(Y[], N)) #max(Y)=42
# priors
beta0 ~ dnorm(0, 0.00001)
for (j in 1 : p) {beta[j] ~ dnorm(0, 0.00001)}
sigma ~ dunif(0.01,100) # vague Gelman prior for sigma
tau<- 1/(sigma*sigma)
} # end of stacks Bayesian p-value model
#Inits:
list(beta0 = 10, beta=c(0,0, 0), sigma = 0.1)
list(beta0 = -10, beta=c(10,10, 10), sigma = 10)
#Data:
list(p = 3, N = 21,
Y = c(42, 37, 37, 28, 18, 18, 19, 20, 15, 14, 14, 13, 11, 12, 8, 7, 8, 8, 9, 15, 15),
x = structure(.Data = c( 80, 27, 89,
80, 27, 88,
75, 25, 90,
62, 24, 87,
62, 22, 87,
62, 23, 87,
62, 24, 93,
62, 24, 93,
58, 23, 87,
58, 18, 80,
58, 18, 89,
58, 17, 88,
58, 18, 82,
58, 19, 93,
50, 18, 89,
50, 18, 86,
50, 19, 72,
50, 19, 79,
50, 20, 80,
56, 20, 82,
70, 20, 91), .Dim = c(21, 3)))
# Results (posterior mean, after 1k burn in):
Dp=0.544
Sp=0.955
Bp=0.354
p_val[1] = .395
p_val[3] = .223
p_val[4] = .126
p_val[21] = .074
# Conclusions:
# Using the omnibus-goodness-of-fit measure, Bayes p-value=0.53.
# Thus the model fits data fairly well. From our previous look at
# the stacks example, we know the Bayesian p-values go in "right direction",
# but are not significant, possibly reflecting the usual "using the
# data twice" problem with the Bayesian p-value.