Covariate-Adjusted Burden Scores By Residualization

Why Residualize Both Quantities?

A burden test asks whether a person-level genetic burden is associated with a quantitative phenotype. The phenotype may also vary with covariates such as age, sex, or ancestry measures. The adjusted question is whether the burden contributes information about the phenotype after accounting for those covariates. Residualization makes that adjustment visible: remove the covariate-predicted component from the phenotype, remove the same covariate component from the burden, and compare what remains.

The burden $b_i$ for individual $i$ is treated here as a single numerical predictor, constructed from that individual's variant data and any chosen weights. The focus is the linear score test for its association with a quantitative phenotype $y_i$, conditional on a specified set of covariates. The calculation is a test of association, not evidence that the burden causes a phenotype difference.

Definition: Let $y$ be the $n$-vector of phenotypes, $b$ the $n$-vector of burdens, and $C$ the $n \times q$ design matrix of covariates, including an intercept. Assume $C$ has full column rank. Define the residual-maker matrix $$ M_C = I - C(C^{\mathsf T}C)^{-1}C^{\mathsf T}. $$ The residualized phenotype and burden are $y_{\perp}=M_Cy$ and $b_{\perp}=M_Cb$. Each is what remains after its least-squares linear projection onto the columns of $C$ has been removed.

The matrix $M_C$ is symmetric and idempotent: $M_C^{\mathsf T}=M_C$ and $M_C^2=M_C$. It also satisfies $M_CC=0$, so the residualized quantities are orthogonal to every included covariate. If $C$ contains an intercept, both residual vectors have mean zero. These properties are what make the adjusted score calculation work.

The Adjusted Linear Score Test

The model for the quantitative trait is $$ y=C\gamma+\beta b+\varepsilon, $$ where $\gamma$ is a vector of covariate coefficients, $\beta$ is the burden coefficient, and $\varepsilon$ is an error vector. The null hypothesis is $H_0:\beta=0$; the alternative for a two-sided test is $H_a:\beta\ne 0$. For the standard linear-model score calculation, assume independent errors with mean zero and common variance $\sigma^2$; normal errors give the usual Gaussian likelihood derivation. The burden and covariates are treated as fixed in the test, conditional on the observed design.

Under $H_0$, fit the covariate-only model. The residual vector from that fit is $M_Cy$. The score for adding $b$, after accounting for the nuisance covariate coefficients, is proportional to the inner product of the residualized burden and residualized phenotype:

U=b_{\perp}^{\mathsf T}y_{\perp} =b^{\mathsf T}M_Cy.

The equality follows from symmetry and idempotence: $b_{\perp}^{\mathsf T}y_{\perp}=b^{\mathsf T}M_C^{\mathsf T}M_Cy=b^{\mathsf T}M_Cy$. The score is large in magnitude when the part of the burden not explained by covariates aligns strongly, positively or negatively, with the part of the phenotype not explained by covariates.

Theorem: Under $H_0$, with $\operatorname{Var}(\varepsilon)=\sigma^2I$, the score $U=b^{\mathsf T}M_Cy$ has mean zero and variance $\sigma^2 b^{\mathsf T}M_Cb=\sigma^2 b_{\perp}^{\mathsf T}b_{\perp}$. If $\sigma^2$ is known and the errors are normal, $U^2/(\sigma^2 b_{\perp}^{\mathsf T}b_{\perp})$ has a $\chi^2$ distribution with 1 degree of freedom under $H_0$. With unknown variance, the usual score statistic replaces $\sigma^2$ by its estimate from the covariate-only model and is compared with a $\chi^2_1$ distribution as a large-sample approximation.

The variance result follows directly from the model under $H_0$. Because $M_CC=0$, the covariate mean contributes nothing to $U$, leaving $U=b^{\mathsf T}M_C\varepsilon$. Its variance is $$ \operatorname{Var}(U)=b^{\mathsf T}M_C(\sigma^2I)M_Cb =\sigma^2 b^{\mathsf T}M_C^2b =\sigma^2 b^{\mathsf T}M_Cb. $$ Thus the standardized squared score is

