# Predict credible interval of Poisson-distributed response based on Lambda credible interval

I am approaching Bayesian inference. Could you review my steps and give me a hand with my model predictions?

I am using a N-mixture model to predict how many individuals N of a rare species inhabit a site i based on the number of individual y observed at multiple visits k of the same site i (probability of detection p < 1):

Ni ∼ Poisson(λ)
yi,k ∼ binomial(Ni,p)


βs and γs are the estimated effects of covariates that affect species abundance and detection probability:

log(λ) = β0 + β1*x1 + β2*x2
logit(p) = γ0 + γ1*x1 + γ2*x2


So far:

1. I sampled the posterior distribution of each β and γ (mean and standard deviation) with MCMC;

2. since the posterior distributions of βs and γs were normal, I simulated 1000 times log(λ) and logit(p) using functions such as

fun1 <- function(beta1, beta2) {
((rnorm(1000, mu0, sd0)) +
(beta1*rnorm(1000, mu1, sd1)) +
(beta2*rnorm(1000, mu2, sd2)))
}

3. I used the link function on the simulated values (e.g., exponentiated the predicted log(λ)) and created 95% credible intervals by selecting the interval between the 2.5% and 97.5% quantiles of my simulations.

Now, what I need to do is to find a credible interval for $N_i$. I explored two approaches, but I wonder which one is correct (if any of the two $i$s) since results are different:

a) simulate for each predicted lambda a discrete value using rpois. This returns credible intervals between discrete values, such as an average prediction abundance of 0.3 with 95% of predictions between 0 and 2 (Increasing the simulation number did not change the interval too much).

b) use a loop to simulate for each lambda 1000 predictions with rpois, and then average them to obtain 1000 average predicted abundance values. This returns credible intervals between continuous values, such as an average prediction abundance of 0.3 with 95% of predictions between 0.15 and 0.45 (Increasing the simulation number tend to narrow the credible interval around the predicted mean).

Any insight on which approach is correct (or on alternative approaches, if these approaches are not valid) will be much appreciated.

• Please merge your accounts; you will then be able to edit your question directly. – mkt - Reinstate Monica Sep 19 '18 at 19:57
• Please register &/or merge your accounts (you can find information on how to do this in the My Account section of our help center), then you will be able to edit & comment on your own question. – gung - Reinstate Monica Sep 19 '18 at 20:37

The difference between the two approaches is that the former generates a credible interval for $N_i$ based on the posterior predictive distribution (or your approximation thereto), but the latter generates a credible interval for $\sum_{i=1}^{1000}N_i / 1000$ (based on the same approximation), which is a very different thing, and doesn't seem particularly meaningful. As the standard deviation of the sample mean tends to be proportional to $1/\sqrt{n}$, it's not surprising that increasing the simulation number tends to narrow the credible interval. (Of course, the credible interval will never get narrower than the C.I. for $\lambda$ itself.)
1. Sample the posterior distribution of each $\beta$ and $\gamma$, as in your step 1,
2. Use those samples to calculate sample values of $\lambda$ based on observed values of $x_1$ and $x_2$,
3. Use the sample values of $\lambda$ to generate sample values of $N$.
This gets you away from relying on the Normality assumption. Note that in step 2 you should probably select which $x_1$ and $x_2$ to use, since obviously $\lambda$ is conditional upon them and the distribution of $N$ will therefore vary as $x_1$ and $x_2$ vary. If you don't, it's not clear just what $N$ represents; it might (in effect) cross species boundaries or sites, making its interpretability doubtful.