Three well known methods for constructing prediction intervals in a generalized linear mixed model (GLMM) are the methods based on pseudo-likelihood, Laplace, and Quadrature approximations. All three of these methods are available in the SAS procedure GLIMMIX. We propose a new method based on a mean squared error (MSE) approximation of the empirical best predictor. Following the approach by Harville and Kackar (1984) for a linear mixed model (LMM), we decompose the prediction error into two terms for the purpose of deriving the MSE approximation. Unlike in the LMM case, however, closed form expressions for the two terms in the subsequent MSE approximation are not available. We confront the computational challenge by proposing a Monte Carlo algorithm for evaluating the plug-in estimators of these two terms. For three illustrated examples, the negative binomial, the Poisson and the Bernoulli GLMMs, numerical results show our prediction interval methodology improves the coverage probability over the three methods available in GLIMMIX. Moreover, our results show that with bootstrap adjustments, our method achieves coverage probabilities satisfactorily close to the nominal level.