# Numerical Approximation of Integrals in GLMM

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?

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.