  While generalized linear mixed models (GLMMs) are a fundamental tool in
applied statistics, many specifications -- such as those involving categorical
factors with many levels or interaction terms -- can be computationally
challenging to estimate due to the need to compute or approximate
high-dimensional integrals. Variational inference (VI) methods are a popular
way to perform such computations, especially in the Bayesian context. However,
naive VI methods can provide unreliable uncertainty quantification. We show
that this is indeed the case in the GLMM context, proving that standard VI
(i.e. mean-field) dramatically underestimates posterior uncertainty in
high-dimensions. We then show how appropriately relaxing the mean-field
assumption leads to VI methods whose uncertainty quantification does not
deteriorate in high-dimensions, and whose total computational cost scales
linearly with the number of parameters and observations. Our theoretical and
numerical results focus on GLMMs with Gaussian or binomial likelihoods, and
rely on connections to random graph theory to obtain sharp high-dimensional
asymptotic analysis. We also provide generic results, which are of independent
interest, relating the accuracy of variational inference to the convergence
rate of the corresponding coordinate ascent variational inference (CAVI)
algorithm for Gaussian targets. Our proposed partially-factorized VI (PF-VI)
methodology for GLMMs is implemented in the R package vglmer, see
https://github.com/mgoplerud/vglmer . Numerical results with simulated and real
data examples illustrate the favourable computation cost versus accuracy
trade-off of PF-VI.
