We present a variational maximization–maximization algorithm for approximate maximum likelihood estimation of generalized linear mixed models with crossed random effects (e.g., item response models with random items, random raters, or random occasion-specific effects). The method is based on a factorized variational approximation of the latent variable distribution given observed variables, which creates a lower bound of the log marginal likelihood. The lower bound is maximized with respect to the factorized distributions as well as model parameters. With the proposed algorithm, a high-dimensional intractable integration is translated into a two-dimensional integration problem. We incorporate an adaptive Gauss–Hermite quadrature method in conjunction with the variational method in order to increase computational efficiency. Numerical studies show that under the small sample size conditions that are considered the proposed algorithm outperforms the Laplace approximation.
