Numerical Integration for Competing-Risks Likelihoods

Why the Likelihood Requires Integration

In a competing-risks joint model, a subject’s event hazards depend on random effects that also connect the longitudinal and event processes. Conditional on those effects, the subject’s likelihood contribution can often be evaluated directly. For maximum likelihood estimation, however, the random effects are unobserved. The contribution used for estimation is therefore an average of the conditional likelihood over their distribution. That average is an integral, and for realistic hazard functions it usually has no convenient closed form.

Gaussian quadrature approximates this integral by evaluating the conditional contribution at a finite set of carefully chosen points and taking a weighted sum. The approximation is deterministic: unlike Monte Carlo integration, it does not rely on random draws. Its accuracy depends on the number and placement of nodes, the dimension and distribution of the random effects, and how sharply the likelihood changes over their range.

Definition: Let $b_i$ be the vector of subject $i$’s shared random effects, with $b_i\sim N_q(0,D)$. If $L_i(b_i;\theta)$ is the subject’s likelihood contribution conditional on $b_i$, then the marginal contribution is $$ L_i(\theta)=\int_{\mathbb{R}^q} L_i(b;\theta)\,\phi_q(b;0,D)\,db, $$ where $\theta$ denotes the model parameters and $\phi_q(\cdot;0,D)$ is the $q$-dimensional normal density with mean zero and covariance matrix $D$. The marginal log likelihood for independent subjects is $\ell(\theta)=\sum_i\log L_i(\theta)$.

For example, conditional on $b_i$, a subject who experiences cause $k$ at time $T_i$ contributes an event hazard for cause $k$ at $T_i$, multiplied by survival from every cause up to that time. A right-censored subject contributes survival from every cause through the censoring time, without an event-hazard factor. In a joint model, the conditional contribution also includes the longitudinal-data likelihood. The integral averages that full conditional contribution over the same random-effects distribution.

The usual modeling conditions still matter: the specified random-effects distribution and cause-specific hazard models must be appropriate; the conditional likelihood must correctly express the assumed dependence structure; and censoring must be noninformative conditional on the modeled covariates and random effects. Quadrature approximates the integral implied by these assumptions. It cannot make misspecified assumptions correct.

Gauss–Hermite Quadrature for Normal Random Effects

The standard Gauss–Hermite rule approximates integrals weighted by $\exp(-x^2)$. To use it for a standard normal random variable $Z$, first write its density explicitly and change variables from $z$ to $x=z/\sqrt{2}$:

E\{f(Z)\} =\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}f(z)e^{-z^2/2}\,dz =\frac{1}{\sqrt{\pi}}\int_{-\infty}^{\infty}f(\sqrt{2}x)e^{-x^2}\,dx \approx \frac{1}{\sqrt{\pi}}\sum_{r=1}^{m}w_r f(\sqrt{2}x_r).

Here $x_r$ and $w_r$ are the nodes and weights for the standard Gauss–Hermite rule. Equivalently, define standard-normal nodes $a_r=\sqrt{2}x_r$ and normalized weights $\tilde w_r=w_r/\sqrt{\pi}$; then the approximation is $\sum_r\tilde w_r f(a_r)$. For two nodes, the standard-normal nodes are $-1,1$, each with weight $1/2$. For three nodes, they are $-\sqrt{3},0,\sqrt{3}$, with weights $1/6,2/3,1/6$, respectively.

Key result: For $b\sim N_q(0,D)$, choose a matrix $A$ such that $AA^\mathsf{T}=D$, for example a Cholesky factor, and write $b=Az$, where the components of $z$ are independent standard normal variables. A tensor-product rule with $m$ nodes per dimension approximates the marginal likelihood by $$ L_i(\theta)\approx \sum_{r_1=1}^{m}\cdots\sum_{r_q=1}^{m} \left(\prod_{j=1}^{q}\tilde w_{r_j}\right) L_i\!\left(A(a_{r_1},\ldots,a_{r_q})^\mathsf{T};\theta\right). $$

The transformation preserves correlation: $b=Az$ has covariance $AA^\mathsf{T}=D$. The product rule evaluates the conditional contribution at every combination of the one-dimensional nodes. Its number of evaluations is $m^q$, so the computational cost can grow quickly with the number of random effects. Gauss–Hermite rules integrate polynomials of limited degree exactly, but a competing-risks likelihood generally includes exponentials of hazards and is not a polynomial. The rule is then an approximation, not an exact identity.

Worked Example: A Two-Node Likelihood for a Cause 1 Event

Suppose one hypothetical subject has a cause 1 event at time $T=2$, with no longitudinal factor included in this simplified illustration. Let $b\sim N(0,1)$, and suppose the conditional cause-specific hazards are constant over follow-up and given by $\lambda_1(t\mid b)=0.2\exp(0.5b)$ and $\lambda_2(t\mid b)=0.1\exp(-0.3b)$, in inverse time units. The conditional event contribution is $g(b)=\lambda_1(2\mid b)\exp\{-2[\lambda_1(b)+\lambda_2(b)]\}$. This uses the event hazard for cause 1 and survival against both causes. We approximate its expectation using the two standard-normal nodes $-1,1$, each weighted $1/2$.

