Quadrature For Multimarker Random-Effects Integrals

Why the Integral Grows With the Marker Panel

A joint model connects repeated measurements from several markers to a subject’s event time through latent subject-specific effects. To obtain a likelihood for the observed data, those effects must be integrated out. Adding marker-specific random intercepts or slopes can therefore enlarge not only the model but also the dimension of the numerical integral required for each subject. This computational consequence matters: a method that is reliable for one or two random effects may become expensive or inaccurate when the random-effects vector is larger.

Let \(Y_{ik}\) denote the observed measurements for subject \(i\) on marker \(k\), let \(T_i\) be the observed event or censoring time, and let \(\delta_i\) indicate whether the event was observed. Stack the subject’s random effects across all markers into a vector \(b_i\) of dimension \(q\). For example, a marker-specific random intercept and slope for each of \(K\) markers gives \(q=2K\). The dimension is the number of random-effect components being integrated, not the number of measurements.

For illustration, suppose \(b_i\sim N_q(0,D)\), where \(D\) is a positive-definite covariance matrix. Conditional on \(b_i\) and the model covariates, write the longitudinal contribution as \(L_{Y_i}(b_i)\), and the event-time contribution as \(L_{T_i}(b_i)\). Under conditional independence of marker outcomes given the random effects, and a conditionally specified survival model, the subject’s marginal contribution is

Here \(\phi_q\) is the \(q\)-variate normal density, \(h_i\) is the conditional hazard, and \(S_i\) is the conditional survival function. The displayed factorization is an assumption of the specified model, not a universal property of multimarker data: residual dependence between markers not represented by \(b_i\) requires a different longitudinal contribution. For censored observations, \(\delta_i=0\), so the survival factor contributes \(S_i(T_i\mid b_i)\), rather than a hazard term.

The numerical methods below assume that the conditional likelihood and random-effects density define a finite integral, and that the log integrand is sufficiently smooth near its dominant region. Laplace approximation additionally needs a well-defined interior mode with negative-definite Hessian there. Both methods can mislead if the assumed conditional models, random-effects distribution, or censoring assumptions are inappropriate; computation does not repair model misspecification.

How Dimension Affects Computation

A tensor-product quadrature rule with \(m\) nodes in each coordinate evaluates the integrand at \(m^q\) points. This exponential growth is the central challenge. Covariance among random effects does not itself add dimensions, but it can rotate and stretch the region where the integrand has appreciable mass, making an unadapted grid inefficient. The following count illustrates the growth before any optimization or data-management overhead is included.

Adaptive Gaussian Quadrature

Gaussian quadrature approximates an integral by a weighted sum of integrand evaluations. Adaptive Gaussian quadrature (AGQ) centers and scales the nodes using the mode and local curvature of the integrand. In a joint likelihood, define the log integrand

Let \(\widehat b_i\) maximize \(\ell_i(b)\), and let \(H_i=-\nabla^2\ell_i(\widehat b_i)\) be its negative Hessian, assumed positive definite for this local approximation. Choose a matrix \(A_i\) such that \(A_iA_i^\mathsf{T}=H_i^{-1}\). For a standard-normal quadrature node \(z\), the corresponding adaptive node in random-effects space is \(b=\widehat b_i+A_i z\). The local curvature therefore places nodes more closely in directions where the integrand is narrow.

For a tensor rule with nodes \(z_r\) and normalized weights \(w_r\), an AGQ approximation can be written as

The weights here sum to one under the standard-normal rule. The exponential correction measures how the actual log integrand differs from its local quadratic approximation at each node. With more nodes, AGQ can represent some of this nonquadratic shape; it does not make a misspecified model correct, and a finite quadrature rule is still an approximation.

Laplace Approximation and Its Connection to AGQ

Laplace approximation replaces the log integrand near its mode by a second-order Taylor expansion. At an interior maximum \(\widehat b_i\), the gradient is zero, so

Exponentiating gives a Gaussian-shaped approximation to the integrand. Integrating that Gaussian over \(\mathbb{R}^q\) gives the Laplace formula:

The formula follows by setting \(u=H_i^{1/2}(b-\widehat b_i)\). The quadratic term becomes \(-u^\mathsf{T}u/2\), while the change of variables has Jacobian \(|H_i|^{-1/2}\). Thus the remaining integral is \(\int_{\mathbb{R}^q}\exp(-u^\mathsf{T}u/2)\,du=(2\pi)^{q/2}\). In the AGQ expression, keeping only the local quadratic term makes the correction inside the weighted sum equal to one. Since the normalized weights sum to one, AGQ then reduces to the Laplace approximation. This gives a direct interpretation: Laplace is the local Gaussian baseline, while AGQ adds evaluations of departures from that baseline.

Comparing the Approximations in a Joint Likelihood

Laplace approximation is usually cheaper once the mode and Hessian are available: it requires no large grid of evaluations around the mode. AGQ can improve on the local Gaussian approximation when the integrand is skewed or otherwise nonquadratic, but tensor-product cost grows rapidly with \(q\) and node count. Neither approximation is uniformly best. In particular, AGQ’s added flexibility is useful only if the quadrature rule covers the important mass of the integrand.

For either method, useful diagnostics include checking that the optimizer finds a stable interior mode, that the negative Hessian is positive definite, and that results are stable under tighter optimization tolerances. For AGQ, compare estimates or log-likelihood values at increasing node counts where feasible. For Laplace, compare with AGQ or another reliable integration method on a representative subset of subjects if computationally possible. Large changes with modest increases in quadrature order, unstable modes, or nearly singular curvature indicate that the approximation or model parameterization needs further attention.

Common Pitfalls and Practical Limits

• Counting observations instead of effects. The integration dimension is \(q\), the length of the random-effects vector. More repeated measurements increase the cost of evaluating the conditional likelihood but do not, by themselves, add integration dimensions.
• Assuming a good mode guarantees a good Laplace approximation. Laplace depends on the integrand being adequately represented by its local quadratic shape over the region that matters to the integral.
• Assuming more AGQ nodes always solve the problem. A tensor rule can become infeasible as \(q\) grows, and a poorly located adaptive region may not cover important probability mass.
• Using a non-positive-definite Hessian in the stated formula. If \(H_i\) is not positive definite at the proposed mode, the local Gaussian integral is not defined as written. Check optimization, parameterization, and model identifiability rather than forcing a determinant into the formula.
• Confusing numerical accuracy with scientific validity. Both methods approximate the likelihood under the chosen model. Neither addresses an incorrect random-effects distribution, unmodeled residual dependence, or inappropriate censoring assumptions.

Monte Carlo integration is a broad alternative when deterministic tensor grids are too costly in higher dimensions, though it has its own simulation error and computational demands. The practical choice depends on the number and structure of random effects, the shape of the integrand, and the precision needed for the inferential task. These methods approximate the marginal likelihood; they do not change the interpretation of the marker-to-hazard associations specified by the joint model.