Q=\frac{(b_{\perp}^{\mathsf T}y_{\perp})^2} {\hat{\sigma}_0^2\,b_{\perp}^{\mathsf T}b_{\perp}}, \qquad \hat{\sigma}_0^2=\frac{y^{\mathsf T}M_Cy}{n-q}.

Here $\hat{\sigma}_0^2$ is the residual mean square from the covariate-only model, with $n-q$ residual degrees of freedom. A nominal two-sided p-value is the upper-tail probability $P(\chi^2_1\ge Q)$. When the variance is estimated, this reference is generally an asymptotic calibration, not an exact finite-sample guarantee. In particular, a very small sample or very few carriers can make the approximation unreliable.

Worked Calculation

Worked Example: Residualizing a Burden and a Phenotype

Suppose six hypothetical individuals have a quantitative phenotype $y$, a burden $b$, and a continuous covariate $x$. The covariate design includes an intercept and $x$.

Individual	$x$	$b$	$y$
1	-2	0	7
2	-1	1	8
3	0	0	9
4	0	1	9
5	1	1	12
6	2	0	15

The covariate values sum to zero, so the intercept and $x$ are orthogonal. The burden values sum to 3, giving mean $0.5$, and their inner product with $x$ is $(-1)(1)+(1)(1)=0$. Consequently, residualizing $b$ against the intercept and $x$ only subtracts its mean:

b_{\perp}=(-0.5,\ 0.5,\ -0.5,\ 0.5,\ 0.5,\ -0.5)^{\mathsf T}.

The phenotype values follow $y=10+2x+e$, where $e=(1,0,-1,-1,0,1)^{\mathsf T}$. Since $\sum e_i=0$ and $\sum x_ie_i=0$, these are exactly the covariate-only residuals: $y_{\perp}=e$. The adjusted score is

U=b_{\perp}^{\mathsf T}y_{\perp} =(-0.5)(1)+(0.5)(0)+(-0.5)(-1) +(0.5)(-1)+(0.5)(0)+(-0.5)(1) =-1.

The residualized burden sum of squares is $b_{\perp}^{\mathsf T}b_{\perp}=6(0.25)=1.5$. The covariate-only residual sum of squares is $1^2+0^2+(-1)^2+(-1)^2+0^2+1^2=4$. There are $q=2$ covariate columns, so $\hat{\sigma}_0^2=4/(6-2)=1$. Therefore,

Q=\frac{(-1)^2}{1(1.5)}=0.6667.

The nominal $\chi^2_1$ upper-tail p-value is about 0.414. This small hypothetical data set illustrates the arithmetic, not a reliable application of the large-sample calibration: with only six observations, that approximation should not be relied on for inference.

Why Both Residuals Matter

Residualizing both quantities is not an arbitrary recipe. It gives the same burden coefficient as fitting the covariates and burden together by least squares. To see this, minimize the residual sum of squares over $\gamma$ for a fixed $\beta$. The part left after fitting $C$ is $M_C(y-\beta b)$. The remaining sum of squares is therefore

(y-\beta b)^{\mathsf T}M_C(y-\beta b) =y^{\mathsf T}M_Cy-2\beta b^{\mathsf T}M_Cy +\beta^2b^{\mathsf T}M_Cb.

Differentiating with respect to $\beta$ and setting the derivative to zero gives $$ -2b^{\mathsf T}M_Cy+2\beta b^{\mathsf T}M_Cb=0, $$ so, provided $b^{\mathsf T}M_Cb>0$, $$ \hat{\beta}=\frac{b^{\mathsf T}M_Cy}{b^{\mathsf T}M_Cb} =\frac{b_{\perp}^{\mathsf T}y_{\perp}}{b_{\perp}^{\mathsf T}b_{\perp}}. $$ This is the partial-regression result: the association is based on what remains in both variables after the same covariates have been removed. The score test uses the numerator of this expression and its null variance, rather than requiring the burden coefficient to be estimated under the alternative.

Worked Example: The Denominator Must Also Be Adjusted

