Worked Two-Cause Joint Model for Disease Events

Why Build a Two-Cause Joint Model?

Suppose a biomarker is measured repeatedly while participants remain at risk of either of two mutually exclusive disease events. The biomarker may change over time, and participants with different underlying biomarker trajectories may also have different risks of each event. Analyzing the measurements and event times separately can miss that dependence, especially when an event stops future biomarker measurement.

A two-cause joint model links the repeated measurements to both cause-specific hazards through subject-level random effects. The two hazards describe the instantaneous rates of the two event types among participants who have experienced neither event yet. Their combination also determines the probability of experiencing each cause by a specified time. These associations are model-based; they do not, by themselves, establish causal effects.

Model Specification and Assumptions

For participant \(i\), let \(Y_{ij}\) be the biomarker measured at time \(t_{ij}\), and let \(m_i(t)\) be the participant's latent mean biomarker trajectory. A simple linear mixed model is

Y_{ij}=m_i(t_{ij})+\epsilon_{ij}, \qquad m_i(t)=\beta_0+\beta_1t+b_{0i}+b_{1i}t.

Here \(\beta_0\) and \(\beta_1\) describe the population mean trajectory, while \(b_i=(b_{0i},b_{1i})^\mathsf{T}\) describes participant-specific deviations. A common specification assumes \(b_i\sim N(0,D)\) and \(\epsilon_{ij}\sim N(0,\sigma^2)\), with measurement errors conditionally independent given \(b_i\). Other longitudinal mean and error structures are possible, but they must be specified and checked for the application.

Let \(T_i\) be time to the first event, and let \(J_i\in\{1,2\}\) denote its cause. The cause-specific hazard for cause \(k\), conditional on the random effects and covariates, is

\lambda_{ik}(t\mid b_i) = \lambda_{k0}(t)\exp\left\{w_i^\mathsf{T}\gamma_k+\alpha_k m_i(t)\right\}, \qquad k=1,2.

In this expression, \(\lambda_{k0}(t)\) is the cause-specific baseline hazard, \(w_i\) is a vector of event-model covariates, and \(\gamma_k\) is its cause-specific coefficient vector. The coefficient \(\alpha_k\) links the latent biomarker level to cause \(k\). Because \(m_i(t)\) contains \(b_i\), the same random effects that govern the biomarker trajectory also enter both event hazards. Thus the model induces dependence between the repeated measurements and both event processes.

Definition: A two-cause joint model combines (1) a probability model for repeated biomarker measurements, (2) cause-specific hazards conditional on the longitudinal random effects, and (3) a shared distribution for those random effects. A cause-specific hazard is the instantaneous event rate for one cause among people still event-free immediately before time \(t\).

The interpretation and likelihood require several assumptions. The two causes are mutually exclusive first events, and the event type is observed correctly. Conditional on the modeled covariates and random effects, the longitudinal errors and event process follow their stated distributions. Censoring is noninformative conditional on the modeled information, so its mechanism does not add a parameter-dependent contribution to the event likelihood. The chosen trajectory, baseline hazards, and association structure must be adequate; for example, the model above links hazards to current latent biomarker level, not to its slope or cumulative history.

Worked Example: Interpreting the Two Biomarker Associations

Suppose biomarker values are recorded in standardized units and the fitted cause-specific association estimates are \(\alpha_1=0.40\) and \(\alpha_2=-0.20\). Interpret a one-unit higher latent biomarker level, holding the other event covariates fixed.

For cause 1, the hazard ratio for a one-unit increase in \(m_i(t)\) is \( \exp(\alpha_1)=\exp(0.40)=1.492 \). Under the model, this is about a 49.2% higher instantaneous cause-1 hazard, conditional on being event-free and on the modeled covariates and random effects.

For cause 2, the corresponding ratio is \( \exp(\alpha_2)=\exp(-0.20)=0.8187 \). This is about an 18.1% lower instantaneous cause-2 hazard under the same conditioning. These are cause-specific associations, not statements that changing the biomarker would cause either hazard to change.

From Cause-Specific Hazards to Cumulative Incidence

Conditional on \(b_i\), the overall hazard of either event is the sum of the two cause-specific hazards. Therefore, the probability of remaining free of both causes through time \(t\), conditional on \(b_i\), is

