Worked Derivation of a Random-Intercept Laplace Likelihood

From the Random-Effects Integrand to a Marginal Likelihood

A binary random-intercept model assigns each cluster a latent intercept that shifts the log odds for every observation in that cluster. The marginal likelihood for the observed responses averages their conditional likelihood over the assumed distribution of that intercept. This tutorial derives a Laplace approximation to that integral, keeping track of the mode, curvature, and normal-density normalization that determine the resulting log likelihood.

Suppose there are \(J\) independent clusters. In cluster \(j\), response \(Y_{ij}\) is binary, and, conditional on the cluster random intercept \(b_j\), responses are independent with \(Y_{ij}\mid b_j\sim\operatorname{Bernoulli}(p_{ij})\), where \(\operatorname{logit}(p_{ij})=\mathbf{x}_{ij}^{\mathsf T}\boldsymbol{\beta}+b_j\). Assume \(b_j\sim N(0,\sigma_b^2)\), independently across clusters, with \(\sigma_b^2>0\). The fixed effects \(\boldsymbol{\beta}\) and variance \(\sigma_b^2\) are held fixed when deriving each integral.

For observed responses \(y_{ij}\), the conditional log likelihood is \(\ell_j(b)=\sum_i\{y_{ij}\eta_{ij}(b)-\log[1+\exp\{\eta_{ij}(b)\}]\}\), where \(\eta_{ij}(b)=\mathbf{x}_{ij}^{\mathsf T}\boldsymbol{\beta}+b\). This expression uses the individual Bernoulli observations; there is no binomial coefficient to add. If the data are instead recorded as a binomial count, the corresponding count likelihood includes its binomial coefficient, which must be retained when calculating the marginal likelihood of that count.

Mode, Curvature, and the Laplace Formula

Write \(p_{ij}(b)=\{1+\exp[-\eta_{ij}(b)]\}^{-1}\). Differentiating the conditional log likelihood and the normal-kernel term gives

The conditional mode \(\hat b_j\) solves \(q'_j(\hat b_j)=0\). Define the negative second derivative at that mode as the scalar curvature \(H_j=-q_j''(\hat b_j)\). Thus \(H_j=\sum_i p_{ij}(\hat b_j)\{1-p_{ij}(\hat b_j)\}+1/\sigma_b^2\), which is strictly positive. In this one-dimensional setting, the Hessian determinant is simply \(H_j\).

These derivatives also establish that there is exactly one finite mode for each cluster at fixed finite \(\boldsymbol{\beta}\) and positive \(\sigma_b^2\): \(q_j''(b)<0\) everywhere, so \(q'_j(b)\) is strictly decreasing. Moreover, the normal penalty makes \(q_j(b)\) tend to negative infinity as \(|b|\) grows. The unique maximum can therefore be found by solving the score equation.

The factor \(\sigma_b^2H_j\) in the final expression is dimensionless: \(H_j\) has inverse-variance units on the random-intercept scale. The normalization is consequential. Keeping the unnormalized kernel but omitting the normal-density constant would add \(\log(\sqrt{2\pi}\sigma_b)\) to the cluster result and would not be the likelihood under the stated normal random-effects distribution.

Combining Clusters and Checking the Calculation

Because clusters are independent and their random intercepts are independent, the joint marginal likelihood factors into cluster integrals: \(L=\prod_{j=1}^J L_j\). Taking logarithms turns that product into a sum. The approximation is applied separately to each cluster integral, using its own mode and curvature, and the resulting log contributions are then added.

A useful calculation check is to verify the score and curvature with the same unrounded mode and probabilities, then round reported results. The score should be close to zero; \(H_j\) must be positive; and the determinant correction must use the product \(\sigma_b^2H_j\). Small rounding differences in intermediate displayed values are expected, but they should not materially alter the final reported precision.

Assumptions, Diagnostics, and Alternatives

The derivation assumes conditionally independent Bernoulli responses given the random intercept, a normal random-intercept distribution with positive variance, independent clusters, and finite fixed predictors. The Laplace step additionally assumes that a local quadratic approximation around the mode adequately represents the integrand’s important mass. A unique mode and positive curvature are guaranteed here by strict concavity, but these facts alone do not ensure that the quadratic approximation is accurate.

Check numerical and modeling warning signs before relying on the approximation:

• Mode check: evaluate \(q'_j(\hat b_j)\) numerically. A visibly nonzero score can indicate failed optimization, inadequate convergence tolerance, or a coding error.
• Curvature check: confirm \(H_j>0\), compute it from the same mode used in the score, and inspect whether the integrand is sharply or weakly concentrated. A broad or strongly non-quadratic integrand can make the local approximation less dependable.
• Scaling check: use the normalizing constant and determinant in the same parameterization as the integral. A change of integration variable requires its Jacobian; dropping it changes the integral.
• Approximation check: compare Laplace results with numerical quadrature where feasible, and examine changes under increased quadrature accuracy. Material disagreement indicates that Laplace error may matter for the likelihood or subsequent estimates.
• Model check: integration accuracy cannot repair a misspecified conditional response model, an unsuitable random-effects distribution, or dependence left unmodeled within clusters.

Numerical quadrature approximates the integral by weighted evaluations across the random-effects distribution and can be more accurate when the integrand is not well represented by one quadratic neighborhood, at additional computational cost. Laplace is a local, mode-and-curvature approximation and is often computationally efficient. Simply substituting the conditional mode is not a competing marginal-likelihood calculation: it ignores both the random-effects penalty and the integral-width correction. The approximation can be limited when cluster information is weak or the integrand is asymmetric; a positive curvature alone does not rule out those cases.