At $b=1$, the hazards are $0.2e^{0.5}\approx0.3297$ and $0.1e^{-0.3}\approx0.07408$. Thus $g(1)=0.2e^{0.5}\exp\{-2[0.2e^{0.5}+0.1e^{-0.3}]\}\approx0.147034$, or $0.1470$ to four significant figures.

At $b=-1$, the hazards are $0.2e^{-0.5}\approx0.1213$ and $0.1e^{0.3}\approx0.1350$. Therefore $g(-1)=0.2e^{-0.5}\exp\{-2[0.2e^{-0.5}+0.1e^{0.3}]\}\approx0.07266$. The two-node approximation to the marginal contribution is $\tfrac12[g(-1)+g(1)]\approx\tfrac12(0.07266+0.147034)=0.1098$. This is an approximation to the event contribution averaged over the random-effects distribution, not a conditional contribution at an “average” random effect.

Correlated Effects and Tensor-Product Nodes

When the random effects are correlated, applying independent quadrature directly to their components would generally use the wrong distribution. Transforming independent standard-normal coordinates by a factor of $D$ handles this correlation. The next example shows the transformation and the weighted sum for a censored subject.

Worked Example: Integrating a Correlated Two-Effect Survival Contribution

Suppose a hypothetical subject is event-free through time 1. Let $b=(b_1,b_2)^\mathsf{T}\sim N_2(0,D)$, where $D=\begin{pmatrix}1&0.6\\0.6&1\end{pmatrix}$. For this calculation, let the constant cause-specific hazards be $\lambda_1=0.2e^{b_1}$ and $\lambda_2=0.1e^{b_2}$. The conditional likelihood contribution is $s(b)=\exp[-0.2e^{b_1}-0.1e^{b_2}]$. A Cholesky factor is $A=\begin{pmatrix}1&0\\0.6&0.8\end{pmatrix}$, so $b_1=z_1$ and $b_2=0.6z_1+0.8z_2$. Both $z_1$ and $z_2$ use two-node quadrature, giving four equally weighted combinations.

At $(z_1,z_2)=(1,1)$, $b=(1,1)$, so $s(b)=\exp[-0.2e-0.1e]\approx0.4425$. At $(1,-1)$, $b=(1,-0.2)$, giving $s(b)=\exp[-0.2e-0.1e^{-0.2}]\approx0.5350$. At $(-1,1)$, $b=(-1,0.2)$, giving $s(b)=\exp[-0.2e^{-1}-0.1e^{0.2}]\approx0.8222$. At $(-1,-1)$, $b=(-1,-1)$, giving $s(b)=\exp[-0.2e^{-1}-0.1e^{-1}]\approx0.8955$.

Each combination has weight $1/4$, so the tensor-product approximation is $\tfrac14(0.4425+0.5350+0.8222+0.8955)\approx0.6738$. The correlation is reflected in the transformed values of $b_2$; using $b_1=z_1,\ b_2=z_2$ instead would impose independent effects and would approximate a different model.

Checking Accuracy Through Parameter Estimates

Accuracy matters not only for individual likelihood contributions but also for the maximum likelihood estimates. If the quadrature approximation is too coarse, its error can alter the shape of the approximated log likelihood, its maximizer, and the numerical derivatives used for standard errors. Comparing estimates and log likelihoods across increasing numbers of nodes is a practical sensitivity check. Agreement is reassuring, though it is not a proof of exactness; disagreement means the lower-resolution fit should not be treated as reliable.

Worked Example: A Two-Cause Model and Quadrature-Dependent Estimates

Consider a hypothetical discrete-time competing-risks model with one event in the interval and two possible causes. Conditional on a standard-normal subject effect $b$, let the probability of cause 1 be $p(b;\alpha)=\operatorname{logistic}(\alpha+b)$, where $\operatorname{logistic}(u)=1/(1+e^{-u})$; the probability of cause 2 is $1-p(b;\alpha)$. Suppose 10 independent subjects yield 6 cause 1 events and 4 cause 2 events, with each subject’s effect integrated out separately. Under a chosen quadrature rule, let $\bar p_m(\alpha)$ denote the approximation to $E[p(b;\alpha)]$. The approximated likelihood is proportional to $\bar p_m(\alpha)^6[1-\bar p_m(\alpha)]^4$, which is maximized when $\bar p_m(\alpha)=6/10=0.6$.

With two nodes, $\bar p_2(\alpha)=\tfrac12[\operatorname{logistic}(\alpha-1)+\operatorname{logistic}(\alpha+1)]$. At $\alpha=0.5127$, the two node probabilities are approximately $\operatorname{logistic}(-0.4873)=0.3805$ and $\operatorname{logistic}(1.5127)=0.8195$. Their average is approximately $0.6000$, giving a two-node estimate of $\alpha\approx0.5127$.