S_i(t\mid b_i) = \exp\left\{ -\int_0^t \left[\lambda_{i1}(u\mid b_i)+\lambda_{i2}(u\mid b_i)\right]du \right\}.

A cause-specific cumulative incidence is the probability of having experienced cause \(k\) by time \(t\), before experiencing the other cause. Conditional on the random effects it is

F_{ik}(t\mid b_i) = \int_0^t S_i(u\mid b_i)\lambda_{ik}(u\mid b_i)\,du.

This expression follows by considering a short interval of length \(du\). The probability of being event-free just before \(u\) is \(S_i(u\mid b_i)\); given that status, the probability of cause \(k\) in that interval is approximately \(\lambda_{ik}(u\mid b_i)du\). Multiplying gives the cause-\(k\) event probability in the interval, and integrating adds these probabilities from time 0 to \(t\). The survival term includes both hazards: a participant who has already experienced the competing cause cannot subsequently experience cause \(k\) as the first event.

For a population-level prediction, random effects are not fixed at zero. The cumulative incidence averages over their distribution:

F_k(t\mid w) = \int F_k(t\mid b,w) f(b)\,db.

In practice this average, and the time integral in the conditional cumulative incidence, may require numerical integration. A fitted conditional association does not directly give a marginal risk difference: the competing hazard and the distribution of random effects both contribute to the cumulative incidence.

Worked Example: Calculating Both Cumulative Incidences

Consider a hypothetical participant whose biomarker and covariates make the two conditional hazards constant over a two-year period: \(\lambda_1=0.12\) per year and \(\lambda_2=0.06\) per year. Assume there are no other event causes and no delayed entry. Calculate the two-year survival probability and cumulative incidences.

The total hazard is \(0.12+0.06=0.18\) per year, so the probability of remaining free of either event at two years is \( S(2)=\exp[-(0.18)(2)]=\exp(-0.36)=0.6977 \).

Because the hazards are constant, for cause 1, \( F_1(2)=\int_0^2 \exp(-0.18u)(0.12)\,du \). Evaluating the integral gives \( F_1(2)=\frac{0.12}{0.18}[1-\exp(-0.36)]=0.2015 \). For cause 2, \( F_2(2)=\frac{0.06}{0.18}[1-\exp(-0.36)]=0.1008 \).

Thus the conditional two-year probabilities are about 20.2% for cause 1, 10.1% for cause 2, and 69.8% for remaining event-free. They sum to 1, up to rounding. Treating the competing event as if it did not exist would not give these cause-specific probabilities, because the competing hazard reduces the time available for either first event.

Likelihood Contribution and Shared-Effect Integration

Let \(\tilde T_i=\min(T_i,C_i)\) be the observed follow-up time, where \(C_i\) is censoring time. Let \(\delta_i=1\) if an event is observed and 0 if the participant is censored; when \(\delta_i=1\), \(J_i\) records the event cause. Conditional on \(b_i\), the event contribution is

L_{Ti}(b_i) = \left\{\prod_{k=1}^2 \lambda_{ik}(\tilde T_i\mid b_i)^{I(J_i=k)}\right\}^{\delta_i} S_i(\tilde T_i\mid b_i).

For an observed cause-\(k\) event, this is survival without either event up to the observed time multiplied by the cause-\(k\) hazard at that time. For right-censoring, it is just the probability of surviving event-free to the censoring time. If \(f(Y_i\mid b_i)\) is the longitudinal-data density and \(f(b_i)\) is the random-effects density, the participant's likelihood contribution is

L_i = \int f(Y_i\mid b_i)L_{Ti}(b_i)f(b_i)\,db_i.

The integration is essential: the random effects are unobserved, and the longitudinal observations and event outcome jointly provide information about them. Substituting an estimated random effect as if it were known generally does not produce this marginal likelihood.

Worked Example: A Conditional Event Contribution

For a hypothetical participant, suppose the latent biomarker is constant at one standardized unit during follow-up, and the cause-specific hazards at that level are constant at \(0.12\) and \(0.06\) per year. A cause-1 event occurs at year 2. At a fixed value of the random effects, calculate the event contribution.

