Covariate-Dependent Association Links

Why Let the Biomarker Link Vary by Covariate?

A biomarker can have a different association with event risk in different baseline groups. For example, a one-unit increase in a biomarker might correspond to a modest change in the hazard among one group and a larger change among another. A model with a single biomarker association coefficient cannot express that difference. Adding an interaction between the biomarker association term and a baseline covariate allows the link strength to vary across groups.

The interaction is about how the biomarker relates to the hazard, not merely whether the groups have different hazards. A baseline covariate can have a direct effect on the hazard even when the biomarker association is the same in each group. Keeping these two roles separate is essential for interpretation.

Model and Interpretation

Let \(m_i(t)\) denote subject \(i\)'s underlying biomarker value at time \(t\), as represented by the longitudinal component of the joint model. Let \(x_i\) be a binary baseline covariate coded 0 or 1, and let \(w_i\) denote any other baseline covariates in the event model. A covariate-dependent current-value link can be written

h_i(t\mid b_i,x_i,w_i) = h_0(t)\exp\left\{ \gamma x_i+\beta^\top w_i+ \left(\alpha+\delta x_i\right)m_i(t) \right\}.

Here \(h_0(t)\) is the baseline hazard, \(b_i\) denotes the subject-specific random effects that help define the latent trajectory, \(\gamma\) is the direct effect of \(x_i\), and \(\delta\) is the interaction coefficient. The model is conditional on the latent trajectory and random effects. In a shared-random-effects joint model, the longitudinal and event processes are linked through those subject-specific effects and the latent biomarker trajectory.

Definition: In this model, the biomarker association coefficient is \(\alpha\) when \(x_i=0\), and \(\alpha+\delta\) when \(x_i=1\). The interaction coefficient \(\delta\) is the difference between the two subgroup-specific log-hazard associations per one-unit increase in \(m_i(t)\). The direct covariate coefficient \(\gamma\) is a separate contribution to the log hazard.

For two otherwise comparable latent trajectories in subgroup \(x\), whose biomarker values at time \(t\) differ by \(\Delta m\), the conditional hazard ratio is

\mathrm{HR}_x(\Delta m) = \exp\left\{\left(\alpha+\delta x\right)\Delta m\right\}.

This contrast holds at the same time and with the same values of other covariates and random effects. The direct effect \(\gamma x\) cancels because both trajectories are in the same subgroup. In contrast, the hazard ratio comparing \(x=1\) with \(x=0\) at the same biomarker value \(m\) is \(\exp(\gamma+\delta m)\), holding other terms fixed. Thus, the group hazard contrast can depend on the biomarker level when \(\delta\ne0\).

The usual no-interaction null hypothesis is \(H_0:\delta=0\). Under that null, the two groups have the same log-hazard association per biomarker unit, though they may still differ through \(\gamma\). A nonzero \(\delta\) indicates different association strengths on the log-hazard scale; it does not by itself establish a causal effect or show that the association is clinically important.

Deriving the Subgroup Contrasts

The subgroup expressions follow by substituting each covariate value into the hazard model. Taking the log hazard first makes the role of the interaction explicit.

Set the subgroup value.
For \(x=0\), the biomarker term is \((\alpha+\delta\cdot0)m_i(t)=\alpha m_i(t)\). For \(x=1\), it is \((\alpha+\delta)m_i(t)\).

Compare two biomarker values within that subgroup.
Subtract the log hazards for values differing by \(\Delta m\). The direct covariate and other shared terms cancel, leaving \(\alpha\Delta m\) for \(x=0\) and \((\alpha+\delta)\Delta m\) for \(x=1\).

Exponentiate the log-hazard difference.
The resulting hazard ratios are \(\exp(\alpha\Delta m)\) and \(\exp\{(\alpha+\delta)\Delta m\}\), respectively. Their ratio is \(\exp(\delta\Delta m)\).

This derivation also shows why testing \(\delta=0\) addresses whether the link differs between groups. Testing \(\alpha=0\) addresses the association in the reference group only. Testing \(\alpha+\delta=0\) addresses the association in the other group, and its uncertainty must account for the covariance between the coefficient estimates.

Worked Example: Calculating Hazard Ratios by Subgroup

Suppose a hypothetical joint model estimates \(\alpha=0.20\) and \(\delta=0.15\), with the biomarker measured in mmol/L. For a two-unit higher latent biomarker value, calculate the conditional hazard ratio in each subgroup. Assume the proportional-hazards model and the specified linear current-value link are adequate for this calculation.

For \(x=0\), the association per mmol/L is \(\alpha=0.20\). A two-unit contrast gives

\mathrm{HR}_0(2) = \exp(0.20\times2) = \exp(0.40) = 1.4918.