With three nodes, $\bar p_3(\alpha)=\tfrac16\operatorname{logistic}(\alpha-\sqrt{3})+\tfrac23\operatorname{logistic}(\alpha)+\tfrac16\operatorname{logistic}(\alpha+\sqrt{3})$. At $\alpha=0.485$, the three probabilities are approximately $0.2232,\ 0.6189,\ 0.9018$. Hence $\bar p_3(0.485)\approx\tfrac16(0.2232)+\tfrac23(0.6189)+\tfrac16(0.9018)=0.6001$, so the three-node estimate is approximately $0.485$. The estimates differ appreciably in this illustration: using only two nodes would not adequately reproduce the three-node result. In an actual fit, one should continue increasing the node count and check whether estimates and relevant likelihood calculations stabilize.

Implementation, Diagnostics, and Alternatives

In practice, quadrature is repeated many times as the optimizer explores different parameter values. The marginal likelihood is a sum of positive weighted contributions; computing its logarithm directly from very small contributions can cause underflow. A stable calculation works with log contributions and uses the log-sum-exp identity: for values $c_r$, $\log\sum_r e^{c_r}=c_{\max}+\log\sum_r e^{c_r-c_{\max}}$, where $c_{\max}=\max_r c_r$. Quadrature weights can be included in $c_r$ by adding their logarithms.

Useful checks include refitting with successively larger node counts, comparing parameter estimates and maximized log likelihoods, and checking whether the optimizer converges to the same solution. Large changes in estimates, standard errors, or likelihood across refinements suggest that integration accuracy may be inadequate. A failure to stabilize can also signal difficult curvature, a poorly scaled parameterization, or optimization problems; it is not automatically evidence that the model itself is wrong. The number of nodes needed can vary across parameters and data sets, especially when random effects are strongly associated with hazards or the likelihood is concentrated in a region that the fixed nodes represent poorly.

Adaptive Gauss–Hermite quadrature centers and scales nodes around a subject-specific region of high conditional likelihood. It can improve efficiency when ordinary nodes are poorly placed, but requires reliable centering and curvature calculations. The Laplace approximation instead uses a local quadratic approximation around an integral’s mode; it can be faster in higher dimensions, but may be inaccurate when the integrand is markedly non-Gaussian. Monte Carlo integration uses random draws and is flexible in higher dimensions, but introduces simulation error that must itself be controlled. These are numerical alternatives for evaluating integrals, not substitutes for specifying the correct conditional likelihood and random-effects distribution.

Approach	Strength	Accuracy consideration
Fixed Gauss–Hermite quadrature	Deterministic weighted sums; effective in low dimensions	Tensor-product cost grows rapidly with dimension; check node stability
Adaptive Gauss–Hermite quadrature	Places nodes near a subject-specific likelihood region	Depends on accurate mode and curvature calculations
Laplace approximation	Often computationally efficient	Relies on a local quadratic approximation to the integrand
Monte Carlo integration	More readily extended to higher dimensions	Has simulation error and may require many draws

A common mistake is to assume that adding nodes must move an estimate steadily toward a known answer, or that a small change between two resolutions proves accuracy. Quadrature error need not decrease monotonically for every likelihood, and two nearby results can still share appreciable error. The correction is to examine several resolutions where feasible, assess stability of the quantities that matter, and report the integration settings. Numerical refinement addresses integration error only; it does not eliminate sampling uncertainty or model misspecification.

Key takeaway: Marginal competing-risks likelihoods average conditional contributions over shared random effects. Gaussian quadrature replaces that average with a weighted sum, using a covariance-preserving transformation for correlated effects. Reliable estimation requires checking that the numerical approximation is sufficiently stable.

Check Your Understanding

Answer each question, then compare with the brief solutions below.

For a two-node standard-normal rule with nodes $-1,1$, what is the quadrature approximation to $E[f(Z)]$?
Why must a correlated random-effects vector be transformed using a factor $A$ satisfying $AA^\mathsf{T}=D$?
A subject is censored at time $C$ with two cause-specific hazards. Which conditional likelihood factors are included?
What does a substantial change in parameter estimates when the number of quadrature nodes increases indicate, and what does it not establish?
Name one numerical alternative to fixed Gauss–Hermite quadrature and one limitation of that alternative.

Solutions. (1) The approximation is $\tfrac12[f(-1)+f(1)]$. (2) If $z$ has independent standard-normal components, then $Az$ has covariance $AA^\mathsf{T}=D$, preserving the specified joint normal distribution. (3) Include survival against both causes through $C$, with no event-hazard factor because the subject had no observed event. (4) It indicates that the lower-resolution integration may be inadequate for the fitted quantities; by itself, it does not establish that the statistical model is misspecified. (5) For example, the Laplace approximation can be faster, but may be inaccurate when the integrand is not well approximated by a local quadratic form.