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: 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.
