Consider a GLMM setting. Let $i=1,...,m$ denote cluster. Let $l_i(y_i|x_i,b_i;\beta)$ denote the log-likelihood (viewed as a function of the parameters) for each cluster $i$, which includes notably conditioning on random intercept $b_i$. Let $f(b_i;\sigma^2) \sim N(0,\sigma^2)$ for all clusters $i$. Thus, the observed data likelihood for all the data (which are independent across cluster) is $$L(y|x,\beta,\sigma^2)=(2\pi \sigma^2)^{-m/2}\prod_{i=1}^m \int \exp\left\{l_i(y_i|x_i,b_i;\beta)\right\}\exp\left\{-\frac{b_i^2}{2\sigma^2}\right\}db_i$$

The integral expressions are approximated using numerical methods, e.g. Laplace or Gauss-Hermite quadrature. How are these not completely unwieldy integrals in certain very generic cases? For example, suppose $y^T_i=(y_{ij})_{j=1}^{n_i}$ is a vector of binary data, where $n_i$ is cluster size. $l_i$ is as in logistic modeling; then

$$l_i(y_i|x_i,b_i;\beta) = \sum_{j=1}^{n_i} \left[y_{ij}(x_{ij}^T\beta+b_i)-\log(1+\exp(x_{ij}^T\beta+b_i))\right]$$

where $n_i$ is cluster size. Regardless of $\beta$ and $\sigma^2$, doesn't this sum (which is effectively linear in $b_i$), become like pretty negative for $n_i$ moderately large and thus the integrand will be mostly zero across values of $b_i$ (we'd expect this as a product of values between 0 and 1). It seems as though numerical methods would just say the above integral(s) are basically the 0 function, which feels bad for precision. So in short, why do these numerical methods work?


1 Answer 1


The trick that makes these work is that we can estimate where the integrand is not close to zero and set up the numerical integration around those points. The technique is called adaptive Gaussian quadrature, and the Laplace approximation is a special case where the integrand is evaluated only at an estimate of the mode of the integrand. If the likelihood conditional on the random effects were Normal, we could write down the posterior modes of the $b_i$, and there are reasonable approximations to the modes or means for the non-Normal case

Using a blindly constructed grid and Gauss-Hermite quadrature can be quite slow (though it's ok for a one-dimensional problem such as random intercepts). This paper has some comparison of Stata implementations using ordinary Gauss-Hermite quadrature (-xtlogit-) and adaptive quadrature (-gllamm-) and other approximations.

  • $\begingroup$ I hadn't realized that there was a distinction between GH and Adaptive GH quadrature. This is a helpful paper. I see that the intuition is the former is essentially simple MC integration and the latter is importance sampling based MC integration. $\endgroup$
    – Winston
    Oct 24, 2023 at 23:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.