Why Joint-Model Likelihoods Need Numerical Integration

Why Marginalization Becomes a Numerical Problem

In a shared random-effects model, a subject’s longitudinal measurements and event time are conditionally related through latent subject-specific effects. The likelihood for that subject is obtained by integrating the conditional likelihood over the random-effects distribution. This step is essential: the observed-data likelihood must account for uncertainty about the unobserved effects, rather than treating them as if they were known.

The longitudinal contribution may have a familiar density, and the event contribution may have a simple conditional form. Yet their product, multiplied by a Gaussian random-effects density, is generally not a Gaussian kernel. The main obstacle is often the survival contribution: a random effect enters a hazard, the hazard is integrated over time to form a cumulative hazard, and the resulting survival term is an exponential of that cumulative hazard. Integrating this nonlinear expression against a Gaussian density usually has no closed-form solution.

Definition: Let $b_i$ be the $q$-dimensional random-effects vector for subject $i$, with $b_i\sim N_q(0,D)$, where $D$ is a positive-definite covariance matrix. Let $Y_i$ denote that subject’s longitudinal data, $T_i$ the observed event or censoring time, and $\delta_i$ the event indicator. If the longitudinal and event processes are conditionally independent given $b_i$ and the modeled covariates, the marginal subject likelihood is $$ L_i(\theta)=\int_{\mathbb{R}^q} f_Y(Y_i\mid b_i;\theta)\, h_i(T_i\mid b_i;\theta)^{\delta_i} S_i(T_i\mid b_i;\theta)\, \phi_q(b_i;0,D)\,db_i . $$ Here $h_i$ is the conditional hazard, $S_i(t\mid b_i)=\exp\{-H_i(t\mid b_i)\}$ is the conditional survival function, $H_i(t\mid b_i)=\int_0^t h_i(u\mid b_i)\,du$ is the cumulative hazard, and $\phi_q$ is the multivariate Gaussian density. The integral has $q$ dimensions.

This expression assumes the conditional likelihood components are correctly specified, censoring is noninformative given the modeled information, and subjects’ random effects are independent across subjects. If entry is delayed, the event contribution must also reflect conditioning on survival to the entry time. These conditions concern the model being integrated; numerical integration cannot repair a misspecified hazard, an incorrect censoring assumption, or an inappropriate random-effects distribution.

For example, a common hazard specification is $h_i(t\mid b_i)=h_0(t)\exp\{\gamma^\mathsf{T}w_i+\alpha m_i(t;b_i)\}$, where $m_i(t;b_i)$ is a longitudinal trajectory depending on the random effects. Even if $m_i(t;b_i)$ is linear in $b_i$, the hazard is exponential in those effects. Its cumulative hazard involves integrating that exponential over time, and the survival term then exponentiates the negative cumulative hazard. The full integrand is therefore not generally a Gaussian density times a linear or quadratic exponential—the special structure that often makes Gaussian integrals analytic.

What the Integration Dimension Counts

The integration dimension is the number of components in one subject’s latent random-effects vector, not the number of visits, event times, or subjects. A random intercept and a random slope give $q=2$; adding a shared frailty effect gives $q=3$. The covariance matrix $D$ may contain correlations among those components, but correlation does not add dimensions. With independent subject-specific random effects, the sample likelihood factors into a product of subject likelihoods, so each subject contributes a $q$-dimensional integral rather than one integral over all subjects’ random effects at once.

The longitudinal record can contain many visits while the integral remains low-dimensional if all those measurements depend on the same few random effects. Conversely, adding another random-effect component increases the dimension even if the subject has only a few observations. This distinction matters computationally: many integration procedures become substantially more expensive as $q$ grows.

Worked Example: A Nonlinear Survival Contribution

Consider a hypothetical subject with one standard-normal random effect, $b\sim N(0,1)$. Suppose the conditional hazard is constant in time given $b$, with $h(t\mid b)=0.4\exp(b)$ per year. We examine the probability of surviving one year, isolating the event contribution to show why Gaussian averaging is difficult.

Conditional survival. Since the hazard is constant over time for fixed $b$, $H(1\mid b)=0.4\exp(b)$, and therefore $$ S(1\mid b)=\exp\{-0.4\exp(b)\}. $$

Marginal survival. Averaging over $b$ requires $$ P(T>1)=\int_{-\infty}^{\infty} \exp\{-0.4\exp(b)\}\, \frac{1}{\sqrt{2\pi}}\exp(-b^2/2)\,db. $$ The factor $\exp\{-0.4\exp(b)\}$ is not a Gaussian kernel. In particular, the integral is not obtained by replacing $b$ with its mean, zero: that would give $\exp(-0.4)\approx0.6703$, which is a conditional survival value at $b=0$, not the marginal survival probability.

