Covariate-Adjusted Multi-Omics Factor Models

Why Put Covariates in the Observation Mean?

In a multi-omics study, measured characteristics such as age, treatment group, collection site, or processing batch can affect many features at once. A factor model may also describe coordinated variation across features and omics views. If a measured effect is left out of the model, that broad pattern can be attributed to a latent factor even when the recorded covariate explains it. The purpose of covariate adjustment here is to represent measured mean differences explicitly, so that estimated factors describe remaining structure.

This distinction is about the model, not a claim that covariate adjustment identifies biological causes. A factor estimated conditional on recorded covariates summarizes residual association under the chosen model. Its interpretation depends on which covariates were measured and included, as well as on the factor model’s assumptions.

Model and Interpretation

Suppose the same $n$ subjects have measurements from $M$ omics views. For view $m$, let $Y^{(m)}$ be the $n \times p_m$ matrix of measurements, with subjects in rows and features in columns. Let $C$ be the $n \times q$ design matrix of measured covariates, including an intercept if the model has feature-specific intercepts. A covariate may have several design columns, for example, indicators for the levels of a categorical variable.

Definition: A covariate-adjusted multi-omics factor model for view $m$ is $$ Y^{(m)} = C B^{(m)} + Z \bigl(L^{(m)}\bigr)^T + E^{(m)}. $$ Here $B^{(m)}$ is a $q \times p_m$ matrix of feature-specific covariate coefficients; $Z$ is an $n \times K$ matrix of subject factor scores shared across views; $L^{(m)}$ is a $p_m \times K$ matrix of view-specific feature loadings; and $E^{(m)}$ is residual variation. The term $C B^{(m)}$ is the modeled covariate contribution to the observation mean. The factor term describes coordinated variation remaining beyond that mean.

The model is commonly fitted by minimizing the summed squared residuals, or through a likelihood with a specified error distribution. For the least-squares interpretation, assume the observations are appropriately scaled, the same subjects are matched across views, and residuals have mean zero conditional on the covariates and factors. A simple likelihood takes residual entries to be independent Gaussian variables with mean zero and common variance. Other error structures can be used, but then the fitting criterion and the justification for preprocessing must reflect those structures.

As in factor models generally, scores and loadings have scale and rotation ambiguities unless constraints or conventions are imposed. The covariate model does not by itself resolve those ambiguities. The central interpretation is that a loading describes a feature’s association with a factor direction, while a score describes a subject’s position along that direction, conditional on the modeled covariate mean. Neither coefficient is automatically causal.

Joint Estimation and Residualization

One can estimate the covariate coefficients and factor structure together. A simpler alternative is to remove the fitted covariate component from each view and estimate factors from the resulting residual matrices. These procedures coincide under particular conditions; residualization is not a universal substitute for joint modeling.

For a derivation, consider least squares with full-column-rank $C$. Define the projection onto the column space of $C$ as $P_C=C(C^TC)^{-1}C^T$, and define the residual projection $M_C=I-P_C$. For fixed $Z$ and $L^{(m)}$, minimizing the squared error over $B^{(m)}$ gives

\widehat{B}^{(m)} =(C^TC)^{-1}C^T\left(Y^{(m)}-Z(L^{(m)})^T\right).

If the factor scores are constrained to be orthogonal to the covariate design, $C^TZ=0$, then $C^TZ(L^{(m)})^T=0$, and therefore

\widehat{B}^{(m)}=(C^TC)^{-1}C^TY^{(m)}.

Substituting this estimate leaves $M_CY^{(m)}-Z(L^{(m)})^T$ as the residual matrix. Thus, under this constraint and a common unweighted least-squares criterion, fitting factors to covariate-residualized data gives the same factor-fitting objective as joint least squares with $C^TZ=0$. This argument also requires the same design and applicable observations for each feature; missingness or feature-specific weights can change the projection.

Worked Example: A Shared Residual Pattern After Covariate Adjustment

