Worked Joint Model for an Ordinal Biomarker

Why Share an Ordinal Trajectory With the Event Model?

An ordinal biomarker records ordered categories rather than measurements on a continuous scale. A participant’s sequence of categories may nevertheless reveal a changing underlying trajectory: movement toward higher categories can be associated with a different event rate. A joint model represents that trajectory while accounting for subject-to-subject variation, then connects that variation to the time-to-event outcome.

The model below uses cumulative-logit probabilities for the ordinal responses and a survival hazard that depends on the same subject-specific random effects. This shared-effects construction allows the marker and event outcomes to be associated after conditioning on measured covariates. It does not, by itself, establish that changes in the marker cause changes in event risk.

Ordinal Trajectory and Survival Submodels

Use an increasing cumulative-logit model. Let \(F(u)=1/(1+\exp(-u))\) be the logistic distribution function, and let \(\alpha_1<\cdots<\alpha_{K-1}\) be cutpoints. Define the subject-specific linear predictor \(\eta_i(t)\) using fixed covariates and a flexible time trajectory, together with random effects:

Here \(x_i\) contains baseline covariates, \(B(t)\) is a chosen spline basis, and \(z_i(t)\) specifies which trajectory components vary by subject. For example, a random intercept and slope use \(z_i(t)=(1,t)^\mathsf{T}\). The vector \(b_i\) has a multivariate normal distribution with mean zero and covariance matrix \(D\). The cutpoint ordering ensures valid, non-overlapping ordinal categories. Because of the minus sign, larger \(\eta_i(t)\) shifts probability toward higher categories.

The category probabilities follow by differencing adjacent cumulative probabilities. Set \(\alpha_0=-\infty\) and \(\alpha_K=+\infty\); then

For the event model, let \(w_i\) contain baseline survival covariates and let \(\lambda\) link the marker’s random effects to the hazard:

The baseline hazard \(h_0(t)\) may be specified parametrically or left flexible. Conditional on \(b_i\), the random effects act multiplicatively on the hazard. A one-unit increase in random-effect component \(b_{ir}\), holding other components and covariates fixed, multiplies the hazard by \(\exp(\lambda_r)\). Its unit depends on the component: a random slope, for instance, is measured in log-odds per unit time.

This specification assumes that, conditional on \(b_i\) and included covariates, the repeated ordinal responses and event process are independent. It also requires an appropriate model for the distribution of \(b_i\), a suitable baseline hazard, and censoring that is non-informative conditional on the modeled information. The proportional-odds assumption is that the same \(\eta_i(t)\) shifts every cumulative logit; the effect is not category-specific.

Worked Category Probabilities

How the Joint Likelihood Is Formed

Conditional on \(b_i\), multiply the ordinal probabilities across visits and the event-time contribution. Write \(q_i(t)=\gamma^\mathsf{T}w_i+\lambda^\mathsf{T}b_i\). The event-time contribution for a right-censored observation is expressed using the hazard and cumulative hazard; it covers both observed events and censoring:

In this shared-effects model \(q_i(t)\) does not vary with time, but writing the hazard contribution in this form makes its event and censoring parts explicit. When \(\delta_i=1\), the hazard at the event time contributes; when \(\delta_i=0\), that factor is absent. Both cases include the probability of remaining event-free through \(T_i\), represented by the exponential of the negative cumulative hazard.

The subject-level likelihood is not obtained by simply substituting \(b_i=0\). The random effects are unobserved, so the likelihood integrates the entire conditional contribution against their density:

This integration is the key justification for joint estimation: each possible random trajectory contributes according to its plausibility under \(N(0,D)\), and the observed marker and event data together inform that plausibility. In practice the integral usually requires numerical methods, particularly when there are several random effects.

Worked Hazard Association

Worked Conditional Likelihood Contribution

Checking the Model and Comparing Alternatives

The model’s interpretation depends on its assumptions, so diagnostics should target both submodels and their connection. Useful checks include:

• Ordinal fit: compare observed and fitted category frequencies across time or predictor ranges. Patterns of systematic misfit can indicate that a single cumulative-logit shift is inadequate.
• Proportional odds: assess whether the predictor’s effect is plausibly common across cumulative cutpoints. Strongly different effects across cutpoints can undermine the shared-slope assumption.
• Trajectory shape: inspect fitted category probabilities over time and look for systematic temporal patterns in residual or predictive checks. A spline that is too restrictive can miss curvature; one that is too flexible can produce unstable subject trajectories.
• Survival fit: assess the baseline-hazard specification and whether the random-effects association adequately captures the intended dependence. A poorly fitting hazard or an overly restrictive linear association can distort \(\lambda\).
• Computation: check optimizer convergence, sensitivity to starting values, and stability when numerical integration accuracy is increased. A convergence message alone does not establish that the likelihood has been accurately evaluated.

A main alternative is to fit a longitudinal ordinal model and a survival model separately. Separate fits can be simpler, but they do not use the event information to help estimate the shared latent trajectory and do not represent the dependence through common random effects. Another joint-model choice links the hazard to the current latent trajectory value \(\eta_i(t)\), rather than to selected random-effect components. That alternative answers a different modeling question; the shared-random-effects link used here directly associates event risk with subject-specific departures from the population trajectory.

The model also has limits. Its conclusions depend on the random-effects distribution, chosen spline basis, proportional-odds structure, and survival link. Sparse categories or few repeated observations can make cutpoints and random-effect variation difficult to estimate. A shared random effect does not automatically resolve unmeasured confounding or guarantee valid causal conclusions. The fitted model is a structured description of association under stated assumptions.

Common Pitfalls

• Reversing the cumulative-logit sign. With \(\operatorname{logit}P(Y\leq k)=\alpha_k-\eta\), a larger \(\eta\) moves probability toward higher categories. Reversing the sign without changing interpretation reverses the direction of the model.
• Using cumulative probabilities as category probabilities. The probability of category \(k\) is a difference of cumulative probabilities, not \(P(Y\leq k)\) itself.
• Interpreting \(\lambda\) without its random effect. A hazard ratio for a one-unit random-effect difference is \(\exp(\lambda_r)\); \(\lambda_r\) alone is a log-hazard contrast, and its units depend on \(b_{ir}\).
• Ignoring the latent-effect integral. A conditional likelihood at one chosen \(b_i\) is not the marginal likelihood. The latter integrates the joint marker-event contribution over the random-effects distribution.