Logistic Score Tests For Binary-Trait Burdens

Why Use a Logistic Score Test for a Burden?

A burden combines selected variant dosages into one predictor, with weights reflecting the analysis choices for the variants. When the trait is binary, the association question is whether that predictor is related to the probability of being a case, after accounting for covariates such as age, sex, or ancestry variables. Logistic regression provides a model for that probability while keeping the test focused on a single burden coefficient.

A score test evaluates the burden coefficient at its null value, zero, after fitting the model that includes the covariates but not the burden. Thus it does not require fitting the alternative model to calculate the test statistic. Covariate adjustment is essential: the part of burden variation predictable from the covariates should not count as independent evidence about the trait.

Definition: For individual $i$, let $Y_i$ equal 1 for a case and 0 for a control. Let $B_i$ be that individual’s fixed, weighted burden, and let $C_i$ be a vector of covariates that includes an intercept. The covariate-adjusted logistic model is $$ P(Y_i=1\mid B_i,C_i)=\pi_i,\qquad \operatorname{logit}(\pi_i)=C_i^\mathsf{T}\gamma+\beta B_i, $$ where $\operatorname{logit}(\pi)=\log\{\pi/(1-\pi)\}$, $\gamma$ is the covariate coefficient vector, and $\beta$ is the burden coefficient. The null hypothesis is $H_0:\beta=0$; a two-sided association test uses $H_a:\beta\ne 0$.

The burden’s variant set and weights are treated as fixed for this test. For example, $B_i$ may be a sum of allele dosages multiplied by predetermined weights. The units of $B_i$ depend on that construction; the coefficient is consequently interpreted per one unit of this particular burden.

Deriving the Adjusted Score Statistic

Under $H_0$, fit the logistic model using $C_i$ alone. Denote its fitted case probability by $\hat{\pi}_i$, and set $\hat{w}_i=\hat{\pi}_i(1-\hat{\pi}_i)$. These fitted quantities are evaluated under the null model, not under a model that includes the burden.

At the null fit, the score for $\beta$, before accounting for nuisance parameters, is the sum of burden times the outcome residual. The information matrix contains both the burden and covariate components. Removing the information attributable to estimating $\gamma$ gives the efficient score and its variance.

U_\beta=\sum_{i=1}^n B_i(Y_i-\hat{\pi}_i),\qquad I_{\beta\beta}=\sum_{i=1}^n \hat{w}_i B_i^2,\qquad I_{\beta\gamma}=\sum_{i=1}^n \hat{w}_i B_iC_i^\mathsf{T},\qquad I_{\gamma\gamma}=\sum_{i=1}^n \hat{w}_i C_iC_i^\mathsf{T}.

Theorem: If $I_{\gamma\gamma}$ is invertible, the score and information for the burden after adjusting for the nuisance covariates are $$ U_{\mathrm{adj}}=U_\beta-I_{\beta\gamma}I_{\gamma\gamma}^{-1}U_\gamma,\qquad I_{\mathrm{adj}}=I_{\beta\beta}-I_{\beta\gamma}I_{\gamma\gamma}^{-1}I_{\gamma\beta}. $$ At the maximum-likelihood fit of the null model, $U_\gamma=0$, so $U_{\mathrm{adj}}=U_\beta$. The one-degree-of-freedom score statistic is $T=U_{\mathrm{adj}}^2/I_{\mathrm{adj}}$. Under $H_0$ and the usual large-sample conditions, $T$ is approximately distributed as chi-square with 1 degree of freedom.

