Joint Models For Multiple Omics Trajectories

Why Model Several Trajectories Together?

A patient’s molecular measurements may evolve together: a shared biological process can affect several omics features, while each feature may also have its own trajectory. If event risk is related to these evolving measurements, fitting separate longitudinal models and treating their fitted values as known predictors can discard uncertainty and cross-feature dependence. A joint model represents the trajectories and event time in one probability model.

The central distinction is between how the omics features depend on one another and how each feature is associated with event risk. Shared latent random effects describe dependence among trajectories. Modality-specific association parameters describe how the fitted trajectory for each feature enters the hazard. Correlated trajectories do not, by themselves, imply equal risk associations.

A Joint Model for Multiple Omics Trajectories

Let \(i\) index subjects, \(k=1,\ldots,K\) index omics features, and \(j\) index measurements. Let \(Y_{ikj}\) be feature \(k\)'s measured value for subject \(i\) at time \(t_{ikj}\). A useful model separates a latent effect shared across features from feature-specific random effects:

Here \(x_{ikj}\) and \(z_{ikj}\) are specified covariate and time-basis vectors; \(\beta_k\) is a feature-specific fixed-effect vector; \(U_i\) is a vector of shared latent effects; \(\lambda_k\) gives feature \(k\)'s loading on those effects; and \(V_{ik}\) is a feature-specific random effect. The subject-specific fitted trajectory, conditional on the latent effects, is

A basic specification assumes \(U_i\sim N(0,G)\) and the stacked vector \(V_i=(V_{i1}^{T},\ldots,V_{iK}^{T})^{T}\sim N(0,D)\), independently of \(U_i\). It may allow dependence between the components of \(V_i\) through \(D\). For a simple Gaussian measurement model, \(\epsilon_{ikj}\sim N(0,\sigma_k^2)\), independently across measurements conditional on the random effects. More general residual covariance structures are possible when justified by the measurement design.

For a single event time \(T_i\), let \(\Delta_i=1\) indicate an observed event and \(\Delta_i=0\) indicate right censoring. A current-value association model is

Here \(h_0(t)\) is the baseline hazard, \(w_i\) contains event-model covariates, \(\gamma\) is their coefficient vector, and \(\alpha_k\) is the association parameter for feature \(k\). The model assumes that, conditional on the latent effects and included covariates, the hazard depends on the specified current fitted trajectory values. It does not automatically establish a causal effect of an omics feature.

The shared effects influence the hazard through the trajectories they help determine. For example, if \(U_i\) is scalar, its contribution to the log hazard through these trajectories is \(\sum_k\alpha_k\lambda_k U_i\). Thus a shared latent effect can be associated with event risk even though the hazard has a separate \(\alpha_k\) for each feature. Adding an additional direct coefficient for \(U_i\) is a different, more elaborate model; it can be difficult to distinguish from the trajectory-mediated contribution and needs a clear scientific rationale and adequate information.

Interpreting the Association Parameters

Conditional on the latent effects and other covariates, increasing feature \(k\)'s fitted value by one unit while holding the other fitted values fixed multiplies the hazard by \(\exp(\alpha_k)\). The unit is the unit in which that feature enters the model. If a feature is standardized, the unit is one standard deviation; if it is log-transformed, it is one unit on the log scale. The interpretation is conditional association, not a comparison of two otherwise interchangeable patients.

When features are correlated, “holding the other fitted values fixed” may describe a contrast that is unusual or poorly supported by the data. The \(\alpha_k\) values are still parameters of the specified model, but their separate estimates can be imprecise. It is useful to report uncertainty intervals and to assess whether the data contain enough independent trajectory variation to distinguish the associations.

Shared Dependence Is Not the Same as Equal Risk Association

The shared random effect can induce correlation between measurements from different features. That dependence is governed by the latent-effect distribution and the feature loadings, not by the \(\alpha_k\) parameters. Conversely, the \(\alpha_k\) values can differ even when the trajectories are strongly correlated. This distinction is central to interpreting a multi-trajectory joint model.

Likelihood and a Justification for Joint Estimation

For a subject with event or censoring time \(T_i\), the event contribution conditional on the random effects is the hazard at \(T_i\), raised to \(\Delta_i\), multiplied by the survival probability through \(T_i\). If \(f_i(Y_i\mid U_i,V_i)\) denotes the joint density of that subject’s longitudinal measurements, the subject-level marginal likelihood contribution is

The integral averages over the unobserved random effects. Its construction relies on a specified longitudinal measurement model, a specified conditional hazard, and a specified random-effects distribution. It also assumes that censoring is noninformative conditional on the modeled quantities, so the censoring mechanism does not require an additional likelihood component. The earlier material on likelihood construction for joint longitudinal-survival models gives the general event-time factor; here the longitudinal density and latent effects must cover all included features.

The reason for integrating rather than plugging in estimated trajectories can be seen directly from the conditional model. Given the random effects, the measurement density and event contribution are modeled together. But those effects are not observed, so the likelihood for the observed data is their probability-weighted average, not the likelihood evaluated at a single estimated value. In general, integrating a likelihood over random effects is not equivalent to substituting their estimated means. Joint estimation therefore carries uncertainty about the latent trajectories into estimation of the event associations.

Assumptions, Diagnostics, and Failure Modes

The interpretation and estimation depend on assumptions that should be checked against the study design and data:

• Trajectory form: the chosen fixed and random time functions adequately describe the mean trajectories over observed times. If curvature remains in residuals, a linear time term may be inadequate; spline models are one possible choice.
• Random-effects structure: the distribution and covariance restrictions for shared and feature-specific effects are plausible. Too few shared factors can miss dependence; too many can be weakly identified.
• Measurement model: the conditional error distributions and residual dependence are adequate for each assay. Heavy tails, outliers, batch effects, or assay-specific heteroscedasticity can distort trajectory estimates.
• Event model: the baseline hazard and chosen trajectory association form are adequate. A current-value link is not interchangeable with a slope or cumulative-exposure link.
• Conditional censoring: censoring is noninformative given the modeled information. If dropout or visit timing depends on unmodeled health status, that assumption may fail.
• Identifiability: enough subjects, repeated measurements, and independent variation exist to distinguish shared effects, feature-specific effects, and \(\alpha_k\) values.

Useful diagnostics include plotting observed measurements and fitted trajectories by feature, inspecting residual patterns against time and fitted values, checking convergence and sensitivity to starting values, and examining estimated random-effect variances and correlations. Very large standard errors, correlations near their parameter boundaries, unstable estimates across reasonable model specifications, or a nearly singular covariance matrix can indicate weak identification. For event associations, compare estimates and uncertainty intervals across defensible trajectory and hazard specifications rather than relying only on whether a coefficient crosses a significance threshold.

A frequent mistake is to interpret a strong correlation between two omics trajectories as evidence that their hazard associations must be the same. Correlation is represented in the random-effects and residual structure; equality of associations is a separate restriction such as \(\alpha_1=\alpha_2\), which must be justified and assessed rather than assumed. Another mistake is to treat fitted trajectories as error-free observed covariates in a survival model. They are estimated from noisy, incomplete measurements, and ignoring that estimation uncertainty can produce misleading association estimates.

Alternatives and Limitations

The joint model is not automatically superior: it trades a more explicit account of measurement and latent dependence for stronger distributional and structural assumptions. With sparse measurements or highly correlated features, separate \(\alpha_k\) estimates may be unstable. A common-factor specification also requires enough information to identify its scale and loadings; without suitable constraints, different parameter values can describe nearly the same trajectory distribution. The model’s conclusions should therefore be framed as conditional on its stated assumptions and measurement scale.