Solved – How to calculate the posterior predictive distribution in WinBUGS

mixed modelposteriorpredictive-modelsstanwinbugs

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)
bloginc   ~ 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

Best Answer

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.

Related Question