The expression for $I_{\mathrm{adj}}$ also shows directly how adjustment works. Define the weighted residual of burden after projection on the covariates as $$ \widetilde B_i=B_i-C_i^\mathsf{T}(C^\mathsf{T}WC)^{-1}C^\mathsf{T}WB, $$ where $W$ is the diagonal matrix with entries $\hat w_i$. Then $$ I_{\mathrm{adj}}=\sum_{i=1}^n \hat w_i\widetilde B_i^2. $$ Because the null-model score equations give $C^\mathsf{T}(Y-\hat\pi)=0$, the adjusted score can also be written as $$ U_{\mathrm{adj}}=\sum_{i=1}^n \widetilde B_i(Y_i-\hat\pi_i). $$ So the test uses the burden variation not explained by the covariates, paired with the outcome residuals from the null logistic model. The weighted projection is important: when fitted probabilities differ between individuals, their contributions to information are not equal.

Worked Score-Test Calculations

Worked Example: Testing a Burden With a Covariate-Adjusted Null Fit

Suppose a hypothetical binary-trait analysis has already fitted its prespecified covariate-only logistic model. The burden score from this null fit is $U_{\mathrm{adj}}=2.50$, and the corresponding covariate-adjusted information is $I_{\mathrm{adj}}=1.21875$. Calculate and interpret the score test. These are hypothetical null-fit summaries, stated to show the calculation.

The statistic is $$ T=\frac{U_{\mathrm{adj}}^2}{I_{\mathrm{adj}}} =\frac{(2.50)^2}{1.21875} =\frac{6.25}{1.21875} =5.1282. $$

For a chi-square distribution with 1 degree of freedom, the nominal two-sided asymptotic p-value is $P(\chi^2_1\ge 5.1282)=0.02354$. At significance level $\alpha=0.05$, this is evidence against $H_0:\beta=0$ in favor of an association between this weighted burden and the binary trait, conditional on the covariates in the null model. It is not evidence that the burden causes the trait, nor is the p-value the probability that $H_0$ is true.

Worked Example: Seeing Covariate Adjustment Within Strata

Consider a hypothetical study with eight participants in two covariate strata, A and B. In each stratum there are two cases and two controls, so the covariate-only logistic fit assigns $\hat\pi_i=0.5$ to all participants. In stratum A, the cases have burdens 2 and 2 and the controls have burdens 0 and 0. In stratum B, cases have burdens 1 and 1 and controls have burdens 0 and 0. The covariate is the stratum indicator, and the model includes an intercept.

Within stratum A, the weighted-mean burden is 1, giving residualized burdens $1,1,-1,-1$ for the two cases and two controls. Within stratum B, the weighted-mean burden is 0.5, giving residualized burdens $0.5,0.5,-0.5,-0.5$. Since all $\hat w_i=0.5(1-0.5)=0.25$, the adjusted score contributions from stratum A sum to $2$: each case contributes $1(1-0.5)=0.5$, and each control contributes $(-1)(0-0.5)=0.5$. In stratum B the corresponding sum is $1$. Thus $U_{\mathrm{adj}}=3$.

The information is the weighted sum of squared residualized burdens: $$ I_{\mathrm{adj}}=0.25\{4(1^2)+4(0.5^2)\} =0.25(4+1)=1.25. $$ Therefore $T=3^2/1.25=7.20$, with a one-degree-of-freedom asymptotic p-value of about 0.0073. This calculation compares burden differences within the covariate strata. Simply treating the raw burdens as though the stratum covariate were absent would not give the adjusted score test.

What the Burden Coefficient Measures

In the stated model, $\beta$ is the change in log odds of being a case associated with a one-unit increase in burden, holding the covariates fixed. For a burden change of $\Delta$, the model’s odds ratio is $\exp(\beta\Delta)$. This interpretation is conditional on the logistic model and on the defined burden scale; it does not describe the effect of changing one specific variant while holding all other components of the burden fixed.

The score and information also give a local, null-based estimate of the coefficient: $$ \hat\beta_{\mathrm{score}}\approx\frac{U_{\mathrm{adj}}}{I_{\mathrm{adj}}}. $$ This follows from a first-order expansion of the score around $\beta=0$: near zero, the score changes approximately as $U_{\mathrm{adj}}(\beta)\approx U_{\mathrm{adj}}(0)-I_{\mathrm{adj}}\beta$. Setting this approximation to zero gives the displayed estimate. It is a score-based approximation, not the maximum-likelihood estimate obtained by fitting the alternative model.

