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.