# Parametric Bootstrap without model refitting?

I have a GLM poisson model that makes predictions of the probability of death for every observation in my dataset after fitting the model. I then take the sum or mean of these probabilities to obtain an overall probability of death or expected number of deaths. Now, I'd like to obtain confidence intervals for the mean/sum of these individual predicted probabilities. Initially I had thought of simply doing a non-parametric bootstrap to obtain these confidence intervals, but after speaking with a college of mine, he advised that I perform a parametric bootstrap, which involved using the covariance matrix of the $\hat{\beta}$'s in some manner. When he explained the method to me, he advised that using this method didn't require refitting the model, which is important to me, since the model takes about 5 minutes to run (and I don't want to refit it for 10,000 bootstraps since this will take too long).

Everything I've read about parametric bootstrap with regression models, requires that I simply generate synthetic bootstrap samples using bootstraped residuals and the original predicted values to perform another regression. But what I'm trying to determine is if there is some way to bootstrap confidence intervals without having to refit models as my colleague (who is not unavailable) may have suggested?

Thanks in advance for any suggestions and/or references you might provide.

## 1 Answer

If you are willing to assume that the sampling distribution of the coefficients is multivariate normal you can derive a sampling distribution of predictions (and e.g. quantiles of this distribution) by sampling $\beta^* \sim \textrm{MVN}(\hat \beta, \Sigma)$ and then computing the predicted values on the basis of each $\beta^*_i$ value. I guess I can understand why this might be called a parametric bootstrap, but it makes much stronger assumptions than the usual PB.

If you were doing this in R with any fitting approach that provides coef() (or fixef(), for the nlme/lme4 family) and vcov() you could do:

betastar <- MASS::mvrnorm(1000,coef(fit),vcov(fit))
pred <- apply(betastar,1,predfun)


In section 7.5.3 of Ecological Models and Data in R (Princeton University Press, 2008: an old draft is online here) I call these population prediction intervals, after Lande et al. 2003.

• Thanks, @Ben Bolker! Much appreciated. I have your book right here on my bookshelf, so I'll take a look at Section 7.5.3. Can you speak to the advantages/disadvantages of this approach compared to traditional bootstrap resampling of cases? My colleague mentioned that this approach would be more appropriate for predicting future Population Prediction Intervals than traditional bootstrap resampling of cases. I believe he made mention that the traditional bootstrap would only be applicable to my current population whereas this method was better suited for future predictions. Can you comment? – StatsStudent Oct 10 '15 at 18:40
• I was also wondering if the intervals calculated in this way would be confidence intervals or prediction intervals? I.e. would they be confidence intervals, i.e. intervals on the predicted mean of the model, or would it be prediction intervals, and give an estimate of an interval in which future observations will fall, taking into account the expected noise on individual measurements? Would it be possible to adapt the method to calculate both for say GLMs with arbitrary distributions? – Tom Wenseleers Mar 14 '18 at 19:58