########################### # 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.