Why a Burden Score Has a One-Degree-of-Freedom Test
A gene may contain many rare variants, but a burden test combines their genotypes into one prespecified predictor. The linear-model score test asks whether that single weighted burden is associated with a quantitative trait. Its central calculation compares the burden’s alignment with the phenotype to the variation expected under no association.
This tutorial considers an intercept-only model: it does not include additional clinical, ancestry, or technical covariates. The variant set and weights are treated as fixed for the test. A later treatment of covariate adjustment changes how the phenotype and burden are handled, so the formulas here should not be read as already adjusting for such variables.
For example, if two variants have dosages \(G_{i1}\) and \(G_{i2}\), and their weights are \(2\) and \(0.5\), then \(B_i=2G_{i1}+0.5G_{i2}\). A larger value means a larger value of this particular weighted sum; it does not, by itself, establish that every included variant has the same biological effect.
Model and Hypotheses
Let \(Y_i\) be a quantitative phenotype for individual \(i\), measured in its stated units. Conditional on the observed burdens, assume the model
Here \(\alpha\) is the intercept, \(\beta_B\) is the change in the conditional mean phenotype per one-unit increase in burden, and \(\sigma^2\) is the common error variance. The score test is for \(H_0:\beta_B=0\) against \(H_a:\beta_B\ne 0\). The two-sided alternative is appropriate unless a directional alternative was specified before examining the data.
The usual model assumptions are independent individuals, a linear conditional mean in this burden, and normally distributed errors with common variance for the exact normal-theory calculations below. The burden must vary across individuals; if every \(B_i\) is identical, its effect cannot be estimated separately from the intercept. The weights and variant inclusion rule must be specified independently of the phenotype for the ordinary null calibration to apply as stated.
Define the centered sums
The intercept is not assumed known. Under \(H_0\), its fitted value is \(\bar{Y}\), so the null residual for individual \(i\) is \(Y_i-\bar{Y}\). Centering the burden accounts for estimating that intercept: evidence for a burden effect comes from the relationship between deviations in burden and deviations in phenotype, not from their overall levels.
Deriving the Score and Its Null Variance
A score is the derivative of the log-likelihood with respect to the parameter being tested, evaluated at its null value. For the normal model, the log-likelihood, apart from a constant, is
Differentiating with respect to \(\beta_B\) gives the score at a general parameter value. At the null fit, \(\beta_B=0\) and \(\hat{\alpha}_0=\bar{Y}\), so the score is
The last equality holds because \(\sum_i(Y_i-\bar{Y})=0\), so replacing \(B_i\) by \(B_i-\bar{B}\) leaves the sum unchanged. Under the null, \(Y_i-\bar{Y}\) are centered normal errors. The null variance of the score, accounting for the unknown intercept, is therefore
To see why the centered sum \(S_{BB}\) appears, write the score as \(\sigma^{-2}\sum_i(B_i-\bar{B})\varepsilon_i\) under the null. Independence and common variance give
Thus, if \(\sigma^2\) is known, the standardized score is normal under \(H_0\); its square has a chi-square distribution with one degree of freedom. There is one degree of freedom because the test concerns one scalar coefficient, \(\beta_B\), regardless of how many variants were combined to form \(B_i\).
In practice \(\sigma^2\) is usually unknown. For the likelihood score test, estimate it under the null by the maximum-likelihood estimate \(\hat{\sigma}_0^2=RSS_0/n\), where \(RSS_0=\sum_i(Y_i-\bar{Y})^2\). Substituting this null estimate into the score and its variance gives
Here \(r_{BY}^2\) is the squared sample correlation between burden and phenotype, when both vary. The likelihood score statistic is compared with a \(\chi^2_1\) distribution as a large-sample approximation under \(H_0\). Using the null residual variance estimated with divisor \(n-1\), rather than the null maximum-likelihood estimate with divisor \(n\), changes the finite-sample statistic; the variance convention should therefore be stated rather than silently mixed.
Worked Example: Constructing a Burden and Testing It
Suppose five hypothetical individuals have dosages at two variants shown below. The prespecified weights are \(w_1=2\) and \(w_2=1\), and the quantitative phenotype is measured in arbitrary trait units. Calculate the estimated-variance score statistic.
| Individual | \(G_{i1}\) | \(G_{i2}\) | \(B_i=2G_{i1}+G_{i2}\) | \(Y_i\) |
|---|---|---|---|---|
| 1 | 0 | 0 | 0 | 2 |
| 2 | 0 | 0 | 0 | 3 |
| 3 | 0 | 1 | 1 | 3 |
| 4 | 0 | 1 | 1 | 5 |
| 5 | 1 | 0 | 2 | 7 |
The burden mean is \(\bar{B}=(0+0+1+1+2)/5=0.8\), and the phenotype mean is \(\bar{Y}=(2+3+3+5+7)/5=4\). The centered products give
The centered burden squares sum to \(S_{BB}=0.64+0.64+0.04+0.04+1.44=2.8\). The null residual sum of squares is \(RSS_0=(-2)^2+(-1)^2+(-1)^2+1^2+3^2=16\), so \(\hat{\sigma}_0^2=16/5=3.2\) trait-units squared. Therefore
Against the large-sample \(\chi^2_1\) reference distribution, the two-sided p-value is approximately \(0.045\). Under the model and the null, this is the probability of observing a score statistic at least this large; it is not the probability that the null hypothesis is true. The positive \(S_{BY}\) indicates that larger burdens tend to accompany larger phenotype values in this sample, but the test alone does not establish a causal effect.
Known-Variance Calculation and Weight Scale
The known-variance version makes the null variance explicit. It can also clarify that the burden’s units affect the score and its standard error in matching ways.
Worked Example: Using a Known Error Variance
In a hypothetical sample of four individuals, burdens are \(0,0,1,1\), and phenotypes are \(1,2,2,5\) trait units. Suppose the error standard deviation \(\sigma\) is known to be \(1.5\) trait units. The means are \(\bar B=0.5\) and \(\bar Y=2.5\). Thus
The score standard deviation under the null is \(\sqrt{\sigma^2S_{BB}}=\sqrt{(1.5)^2(1)}=1.5\) trait units per burden-scale unit. Hence \(Z=2/1.5=1.333\), and \(Z^2=1.778\). The exact normal-theory two-sided p-value is about \(0.182\), equivalently the upper-tail probability for \(\chi^2_1\) at \(1.778\). This calculation relies on the stipulated known variance; it is not the usual estimated-variance procedure.
Worked Example: Rescaling All Weights
Consider four hypothetical individuals with burden values \(B=(0,1,2,1)\) and phenotype values \(Y=(1,1,4,2)\). The means are \(\bar B=1\) and \(\bar Y=2\), giving \(S_{BB}=2\), \(S_{BY}=3\), and \(RSS_0=6\). The score statistic is
If every weight is doubled, every burden becomes \(B_i^*=2B_i\). Then \(S_{BB}^*=4S_{BB}=8\), \(S_{BY}^*=2S_{BY}=6\), and \(RSS_0\) is unchanged. The recalculated statistic is \(4(6^2)/[6(8)]=3.000\). In general, multiplying all weights by a nonzero constant \(c\) multiplies \(S_{BY}\) by \(c\) and \(S_{BB}\) by \(c^2\), so their squared ratio is unchanged. This invariance applies only to a common rescaling; changing relative weights changes the burden and can change the test result.
Interpretation, Scope, and Common Pitfalls
This test is a one-predictor association test, not a separate test for every variant in the gene. Its result depends on the chosen mask, dosage coding, and relative weights because those choices define \(B_i\). The score statistic has no direction because it is squared; inspect \(S_{BY}\) or the signed standardized score to determine whether higher burdens align with higher or lower phenotypes in the sample.
The large-sample \(\chi^2_1\) calibration should not be treated as automatically reliable in every rare-variant sample. If burden variation is extremely limited, the usual approximation may be poor. This derivation also does not address binary outcomes, additional covariates, or tests designed to accommodate variants with differing effect directions.
- Using uncentered products as if the intercept were known. The fitted null includes an intercept, so the score uses \(S_{BY}\), the centered burden-phenotype cross-product. An overall shift in all burden values cannot supply evidence for a slope.
- Counting variants as degrees of freedom. The test has one degree of freedom because it tests one burden coefficient. It does not test one coefficient per included variant.
- Using the wrong variance convention without saying so. The likelihood score formula above uses \(RSS_0/n\), the null maximum-likelihood estimate. Replacing it by \(RSS_0/(n-1)\) changes the numerical statistic.
- Calling a small p-value the probability that no association exists. A p-value is calculated under the null model; it is not a probability assigned to \(H_0\), and an association does not by itself imply causation.
Check Your Understanding
Use the model assumptions and formulas above to answer the following.
- Why does estimating an intercept lead to a centered burden sum \(S_{BB}\) in the null variance?
- If \(S_{BY}=0\) while \(S_{BB}>0\), what is the score statistic, and what does that say about the observed linear alignment?
- Why does combining 20 variants into one prespecified burden still produce a one-degree-of-freedom test?
- If all burden weights are multiplied by \(-3\), which quantities change, and does the estimated-variance score statistic change?
- State the assumptions used for the exact known-variance normal test, and distinguish them from the large-sample calibration with estimated variance.