Worked Quantitative-Trait Burden Analysis

From Variant Data to a Quantitative-Trait Association

A burden analysis combines information across a selected set of variants into a participant-level score, then asks whether that score is associated with a quantitative phenotype. The combination is useful when the scientific question concerns the aggregate effect of a gene’s included variants rather than the separate effect of each rare variant. But the regression result is interpretable only after the score, adjustment variables, and coefficient units are clear.

This tutorial works through a hypothetical analysis from genotype indicators to an adjusted regression estimate and test. The variant set and weights are treated as already specified. The analysis tests the linear association between the resulting burden and the quantitative trait; it does not establish that the variants cause a change in the trait.

Definition: For participant \(i\), let \(G_{ij}\) be the genotype measure for included variant \(j\), and let \(w_j\) be its prespecified weight. A weighted burden score is \(B_i=\sum_{j=1}^{m} w_jG_{ij}\), where \(m\) is the number of included variants. A quantitative-trait burden regression models the phenotype \(Y_i\) using \(B_i\) and any specified covariates. The burden coefficient is the modeled change in mean phenotype per one-unit increase in burden score, conditional on those covariates.

Calculate the Participant-Level Score

A score depends on both the genotype coding and the weights. In the example below, the genotype entries are carrier indicators: \(G_{ij}=1\) if participant \(i\) carries an included alternate allele at variant \(j\), and \(G_{ij}=0\) otherwise. The two hypothetical variants have weights 1 and 2. Thus, a carrier of both has a score of 3. These weights are simply part of this example’s specified score; the numerical value of a score is not itself a biological unit.

Worked Example: Constructing Weighted Burden Scores

Suppose eight hypothetical participants have the carrier indicators shown below. For each participant, calculate \(B_i=1G_{i1}+2G_{i2}\).

Participant	Variant 1 carrier \(G_{i1}\)	Variant 2 carrier \(G_{i2}\)	Burden \(B_i\)
A	0	0	0
B	0	0	0
C	1	0	1
D	1	0	1
E	0	1	2
F	0	1	2
G	1	1	3
H	1	1	3

For participant G, for example, \(B_G=1(1)+2(1)=3\). For participant E, \(B_E=1(0)+2(1)=2\). The full score vector is \((0,0,1,1,2,2,3,3)\). These scores, rather than the separate carrier indicators, are the predictor in the one-score burden regression below.

Specify the Adjusted Regression

For one covariate \(C_i\), such as a measured technical or demographic variable, use the linear model

Y_i=\alpha+\beta_B B_i+\gamma C_i+\varepsilon_i.

Here \(\alpha\) is the intercept, \(\beta_B\) is the burden coefficient, and \(\gamma\) is the covariate coefficient. The usual exact small-sample \(t\) inference assumes independent observations, a correctly specified linear conditional mean, errors with mean zero and constant variance, and normally distributed errors conditional on the predictors. The design matrix must have full column rank: in particular, the burden score cannot be an exact linear combination of the intercept and covariates. In an actual analysis, the chosen covariates and score should be appropriate to the study design and specified before interpreting the test.

The null hypothesis for the aggregate association is \(H_0:\beta_B=0\); a two-sided alternative is \(H_a:\beta_B\ne 0\). This is a test of the burden’s linear association with the phenotype, conditional on the included covariates. It does not test each variant separately.

Key result: After adjusting both phenotype and burden for the intercept and covariates, the least-squares burden estimate is the cross-product of the residualized burden and residualized phenotype divided by the residualized burden sum of squares. This is equivalent to the burden coefficient in the full regression.

To see why, let \(Z\) be the design matrix containing the intercept and covariates, and let \(M_Z=I-Z(Z^TZ)^{-1}Z^T\) be the residual-maker matrix, assuming \(Z\) has full column rank. Then \(\widetilde B=M_ZB\) and \(\widetilde Y=M_ZY\) are the parts of the burden and phenotype not linearly explained by \(Z\). The least-squares objective, for fixed \(\beta_B\), is minimized over the nuisance coefficients by projecting \(Y-\beta_BB\) onto the columns of \(Z\). The remaining squared error is therefore \((Y-\beta_BB)^TM_Z(Y-\beta_BB)\). Differentiating this expression with respect to \(\beta_B\) and setting the derivative to zero gives

-2B^TM_ZY+2\beta_B B^TM_ZB=0, \qquad \widehat{\beta}_B =\frac{B^TM_ZY}{B^TM_ZB} =\frac{\widetilde B^T\widetilde Y}{\widetilde B^T\widetilde B}.

The adjustment matters: the estimate uses the variation in burden that remains after its linear relationship with the covariates has been removed, and the corresponding part of the phenotype. A raw burden–phenotype comparison can differ if burden and covariates are associated.

Fit the Model and Test the Burden

Worked Example: Covariate-Adjusted Quantitative-Trait Burden Regression

Use the burden scores above. Suppose \(C_i\) is a hypothetical centered covariate, and the quantitative phenotype \(Y_i\) is measured in assay units. All eight observations below are hypothetical.

Participant	Burden \(B_i\)	Covariate \(C_i\)	Phenotype \(Y_i\) (assay units)
A	0	-1	10.5
B	0	1	9.5
C	1	-1	10.5
D	1	1	13.5
E	2	-1	14.5
F	2	1	13.5
G	3	-1	14.5
H	3	1	17.5

Check the design. The covariate takes both values at each burden level, and the burden varies across participants. Thus, here the intercept, burden, and covariate are not perfectly collinear. Also, \(\bar B=(0+0+1+1+2+2+3+3)/8=1.5\), \(\bar C=0\), and \(\bar Y=13\).

