# What does it mean when a low number of quadrature points gives a very different GLMM fit?

I am interested in a logistic regression model with 10 fixed-effects parameters and random intercepts, which I can fit using the lme4::glmer function in R. The value of the nAGQ parameter determines the number of quadrature points used in approximating the likelihood, with 1 corresponding to a Laplace approximation and 0 using "a faster but less exact form of parameter estimation for GLMMs by optimizing the random effects and the fixed-effects coefficients in the penalized iteratively reweighted least squares step" (glmer documentation).

The fixed-effects parameter estimates appear to stabilise once nAGQ gets above 8-10: (The coloured horizontal lines show the estimates from the same model fitted using mgcv::gam (red), gamm4::gamm4 (blue) and SPSS's GENLINMIXED (green) command)

I appreciate that as we increase the number of quadrature points, we get a closer approximation to the log-likelihood, but I came across this message posted by Jim Lindsey to the r-help list in 2001:

When discussing GLMMs using Gauss-Hermite integration, the question of approximations is essentially a red herring. No matter the number of quadrature points, the likelihood is always exact (as exact as any likelihood can be on a digital computer). It is a finite mixture. The approximation question arises in the sense that this finite mixture is more or less close to a Gaussian mixing distribution, which is a completely artificial choice in the first place. It is quite possible for the model with very few quadrature points to fit better than one with sufficient for a very close approximation to the normal mixing distribution, indicating that normal mixing is not a good choice.

It is false to say that the properties of this approximation, in the latter sense are unknown. GLMMs using Gauss-Hermite go back at least to an unpublished tech report by Don Pierce in the mid 70s and the published paper by John Hinde in 1982. Since then, there is a vast literature on the subject of the approximation in the second sense above, especially for the model in question, mixed logistic regression, including work by Alan Agresti, Murray Aitkin, David Brillinger, Bruce Lindsay, etc. (There is also another literature on approximations replacing Gauss-Hermite, such as Breslow and Clayton.) The most recent published reference that I am aware of is earlier this year, but I have refereed others that are not yet in print. 15 to 20 quadrature points gives an extremely close numerical approximation to the normal mixing distribution for most data sets, for what that is worth.

My question has a few parts:

(1) Does this mean that if the results from using a PIRLS fit or Laplace approximation are very different to using a higher number of quadrature points, it suggests that the Gaussian assumption for the random effects is not well-supported by the data?

(2) If so, would it ever make sense to take the PIRLS / Laplace results as the fitted model? Or is a higher number of quadrature points always 'better' (convergence & computational time issues not withstanding)?

(3) In my example, the log-likelihood of the fitted model takes its maximum at nAGQ = 4 (see below). Does this mean that this is the optimal fit (I appreciate that it is only a tiny difference), or are these fitted log-likelihoods not comparable because they are different approximations? (Background information: I am trying to help someone who is only comfortable using SPSS. I haven't been able to find any documentation of the GENLINMIXED fitting algorithm, but the parameter estimates that it gives are very close to those of glmer with nAGQ = 0)

A few points:

• according to this page, SPSS GENLINMIXED uses (one of) pseudo-likelihood, marginal quasi-likelihood, or penalized quasi-likelihood to estimate GLMMs. I'm not sure about MQL/PL, but PQL (at least standard PQL: there are some "improved"/second-order variants out there) is pretty well known to behave badly in many of the same cases where Laplace approximation fails and Gauss-Hermite quadrature is necessary (Breslow 2004). In particular, GH is necessary where the sampling distributions of the conditional modes are far from Gaussian, which happens when there is little information per cluster - e.g., binary responses with few observations per group (Biswas, 2015). These are the cases where we expect a large number of quadrature points to be necessary, e.g. the (in)famous "toenail" (oncholysis) data set (Lesaffre and Spiessens 2001).
• the failure of PIRLS/Laplace doesn't invalidate the assumption that the random effects have a Gaussian distribution; rather, it suggests that the sampling distributions of the conditional modes are not Gaussian. In particular, adaptive GHQ with order $n$ works perfectly if the sampling distributions are proportional to a Gaussian times an $n-1$st order polynomial; $n=1$ (Laplace) works if the sampling distributions are exactly Gaussian, $n=2$ works if the sampling distributions are $\propto (a+bx) G(x)$, $n=3$ works if the sampling distribution is quadratic $\times$ Gaussian, etc.
• Lindsey's comments (that GHQ with a small number of quadrature points can be considered an estimate of a particular finite mixture model), and your comment (that PIRLS/Laplace might actually give the best answer), is interesting - I hadn't thought of it that way. I'd be a little uncomfortable since I know so little about the model that we're implicitly fitting in that case, though ...
• I tend to think of different nAGQ values as being different approximations to the same model, but, again, Lindsey's point is a good one. AFAICT the models should be of the same complexity (there are the same number of parameters), so we shouldn't have to worry about whether they're nested or not. In this case I don't think I'd take that bump at nAGQ=4 very seriously though ...

Biswas, Keya. “Performances of Different Estimation Methods for Generalized Linear Mixed Models.,” 2015. https://macsphere.mcmaster.ca/handle/11375/17272.

Breslow, N. E. “Whither PQL?” In Proceedings of the Second Seattle Symposium in Biostatistics: Analysis of Correlated Data, edited by Danyu Y. Lin and P. J. Heagerty, 1–22. Springer, 2004. http://www.bepress.com/uwbiostat/paper192/.

Lesaffre, Emmanuel, and Bart Spiessens. 2001. “On the Effect of the Number of Quadrature Points in a Logistic Random Effects Model: An Example.” Journal of the Royal Statistical Society: Series C (Applied Statistics) 50 (3): 325–35. doi:10.1111/1467-9876.00237.

• Thanks for this great answer, it's really cleared a few things up. My outcome is binary and I only have 1 or 2 observations per cluster, so it fits under the "little information per cluster" scenario. So in your opinion, I should trust the GHQ results with nAGQ high enough so that estimates have stabilised? – Mark Sep 7 '18 at 3:45
• Yes, the general advise is to start increasing the number of quadrature points until the results are stabilized. Currently, glmer() only allows you to fit the adaptive Gaussian quadrature when you only have a random intercept for a single grouping factor. If you want to include a random slope, you can give a try in the GLMMadaptive package: drizopoulos.github.io/GLMMadaptive – Dimitris Rizopoulos Sep 7 '18 at 5:03