A numerical illustration. A coarse Simpson approximation on $[-4,4]$ with step size 1 uses the integrand values $g(b)=\phi(b)\exp\{-0.4\exp(b)\}$. At $b=-4,-3,\ldots,4$, these are approximately $0.000133, 0.004344, 0.05115, 0.20886, 0.26742, 0.08155, 0.002809, 0.00000144,$ and approximately $0$, respectively. Simpson’s formula gives $$ \int_{-4}^{4}g(b)\,db\ \approx\ \frac{1}{3}\left[g(-4)+g(4)+4\{g(-3)+g(-1)+g(1)+g(3)\} +2\{g(-2)+g(0)+g(2)\}\right]\approx 0.6073. $$ This is a deliberately coarse numerical illustration, not a precision estimate; a practical calculation would check accuracy with a finer integration scheme and control the omitted tails. The standard-normal probability outside $[-4,4]$ is about $0.0000634$, and because the survival factor is between zero and one, the omitted survival integral is no larger than that tail probability. The key point is that evaluating the Gaussian expectation requires integrating a nonlinear function, not merely computing a Gaussian mean.

Why a Simple Conditional Formula Does Not Settle the Integral

A closed-form conditional survival expression does not imply a closed-form marginal survival expression. In the example, the conditional survival is simple for every fixed $b$; the difficulty appears only after averaging over the Gaussian distribution. Likewise, having a closed form for $H_i(t\mid b_i)$ as a function of $b_i$ does not ensure that its exponential has an analytic Gaussian integral.

This is a general distinction between evaluating the integrand and integrating the integrand. The conditional event contribution may be straightforward to calculate at any proposed $b_i$, while the marginal likelihood still requires evaluating the area under a function over all possible random-effects values. Numerical methods approximate that integral; they do not change the likelihood being defined.

Worked Example: Three Random Effects Give a Three-Dimensional Integral

Suppose, hypothetically, that each subject has a random intercept $b_{0i}$ and random slope $b_{1i}$ for a longitudinal marker, plus a shared frailty $u_i$ in the event hazard. Let $b_i=(b_{0i},b_{1i},u_i)^\mathsf{T}\sim N_3(0,D)$, with a general covariance matrix $D$, and suppose the event hazard depends on both the marker trajectory and $u_i$.

For each subject the marginal likelihood is an integral over $b_{0i},b_{1i},u_i$: $$ L_i(\theta)=\int_{\mathbb{R}^3} f_Y(Y_i\mid b_{0i},b_{1i};\theta) h_i(T_i\mid b_{0i},b_{1i},u_i;\theta)^{\delta_i} S_i(T_i\mid b_{0i},b_{1i},u_i;\theta) \phi_3(b_i;0,D)\,db_{0i}\,db_{1i}\,du_i. $$

The random effects are correlated through $D$, but there are still three integration variables. If the hypothetical data set has 80 independent subjects with this same structure, the likelihood is $L(\theta)=\prod_{i=1}^{80}L_i(\theta)$, or $\ell(\theta)=\sum_{i=1}^{80}\log L_i(\theta)$. Thus there are 80 separate three-dimensional subject integrals, not one 240-dimensional integral, under the stated independence of subjects’ random effects. Each subject’s number of longitudinal visits changes the form of $f_Y$, but does not change the integration dimension unless the random-effects vector itself is enlarged.

When Exact Simplification Is Possible

“Generally lacks a closed form” does not mean that analytic integration is impossible in every model. A useful simplification occurs when a factor in the conditional likelihood does not depend on the random effects: that factor can be taken outside the integral. But this alone does not make the remaining integral equal to one. For example, if the hazard is independent of $b_i$, the event contribution comes outside, while the longitudinal likelihood may still depend on $b_i$ and must still be averaged against its Gaussian distribution.

An exact result is available when the remaining integral has a tractable form. One familiar case is a Gaussian longitudinal model with a Gaussian random intercept and no random-effect dependence in the event process. The marginal longitudinal density is then a Gaussian convolution. This is a special structural property, not a consequence of the event hazard being independent of the random effect by itself.

Worked Example: An Analytic Gaussian Case

Suppose a hypothetical subject has one longitudinal measurement $Y\mid b\sim N(b,0.25)$, where $b\sim N(0,1)$. The event hazard is constant at $0.2$ per year and does not depend on $b$. The subject has $Y=1.2$ and an event at $T=2$ years. Assume conditional independence of the measurement and event given $b$.

Longitudinal marginal density. Since $Y=b+\epsilon$, with independent $b\sim N(0,1)$ and $\epsilon\sim N(0,0.25)$, $Y\sim N(0,1.25)$. Thus $$ f_Y(1.2)=\frac{1}{\sqrt{2\pi(1.25)}}\exp\left\{-\frac{1.2^2}{2(1.25)}\right\} =\frac{1}{\sqrt{7.85398}}\exp(-0.576)\approx 0.2006. $$ The event density contribution is $$ h(2)S(2)=0.2\exp\{-0.2(2)\}=0.2\exp(-0.4)\approx0.1341\ \text{per year}. $$