Use the six-person data above. If the phenotype is residualized but the burden is not, the numerator $b^{\mathsf T}y_{\perp}$ still happens to equal $-1$: $0(1)+1(0)+0(-1)+1(-1)+1(0)+0(1)=-1$. This works because $y_{\perp}$ is orthogonal to the covariates. But the unadjusted burden sum of squares is $b^{\mathsf T}b=3$, not the adjusted value $b_{\perp}^{\mathsf T}b_{\perp}=1.5$.

Using the raw burden denominator would produce $$ Q_{\mathrm{wrong}}=\frac{(-1)^2}{1(3)}=0.3333, $$ rather than the correct $Q=0.6667$. The mismatch is not a harmless alternative calculation: it uses the total burden variation, including the component explained by covariates, to measure the uncertainty in an adjusted association. The score test requires the residualized burden sum of squares in its denominator.

When the Adjusted Test Has No Information

A test requires variation in the burden that is not already explained by the covariates. In matrix terms, $b_{\perp}^{\mathsf T}b_{\perp}$ must be positive. If $b$ lies exactly in the column space of $C$, then $M_Cb=0$. Its adjusted score variance is zero, and the burden coefficient cannot be separated from the covariate effects in this model. This is a lack of identifiability, not evidence for no association.

Worked Example: A Burden Fully Explained by a Covariate

Suppose four hypothetical individuals have a binary covariate $z=(0,0,1,1)^{\mathsf T}$, and the burden is $b=(0,0,1,1)^{\mathsf T}$. Let $C$ contain an intercept and $z$. Since $b=z$, the burden is already a column of $C$. Regressing $b$ on $C$ gives fitted values $(0,0,1,1)^{\mathsf T}$, hence $$ b_{\perp}=b-C\hat{\delta}=(0,0,0,0)^{\mathsf T}, \qquad b_{\perp}^{\mathsf T}b_{\perp}=0. $$

The score statistic would require division by zero, so it is undefined. Every burden difference in these data coincides perfectly with a difference in $z$; the data cannot distinguish their associations with the phenotype. The appropriate conclusion is that this design does not identify a separate burden effect, not that the burden has passed a test of no association.

Assumptions and Practical Checks

Before interpreting a residualized score, check that the covariate matrix includes the intended adjustment variables and an intercept, has full column rank, and has fewer columns than observations. The linear model assumes the mean phenotype is represented by the included covariate terms and the burden term; the standard score variance also assumes independent, constant-variance errors. Normality supports the Gaussian derivation, while the estimated-variance chi-square reference is an approximation whose quality depends on the sampling situation. A poorly calibrated reference distribution cannot be repaired merely by residualizing.

Residualization adjusts for the covariates included in $C$; it does not automatically correct for omitted confounding, an inappropriate mean model, or other violations of the assumptions. It also does not change the burden definition: the mask, variant weights, and coding determine $b$ before this adjustment is performed. The resulting test concerns the one burden vector supplied to the model.

Key takeaway: For a quantitative-trait burden score test with covariates, compute $y_{\perp}=M_Cy$ and $b_{\perp}=M_Cb$, then use $U=b_{\perp}^{\mathsf T}y_{\perp}$ and variance estimate $\hat{\sigma}_0^2 b_{\perp}^{\mathsf T}b_{\perp}$. If the residualized burden has zero variance, the adjusted burden effect is not identifiable.

Check Your Understanding

Answer each question using the residualization framework above.

What does $M_Cb$ represent, and why is it orthogonal to the columns of $C$?
In the score statistic, why is the denominator based on $b_{\perp}^{\mathsf T}b_{\perp}$ rather than $b^{\mathsf T}b$?
State the null hypothesis and the assumptions used to obtain the score's variance.
If the residualized burden is zero, what does that imply about whether the burden effect can be tested separately from the covariates?
Why should the nominal chi-square p-value from the six-person example not be treated as dependable evidence?

Individual	\(x\)	\(b\)	\(y\)
1	-2	0	7
2	-1	1	8
3	0	0	9
4	0	1	9
5	1	1	12
6	2	0	15