Thus, at a given time, a two mmol/L higher latent biomarker value corresponds to a hazard 1.4918 times as large in subgroup 0, conditional on the model terms held fixed. For \(x=1\), the association per mmol/L is \(\alpha+\delta=0.35\), so

\mathrm{HR}_1(2) = \exp\{(0.20+0.15)\times2\} = \exp(0.70) = 2.0138.

The ratio of these two subgroup-specific hazard ratios is \(\exp(\delta\times2)=\exp(0.30)=1.3499\). This indicates a larger modeled hazard contrast per two-unit biomarker difference in subgroup 1. It does not mean that the subgroup 1 hazard itself is 1.3499 times the subgroup 0 hazard: that between-group comparison also depends on \(\gamma\) and the biomarker value.

Assessing Whether the Link Differs

A common large-sample assessment uses a Wald statistic for the interaction coefficient. If \(\widehat{\delta}\) is the estimate and \(SE(\widehat{\delta})\) its standard error, then \(z=\widehat{\delta}/SE(\widehat{\delta})\) is compared with a standard normal distribution under \(H_0:\delta=0\). The standard error must come from the fitted joint model and reflect uncertainty in the estimated parameters. This assessment relies on an adequate model, regular estimation behavior, and a sufficiently accurate large-sample approximation.

A confidence interval for \(\delta\) is useful alongside the test because it shows the range of interaction strengths compatible with the data. If the goal is to describe each subgroup's association, calculate its uncertainty using the coefficient covariance matrix; do not treat \(\widehat{\alpha}\) and \(\widehat{\delta}\) as independent unless their estimated covariance is zero.

Worked Example: Testing the Interaction and Quantifying Subgroup Uncertainty

Consider hypothetical estimates \(\widehat{\alpha}=0.20\), \(\widehat{\delta}=0.12\), \(SE(\widehat{\alpha})=0.08\), and \(SE(\widehat{\delta})=0.10\), with \(\operatorname{Cov}(\widehat{\alpha},\widehat{\delta})=-0.0024\). Test whether the association differs between the two groups at the 5% level. Assume the Wald approximation is reasonable.

For the interaction, \(z=0.12/0.10=1.20\). The two-sided normal-approximation p-value is approximately 0.230. Since \(0.230>0.05\), these estimates do not provide sufficient evidence to reject \(H_0:\delta=0\) at the 5% level. This is not proof that the subgroup associations are identical.

The subgroup 0 coefficient is \(\widehat{\alpha}=0.20\), with a 95% Wald interval \(0.20\pm1.96(0.08)=(0.0432,0.3568)\). Exponentiating gives an interval of approximately \((1.044,1.429)\) for the hazard ratio per one-unit higher biomarker value. For subgroup 1, the coefficient estimate is \(0.20+0.12=0.32\). Its variance is

\operatorname{Var}(\widehat{\alpha}+\widehat{\delta}) = 0.08^2+0.10^2+2(-0.0024) = 0.0116, \qquad SE(\widehat{\alpha}+\widehat{\delta}) = \sqrt{0.0116} = 0.1077.

The subgroup 1 interval on the log-hazard scale is \(0.32\pm1.96(0.1077)=(0.1089,0.5311)\). Exponentiating gives approximately \((1.115,1.701)\) for the hazard ratio per one-unit higher biomarker value. These are model-based 95% confidence intervals: under repeated sampling and the stated approximation, the procedure has 95% coverage for its target parameter. The nonsignificant interaction test and the positive subgroup-specific estimates answer different questions; neither should be substituted for the other.

Worked Example: Separating the Direct Group Effect From Effect Modification

Suppose a hypothetical model has \(\gamma=0.30\) and \(\delta=0.08\), with biomarker value \(m\) measured in mmol/L. Calculate the hazard ratio for subgroup 1 versus subgroup 0 at \(m=2\) and \(m=4\), holding the remaining terms fixed. Assume the fitted model's log-linear interaction is appropriate over this range.

At a fixed biomarker value \(m\), the group hazard ratio is \(\exp(\gamma+\delta m)\). At \(m=2\),

\exp(0.30+0.08\times2) = \exp(0.46) = 1.5841.

At \(m=4\),

\exp(0.30+0.08\times4) = \exp(0.62) = 1.8589.

The group hazard ratio is larger at the higher biomarker value because the model includes a positive interaction. The ratio of the two group hazard ratios is \(\exp\{0.08(4-2)\}=\exp(0.16)=1.1735\). The direct covariate contribution \(\gamma\) is present in both comparisons, whereas the change across biomarker levels is attributable to \(\delta\). If a more interpretable reference biomarker value is desired, the biomarker can be centered at that value; centering changes the interpretation of the direct group coefficient but not the interaction coefficient or the fitted contrasts.