Consider four hypothetical subjects, two omics features, and a continuous covariate $x=(0,1,2,3)^T$. Let $r=(-1,1,1,-1)^T$. Suppose the measured values follow $Y^{(1)}=10+2x+r$ and $Y^{(2)}=4+0.5x-0.5r$, where the first feature is an RNA measurement and the second a protein measurement, in their respective arbitrary assay units.

Subject	$x$	$r$	RNA feature	Protein feature
1	0	-1	9	4.5
2	1	1	13	4.0
3	2	1	15	4.5
4	3	-1	15	6.0

Take $C$ to have columns $1$ and $x$. The sum of the entries of $r$ is $-1+1+1-1=0$, and its inner product with $x$ is $0(-1)+1(1)+2(1)+3(-1)=0$. Thus $r$ is orthogonal to both covariate columns. The fitted covariate values are $10+2x=(10,12,14,16)^T$ for RNA and $4+0.5x=(4,4.5,5,5.5)^T$ for protein. Subtracting these values from the observed data gives residuals $r=(-1,1,1,-1)^T$ and $-0.5r=(0.5,-0.5,-0.5,0.5)^T$, respectively.

The two residual columns form the matrix $[r,\,-0.5r]$, which has rank one: its second column is $-0.5$ times its first. A one-factor representation has scores $r$ and loadings $1$ for RNA and $-0.5$ for protein. The shared factor describes a coordinated residual pattern after adjustment for $x$; it does not describe the covariate trend itself.

What Does “Conditional on Covariates” Mean?

The factor scores are interpreted relative to the fitted mean $C B^{(m)}$. For example, a positive score indicates that a subject lies in the positive direction of the estimated residual pattern, after accounting for the included covariates. A loading indicates how a feature varies along that direction. The same score can have different feature manifestations across views because each view has its own loading matrix.

This is a model-based conditional interpretation, not a guarantee that every estimated score is independent of every covariate in a finite sample. If scores are explicitly constrained by $C^TZ=0$, they are orthogonal to the design in the fitted data. If scores instead have a population distribution independent of $C$, sample correlations can still occur. The modeling choice should be stated because it affects what “adjusted” means.

Worked Example: A Covariate Can Mimic a Shared Factor

Suppose four hypothetical subjects have a binary covariate $x=(0,0,1,1)^T$, and two measured features are RNA values $(2,2,6,6)^T$ and protein values $(5,5,8,8)^T$. With an intercept and $x$ in $C$, the fitted means are $2+4x$ and $5+3x$. Subtracting those means gives residual vectors $(0,0,0,0)^T$ in both views, so these data contain no residual factor variation under this exact model.

If the covariate is omitted, centering each feature gives $(-2,-2,2,2)^T$ and $(-1.5,-1.5,1.5,1.5)^T$. The centered columns are proportional: the second equals $0.75$ times the first. A one-factor analysis of these unadjusted values can therefore find a perfectly shared direction. Here that direction reflects the measured group difference. It would be misleading to interpret it as unexplained latent structure without accounting for $x$.

Overlap, Diagnostics, and Failure Modes

A covariate-adjusted model still needs enough information to distinguish its components. If a latent score direction lies in, or nearly in, the column space of $C$, then variation along that direction can be represented by the covariate term, the factor term, or both. This is the confounding problem between latent factors and covariates. An orthogonality constraint assigns covariate-aligned variation to the covariate component, but that assignment is a modeling convention and may remove scientifically meaningful variation if the covariate is not intended to define the adjustment target.

For example, if $x=(0,1,2,3)^T$ and an observed feature is exactly $2+3x$, the data alone do not distinguish a covariate slope of 3 and no factor from a covariate slope of 1 plus a factor contribution $2x$. Requiring scores to be orthogonal to the intercept and $x$ excludes the latter representation. The constraint makes the decomposition well-defined in that direction; it does not establish which interpretation is biologically correct.

