# How do I calculate the posterior predictive distribution in WinBUGS?

I would like to work out the posterior predictive distribution from the following multilevel (mixed effect) model in WinBUGS. The example is taken from an example I found online for illustrative purposes.

Yi ~ Binomial(pi,1) where

logit(pi) = b0+ b1 log(income) + b2 distance + b3 dropout + b4 college + uj(i)

Non-informative priors are given for the fixed effects, assuming bk ~ Normal(0,0.000001). The second parameter is the precision (the reciprocal of the variance), so the variance is one million. We assume that

uj(i) ~ N(0, t)

where the precision t has a gamma prior with parameters 0.001 and 0.001, so the mean is one and the variance is 1000.

In WinBUGS:

model {
# N observations
for (i in 1:N) {
hospital[i] ~ dbern(p[i])
logit(p[i]) <- bcons + blonginc*loginc[i] + bdistance*distance[i] +
bdropout*dropout[i] + bcollege*college[i] + u[group[i]]
}
# M groups
for (j in 1:M) {
u[j] ~ dnorm(0,tau)
}
# Priors
bcons     ~ dnorm(0.0,1.0E-6)
bdistance ~ dnorm(0.0,1.0E-6)
bdropout  ~ dnorm(0.0,1.0E-6)
bcollege  ~ dnorm(0.0,1.0E-6)
# Hyperprior
tau ~ dgamma(0.001,0.001)
}


I can fit this model using stan_glmer which automatically gives me the sample average of the posterior predictive distribution mean_ppd. I would like to work out the same thing in Winbugs but don't know how

To monitor a PPD, you just write out the distribution of the variable that you want to look at. That variable can be an observed node, or a function of the observed node and other variables in your model.

For example, say you want to see if the node "hospital" is fitted well. You need one additional line of code like this:

hospital_1_ppd ~ dbern(p)

, and add 'hospital_1_ppd' to the list of variables you are monitoring. You can figure out the rest.

Hope it helps.

• Thanks for your reply I do see what you're saying but in this dataset there are 1061 hospitals which means I would have to write the line out hospital_1_ppd ~ dbern(p) 1061 times and then manually out of winbugs take the average. Do you know if there is a faster way by writing a few lines of code to give me the sample average (of all 1061 hospitals) of the posterior predictive distribution? – gemster Jul 11 '17 at 8:06
• I would use the following code: for(i in 1:1061) {hospital_pdd[i] <- dbern(p[i])}. Iterate over all cases. You can manipulate the samples in whatever way you want. – CANZ Jul 17 '17 at 15:45
• Also for all your future BUGS questions, you can try this mailing list: BUGS@jiscmail.ac.uk. – CANZ Jul 17 '17 at 15:51