Worked Example: Coefficient Scale and Odds Ratio

Use the hypothetical score and information in the first example. The local score-based coefficient is $$ \hat\beta_{\mathrm{score}}\approx\frac{2.50}{1.21875}=2.051. $$ If burden is measured in the units used for that analysis, this corresponds locally to an increase of about 2.051 in log odds per burden unit. Exponentiating gives $\exp(2.051)\approx7.78$, a local odds-ratio approximation per one burden unit. It should not be reported as the fitted alternative-model odds ratio without fitting that model.

Now define the same burden in doubled units, $B_i'=2B_i$. Its score doubles to $U_{\mathrm{adj}}'=5.00$, and its information becomes $I_{\mathrm{adj}}'=4(1.21875)=4.875$. The score statistic is unchanged: $$ T'=\frac{(5.00)^2}{4.875}=5.1282. $$ The local coefficient on the doubled scale is $5.00/4.875=1.026$, half the coefficient on the original scale, as required because one doubled unit represents twice as much burden. A coefficient’s numerical value therefore depends on how burden units are defined; the test of a zero coefficient does not change under a nonzero rescaling.

Assumptions and Common Pitfalls

The chi-square calibration is an asymptotic result, not an automatic guarantee for every rare-variant data set. The model assumes independent observations conditional on the covariates, a Bernoulli outcome with a correctly specified null logistic mean, and a burden and covariate specification chosen independently of the observed outcomes for the test at hand. The null fit must be estimable, and the adjusted information must be positive. The approximation also relies on enough information for the score statistic to behave approximately as its large-sample reference distribution. When very few individuals carry informative variants, that approximation may be poor.

The burden test addresses the single predictor $B_i$. Its result depends on the variants included and their weights, and it does not identify which component of the burden is associated. The sign of the score indicates the direction of the score evidence: a positive score points toward larger case probability as burden increases in the specified model. A two-sided score statistic squares the score, so the test statistic itself does not retain that direction.

Common mistake: “The burden score is $\sum_i B_i(Y_i-\hat\pi_i)$, so I can use $\sum_i\hat w_iB_i^2$ as its variance even when covariates predict burden.” This omits the information shared between burden and covariates. The correction is to use $I_{\mathrm{adj}}=I_{\beta\beta}-I_{\beta\gamma}I_{\gamma\gamma}^{-1}I_{\gamma\beta}$, or equivalently the weighted sum of squared burden residuals after projecting on the covariates. Another mistake is to interpret $\beta$ as a log odds change per allele when the predictor is a weighted burden: its unit is one unit of the defined burden, not necessarily one allele.

Specify the burden model.
Define the weighted burden, the binary outcome, and the covariates, including an intercept.

Fit the null model.
Fit the covariate-only logistic model and obtain fitted probabilities and weights.

Calculate the adjusted score and information.
Use the burden score and remove the component of information explained by the covariates.

Test and interpret.
Calculate $T=U_{\mathrm{adj}}^2/I_{\mathrm{adj}}$, compare with the one-degree-of-freedom chi-square reference, and interpret any coefficient on the specified burden scale.

Key takeaway: A covariate-adjusted logistic burden score tests whether burden variation beyond that explained by covariates is associated with the binary trait. The score statistic is the squared adjusted score divided by adjusted information, and the model coefficient represents log odds change per burden unit.

Check Your Understanding

Answer each question using the model and notation above.

Why is the null logistic model fitted without the burden before calculating the score test?
What does $I_{\mathrm{adj}}$ subtract from the unadjusted burden information?
If the adjusted score is $-2$ and adjusted information is 1, what is the score statistic, and what direction does the score suggest?
In a logistic burden model, what does $\exp(\beta)$ represent when the burden is measured in its defined units?
Why can changing the numerical scale of a burden change its coefficient but leave the zero-coefficient score test unchanged?