Laplace Approximation with Cluster-Level Binomial Counts

Why the Grouping Convention Matters

When several binary observations in a cluster share the same covariates, they can be represented by a success count and a trial count. This is a compact data representation, but it is important to distinguish two likelihoods: the probability of one particular sequence of individual outcomes, and the probability of observing a total count. The latter includes a binomial coefficient because many outcome sequences produce the same count.

In a generalized linear mixed model, the shared random intercept also induces dependence between outcomes in the same cluster after that intercept is integrated out. The right starting point is therefore the conditional likelihood given the intercept; the marginal cluster likelihood follows by integrating that conditional likelihood against the random-intercept distribution.

From Individual Outcomes to the Cluster Integrand

Let \(Y_j=\sum_{i=1}^{m_j}Y_{ij}\) be the number of successes in cluster \(j\), and write its observed value as \(y_j\). Set \(\eta_j(b)=x_j^\mathsf{T}\beta+b\) and \(p_j(b)=\{1+\exp[-\eta_j(b)]\}^{-1}\). If the individual outcomes and their labels are retained, their conditional likelihood is \(p_j(b)^{y_j}\{1-p_j(b)\}^{m_j-y_j}\). If only the count is observed, its conditional probability is the binomial probability, which also includes \({m_j\choose y_j}\).

The second equality follows from \(p=\exp(\eta)/(1+\exp(\eta))\) and \(1-p=1/(1+\exp(\eta))\). Integrating over the random intercept gives the marginal likelihood contribution of the count:

For numerical work, separate the binomial coefficient and the normalizing constant of the normal density from the exponent. Define \(h_j(b)=y_j\eta_j(b)-m_j\log\{1+\exp(\eta_j(b))\}-b^2/(2\sigma_b^2)\). Then the likelihood is the integral of \(\exp\{h_j(b)\}\), multiplied by the factors outside that exponent:

What the Binomial Coefficient Does—and Does Not Do

For fixed \(m_j\) and \(y_j\), the factor \({m_j\choose y_j}\) does not depend on \(b\), \(\beta\), or \(\sigma_b^2\). It can therefore be omitted when locating the mode of \(h_j\), or when differentiating \(h_j\) with respect to \(b\). It cannot be omitted from the probability of the count or from the corresponding marginal likelihood value. In particular, likelihood contributions for counts and for one labelled sequence differ by this factor even after integration: it is constant with respect to \(b\), so it remains outside the integral.

If the same observed data are represented consistently and the same model is compared across parameter values, a parameter-independent factor does not change the maximum-likelihood estimates. But likelihood values are not interchangeable across representations unless the constants are handled consistently. The grouped binomial likelihood is the appropriate probability model when the recorded observation is the count.

The Laplace approximation replaces \(h_j(b)\) near its mode by its quadratic Taylor approximation. At an interior mode, \(h_j'(\hat b_j)=0\), so that approximation is \(h_j(\hat b_j)-H_j(b-\hat b_j)^2/2\). Integrating the resulting Gaussian-shaped function gives the usual one-dimensional Laplace formula. Keeping the likelihood’s external factors yields:

This expression shows why the binomial coefficient should be retained in the marginal likelihood even though it has no effect on the conditional mode or curvature. The approximation uses the local height of the integrand and its local width; it is most persuasive when the integrand is concentrated around a smooth, well-defined mode and is reasonably close to a Gaussian there.

When Counts Can Be Grouped

A single binomial count is a valid representation only when the trials being combined have the same conditional success probability given the cluster effects and model parameters. In the random-intercept model above, that condition holds within cluster when all trials share \(x_j\). If covariates vary between individuals, their probabilities differ even after conditioning on the common intercept; one total success count generally cannot preserve the individual likelihood.

Diagnostics, Alternatives, and Limitations

Before applying the approximation, check that the data representation matches the model. Within each proposed group, verify that trials share the relevant covariates and the same conditional probability. Confirm that \(0\le y_j\le m_j\), that the random-intercept variance is positive, and that the integral includes the normal density and the count coefficient. A missing normalizing constant changes the marginal likelihood; an omitted binomial coefficient changes the count likelihood.

The strict concavity displayed above is also a useful computational check for this one-dimensional random-intercept integrand: there should be no competing local modes. A root-finding result that does not satisfy \(y_j-m_jp_j(\hat b_j)-\hat b_j/\sigma_b^2\approx0\), or a nonpositive computed \(H_j\), indicates an implementation or numerical error. Even with a correctly found mode, the quadratic approximation can be inaccurate if the integrand is appreciably skewed or spread out. Comparing against adaptive quadrature is a direct alternative when a more accurate numerical integral is needed; the exact integral is the reference target, while Laplace is a local approximation.

The grouped-count form is efficient and exact under its assumptions, but it does not make outcomes marginally independent: integrating over a shared \(b_j\) couples them. Nor does a correct count likelihood guarantee that a Laplace approximation is accurate for every cluster or parameter value. These are separate issues—one concerns the probability model and data grouping, the other the numerical approximation to its integral.