Residual Correlations After Covariate Adjustment

Why Adjust Expression Before Building a Network?

Gene-expression measurements can co-vary because the genes participate in related biology. They can also co-vary because samples differ in measured features such as processing batch, age, tissue composition, or treatment group. A correlation network built directly from expression values combines these sources. If the scientific aim is to describe co-expression beyond specified measured covariates, one approach is to regress those covariates out of each gene’s expression and correlate the resulting residuals.

This adjustment changes the meaning of an edge. It is no longer simply a summary of how two genes vary together across the observed samples. It describes their linear co-variation after subtracting the parts of their expression fitted by the selected covariates. The interpretation depends on which covariates were included, how they were modeled, and whether the regression assumptions are reasonable.

Definition and Linear Model

Suppose expression is measured for \(n\) samples and \(p\) genes. Let \(Y\) be the \(n \times p\) expression matrix, with one sample per row and one gene per column. Let \(C\) be an \(n \times q\) design matrix containing the covariates to adjust for. It should include a column of ones when an intercept is intended. For gene \(g\), the ordinary least-squares model is \(Y_g=C\beta_g+\varepsilon_g\), where \(Y_g\) is that gene’s expression vector.

For a full-column-rank design, the projection and residual-maker matrices are \(H=C(C^\mathsf{T}C)^{-1}C^\mathsf{T}\) and \(M=I-H\). Thus the matrix of residuals is \(E=MY\). When the same samples and design are used for every gene, the residuals are orthogonal to the covariates in \(C\): \(C^\mathsf{T}E=0\). With an intercept, each residual column also sums to zero, so its Pearson correlation can be calculated directly from cross-products:

A standard interpretation requires more than the ability to compute this number. The covariate effects should be adequately represented by the chosen linear model, residual variation should be meaningful rather than dominated by measurement or model failure, and the samples should be comparable across genes. For classical covariance estimation, independent sample errors with a common covariance matrix across samples give a useful reference model. If residual variances are zero, the correlation is undefined. Residual correlation itself does not require normally distributed errors, though normality and other assumptions matter for particular inferential procedures.

Worked Calculations

Why the Residual Matrix Matters

The residual correlation matrix has a useful mathematical property when every gene uses the same complete set of samples and the same design matrix. The residual cross-product matrix is \(E^\mathsf{T}E\), which is positive semidefinite: for any vector \(v\), \(v^\mathsf{T}E^\mathsf{T}Ev=(Ev)^\mathsf{T}(Ev)\geq 0\). Rescaling its rows and columns by positive residual standard deviations preserves this property, so the resulting correlation matrix is also positive semidefinite, provided no gene has zero residual variance.

This property depends on consistent construction. If different gene pairs use different samples because of pairwise missing-data deletion, the resulting pairwise correlations need not form a positive-semidefinite matrix. Likewise, fitting different covariate designs across genes does not generally produce one shared residual cross-product matrix. These issues matter when a downstream method requires a coherent matrix, rather than a list of pairwise values.

There is also a sample-size constraint. If the design has rank \(q\), residual vectors lie in a subspace of dimension at most \(n-q\). When \(p\) is large relative to \(n-q\), the residual correlation matrix can be singular, and the individual correlations may be noisy. Residualization does not create information or make a small sample sufficient for a large network. Estimation approaches such as shrinkage can address instability in a correlation matrix, but they do not repair a poorly specified covariate model.

What an Adjusted Edge Does and Does Not Mean

An edge based on residual correlation describes linear association remaining after the fitted effects of the included covariates have been removed gene by gene. It is not the same as an unadjusted marginal edge, which summarizes association in the original expression values. It is also not, by itself, evidence of conditional independence or a direct regulatory relationship. Residualizing measured covariates does not adjust one gene for all the other genes.

In the special case of a shared linear adjustment, the calculation is based on correlations between vectors made orthogonal to the covariate design. Its interpretation still relies on the specified model and measured variables. Other approaches to describing gene relationships condition on different quantities and answer different questions; the network’s edge definition must match the intended scientific target.

Covariate selection is therefore a scientific decision, not just a computational one. Adjusting for a known processing batch may help separate technical variation from biological co-expression. Adjusting for a biological feature can instead remove variation that is part of the scientific signal. If a covariate lies on a pathway connecting a biological exposure to gene expression, adjusting for it changes the target association. If covariate inclusion is informed by the observed expression correlations, the resulting network can also be difficult to interpret without accounting for that selection.

Diagnostics, Alternatives, and Common Pitfalls

Before interpreting edges, check that the design matrix has full column rank and that every fitted gene has enough observations to estimate its coefficients. Inspect whether categorical covariates have represented levels, whether continuous covariate effects appear plausibly linear, and whether outliers or high-leverage samples dominate fitted values. Plotting gene expression against important continuous covariates, and examining residual patterns against fitted values or batch, can reveal trends the model failed to remove. If a linear adjustment is inadequate, a more appropriate model may require transformed variables, nonlinear terms, or a different modeling strategy.

The main alternatives answer related but distinct questions. An unadjusted correlation network is simpler and may be appropriate when the target is total observed co-expression, but it retains technical and biological covariate effects. Residual correlation explicitly removes the fitted contribution of selected covariates. Pairwise regression can allow a different adjustment set for each gene pair, but then the resulting correlations may not form a coherent positive-semidefinite matrix. None of these choices is automatically best; state the target and the adjustment set.

• Common mistake: “The batch-adjusted correlation is high, so the genes must be directly biologically connected.” Correction: It is a high residual association after the chosen batch model; unmeasured factors, other biological processes, and indirect relationships may still explain it.
• Common mistake: “The covariate is measured, so removing it can only improve the network.” Correction: Adjustment changes the estimand and can remove meaningful biological variation if the covariate is part of the biology being studied.
• Common mistake: “Residuals are orthogonal to the covariate, so the model has removed every effect of it.” Correction: OLS guarantees sample orthogonality to the fitted columns, not that the chosen functional form captures the true relationship or that residuals are free of all covariate-related structure.