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?
