Solved – Truncating a posterior predictive distribution in JAGS

bayesianjagsposterior

I have run into an error associated with truncating a distribution in JAGS.

In my minimum reproducible example, I have data for 9 observations and would like to find a posterior predictive distribution for the 10th observation. To do this, I include the 10th observation as an NA and estimate its posterior predictive distribution as the variable pi10.

jagsdata <- data.frame(Y = c(47, 126, 68, 43, 67, 80, 61, 9, 26, NA))

model.string <- "
model{
 for (k in 1:10){
    Y[k]  ~ dlnorm(Z[k], tau.sp[k])
    tau.sp[k] ~ dgamma(0.01,0.01)
    Z[k] <- beta.o + beta.sp[k]
 }
  for (g in 1:10) {
    beta.sp[g] ~ dnorm(0, 0.0001)
  }
  beta.o    ~ dgamma (2, 0.04)
  pi10   <- Y[10]
}
"
writeLines(model.string, con = 'jagstest.bug')

library(rjags)
j.model  <- jags.model(file = "jagstest.bug", 
                       data = jagsdata, 
                       n.adapt = 500, 
                       n.chains = 4)
mcmc.object <- coda.samples(model = j.model,
                            variable.names = c('pi10'),
                            n.iter = 5000)

This works, but I would like to truncated the distribution of Y, for example by using the T(1,200). However replacing line 4 above with

    Y[k]  ~ dlnorm(Z[k], tau.sp[k])T(1,200)

gives the error:

Unobserved node inconsistent with unobserved parents at initialization

Although Y with a normal distribution does not give an error.

    Y[k]  ~ dnorm(Z[k], tau.sp[k])T(1,200)

I have read through the JAGS manual section 7 and some examples online, but it is not clear to me how to implement this or why I am getting this error.

Suggestions appreciated.

Best Answer

You can avoid this problem altogether by sampling from the untruncated distribution of Y[k], then (in R) discarding all samples for which Y[k] doesn't lie within the constraint bounds. This is a perfectly valid operation, however, if you have few posterior observations in the feasible region, you'll naturally have a large simulation error associated with your posterior distribution(s).

As a side note, you might want to avoid the Gamma(0.01,0.01) prior on the variance; see for example this presentation by Andrew Gelman and this paper, also by Gelman, for reasons why and alternative suggestions.

Related Solutions

Solved – JAGS: posterior predictive distribution

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?

Answer

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).

Solved – Posterior predictive distribution vs MAP estimate

I often think of it this way. In the fully Bayesian approach, we find the integral

$$p(x^*|X) = \int p(x^*|\theta) p(\theta|X) \text{ d}\theta$$

as integrating over all possible models (infinitely many in fact), and we make a prediction taking all of these models "into consideration". As this is often intractable, we use the MAP estimate of the posterior $p(\theta|X)$, which corresponds to evaluating the same integral but this time using a infinitely small part of $p(\theta|X)$, namely at its maximum. In other words, we multiply $p(x^*|\theta)$ with a new "delta-distribution" located at the max of the posterior distribution and integrate this to obtain the prediction.

The difference is therefore rather obvious: a fully Bayesian treatment corresponds to an infinite ensemble of models, where a given prediction $p(x|\textbf{x},\theta)$ is weighted by the model probability $p(\theta|\textbf{x})$, i.e. more likely models will contribute more to the prediction. The MAP estimate of the parameters will give you the prediction from one specific model, namely the most likely one according to Bayes theorem. Ensemble theory shows us that we often obtain better generalization and more accurate predictions and therefore this will often be "better" than the MAP.

Hope this helps.

Best Answer

Related Solutions

Solved – JAGS: posterior predictive distribution

Solved – Posterior predictive distribution vs MAP estimate

Related Question