Deriving the Burden Score Test

Why a Burden Score Has One Degree of Freedom

A rare-variant set can contain many genotype variables, but a burden test replaces them with one prespecified summary per individual. That summary allows a single regression coefficient to represent the association being tested. The score test asks whether the observed phenotype values align with that burden more strongly than expected under a model with no burden effect.

The derivation here uses an intercept-only Gaussian regression apart from the burden. This keeps the source of the score and its null variance explicit. The model does not adjust for additional covariates; the later subject of covariate-adjusted tests requires accounting for those covariates in the null residuals and information.

Burden and Regression Model

Definition: For individual $i$, let $G_{ij}$ be the genotype dosage at variant $j$, for $j=1,\ldots,m$. Given fixed weights $w_j$, define the weighted burden $$ B_i=\sum_{j=1}^{m}w_jG_{ij}. $$ The weights and genotype coding are treated as fixed for the test. A burden test in this tutorial tests the coefficient $\beta$ in $Y_i=\mu+\beta B_i+\varepsilon_i$.

The dosage $G_{ij}$ may, for example, count alternate alleles as 0, 1, or 2; a carrier indicator is another possible coding. The choice of weights and coding determines what the burden represents. Once chosen, the test treats $B_i$ as a single predictor. In this model, $\beta$ is the mean change in phenotype per one-unit increase in burden, measured in phenotype units per burden unit.

Assume the individuals are independent, the burden values are fixed conditional on the observed genotypes, and $Y_i=\mu+\beta B_i+\varepsilon_i$, where the errors are independent with mean zero and common variance $\sigma^2>0$. For the exact finite-sample reference distribution described below, additionally assume $\varepsilon_i\sim N(0,\sigma^2)$ and that $\sigma^2$ is known. The null hypothesis is $H_0:\beta=0$. The alternative is $H_a:\beta\ne 0$. A burden test also requires variation in the burden values: if all $B_i$ are identical, the data cannot distinguish a burden effect from the intercept.

Deriving the Score and Its Null Variance

Under the Gaussian model with known $\sigma^2$, the log likelihood, up to a constant, is

\ell(\mu,\beta)=-\frac{1}{2\sigma^2}\sum_{i=1}^{n}(Y_i-\mu-\beta B_i)^2.

At $\beta=0$, the maximum-likelihood estimate of the intercept is $\hat{\mu}_0=\bar{Y}$. Differentiating the log likelihood with respect to $\beta$, then evaluating at $\beta=0,\mu=\bar{Y}$, gives the score for $\beta$ under the null:

U_\beta =\left.\frac{\partial\ell}{\partial\beta}\right|_{\beta=0,\,\mu=\bar{Y}} =\frac{1}{\sigma^2}\sum_{i=1}^{n}B_i(Y_i-\bar{Y}) =\frac{1}{\sigma^2}\sum_{i=1}^{n}(B_i-\bar{B})(Y_i-\bar{Y}).

The last equality holds because $\sum_i(Y_i-\bar{Y})=0$. It also shows why the intercept matters: the score is based on centered burdens and centered phenotypes, not on their uncentered product.

To find the null variance, condition on the observed burdens. Under $H_0$, $Y_i=\mu+\varepsilon_i$. Thus the centered phenotype residuals have covariance $\sigma^2(I-\mathbf{1}\mathbf{1}^{T}/n)$, where $I$ is the $n$-by-$n$ identity matrix and $\mathbf{1}$ is the vector of ones. Consequently, the variance of the centered cross-product is $\sigma^2\sum_i(B_i-\bar{B})^2$. Dividing by $\sigma^4$, as in the score, gives the null variance below. Equivalently, this is the information for $\beta$ after accounting for the intercept.

Result: Define $S_{BB}=\sum_i(B_i-\bar{B})^2$. Under $H_0$, the score has mean zero and null variance (the information for $\beta$) $$ I_\beta=\operatorname{Var}_0(U_\beta)=\frac{S_{BB}}{\sigma^2}. $$ The standardized score is $Z=U_\beta/\sqrt{I_\beta}$, and the two-sided score statistic is $Q=U_\beta^2/I_\beta$.

