I used a Monte-Carlo approach to estimate the prediction interval for a new observation from a GLM using a Negative Binomial distribution. I used this method for linear models and got reliable estimations of the prediction interval but I am not enough confident with statistics to know if this method could also be applied for non-Gaussian distributions.
Could you please tell me if it makes sens to use it or not?
Method:
- compute the se of the expected value with a Wald approach (on the link).
- to use this se to randomly simple 1000 values from a Gaussian distribution with the expected value as $\mu$ and se as $\sigma$.
- for each random value, to randomly sample 1000 values from a negative binomial distribution with the exponential of the random value as mu.
- The prediction interval was computed as the 2.5 & 97.5 quantile of the 1e6 random values obtained for a given expected value.
Does this methodology makes sense?
Regards,
Maxime
Best Answer
It makes sense if you're taking the $\theta$ as known. If you want to incorporate the uncertainty in $\theta$ you would either simulate from the joint distribution (asymptotically Gaussian) of the other parameters and $\theta$, or you could simulate one conditionally on the other (since $\theta$ is an input to the GLM, I'd suggest simulating $\theta$ from the normal approximation to its profile likelihood and then $\mu$ from the conditional of $\mu$ on $\theta$, since you get that bit of information out of the GLM). Then proceed as above to simulate the negative binomial.
This is quite similar to an approach Ripley has discussed (for example on the R-help mailing list) -- e.g. here (though he was particularly terse in that one).
[If anyone has a reference for that approach, it would round out the answer to this question nicely.]