Useful checks include verifying the covariate design and coding, inspecting whether covariate columns are redundant, and checking whether fitted scores or candidate factor directions strongly overlap the covariate space. Compare covariate coefficients and residual patterns under reasonable model specifications. Residual association with a known batch variable, for instance, can signal that its effect was omitted, misspecified, or inadequately modeled. Diagnostics cannot prove that all relevant covariates were measured.

Worked Example: Why an Overlapping Factor Cannot Be Separated by the Data Alone

For four hypothetical subjects, let $x=(0,1,2,3)^T$ and suppose one RNA feature has values $y=(2,5,8,11)^T$. The observed values satisfy $y=2+3x$ exactly. One representation uses an intercept of 2, a covariate slope of 3, and no factor contribution.

A different representation uses an intercept of 2, a covariate slope of 1, scores $z=x$, and loading 2. The fitted values are then $2+1x+2z=2+3x=(2,5,8,11)^T$. Both representations reproduce every observation, so the data do not distinguish them without an additional constraint or modeling assumption. Since $z=x$ lies in the covariate space, the constraint $C^TZ=0$, with $C$ containing the intercept and $x$, rules out this factor representation. The resulting factor interpretation is explicitly conditional on that constraint.

Choosing an Adjustment Strategy

Joint estimation models covariate effects and factors in one objective, which is useful when the likelihood, missing-data pattern, or error structure needs to be handled together. Residualizing first can be computationally convenient and is justified by the projection argument when its conditions hold. If views have different missingness, feature-specific weights, non-Gaussian outcomes, or covariance structures that require weighted fitting, ordinary least-squares residualization may not reproduce the joint fit.

A different model treats covariates as predictors of factor scores, for example, $Z=C A+U$. That asks whether covariates explain variation in the scores; it does not simply place covariate effects in each feature’s observation mean. It may be appropriate for a different scientific question, but it should not be confused with adjusting the feature means and then interpreting residual factors.

The main limitations are practical and interpretive. Misspecified covariate effects, omitted covariates, nonlinear relationships represented only by linear terms, and interactions that are left out can leave structured variation in the residuals. Conversely, adjusting for a variable that lies on a pathway of interest can remove variation one intended to study. Covariate adjustment changes the target of factor interpretation; it does not automatically improve every analysis.

Common pitfall: “The estimated factor is associated with treatment, so the factor must represent a treatment effect.” This reverses the model’s roles. If treatment is included in $C$, its modeled mean effect is in $C B^{(m)}$; the factor describes residual structure conditional on that term. Association between estimated scores and treatment may still arise from sampling variation, model misspecification, or a score model that does not impose exact orthogonality. It is not by itself evidence of a causal treatment effect.

Check Your Understanding

Answer each question, then compare with the brief solutions.

In $Y^{(m)}=CB^{(m)}+Z(L^{(m)})^T+E^{(m)}$, which term contains the feature-specific mean effects of measured covariates?
For $x=(0,1,2)^T$ and $r=(-1,2,-1)^T$, check whether $r$ is orthogonal to the intercept and to $x$.
Why can omitting a measured batch variable produce an apparently shared factor across two omics views?
State one condition under which ordinary residualization followed by factor fitting matches joint least squares, and one situation where the equivalence may fail.
Does a factor estimated conditional on treatment establish that the factor is a causal treatment effect? Explain.

Solutions. (1) $CB^{(m)}$ contains the covariate contributions to the modeled feature means. (2) The entries of $r$ sum to $-1+2-1=0$; its inner product with $x$ is $0(-1)+1(2)+2(-1)=0$. Thus $r$ is orthogonal to both columns of the intercept-and-$x$ design. (3) If batch shifts many features in each view, the shared shift remains in the data when batch is omitted and can appear as coordinated latent variation. (4) With full-rank $C$, a common unweighted least-squares criterion, and scores constrained by $C^TZ=0$, residualization yields the same factor-fitting objective. Feature-specific weights or missingness can invalidate that ordinary projection argument. (5) No. The factor summarizes modeled residual association conditional on treatment and other covariates; causal interpretation requires additional design and assumptions and does not follow from the factor model alone.