Substituting the expressions for score and information yields a useful computational form:

Q=\frac{\left[\sum_i(B_i-\bar{B})(Y_i-\bar{Y})\right]^2} {\sigma^2\sum_i(B_i-\bar{B})^2}.

With fixed burdens and Gaussian errors of known variance, $Z\sim N(0,1)$ exactly under $H_0$, so $Q\sim\chi^2_1$ exactly. More generally, when $\sigma^2$ is unknown and replaced by a consistent estimate, the standardized score is asymptotically standard normal under regularity conditions, and $Q$ is asymptotically $\chi^2_1$. For a level-$\alpha$ test, reject $H_0$ when $Q$ exceeds the upper-tail $\chi^2_1$ critical value, or equivalently when its calibrated p-value is below $\alpha$.

Worked Examples

Worked Example: Building a Weighted Burden and Testing It

Consider six hypothetical individuals with two genotype dosages. Use weights $w_1=2$ and $w_2=1$, and suppose the phenotype values and genotypes are as follows. Assume the Gaussian model and known error variance $\sigma^2=2.5$.

Individual	$G_{i1}$	$G_{i2}$	$B_i=2G_{i1}+G_{i2}$	$Y_i$
1	0	0	0	1
2	0	1	1	2
3	1	0	2	2
4	0	1	1	1
5	1	1	3	4
6	0	1	1	2

The burden mean is $\bar{B}=8/6=1.3333$, and the phenotype mean is $\bar{Y}=12/6=2$. The centered cross-product can also be calculated as $\sum_iB_iY_i-n\bar{B}\bar{Y}=21-6(1.3333)(2)=5.000$. The burden sum of squares is $S_{BB}=\sum_iB_i^2-n\bar{B}^2=16-64/6=5.3333$.

Therefore $U_\beta=5.000/2.5=2.000$, and $I_\beta=5.3333/2.5=2.1333$. The score statistic is $Q=2.000^2/2.1333=1.875$, with a $\chi^2_1$ upper-tail p-value of approximately 0.171. At $\alpha=0.05$, this test does not reject $H_0$. The result is not evidence that the burden has no association; it indicates that this test does not provide sufficient evidence against a zero coefficient at that threshold.

Worked Example: Checking the Centered Cross-Product

For a second hypothetical data set, let $B=(0,1,1,3,2)$ and $Y=(1,2,1,4,3)$, with known Gaussian error variance $\sigma^2=1.5$. The means are $\bar{B}=1.4$ and $\bar{Y}=2.2$. Compute each centered product to avoid letting a single arithmetic error alter the test.

The five products are $(-1.4)(-1.2)=1.68$, $(-0.4)(-0.2)=0.08$, $(-0.4)(-1.2)=0.48$, $(1.6)(1.8)=2.88$, and $(0.6)(0.8)=0.48$. Their sum is $5.60$. Also, $S_{BB}=1.96+0.16+0.16+2.56+0.36=5.20$.

It follows that $U_\beta=5.60/1.5=3.7333$, $I_\beta=5.20/1.5=3.4667$, and $Q=3.7333^2/3.4667=4.0205$. Equivalently, $Q=5.60^2/(1.5\times5.20)=4.0205$. The $\chi^2_1$ upper-tail p-value is approximately 0.0449. Under the stated model and known-variance calibration, reject $H_0$ at the 0.05 level. The arithmetic and the decision both depend on using 2.88 for the fourth centered product.

Worked Example: Rescaling the Burden

Use the data in the preceding example, but double every burden value, as would happen if all the weights were multiplied by 2. Then the new burden has mean $2.8$, its centered cross-product with $Y$ is $2(5.60)=11.20$, and its centered sum of squares is $2^2(5.20)=20.80$.