Adjust for the covariate. The least-squares covariate slope is \(0.5\) assay units per covariate unit: within burden levels 0, 1, 2, and 3, the phenotype differences \(Y(C=1)-Y(C=-1)\) are respectively \(-1, 3, -1,\) and \(3\). The overall adjusted-regression covariate slope is \(0.5\), since burden is orthogonal to the centered covariate in this design and \(\widehat\gamma=\sum_i C_iY_i/\sum_i C_i^2=4/8=0.5\). Removing the intercept and this covariate gives \(\widetilde Y_i=Y_i-13-0.5C_i\). In participant order, the adjusted phenotype values are \((-2,-4,-2,0,2,0,2,4)\). Because burden is balanced across the two covariate values at every score level, its residualized values are simply \(B_i-\bar B=(-1.5,-1.5,-0.5,-0.5,0.5,0.5,1.5,1.5)\).

Calculate the burden slope. The adjusted burden sum of squares and cross-product are

\widetilde B^T\widetilde B =2(1.5^2)+2(0.5^2)+2(0.5^2)+2(1.5^2)=10,

and

\widetilde B^T\widetilde Y =3+6+1+0+1+0+3+6=20.

For instance, participant B contributes \((-1.5)(-4)=6\), not 2. Therefore \(\widehat{\beta}_B=20/10=2.00\) assay units per burden-score unit.

Recover the fitted model and uncertainty. The fitted intercept is \(\bar Y-\widehat{\beta}_B\bar B-\widehat\gamma\bar C=13-2(1.5)-0=10\). The fitted equation is \(\widehat Y_i=10+2B_i+0.5C_i\). In participant order, the fitted values are \(9.5,10.5,11.5,12.5,13.5,14.5,15.5,16.5\). The residuals are \(1,-1,-1,1,1,-1,-1,1\), so the residual sum of squares is \(RSS=8(1^2)=8\). There are \(8-3=5\) residual degrees of freedom, because three regression coefficients were estimated. Thus \(MSE=RSS/df=8/5=1.6\) assay-units squared.

The standard error of the burden estimate is

SE(\widehat{\beta}_B) =\sqrt{\frac{MSE}{\widetilde B^T\widetilde B}} =\sqrt{\frac{1.6}{10}} =0.4000\ \text{assay units per burden-score unit}.

For \(H_0:\beta_B=0\), the test statistic is \(t=2.00/0.4000=5.00\), with \(df=5\). The two-sided \(p\)-value is approximately \(0.0041\). Under the stated model assumptions, this is evidence against a zero linear burden coefficient, conditional on \(C\); it is not the probability that the null hypothesis is true.

A 95% confidence interval uses \(t_{0.975,5}\approx2.571\):

2.00\pm 2.571(0.4000) =2.00\pm1.028 \approx[0.972,\ 3.028]\ \text{assay units per burden-score unit}.

In repeated samples under the model, the interval procedure would cover the true conditional coefficient 95% of the time. For this hypothetical analysis, the estimated mean phenotype is 2.00 assay units higher for a one-unit increase in burden score, comparing participants with the same covariate value. The interval quantifies uncertainty in that coefficient under the model.

Translate the Coefficient to a Score Contrast

Worked Example: Comparing Two Burden Scores

Using the fitted coefficient \(2.00\) assay units per score unit, compare score 0 with score 3 while holding \(C\) fixed. The fitted mean difference is

\widehat Y(B=3,C)-\widehat Y(B=0,C) =\widehat{\beta}_B(3-0) =2.00(3)=6.00\ \text{assay units}.

Scaling the coefficient interval by the positive contrast of 3 gives \(3[0.972,3.028]=[2.916,9.084]\) assay units for this fitted mean difference. This comparison concerns the modeled score levels in the hypothetical data. It does not imply that changing a person’s burden score would cause the phenotype to change by this amount.

Interpretation and Common Pitfalls

A burden coefficient refers to the chosen score, including its genotype coding and weights. If all weights were multiplied by a positive constant, the numerical burden scores would scale by that constant and the coefficient per score unit would scale inversely; the fitted phenotype differences would be unchanged. Thus, report how the score was formed as well as the coefficient and its units. The estimate summarizes the aggregate linear association encoded by the score; it does not identify which included variant is associated, nor does it show that all variant effects are equal.

Using the raw instead of adjusted relationship. A raw correlation or slope does not generally equal the coefficient adjusted for covariates. Use the covariate-adjusted regression, or equivalently residualize both phenotype and burden against the same intercept and covariates.
Misreading the score unit. “Two assay units per score unit” is not “two assay units per variant.” The score unit is determined by the variant coding and weights.
Overstating the test. A small \(p\)-value is evidence against a zero aggregate coefficient under the model; it is not a probability that \(H_0\) is true and does not establish causation.
Ignoring model or design problems. A coefficient is not estimable if the burden is perfectly collinear with the covariates. Incorrect functional form, dependent observations, or substantial heteroskedasticity can also make the stated standard error and \(t\)-based inference unreliable.

Key takeaway: A quantitative-trait burden analysis proceeds from a specified weighted score to a covariate-adjusted linear regression. The coefficient estimates the conditional phenotype change per score unit, while its standard error and test quantify uncertainty under the regression assumptions.

Check Your Understanding

Use the definitions and calculations above to answer the following questions.

If a participant carries variant 1 but not variant 2 in the first example, what is that participant’s weighted burden score?
Why does the adjusted slope use residualized burden and phenotype rather than their unadjusted cross-product?
In the regression example, what are the residual degrees of freedom, and why are there three fitted coefficients?
State the estimated burden coefficient with its units and explain what is held fixed in its interpretation.
Does the two-sided \(p\)-value of approximately \(0.0041\) give the probability that the null hypothesis is true? Explain briefly.