Because the event factor is independent of $b$, it multiplies the exact marginal longitudinal density. The joint marginal density is therefore approximately $0.2006\times0.1341=0.02689$, in measurement-units$^{-1}$ year$^{-1}$. Here the Gaussian longitudinal integral is analytic. If the event hazard instead depended on $b$, its survival factor would remain inside the integral and this simple Gaussian convolution would not, in general, evaluate the full joint likelihood.

Numerical Integration: Diagnostics and Alternatives

The practical question is not whether an approximation is exact, but whether its error is small enough for the likelihood calculations and conclusions at hand. A numerical integration result should be checked rather than treated as automatically reliable. Useful checks include:

Stability under increased accuracy. Recompute the likelihood using a more accurate integration setting. Material changes in the log likelihood or parameter estimates indicate that the original calculation was not sufficiently resolved.
Behavior across subjects. Inspect whether a few subjects contribute unusually unstable or extreme log-likelihood values. Their posterior random-effects distributions may be concentrated away from the Gaussian prior’s center, making the integral harder to approximate.
Optimization stability. Check whether changing integration accuracy changes convergence, estimates, or standard errors. Small changes in the likelihood approximation can matter near a flat optimum or when comparing models.
Numerical scale. Products of small densities can underflow. Computing subject contributions on the log scale and using stable summation helps, but does not substitute for accurate integration.

The main alternatives make different compromises. Gaussian quadrature approximates an integral using weighted evaluations and is often effective at modest dimension; its particular construction is the subject of the next step in this sequence. Adaptive quadrature concentrates evaluations in regions important to the integrand. Laplace approximation uses a local expansion near the integrand’s dominant region and can be computationally attractive, but accuracy may suffer when the integrand is skewed or poorly represented by a local quadratic approximation. Monte Carlo integration averages evaluations at random draws; it handles more general integrands but introduces simulation error that must itself be assessed. None makes an incorrect conditional model correct, and no single method is uniformly best for every dimension and posterior shape.

A further limitation is dimensionality. As $q$ increases, covering the regions that matter in a multivariate integral becomes more demanding. A method that is accurate for a random intercept may be inadequate for a random intercept, slope, and several additional effects. This is why the random-effects structure is both a scientific modeling decision and a computational one: extra effects may be meaningful, but they increase the burden of marginal likelihood evaluation.

Common Pitfall

Common mistake: “The hazard does not depend on the random effects, so the random-effects integral disappears.” This is only true if the entire remaining integrand no longer depends on those effects or can be integrated analytically. In a joint model, the longitudinal density may still depend on them. Correctly, one takes the event factor outside the integral and then evaluates the remaining longitudinal expectation; the expectation equals one only when the remaining factor is itself a normalized density integrated over all its data, not generally for the observed data contribution.

The distinction is especially important when explaining why an integral is easy or difficult. State which factors depend on $b_i$, take outside only factors that do not, and then inspect the form of what remains. A Gaussian random-effects distribution does not make every expectation analytic: the function being averaged matters.

Key takeaway: The marginal subject likelihood is a Gaussian expectation of the conditional joint contribution. Nonlinear survival terms usually destroy the Gaussian-kernel structure needed for a closed-form integral. The number of random-effects components determines the integration dimension; exact simplifications require additional structure in every factor that remains inside the integral.

Check Your Understanding

Answer each question, then compare with the brief solutions.

A subject has a random intercept and random slope, with 12 longitudinal visits. What is the dimension of the random-effects integral, assuming no other random effects?
Why is $\exp\{-c\exp(b)\}$, averaged over $b\sim N(0,1)$, generally not evaluated by substituting $b=0$?
If the event hazard is independent of $b_i$, what happens to the event factor in the marginal likelihood? What may remain to be integrated?
A model adds a shared frailty to a random intercept and random slope. How does the per-subject integration dimension change?
Name one diagnostic check for a numerical marginal likelihood and one limitation of Monte Carlo integration.

Solutions. (1) The integral is two-dimensional: visits affect the longitudinal likelihood but not the number of random-effect coordinates. (2) Substitution gives the conditional value at the mean, whereas the marginal value is an average over the full distribution of $b$; a nonlinear function’s average is not generally its value at the average. (3) The event factor can be taken outside, but the longitudinal contribution may still depend on $b_i$ and require integration. (4) It increases from two to three dimensions. (5) One can check stability after increasing integration accuracy; Monte Carlo integration has simulation error, which can be substantial unless enough draws are used and assessed.)