With the same $\sigma^2=1.5$, the new score is $U_\beta=11.20/1.5=7.4667$, and the new information is $I_\beta=20.80/1.5=13.8667$. Thus $Q=7.4667^2/13.8667=4.0205$, the same statistic as before. Rescaling changes the coefficient’s units and the magnitudes of the score and information, but not this one-degree-of-freedom test. This invariance does not mean that changing the relative weights leaves the test unchanged.

Assumptions, Diagnostics, and Alternatives

The test’s interpretation depends on both the model and the burden construction. Individuals must be independent under the model, the error variance must be positive, and the burdens must have nonzero variation. The exact $\chi^2_1$ calibration additionally requires Gaussian errors and known $\sigma^2$. If the variance is estimated, its estimate must be suitable and consistent for the asymptotic calibration; a small sample alone does not guarantee a reliable approximation. For non-Gaussian outcomes, this derivation does not establish a valid reference distribution.

Before trusting a result, check that genotype coding and weights produce the intended burden, that no variant or individual has been inadvertently misaligned, and that $S_{BB}>0$. Inspect the phenotype-versus-burden relationship and residuals for strong departures from the assumed mean and constant-variance model. A few individuals with large burdens or extreme phenotypes can dominate the cross-product. Such features warn that the model or calibration may be inadequate; they should not be treated as proof of an association or as a reason to silently alter the analysis.

The burden test’s main alternative at this level is testing variant effects separately or fitting a multi-parameter regression with distinct coefficients for the variants. Those approaches do not collapse the genotype vector to one predictor, but they test a different parameter structure and require more parameters. The burden test is simpler and has one degree of freedom, but its conclusion concerns the chosen weighted aggregate; it does not identify which variant is associated or estimate each variant’s effect. If distinct variant contributions do not align with the aggregate, the single burden coefficient may not summarize them well. Choice among these approaches depends on the scientific target and model assumptions.

Common Pitfall

Common mistake: Use $\sum_i B_iY_i$ as the score numerator even though the model includes an intercept. Correction: Center both variables. Since $\sum_i(Y_i-\bar{Y})=0$, the correct numerator is $\sum_i(B_i-\bar{B})(Y_i-\bar{Y})$. Also verify each product and sum before calculating the statistic: one incorrect term changes both the score and the p-value.

Key takeaway: For the intercept-only Gaussian burden model, the score test compares the squared centered burden–phenotype cross-product with its null variance. The resulting statistic has one degree of freedom, with exact $\chi^2_1$ calibration for Gaussian errors and known variance, and asymptotic $\chi^2_1$ calibration under suitable large-sample conditions.

Check Your Understanding

Answer each question, then compare with the solutions.

If every individual has the same burden, can this model estimate a distinct burden effect and intercept? Why?
For fixed burdens, $S_{BB}=8$, centered cross-product $=4$, and known $\sigma^2=2$. Calculate $U_\beta$, $I_\beta$, and $Q$.
What null reference distribution applies exactly when errors are Gaussian and $\sigma^2$ is known?
If all burden weights are multiplied by a nonzero constant $c$, what happens to $Q$?
Does a significant burden test identify which variant in the set is associated?

Solutions. (1) No. With constant burden, the burden is indistinguishable from a change in the intercept, and $S_{BB}=0$. (2) $U_\beta=4/2=2$, $I_\beta=8/2=4$, and $Q=2^2/4=1$. (3) $Q\sim\chi^2_1$ exactly under the stated Gaussian known-variance assumptions. (4) The score is multiplied by $c$, information by $c^2$, and their ratio $U_\beta^2/I_\beta$ is unchanged. (5) No. It tests the coefficient of the aggregate burden, not separate variant coefficients.

Individual	\(G_{i1}\)	\(G_{i2}\)	\(B_i=2G_{i1}+G_{i2}\)	\(Y_i\)
1	0	0	0	1
2	0	1	1	2
3	1	0	2	2
4	0	1	1	1
5	1	1	3	4
6	0	1	1	2

What You'll Learn