Suppose I use R to fit a Generalized Linear Mixed Model from the binomial family and with a logit link. How do I obtain the prediction intervals (as opposed to the confidence intervals) for the fixed effects (as opposed to the random effects) in a way that incorporates the variability described by the random effects? Thank you.
We figured out a good way to do this that is good enough for our purposes (which is utilizing the predictive distribution to draw samples of empirically estimated fixed effects to calibrate an agent-based model). Here's the step-by-step process:
Simulate the $\beta_j$'s from the variance-covariance matrix output of the GLMM.
Sample the standard error of the standard deviation of the random intercepts from a chi-square distribution with $n - 1$ degrees of freedom. See this R-sig-ME post for details.
From the $\beta_j$'s estimated in step one, sample the $\beta_j$'s again from a normal distribution with mean equal to the $\beta_j$ and standard deviation equal to that sampled in step two.
Repeat steps 1-3 many times, saving the results at each iteration, to obtain an empirical predictive distribution of the $\beta_j$'s.