Laplace Approximation for Poisson Random-Slope Models

Why the Conditional Mode Matters

In a Poisson generalized linear mixed model, the random effects are not observed. The marginal likelihood for a cluster therefore requires integrating its conditional likelihood over the random-effects distribution. With random intercepts and slopes, this is generally a multidimensional integral whose integrand need not have a closed-form integral. The Laplace approximation replaces the log integrand near its maximum by a quadratic function, making the integral locally Gaussian.

The approximation depends on two quantities evaluated at the conditional mode: the log integrand itself and its curvature there. The mode locates the largest contribution to the integral; the curvature determines the volume around that mode that contributes appreciably. A correct mode with an incorrect Hessian, or a correct Hessian evaluated at the wrong point, gives an incorrect approximation.

Model, Log Integrand, and Laplace Formula

For cluster \(j\), let \(Y_{ij}\) be the count for observation \(i\), with fixed-effect design vector \(x_{ij}\) and random-effect design vector \(z_{ij}=(1,t_{ij})^\top\). The random effect \(b_j=(b_{0j},b_{1j})^\top\) contains an intercept and a slope for covariate \(t_{ij}\). The model is

Here \(D\) is a positive-definite covariance matrix. The conditional independence assumption applies within a cluster given its random effect; random effects are also assumed independent across clusters. The marginal likelihood integrates over the Gaussian distribution of \(b_j\), rather than treating its realized value as known.

More explicitly, up to terms that do not depend on \(b\),

The factorial terms do not affect the mode or Hessian, but they belong in the likelihood value \(h_j(\widehat b_j)\). The Gaussian normalizing constant and covariance determinant also belong in that value. Omitting these constants may leave the mode unchanged while shifting the marginal likelihood.

Score and Curvature

Differentiating with respect to \(b\) gives the score equation

Thus the mode solves \(\nabla h_j(\widehat b_j)=0\), and the curvature used in the approximation is \(K_j=K_j(\widehat b_j)\). The Poisson contribution is a weighted sum of outer products, with positive weights \(\mu_{ij}\). Since \(D^{-1}\) is positive definite, \(K_j(b)\) is positive definite for every finite \(b\) under this model. This gives strict concavity of \(h_j\), so any solution to the score equation is its unique global maximum. In practice, numerical convergence still needs checking: a computation can stop before reaching that solution.

For a random intercept and slope, write \(z_{ij}=(1,t_{ij})^\top\). The entries of the negative Hessian are

The off-diagonal term matters: intercept and slope curvature are usually coupled, and the determinant is not generally the product of the diagonal entries. If the random effects are correlated, the off-diagonal entries of \(D^{-1}\) also enter the Hessian.

Why the Formula Has This Form

The local quadratic approximation follows from a Taylor expansion about the mode. Because the score is zero there, the linear term vanishes:

Exponentiating this expression gives \( \exp\{h_j(\widehat b_j)\} \) times the kernel of a \(q\)-dimensional Gaussian function with precision matrix \(K_j\). Integrating that Gaussian function over \(b\) gives \((2\pi)^{q/2}|K_j|^{-1/2}\). Taking the logarithm yields the stated Laplace formula. The approximation is exact for a quadratic log integrand; otherwise its accuracy depends on how well that quadratic describes the integrand over the region contributing to the integral.

Diagnostics, Alternatives, and Limitations

A reliable calculation should check the mode, not merely accept an optimizer’s returned value. Evaluate the score at the proposed mode and confirm that its components are close to zero on a numerically meaningful scale. Confirm that the Hessian is symmetric and that its eigenvalues are positive, or use a Cholesky factorization to verify positive definiteness. Recompute the log likelihood with all factorial and Gaussian density terms included. For multiple clusters, form each cluster’s contribution using its own mode and Hessian, then sum the contributions.

A random slope can be weakly informed when a cluster has few observations or little variation in \(t_{ij}\). The Gaussian prior still contributes curvature, but it does not create data information about the slope. A very broad or asymmetric integrand can make the local quadratic approximation inaccurate even when the mode and Hessian are computed correctly. Sparse counts, influential observations, extreme fitted means, and near-confounding of intercept and slope are reasons to compare the result with more direct integration where feasible.

Adaptive Gauss-Hermite quadrature approximates the same random-effects integral by evaluating the integrand at weighted nodes centered and scaled around the conditional mode. It can represent departures from a local Gaussian shape more accurately as the quadrature order increases, at greater computational cost. Laplace uses only the mode and local curvature, so it is usually cheaper but can be less accurate when the integrand is not well described by its local quadratic approximation. Neither method repairs a misspecified Poisson mean-variance relationship, an unsuitable random-effects distribution, or a poorly chosen random-effects structure.

One common mistake is to use only the diagonal Hessian entries and omit the off-diagonal term when computing \(|K|\). That treats intercept and slope curvature as independent even when the data couple them; for a two-dimensional matrix, the determinant is \(K_{00}K_{11}-K_{01}K_{10}\). Another mistake is to use \(+\nabla^2h\) as the curvature matrix. At a maximum, \(\nabla^2h\) is negative definite, whereas the Laplace formula uses \(K=-\nabla^2h\), which is positive definite.