The conditional survival probability is \(S(2\mid b)=\exp[-(0.12+0.06)(2)]=0.6977\). Since the observed event is cause 1, the event contribution is \(S(2\mid b)\lambda_1(2\mid b)=0.6977(0.12)=0.08372\) per year.

For illustration, suppose three biomarker observations have conditional residuals \(0, 0.2,-0.1\) standardized units, with independent normal errors and \(\sigma=0.5\). Their conditional density product is \( \prod_{e\in\{0,0.2,-0.1\}} \frac{1}{\sqrt{2\pi}(0.5)}\exp[-e^2/(2(0.5)^2)] \) \(=(0.7979)(0.7365)(0.7821)=0.4596 \). The complete-data contribution at this fixed \(b\), before multiplying by the random-effects density, is therefore \(0.4596(0.08372)=0.03848\), in the corresponding measurement-density units times per year. The model likelihood contribution for the participant instead integrates this product over possible \(b\), weighted by \(f(b)\).

Estimation, Diagnostics, and Alternatives

The likelihood commonly requires numerical integration over the random effects; the earlier material on numerical integration for competing-risks likelihoods is relevant to this computational step. Baseline hazards may be represented parametrically or piecewise, and their form should be selected with the time scale and event data in mind. After fitting, evaluate whether the model is describing both parts of the data adequately, rather than judging it only by convergence or by whether the association estimates appear plausible.

Check the longitudinal component. Examine residual patterns over time and across fitted values, and assess whether the assumed trajectory and error distribution are credible. Systematic curvature may indicate that a linear mean trajectory is inadequate.
Check event-model fit by cause. Compare predicted and observed event patterns separately for each cause. A poor fit for a sparse cause can make its association estimate unstable.
Check the association structure. Consider whether current latent level is scientifically and empirically defensible. A model using current level can be misspecified if risk is better related to another feature of the trajectory.
Check computation and sensitivity. Monitor convergence, integration accuracy, and sensitivity to baseline-hazard and random-effects specifications. Large changes under reasonable alternatives signal that results are model-dependent.

Separate cause-specific models can describe associations with each cause, but fitting them apart from the longitudinal model does not represent the shared dependence through the latent trajectory in the same way. A subdistribution-hazard approach targets a different summary of competing-risk incidence; it should not be interpreted as the same hazard or as a substitute for the joint likelihood when the target is a biomarker-linked longitudinal and cause-specific model. The joint model itself has limitations: it depends on distributional and structural assumptions, may be difficult to estimate with few events for one cause, and does not turn an observational association into a causal effect.

Common pitfall: It is incorrect to calculate cause-1 cumulative incidence by integrating its hazard as though cause 2 did not exist. That calculation omits the chance of experiencing cause 2 first. Use \(S_i(u\mid b_i)\), which is based on the sum of both cause-specific hazards, inside the cause-specific cumulative-incidence integral.

Check Your Understanding

Try each problem, then compare your reasoning with the solutions below.

If \(\alpha_1=0.30\), what is the hazard ratio for a one-unit higher latent biomarker level for cause 1?
With constant hazards \(\lambda_1=0.10\) and \(\lambda_2=0.10\) per year, write the expression for the cause-1 cumulative incidence at time \(t\).
Why does the likelihood contribution for a cause-2 event include \(\lambda_2(\tilde T_i\mid b_i)\) rather than \(\lambda_1(\tilde T_i\mid b_i)\)?
What does integrating over \(b_i\) accomplish in the participant likelihood?
Give one reason a fitted model could yield an unstable estimate of \(\alpha_2\).

Solutions. (1) The ratio is \(\exp(0.30)=1.350\), so the modeled cause-1 hazard is 35.0% higher per one-unit increase, conditional on the other modeled quantities. (2) The total hazard is \(0.20\) per year, so \(F_1(t)=\int_0^t e^{-0.20u}(0.10)\,du=\frac{0.10}{0.20}[1-e^{-0.20t}]\). (3) An observed cause-2 event contributes its own cause-specific hazard; the survival factor accounts for remaining free of both causes before that event. (4) The random effects are latent, so integration averages the joint longitudinal and event contribution over their modeled distribution instead of treating an estimated effect as known. (5) Few cause-2 events provide limited information about its association, making \(\alpha_2\) poorly determined and potentially sensitive to modeling choices.