Why use a burden test for a quantitative phenotype?
A rare-variant set may contain many variants, each too uncommon to estimate its association reliably on its own. A burden test replaces their separate genotype columns with one prespecified weighted score for each participant, then tests whether that score is associated with a quantitative outcome. The reduction to one predictor is useful, but it also imposes a substantive assumption: the variants contribute in a common direction after the chosen allele coding and weights are applied.
For a quantitative phenotype, linear regression gives the burden coefficient a direct scale: the expected change in the phenotype for a one-unit increase in the weighted burden, conditional on the covariates. The covariate-adjusted score test asks whether that coefficient is zero while accounting for the association between the burden and the covariates.
Model, assumptions, and adjusted score
The genotype dosage coding and variant set define what enters \(B\); the weights determine how strongly each variant contributes. For example, if genotypes count copies of the allele of interest, a positive weight makes an additional copy increase the burden. Weights should be chosen without using the phenotype association results being tested. A rescaling of all weights rescales \(B\) and inversely rescales \(\beta\), while leaving the fitted values and test statistic unchanged.
The usual model treats \(C\), \(B\), and the weights as fixed conditional on the observed design. It requires a correctly specified conditional mean, independent errors with mean zero and common variance \(\sigma^2\), and no unmodelled dependence that invalidates that error structure. Normal errors, \(\varepsilon\sim N(0,\sigma^2 I)\), support the usual likelihood derivation; the score statistic with estimated variance is generally calibrated by a large-sample chi-square approximation. For exact small-sample inference under normality, the full-model coefficient estimate divided by its full-model standard error has a t distribution under the null, as described below.
Assume \(C\) has full column rank and includes an intercept. Define the residual-maker matrix \(M_C=I-C(C'C)^{-1}C'\). It projects a vector onto the part orthogonal to the columns of \(C\). The burden must have some variation left after covariate adjustment: \(Q=B'M_CB>0\). Otherwise, the burden is fully explained by the covariates and its separate coefficient cannot be estimated.
This statistic tests \(H_0:\beta=0\) against a two-sided alternative. A large value of \(S\) indicates that the adjusted burden score is large relative to its null-model standard error. The resulting p-value measures how unusual a statistic at least this large would be under the null model and its calibration; it is not the probability that the null hypothesis is true.
Why residualizing the burden gives the adjusted estimate
The adjustment can be justified by separating the burden and outcome from their linear relationships with the covariates. For any fixed \(\beta\), the least-squares choice of \(\gamma\) fits the covariate component of \(Y-B\beta\). The remaining residual vector is therefore \(M_C(Y-B\beta)\). The part of the residual sum of squares that depends on \(\beta\) is
Differentiating this quadratic with respect to \(\beta\) gives \(-2U+2\beta Q\). Setting the derivative to zero yields \(\hat{\beta}=U/Q\). This is the residualization principle: regress the outcome and burden on the covariates, then compare their covariate-residualized parts. The score \(U\) is the slope’s numerator, while \(Q\) measures the adjusted burden variation available to estimate it.
Worked Example: A score test with an orthogonal covariate
Suppose six hypothetical individuals have quantitative outcomes \(Y=(7,2.5,4,2,7.5,7)'\), a continuous covariate \(x=(-2,-1,0,0,1,2)'\), and a weighted burden \(B=(1,0,-1,-1,0,1)'\). Use an intercept and \(x\) as the columns of \(C\). The values were constructed from \(Y=5+0.5x+2B+e\), with \(e=(1,-2,1,-1,2,-1)'\); direct substitution gives the stated outcomes. Both \(B\) and \(e\) have zero sum, \(x'B=0\), and \(x'e=B'e=0\), so the covariate and error are orthogonal to the burden in this example.
The burden is already orthogonal to the intercept and \(x\), so \(M_CB=B\), and \(Q=B'B=4\). Removing the intercept and covariate from \(Y\) leaves \(2B+e\). Thus \(U=B'(2B+e)=2(4)+0=8\). The null-model residual sum of squares is \((2B+e)'(2B+e)=4(4)+12=28\), because \(B'e=0\) and \(e'e=12\). The null model has \(6-2=4\) residual degrees of freedom, so \(s_0^2=28/4=7\).
The score statistic is \(S=8^2/(7\cdot4)=2.286\), giving an approximate chi-square p-value of \(0.131\) on 1 degree of freedom. At a 0.05 significance level, this score test does not reject \(H_0\). This small illustrative sample is not a setting in which the large-sample calibration should be relied on as an exact one.
The fitted burden coefficient is \(\hat{\beta}=8/4=2\) outcome units per burden unit. In the full regression, the residual sum of squares is \(e'e=12\), with \(6-2-1=3\) residual degrees of freedom. Its coefficient standard error is \(\sqrt{(12/3)/4}=1\). Under normal errors, a 95% t interval with 3 degrees of freedom is \(2\pm3.182(1)\), or approximately \((-1.182,5.182)\) outcome units per burden unit. This wide interval reflects substantial uncertainty; the fitted coefficient is an estimate, not proof of an effect.
Covariate adjustment can change the burden estimate
If burden and covariates are associated, an unadjusted regression generally answers a different question from the adjusted model. The following example shows the calculation explicitly. Its small sample is for arithmetic illustration, not a recommendation to use asymptotic inference with so few individuals.
Worked Example: Adjustment changes the estimated effect
For five hypothetical individuals, let \(x=(0,1,2,3,4)'\), \(B=(0,0,1,1,0)'\), and \(Y=(11,11,16,20,18)'\). Use an intercept and \(x\) as covariates. These outcomes satisfy \(Y=10+2x+3B+e\), where \(e=(1,-1,-1,1,0)'\). The residual vector has zero sum and \(x'e=B'e=0\), so it is orthogonal to the intercept, covariate, and burden.
The unadjusted centered sums are \(B'M_1B=2-5(0.4)^2=1.2\), \(B'M_1x=5-5(0.4)(2)=1\), and \(B'M_1Y=36-5(0.4)(15.2)=5.6\), where \(M_1\) removes the intercept. The unadjusted coefficient is \(5.6/1.2=4.667\) outcome units per burden unit.
For adjustment, \(x'M_1x=30-5(2)^2=10\). Consequently, \(Q=B'M_CB=1.2-1^2/10=1.1\). Also \(x'M_1Y=175-5(2)(15.2)=23\), so \(U=B'M_CY=5.6-(1)(23)/10=3.3\). The adjusted coefficient is \(3.3/1.1=3.000\) outcome units per burden unit. It recovers the coefficient in the specified data-generating equation because the residual is orthogonal to the predictors.
The null residual sum of squares is \(3^2(1.1)+e'e=9.9+4=13.9\), with \(5-2=3\) degrees of freedom. Therefore \(s_0^2=13.9/3=4.633\) and \(S=3.3^2/(4.633\cdot1.1)=2.136\). The approximate chi-square p-value is \(0.144\). The large difference between the unadjusted estimate, 4.667, and adjusted estimate, 3.000, illustrates why covariate adjustment matters; neither estimate alone establishes causation.
From a weighted score to its interpretation
The coefficient’s units depend on both the phenotype and the burden construction. If the outcome is measured in milligrams per deciliter and a burden unit is one weighted dosage unit, \(\beta\) is measured in milligrams per deciliter per weighted dosage unit. It is not automatically the effect of any single variant. The score test addresses association with the aggregate burden defined by the set, dosage coding, and weights.
Worked Example: Constructing and interpreting a weighted burden
Consider four hypothetical individuals with two dosage columns \(G_1=(0,1,0,1)'\) and \(G_2=(0,0,1,1)'\). Choose fixed weights \(w_1=2\) and \(w_2=1\), giving \(B=2G_1+G_2=(0,2,1,3)'\). Let the quantitative outcomes be \(Y=(1,4,2,7)'\), and use an intercept only.
The means are \(\bar B=1.5\) and \(\bar Y=3.5\). The adjusted burden sum of squares is \(Q=\sum_i(B_i-\bar B)^2=2.25+0.25+0.25+2.25=5\). The score is \(U=\sum_i(B_i-\bar B)(Y_i-\bar Y)=3.75+0.25+0.75+5.25=10\). Hence \(\hat{\beta}=10/5=2\) outcome units per weighted burden unit.
The null residual sum of squares is \(21\), so \(s_0^2=21/(4-1)=7\). The score statistic is \(S=100/(7\cdot5)=2.857\), with approximate chi-square p-value \(0.091\). The full-model residuals are \((0.5,-0.5,-0.5,0.5)'\), giving residual sum of squares \(1\), residual variance estimate \(1/(4-2)=0.5\), and coefficient standard error \(\sqrt{0.5/5}=0.3162\). The 95% t interval with 2 degrees of freedom is \(2\pm4.303(0.3162)\), approximately \((0.639,3.361)\) outcome units per weighted burden unit. The weight of 2 for \(G_1\) is part of the definition of that unit; it does not mean the estimated coefficient is an individual-variant effect.
Diagnostics, alternatives, and limitations
Before interpreting a result, check that the burden was constructed as intended and that its weights and allele directions are explicit. Inspect \(Q\): a value near zero means little burden variation remains after covariate adjustment, making the coefficient unstable. Check for influential observations, nonlinearity in the conditional mean, unequal residual spread, and dependence among individuals. These issues can undermine the standard errors or test calibration even when the algebra is correct. If residual spread or dependence is not adequately represented by the model, the simple \(s_0^2Q\) variance calculation is not justified.
A burden test is most naturally suited to a set for which the chosen coding and weights produce a meaningful common-direction aggregate signal. Opposing variant effects can cancel in the weighted sum. The test also compresses the set to one dimension: it cannot identify which variant drives an association, and a significant burden result does not demonstrate that every included variant has an effect.
| Approach | What it tests or estimates | Main distinction |
|---|---|---|
| Unadjusted burden regression | Association between outcome and burden without covariates | Can differ from the conditional association when burden is related to covariates. |
| Covariate-adjusted burden score test | Whether one prespecified burden coefficient is zero, conditional on covariates | Uses one aggregate predictor and a score statistic evaluated under the null. |
| Full-model coefficient test | Whether the fitted burden coefficient differs from zero, using the full-model standard error | Provides a direct coefficient estimate and interval; its finite-sample t calibration differs from the plug-in score calibration. |
| Variant-by-variant regression | Separate coefficients for individual variants | Does not impose the same single aggregate effect, but rare alleles may give imprecise estimates. |
The score statistic and the full-model coefficient test target the same null in this linear model, but their variance estimates and finite-sample calibrations are not interchangeable. Under normal errors, the full-model coefficient divided by its full-model standard error has a t distribution with \(n-\operatorname{rank}(C)-1\) degrees of freedom under the null. The plug-in score statistic above instead uses the null-model variance estimate and a large-sample chi-square reference. In small samples, report the calibration used and avoid treating an asymptotic p-value as exact.
Common mistake
Wrong reasoning: “The burden coefficient from a regression without covariates is the adjusted gene-set effect, because the same burden vector was used.” The fact that the predictor is unchanged does not make the estimand unchanged. If \(B\) is associated with a covariate that predicts \(Y\), the unadjusted slope includes that covariate-related pattern.
Correction: Residualize both outcome and burden against the covariate design, or fit the full model directly. The adjusted coefficient is \(B'M_CY/(B'M_CB)\), provided \(B'M_CB>0\). State the burden’s construction and report the effect in outcome units per burden unit; do not describe it as a per-variant effect.
Check Your Understanding
Try each question, then check the solution immediately below it.
- If \(B'M_CB=0\), can the adjusted burden coefficient be estimated separately from the covariate effects? Solution: No. The burden has no residual variation after projection onto the covariates, so the denominator of \(\hat{\beta}=B'M_CY/(B'M_CB)\) is zero and the coefficient is not separately identifiable.
- A fitted coefficient is 1.8 millimeters of phenotype per burden unit. What does it describe? Solution: Conditional on the covariates, a one-unit increase in the specifically defined weighted burden corresponds to an estimated 1.8-millimeter increase in the conditional mean phenotype. It is not necessarily the effect of one variant.
- Why might an unadjusted burden estimate differ from an adjusted estimate? Solution: If burden is associated with covariates that predict the outcome, the unadjusted regression includes variation attributable to those covariates. Adjustment uses the parts of burden and outcome left after accounting for the covariates.
- Does a small score-test p-value show that every variant in the set affects the phenotype in the same direction? Solution: No. It is evidence against a zero coefficient for the prespecified aggregate burden under the model. It neither identifies individual effects nor establishes that every variant contributes.