# JAGS: posterior predictive distribution

I am fitting a simple linear regression model with R & JAGS:

model {
for (i in 1:length(ri)){
ri[i] ~ dnorm(mu[i],tau)
mu[i] <- alpha + b.vr*vr[i] + b.ir*ir[i]
}

#posterior predictive distribution
for (j in 1:length(pvr)){
pmu[j] <- alpha + b.vr*pvr[j] + b.ir*pir[j]
pri[j] ~ dnorm(pmu[j],tau)
}

#priors for regression
alpha ~ dnorm(0,0.01)
b.vr ~ dunif(-10,0)
b.ir ~ dunif(0,10)

#prior for precision
tau ~ dgamma(0.001,0.001)
sigma <- 1/sqrt(tau)
}


I am trying to calculate the posterior predictive distribution with new data (pvr and pir in the model). But somehow the reported results (pri) do not make sense (the means of pri are smaller than expected).

Could someone explain me, is something wrong with the model specification?

• could you at a minimum provide the input data that you are using? – David LeBauer Feb 15 '11 at 20:24

disclaimer I still don't fully understand your model; but without at least a reproducible example, this is the best I can offer. It is not clear exactly what you are doing here. For example, how are pvr and pir calculated? Would it make sense to calculate them inside the same model?

I am assuming that your data includes observations for mu[] but not pmu[] and you want to estimate pmu[j] given j values of pvr and pir.

Append the pir and pvr to the ir and vr columns, get rid of the second for loop, and then consider the values of mu[] estimated using pir and pvr to be the posterior predictive estimates of mu. Then replace the two for loops with this:

for (i in 1:length(ri)+length(pri)){
ri[i] ~ dnorm(mu[i],tau)
mu[i] <- alpha + b.vr*vr[i] + b.ir*ir[i]
}


I have done something similar, but without predicted regressors, similar to the example given by Gelman et al in 'Bayesian Data Analysis' (pp 598-599 starting under posterior predictive simulations).

• @David pvr and pri are user supplied data: they are the predicted regressors... – teucer Feb 16 '11 at 8:10
• @teucer are they predicted from user-supplied values of ir and vr? if so, how? and if so, why not calculate pri and pvr yourself inside the JAGS model? – David LeBauer Feb 16 '11 at 16:06
• @teucer, this is a shot in the dark since I still don't understand your model, and this might be the same as what you already have, but would it make sense to append the pir and pvr to the ir and vr columns, get rid of the second for loop, and then consider the values of mu[] estimated using pir and pvr to be the posterior predictive estimates of mu? – David LeBauer Feb 16 '11 at 16:07
• @David I want to have some credible interval for the predictions (pri), this is why having mu is not enough. If I append pvr and pir I would get alpha, b.vr and b.ir for the whole ri, right? Or can I have ri with missing values e.g. (1,2,3,NA,NA) and append the regressors? – teucer Feb 16 '11 at 22:18
• @teucer yes, you can have ri with missing values; this is a neat feature of bugs that simplifies calculation of a predictive interval. – David LeBauer Feb 17 '11 at 15:46