  A method is introduced for approximate marginal likelihood inference via
adaptive Gaussian quadrature in mixed models with a single grouping factor. The
core technical contribution is an algorithm for computing the exact gradient of
the approximate log marginal likelihood. This leads to efficient maximum
likelihood via quasi-Newton optimization that is demonstrated to be faster than
existing approaches based on finite-differenced gradients or derivative-free
optimization. The method is specialized to Bernoulli mixed models with
multivariate, correlated Gaussian random effects; here computations are
performed using an inverse log-Cholesky parameterization of the Gaussian
density that involves no matrix decomposition during model fitting, while Wald
confidence intervals are provided for variance parameters on the original
scale. Simulations give evidence of these intervals attaining nominal coverage
if enough quadrature points are used, for data comprised of a large number of
very small groups exhibiting large between-group heterogeneity. In contrast,
the Laplace approximation is shown to give especially poor coverage and high
bias for data comprised of a large number of small groups. Adaptive quadrature
mitigates this, and the methods in this paper improve the computational
feasibility of this more accurate method. All results may be reproduced using
code available at \url{https://github.com/awstringer1/aghmm-paper-code}.
