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

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