Why Model Omics Layers Together?
Repeated measurements from different omics modalities can move together within a person over time. A subject whose expression measurements tend to remain above the population mean might also have protein measurements that tend to remain above their mean. A multivariate mixed model represents this pattern through correlated subject-specific random effects. Its purpose is to describe and estimate shared variation across layers while accounting for repeated observations; a fitted covariance is not, by itself, evidence that one layer causes changes in another.
The earlier material on mixed models for repeated omics measurements establishes the basic idea of subject-specific random effects. Here the model is extended so that random effects from different modalities can covary. The focus is the covariance structure and its interpretation, not time-to-event modeling or flexible trajectory functions.
Model and Assumptions
Suppose subject \(i\) has measurements on \(K\) omics modalities. Let \(Y_{ij}(t)\) be the response for subject \(i\), modality \(j\), at time \(t\). Responses may be transformed or normalized assay values; their units and scale must be stated. A common model is
Here \(\mathbf{x}_{ij}(t)\) contains fixed-effect predictors, \(\boldsymbol{\beta}_j\) contains modality-specific population coefficients, and \(\mathbf{z}_{ij}(t)\) contains predictors for random effects. The vector \(\mathbf{b}_{ij}\) contains the subject-specific random effects for modality \(j\), such as a random intercept and possibly a random slope. Stack these vectors across modalities as \(\mathbf{b}_i=(\mathbf{b}_{i1}^\mathsf{T},\ldots,\mathbf{b}_{iK}^\mathsf{T})^\mathsf{T}\).
A standard specification assumes \(\mathbf{b}_i\sim N(\mathbf{0},G)\), independent across subjects. Conditional on these random effects and the modeled covariates, the errors have mean zero and covariance \(R_i\), and are independent of the random effects. The errors may be correlated across modalities or measurement times if that dependence is explicitly represented in \(R_i\). Conditional errors from different subjects are usually taken to be independent. The fixed-effect mean structure must be appropriate, and the visit schedule must not be informative in a way the model ignores. For the usual likelihood analysis, the random effects and errors are commonly assumed Gaussian; departures can affect estimates and especially small-sample inference.
The model need not use a separate random effect for every assay feature. It is often fit to a small set of prespecified or summarized responses. A covariance matrix with too many random effects relative to the number of subjects or repeated observations may be weakly identified.
What the Cross-Layer Covariance Means
Stack all observations for subject \(i\) into a vector \(\mathbf{Y}_i\). Let \(X_i\) and \(Z_i\) be the corresponding fixed- and random-effect design matrices, and let \(\boldsymbol{\beta}\) stack the fixed effects. The model can then be written as
The marginal covariance follows from the law of total covariance. Conditional on \(\mathbf{b}_i\), the covariance of \(\mathbf{Y}_i\) is \(R_i\), while its conditional mean is \(X_i\boldsymbol{\beta}+Z_i\mathbf{b}_i\). Therefore
This separates two sources of dependence. \(Z_iGZ_i^\mathsf{T}\) describes dependence induced by persistent subject-specific deviations in trajectories. \(R_i\) describes residual dependence after conditioning on those effects. A positive cross-layer random-intercept covariance means subjects with higher-than-population intercepts in one modality tend, on average, to have higher-than-population intercepts in the other. It does not necessarily imply that their measurements are strongly correlated at every pair of times: the design matrices and any residual covariance also matter.
Worked Example: Interpreting Cross-Layer Random Intercepts
Suppose two assay responses have random-intercept variances \(0.52\) and \(0.34\), in squared standardized assay units, and their random-intercept covariance is \(0.18\). The correlation of the subject-specific intercepts is covariance divided by the product of their standard deviations:
The denominator is \(\sqrt{0.1768}\approx0.420476\), which rounds to \(0.4205\). The estimated correlation is about \(0.428\), a moderate positive association between the subjects' latent assay-specific intercept deviations. It is not the correlation between all observed assay values: those also include fixed effects, time patterns, and residual variation. For these two variances and covariance, the implied \(2\)-by-\(2\) covariance block is admissible because \(0.18^2=0.0324\) is less than \(0.52(0.34)=0.1768\).
With random slopes, cross-layer covariance can depend on both measurement times. Consider two modalities with random intercepts and slopes, denoted \(b_{10},b_{11}\) and \(b_{20},b_{21}\). The subject-specific deviations at times \(t\) and \(s\) are \(b_{10}+t b_{11}\) and \(b_{20}+s b_{21}\). Expanding their covariance gives
Thus the cross-layer random-intercept covariance is only one part of shared trajectory variation. Cross-covariances involving slopes determine whether a subject's deviation in the rate of change in one modality tends to align with an intercept or slope deviation in another.
Worked Example: Shared Trajectory Covariance at Two Times
Suppose two standardized assay trajectories are modeled with random intercepts and slopes, with time measured in months. Let the cross-layer covariance terms be \(G_{10,20}=0.18\), \(G_{10,21}=0\), \(G_{11,20}=0\), and \(G_{11,21}=0.04\). The covariance between the first modality's subject deviation at month \(t=1\) and the second modality's deviation at month \(s=2\) is
The covariance is \(0.26\) squared standardized assay units. The \(0.18\) term is shared intercept variation; the \(0.08\) contribution comes from correlated slopes over the two-month product \(ts\). This calculation is valid only as part of a full admissible \(G\): all its variances and covariances together must form a positive semidefinite matrix, not merely yield a plausible selected covariance.
Fitting and Comparing Covariance Structures
Parameters are commonly estimated by maximum likelihood or restricted maximum likelihood (REML). Estimation chooses fixed effects and covariance parameters to fit the observed repeated measurements under the specified model, subject to covariance constraints such as positive semidefiniteness. REML accounts for estimation of fixed effects when estimating covariance parameters; comparisons of models with different fixed-effect structures generally use maximum likelihood rather than REML likelihoods. Standard errors and intervals for covariance parameters can be less reliable than those for well-supported fixed effects, particularly with few subjects.
The most flexible choice is an unstructured \(G\), which estimates every allowed variance and covariance. This can represent varied cross-layer patterns, but the number of parameters increases quickly. More restrictive structures can be appropriate when scientifically justified, for example allowing correlated random intercepts across modalities while setting cross-layer slope covariances to zero. Constraints should reflect a defensible model, not be introduced only because they make results look cleaner.
| Approach | What it represents | Main trade-off |
|---|---|---|
| Separate univariate mixed models | Repeated dependence within each modality | Do not estimate cross-layer random-effect covariance; simpler to fit |
| Multivariate mixed model | Joint fixed effects and cross-layer random-effect covariance | Can describe shared trajectories, but requires more covariance parameters |
| Restricted cross-layer structure | Selected shared effects, with other covariances fixed to zero or constrained | More stable when justified, but can miss real dependence if too restrictive |
Worked Example: Cross-Layer Dependence and a Treatment Contrast
In a hypothetical fitted model, the estimated treatment effects are \(0.30\) standardized units for one modality and \(0.10\) standardized units for another. Suppose their estimated fixed-effect variances are \(0.04\) and \(0.05\), and their covariance is \(0.015\). The estimated difference is \(0.30-0.10=0.20\) standardized units. Its standard error uses the covariance:
Using a normal approximation, a \(95\%\) interval is \(0.20\pm1.96(0.245)\), or approximately \((-0.280,0.680)\) standardized units. This interval includes zero, so these hypothetical results do not establish a difference between the treatment effects at the \(5\%\) level under the approximation. The fixed-effect covariance in this calculation is part of the fitted model's coefficient uncertainty; it should not be confused with a single entry of \(G\), the random-effect covariance matrix.
Diagnostics, Failure Modes, and Limits
A model's cross-layer covariance is informative only if its assumptions and data support it. Inspect residuals by modality and time for mean patterns, unequal spread, and remaining temporal or cross-layer dependence. Compare observed within-subject trajectories with fitted conditional trajectories, while remembering that estimated random effects are shrunken predictions rather than directly observed subject traits. Check convergence warnings, boundary estimates, near-zero variance components, and very high estimated correlations. A singular fit or a covariance matrix with eigenvalues near zero can indicate that the data do not support the requested random-effect structure.
Also verify that modalities are on interpretable and reasonably modeled scales. A few extreme values can dominate Gaussian fits; a nonlinear mean trend or changing variability can be mistaken for covariance structure. If random effects and residuals are both allowed to explain similar cross-layer dependence, they may be difficult to distinguish without adequate repeated observations or design variation. Sensitivity to a simpler \(G\), a different defensible \(R_i\), or removal of a weakly supported random slope is useful evidence about stability, not a license to select whichever result is preferred.
- Too many covariance parameters. A large unstructured matrix may fail to converge or yield unstable estimates when the number of subjects is limited.
- Confusing covariance sources. Cross-layer residual correlation at a visit is not the same as persistent covariance between subject-level trajectories.
- Overinterpreting correlation. A covariance summarizes modeled association after specified adjustments; it does not establish biological direction or causation.
- Ignoring observation patterns. If measurement or dropout depends on unobserved outcomes beyond the model's assumptions, ordinary mixed-model likelihood results may be biased.
A multivariate mixed model also has limits. Its conclusions depend on the selected mean and covariance structures, and it may be computationally demanding for many layers, features, or random effects. It describes repeated outcomes jointly but does not automatically model a survival process or correct informative dropout. If a simpler analysis is adequate, separate univariate models can be easier to interpret, though they discard cross-layer covariance information and cannot directly estimate shared subject-level effects.
Check Your Understanding
Answer each question, then compare with the brief solutions.
- What does a positive covariance between two modalities' random intercepts describe?
- In a multivariate mixed model, what two terms contribute to the marginal covariance of the observations?
- Why does a valid cross-layer covariance require checking the full \(G\), rather than only checking whether one covariance is positive?
- A model estimates cross-layer residual covariance but sets cross-layer random-effect covariance to zero. Does it describe shared subject-level trajectories across modalities?
Solutions. (1) Subjects with higher-than-population intercept deviations in one modality tend, on average, to have higher-than-population intercept deviations in the other, conditional on the model. (2) The terms are \(Z_iGZ_i^\mathsf{T}\), induced by random effects, and \(R_i\), the conditional residual covariance. (3) Every variance and covariance must combine into a positive semidefinite matrix; individually plausible entries can still produce an inadmissible matrix. (4) No. It allows residual dependence, but the model does not attribute cross-layer dependence to correlated subject-level random effects.