Why Model Two Biomarkers Together?
A clinical event may be related to more than one evolving biomarker. For example, a patient’s inflammatory burden and cardiac function may both be relevant to event risk, and their underlying levels may be associated. Fitting separate survival models can obscure their conditional contributions: each model may attribute to one marker information that is also carried by the other. A joint model represents both longitudinal processes and relates their estimated current values to the event hazard in one framework.
The focus here is a two-marker model with correlated subject-specific random intercepts and marker-specific current-value hazard effects. The random-intercept correlation describes how subjects’ underlying marker levels vary together. The two hazard association parameters describe how each marker’s latent current value relates to event risk, conditional on the other marker and included covariates. These are different quantities.
Model Definition and Assumptions
The design vector \(X_{ik}(t)\) can contain time and other longitudinal predictors. In the random-intercept version above, subjects differ in their underlying marker levels, while the fixed-effects part determines the mean trajectory. The covariance matrix \(D\) must be positive semidefinite; in particular, \(|\tau_{12}|\leq \tau_1\tau_2\). The random-intercept correlation is \(\rho_b=\tau_{12}/(\tau_1\tau_2)\) when both variances are positive.
A standard specification assumes the random-effect vectors are independent across subjects and that, conditional on them and the modeled covariates, the longitudinal errors have mean zero and the stated residual covariance structure. A simple version uses independent errors across observations and markers, with marker-specific residual variances. This is an assumption to check, not a consequence of allowing correlated random intercepts: correlation in the subject-level intercepts does not automatically account for residual dependence between measurements taken at the same visit.
The survival submodel assumes the event hazard at time \(t\) depends on the two latent marker values at that time, not directly on their observed noisy measurements. It also assumes the event and longitudinal processes are conditionally independent given the random effects, covariates, and latent trajectories. Censoring and the observation process must be handled under appropriate noninformative assumptions, conditional on modeled information. The model does not by itself make the marker effects causal.
How to Read the Association Parameters
Holding the other marker, covariates, and latent subject effects fixed, a one-unit increase in marker 1 multiplies the hazard by \(\exp(\alpha_1)\); a one-unit increase in marker 2 multiplies it by \(\exp(\alpha_2)\). The units are those of the markers as modeled. Thus, changing marker scaling changes the numerical coefficient and its unit, although the implied hazard contrast for the same physical change remains the same.
If both latent marker values change, the hazard ratio follows directly from the linear predictor. Comparing two otherwise identical latent profiles at time \(t\), with differences \(\Delta m_1\) and \(\Delta m_2\), gives
This is a conditional comparison under the fitted model. It is not generally a marginal population hazard ratio, because averaging over correlated random effects does not in general preserve the same hazard-ratio form.
Worked Example: Specifying a Two-Marker Clinical Model
Consider a hypothetical cohort in which inflammatory marker \(Y_1\) is measured in mg/L and cardiac function \(Y_2\) is measured in percentage points. A few hypothetical observations are shown below. The time unit is months, and “event” indicates whether the event occurred during follow-up.
| Subject | Month | Inflammatory marker (mg/L) | Cardiac function (percentage points) | Follow-up status |
|---|---|---|---|---|
| A | 0 | 8.0 | 42 | Event at month 9 |
| A | 3 | 10.0 | 39 | Event at month 9 |
| B | 0 | 4.0 | 55 | Censored at month 12 |
| B | 3 | 3.5 | 56 | Censored at month 12 |
| C | 0 | 6.0 | 48 | Event at month 7 |
| C | 3 | 7.0 | 46 | Event at month 7 |
One possible pair of longitudinal mean structures is \(m_{i1}(t)=\beta_{10}+\beta_{11}t+b_{i1}\) and \(m_{i2}(t)=\beta_{20}+\beta_{21}t+b_{i2}\). The random intercepts can be correlated through \(\tau_{12}\), so a subject whose underlying inflammatory level is high may tend, on average, to have either higher or lower underlying cardiac function. The direction is estimated from the data; the table alone is not enough to establish it reliably.
The paired data contribute to estimating both marker trajectories and their random-effect covariance, while event and censoring information contributes through the survival likelihood. For subject A, the observed values rise in the first marker and fall in the second. The model does not plug those two noisy measurements directly into the hazard; it estimates the latent trajectories and evaluates their current values at event times. The table illustrates the structure, not a basis for claiming stable parameter estimates from only three subjects.
Deriving the Joint Contribution
The essential construction is to condition on the random effects first, then integrate them out. Let \(T_i\) be the observed follow-up time and \(\delta_i\) indicate an event (\(\delta_i=1\)) rather than right censoring (\(\delta_i=0\)). Write \(L_i^{\mathrm{long}}(b_i)\) for the conditional longitudinal-data likelihood and \(S_i(t\mid b_i)\) for survival under the conditional hazard. The subject’s joint likelihood contribution is
Here \(\phi_2(b_i;0,D)\) is the bivariate normal density of the random intercepts. The event factor is the conditional hazard at the observed time when an event occurs; the survival factor accounts for remaining event-free until that time, including for a censored subject. Integrating the product, rather than fitting each marker and survival outcome as unrelated pieces, preserves the dependence induced by the shared subject-specific effects.
The integral is generally not available in closed form for a flexible baseline hazard and nonlinear likelihood contributions. It can be evaluated by numerical integration or approximated using established likelihood or Bayesian methods. The key conceptual point is that both markers’ random intercepts are integrated jointly under \(D\), not as two unrelated one-dimensional distributions.
Worked Example: Interpreting Two Fitted Hazard Effects
Suppose a hypothetical joint-model fit gives \(\hat\alpha_1=0.08\) per mg/L for the inflammatory marker and \(\hat\alpha_2=-0.05\) per percentage point for cardiac function. These values are illustrative. Consider two latent profiles at the same time and with equal covariates, where the second profile has inflammatory marker value 2 mg/L higher and cardiac function 4 percentage points lower.
The conditional log hazard ratio is
Therefore the conditional hazard ratio is
Under this fitted model, the profile with both the higher inflammatory marker and the lower cardiac function has an estimated instantaneous event hazard about 1.43 times that of the comparison profile, conditional on the other modeled quantities. This is a hazard comparison, not a statement that the event probability over a specified period is 43.3% higher.
For a one-unit increase in the inflammatory marker alone, holding the other marker fixed, the hazard ratio is \(\exp(0.08)=1.083\), or about an 8.3% higher conditional hazard per mg/L. For a one-percentage-point increase in cardiac function alone, it is \(\exp(-0.05)=0.9512\), or about a 4.88% lower conditional hazard. Interpretation depends on the specified units and the model’s linear current-value association.
What the Correlation Does—and Does Not—Mean
The parameter \(\rho_b\) describes the association between subjects’ latent baseline levels after accounting for the fixed-effects structures. It is not the correlation between the two observed marker columns, which can also reflect time trends, residual variation, and covariates. Nor is it the correlation between \(\hat\alpha_1\) and \(\hat\alpha_2\). The association parameters concern the survival submodel; their sampling covariance is part of the uncertainty in estimation.
A correlated random-intercept distribution can improve the model when subjects’ underlying levels co-vary. However, if the two latent trajectories are highly collinear, the data may not clearly separate their hazard associations. In that case, estimates of \(\alpha_1\) and \(\alpha_2\) can be imprecise or sensitive to model choices even when a combined risk contrast is more stable. Marker scaling, considered in the earlier material on scaling markers before estimating association effects, makes coefficients easier to compare numerically but does not create information that the data lack.
Worked Example: Correlated Intercepts Versus Independent Intercepts
Suppose a hypothetical fitted model estimates \(\hat\tau_1=2.0\) mg/L, \(\hat\tau_2=5.0\) percentage points, and \(\hat\tau_{12}=-6.0\) in the corresponding product units. Then
The estimated random-intercept correlation is \(-0.60\): subjects with higher underlying inflammatory levels tend to have lower underlying cardiac function in this model. The covariance is admissible because \(|-6.0|\leq(2.0)(5.0)=10.0\). This describes between-subject dependence in latent levels; it does not say that raising inflammation causes cardiac function to decline.
An independent-intercept model would set \(\tau_{12}=0\). That may be a useful simpler alternative if the data support it, but it imposes a substantive restriction on the latent marker distribution. A shared random-effect or shared-marker approach imposes a different structure: it can represent common variation, but it need not preserve the two distinct marker-specific current-value effects specified here. The choice should follow the scientific question and model diagnostics rather than convenience alone.
Diagnostics, Failure Modes, and Limitations
Check each marker’s longitudinal fit and the survival component, as well as their joint behavior. Useful checks include comparing observed trajectories with fitted subject-specific and population-average trajectories; inspecting marker-specific residual patterns over time; checking whether residual spread changes with fitted values or time; and assessing whether event patterns are compatible with the assumed baseline hazard and current-value relationship. Cross-marker residual checks are particularly important: correlated random intercepts do not excuse unexplained residual correlation. The earlier material on cross-marker residual checks in joint models treats that diagnostic question directly.
Other warning signs include estimates that change substantially under reasonable scaling or baseline-hazard choices, very wide intervals for the two \(\alpha\) parameters, unstable covariance estimates near the boundary of the positive-semidefinite region, and poor numerical integration or convergence. A linear current-value effect can also be inadequate if risk changes nonlinearly with a marker. A random-intercept model may miss subject-specific trends if slopes vary meaningfully. These are model-fit problems, not reasons to interpret unstable coefficients as precise biological effects.
Compared with fitting separate one-marker joint models, the two-marker model estimates conditional associations while accounting for both markers together. Separate fits can be useful descriptive comparisons, but their marker effects do not answer the same conditional question and can be affected by omission of the other marker. Compared with a model imposing independent random intercepts, the correlated version is more flexible but estimates an additional covariance parameter. Compared with a single shared latent marker, it retains marker-specific hazard coefficients, at the cost of potentially greater collinearity and estimation complexity.
The interpretation remains conditional and model-based. Biomarker measurements may be affected by treatment changes, informative visit timing, unmeasured confounding, or measurement error beyond the specified residual model. A joint association does not establish a causal mechanism, and a fitted conditional hazard ratio does not directly give an individual’s absolute event probability. Results should be reported with units, uncertainty intervals, model assumptions, and appropriate cautions about identification.
Common Pitfall
Correction: The correlation describes an estimated between-subject association between latent marker levels. It does not establish direction or causation. The survival coefficient is a conditional association with the hazard under the specified joint model, not automatically a causal effect.
Check Your Understanding
Answer each question, then compare with the solutions below.
- If \(\alpha_1=0.10\) per mg/L, what is the conditional hazard ratio for a 3 mg/L increase in marker 1 with marker 2 fixed?
- If marker 1 increases by 1 mg/L and marker 2 decreases by 2 percentage points, with \(\alpha_1=0.10\) and \(\alpha_2=-0.04\), what is the combined conditional hazard ratio?
- What does a negative random-intercept correlation mean, and what does it not establish?
- Why does allowing correlated random intercepts not automatically account for cross-marker residual dependence?
Solutions. (1) The log hazard contrast is \(0.10(3)=0.30\), so the hazard ratio is \(\exp(0.30)=1.350\). (2) The log contrast is \(0.10(1)+(-0.04)(-2)=0.18\); thus the ratio is \(\exp(0.18)=1.197\). (3) It indicates that subjects with higher latent levels of one marker tend to have lower latent levels of the other, under the model; it does not establish causation or a temporal effect of one marker on the other. (4) Random-intercept covariance describes between-subject variation in latent levels, whereas residual dependence concerns within-subject deviations around those levels. They are distinct parts of the model and require separate assessment.