I have a relatively complex logistic model I built using generalized estimating equations (GEE). The model is used to predict population mortality. I have gone about estimating the mortality in the population by (1) fitting a model to my data, (2) making individual-level predictions of mortality from the data, and finally by (3) taking the mean of these predictions. I am particularly interested in examining the shape of the model as I vary one predictor, holding all other variables constant. To this end, I need to calculate a confidence band around the population prediction function from my model. I can calculate population prediction intervals just fine (using the advice given by @Ben Bolker here Parametric Bootstrap without model refitting?), but now I'm trying to figure out how I might be able to bootstrap confidence bands/simultaneous confidence intervals instead of just bootstrap point-wise confidence intervals. Can anyone point me in the right direction?