# 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[1]" is fitted well. You need one additional line of code like this:

hospital_1_ppd ~ dbern(p[1])

, 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[1]) 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? Jul 11, 2017 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. Jul 17, 2017 at 15:45
• Also for all your future BUGS questions, you can try this mailing list: [email protected]. Jul 17, 2017 at 15:51