Why Residualize Both Variables?
In a quantitative-trait GWAS, the SNP coefficient is usually estimated while accounting for covariates such as age, sex, or ancestry-related variables. A multiple regression does this in one step. The Frisch–Waugh–Lovell theorem gives an equivalent view: remove the part explained by the covariates from both the phenotype and the genotype, then regress one remainder on the other. This view clarifies what “adjusted for covariates” means algebraically.
The equivalence is exact for the coefficient when the same observations and covariates are used in both regressions. It is not a claim that covariates can be ignored, nor that residualization changes the scientific model. The full model remains the reference for interpretation and standard errors.
Here, \(P_C\) projects a vector onto the column space of the covariates, and \(M_C\) subtracts that projection. Thus \(M_Cy\) contains the phenotype variation not linearly explained by \(C\); \(M_Cg\) contains the genotype variation not linearly explained by \(C\). Both residual vectors are orthogonal to every column of \(C\). The intercept in \(C\) ensures that each residual vector also has mean zero.
The model’s algebraic coefficient equivalence requires that \(C\) have full column rank and that \(g\) not lie entirely in the column space of \(C\). Equivalently, \(g'M_Cg>0\). These conditions ensure that the covariate coefficients and the SNP coefficient are identifiable. They do not, by themselves, guarantee valid uncertainty estimates. For conventional model-based inference, one commonly assumes errors with conditional mean zero given \(C\) and \(g\), and conditional covariance \(\sigma^2 I\); exact small-sample \(t\)-inference additionally assumes conditionally multivariate normally distributed errors. The coefficient identity itself does not require normally distributed errors.
The Frisch–Waugh–Lovell Derivation
The theorem follows directly from the least-squares normal equations. For any fixed value of \(\beta\), the best-fitting covariate contribution is the projection of \(y-g\beta\) onto the columns of \(C\). The part left after that projection is \(M_C(y-g\beta)\). Therefore, choosing \(\beta\) in the full model is equivalent to choosing it to minimize the squared length of this remaining vector.
Differentiate the last squared length with respect to \(\beta\), set the derivative to zero, and solve. Since \(M_C\) is symmetric and idempotent, the resulting expression is:
The middle expression is precisely the no-intercept least-squares slope from regressing the residualized phenotype on the residualized genotype. The denominator must be positive; if it is zero, the genotype is completely explained by the covariates and the SNP coefficient cannot be separately estimated. With an intercept in \(C\), both residual vectors have mean zero, so the no-intercept residual regression is appropriate. Adding an intercept to that residual regression leaves its slope unchanged.
Worked Example: Residualization With an Intercept Only
Suppose four hypothetical individuals have dosage values \(g=(2,2,0,0)\) and phenotype values \(y=(1,0,0,-1)\). Adjust only for an intercept. The genotype mean is \(\bar g=1\), and the phenotype mean is \(\bar y=0\), so the residuals after removing the intercept are \(M_Cg=(1,1,-1,-1)\) and \(M_Cy=(1,0,0,-1)\).
The adjusted slope is the cross-product of these residuals divided by the squared length of the residualized genotype:
In this hypothetical sample, the fitted phenotype difference associated with a one-unit increase in dosage is \(0.5\) phenotype units. The example illustrates the arithmetic of the adjusted coefficient; with an intercept as the only covariate, it is also the ordinary unadjusted slope.
Worked Example: Adjusting for a Covariate
Consider four hypothetical observations with covariate values \(x=(-1,-1,1,1)\), genotype dosages \(g=(0,1,1,2)\), and phenotypes \(y=(2,3,5,8)\). The covariate matrix contains an intercept and \(x\). The covariate-only fitted genotype values are the group means \(0.5,0.5,1.5,1.5\). Thus the residualized genotype is \(M_Cg=(-0.5,0.5,-0.5,0.5)\).
The covariate-only fitted phenotype values are \(2.5,2.5,6.5,6.5\), giving \(M_Cy=(-0.5,0.5,-1.5,1.5)\). The FWL coefficient is:
The dosage variation remaining after accounting for \(x\) is associated with a fitted increase of \(2\) phenotype units per dosage unit. The raw, intercept-only slope would be \(3\): \(\bar g=1\), \(\bar y=4.5\), so its numerator is \(2.5+0+0+3.5=6\), and its denominator is \(1+0+0+1=2\). The difference between \(3\) and \(2\) shows why the covariate adjustment matters in these data.
What the Residual Regression Does—and Does Not Do
The residual-regression interpretation is a partial association. It compares phenotype and genotype variation that remains after their linear relationships with the specified covariates have been removed. It does not establish that changing genotype causes a phenotype change. Nor does residualization repair an inappropriate covariate set, nonlinear covariate relationships, or problems in the underlying data.
The theorem also explains why residualizing only one variable is generally insufficient. If only \(y\) is residualized while the original \(g\) is retained, the resulting regression does not, in general, use the genotype variation independent of \(C\). The adjusted coefficient requires both \(M_Cy\) and \(M_Cg\). The same covariate matrix and analysis sample must be used in both residualizations.
For the coefficient, direct multiple regression and residualization are equivalent. Direct regression is usually the clearest way to obtain the full model output. Residualization is useful for understanding the coefficient, visualizing the partial relationship, or implementing repeated calculations efficiently. For standard errors and tests, use the full model or an equivalent method that accounts for the covariate degrees of freedom and the assumed error structure. A residual-on-residual slope alone does not justify treating its residual regression as a new independent-sample model.
Worked Example: A Rank Failure
Suppose a hypothetical SNP dosage is \(g=(0,1,2,1)\), and \(C\) contains an intercept and a covariate equal to that dosage. Because \(g\) is itself a column of \(C\), its projection onto \(C\) is exactly \(g\), so \(M_Cg=(0,0,0,0)\). Consequently:
The FWL coefficient formula has a zero denominator. In the full model, the SNP column duplicates a covariate column, so infinitely many combinations of their coefficients produce the same fitted values; the separate SNP effect is not identifiable. The right response is not to report a residualized slope, but to recognize that the model cannot distinguish the two terms.
Diagnostics and Common Pitfalls
A practical check is to verify that each residualized vector is orthogonal to the columns of \(C\), up to numerical precision: \(C'M_Cy=0\) and \(C'M_Cg=0\). Inspect the residualized genotype’s spread as well. A near-zero \(g'M_Cg\) indicates little genotype variation remains after adjustment; the coefficient may be numerically unstable even if the denominator is not exactly zero. Also check that both residualizations use identical rows and covariate coding.
- Residualizing only the phenotype: “The covariate has been removed from \(y\), so I can use the original genotype.” This is generally wrong because the covariate-related component of \(g\) remains. Residualize both vectors against the same \(C\).
- Using the wrong denominator: The slope denominator is \((M_Cg)'(M_Cg)\), not the squared length of the phenotype residuals.
- Ignoring rank deficiency: If \(g\) is completely explained by \(C\), the adjusted SNP coefficient is not estimable.
- Using residual regression output for inference without checking it: The coefficient identity does not automatically make every standard error or degrees-of-freedom calculation from a separate residual regression valid.
A useful comparison is with unadjusted regression: it regresses the original phenotype on the original genotype and does not remove covariate-related variation. Multiple regression adjusts jointly; residualization expresses that same adjustment as two projections followed by a slope calculation. Neither approach makes the result causal by itself. Both depend on the model specification, the observations included, and sufficient genotype variation after adjustment.
Check Your Understanding
Try each question, then compare your reasoning with the solutions.
- In one sentence, state what each residualized vector represents.
- Suppose \(M_Cg=(1,-1,0)\) and \(M_Cy=(2,0,1)\). Calculate the FWL slope.
- Why can a SNP coefficient not be estimated if \(M_Cg\) is the zero vector?
- True or false: residualizing \(y\) against \(C\) and regressing it on the original \(g\) always gives the adjusted coefficient. Explain.
- What is the role of the full model when obtaining standard errors after using residualization to understand the coefficient?
Solutions. (1) \(M_Cy\) is the phenotype after removing its linear projection on the covariates; \(M_Cg\) is the genotype after removing its linear projection on those same covariates. (2) The numerator is \(1(2)+(-1)(0)+0(1)=2\), and the denominator is \(1^2+(-1)^2+0^2=2\), so \(\hat{\beta}=2/2=1\). (3) Its squared length \(g'M_Cg\) is zero, so there is no remaining genotype variation with which to distinguish the SNP coefficient from the covariate effects. (4) False. The original genotype still contains covariate-related variation; both \(y\) and \(g\) must be residualized against the same \(C\). (5) The full model, or a valid algebraically equivalent calculation, accounts for the covariate degrees of freedom and the assumed error structure; the residual-regression coefficient identity alone does not guarantee valid uncertainty estimates.