Assumptions, Diagnostics, and Alternatives

The interaction has a useful interpretation only if the underlying model is credible. The hazard model above assumes that the conditional hazard is proportional over time, that the biomarker association is linear on the log-hazard scale within each subgroup, and that the coefficient difference represented by \(\delta\) is constant over time. It also assumes that the latent trajectory is adequately specified and that the shared-random-effects structure captures the relevant dependence between repeated biomarker measurements and the event process. Inference additionally depends on appropriate handling of censoring, commonly requiring censoring to be conditionally non-informative given the modeled information.

Check whether both groups have enough biomarker and event information to estimate their associations, and whether their biomarker ranges overlap. Sparse events in one group can make the interaction estimate unstable. Large standard errors, sensitivity to a few subjects, convergence warnings, or highly correlated coefficient estimates are reasons to qualify the result. Examine whether the proportional-hazards and trajectory assumptions are plausible, and consider whether the chosen linear biomarker scale is reasonable. A statistically clear interaction cannot repair a misspecified trajectory or an inappropriate functional form.

A model without an interaction imposes one common biomarker association and is more parsimonious, but it cannot represent genuine subgroup differences. Fitting separate models by subgroup allows more freedom but does not directly estimate a single coefficient difference and can make comparisons difficult, especially when estimates are imprecise. A stratified baseline hazard can allow groups to have different baseline hazard shapes; that is distinct from allowing their biomarker association coefficients to differ. When a linear interaction is not credible, a more flexible interaction may be considered, but it requires enough data and a clearly specified target contrast.

For a continuous baseline covariate \(z_i\), choose a reference value \(z_{\mathrm{ref}}\) and use \(\{\alpha+\delta(z_i-z_{\mathrm{ref}})\}m_i(t)\). The coefficient \(\delta\) then describes the change in log-hazard association per one-unit increase in \(z_i\). Centering clarifies the reference association, while the same modeling cautions apply: the assumed change in association must be reasonable over the covariate range.

Common Mistakes

Confusing interaction and direct effect. A nonzero \(\gamma\) means the groups differ in log hazard at the reference biomarker value; it does not show that the biomarker association differs. That question concerns \(\delta\).
Reporting only the reference-group association. In subgroup 1 the association is \(\alpha+\delta\), not \(\delta\) alone.
Ignoring the coefficient covariance. The variance of \(\widehat{\alpha}+\widehat{\delta}\) includes \(2\operatorname{Cov}(\widehat{\alpha},\widehat{\delta})\).
Calling a nonsignificant interaction proof of equality. A large p-value indicates insufficient evidence against the tested null at the chosen level, not evidence that the groups are identical.
Interpreting a subgroup hazard ratio as a causal effect. The model describes conditional associations unless additional causal assumptions and a suitable design justify a causal interpretation.

Key takeaway: With a binary baseline covariate, the subgroup-specific biomarker associations are \(\alpha\) and \(\alpha+\delta\). Calculate contrasts from those coefficients, assess their uncertainty with the fitted covariance matrix, and use the interaction test to evaluate whether the modeled association differs across groups.

Check Your Understanding

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

If \(\alpha=0.10\) and \(\delta=0.05\), what are the biomarker association coefficients in the two groups? What is the hazard ratio for a three-unit increase in the biomarker in group 1?
In the model, what does \(\gamma\) describe, and why does it not test whether the biomarker link differs across groups?
Suppose \(\widehat{\alpha}=0.30\), \(\widehat{\delta}=-0.20\), and \(SE(\widehat{\delta})=0.08\). Calculate the Wald statistic for testing \(H_0:\delta=0\), and state the two-sided conclusion at the 5% level using a normal approximation.
Why is the variance of \(\widehat{\alpha}+\widehat{\delta}\) not generally the sum of the two coefficient variances?

Solutions. (1) The coefficients are 0.10 for group 0 and \(0.10+0.05=0.15\) for group 1. The three-unit hazard ratio in group 1 is \(\exp(0.15\times3)=\exp(0.45)=1.5683\). (2) \(\gamma\) is the direct group contribution to log hazard at the biomarker reference value (zero as written); it does not compare the biomarker coefficients. The interaction coefficient \(\delta\) is their difference. (3) \(z=-0.20/0.08=-2.50\), giving a two-sided normal-approximation p-value of about 0.0124. Reject \(H_0\) at the 5% level; the data provide evidence of a difference in modeled association strengths. (4) The estimates may be correlated, so the variance is \(\operatorname{Var}(\widehat{\alpha})+\operatorname{Var}(\widehat{\delta})+2\operatorname{Cov}(\widehat{\alpha},\widehat{\delta})\). Omitting the covariance is justified only when it is zero.