Tutorials › High-Dimensional Genomic Data Analysis › Logistic Burden Tests for Binary Phenotypes

Rare-Variant Burden and Variance-Component Tests · Tutorial 328 of 795

Logistic Burden Tests for Binary Phenotypes

Build and interpret covariate-adjusted logistic burden score tests, including how allele weights define the tested burden.

Advanced 10 min read

What You'll Learn

  • Define a weighted variant burden from genotype dosages and prespecified allele weights.
  • Specify the logistic null model for a binary phenotype and identify the conditions for score-test calibration.
  • Derive the efficient burden score and its variance after nuisance-covariate adjustment.
  • Calculate a one-degree-of-freedom test statistic and p-value from hypothetical case-control data.
  • Recognize how weighting, confounding, sparse data, and model misspecification can affect conclusions.

Why Use a Logistic Burden Test?

For a binary phenotype, such as case status, a rare-variant burden test asks whether a weighted summary of variants in a defined set is associated with the outcome. It is useful when the scientific alternative is that variants in the set tend to affect risk in a common direction. Instead of estimating a separate effect for every rare variant, the test evaluates one coefficient for the combined burden.

The outcome requires a logistic model rather than a linear model for a quantitative response. Covariate adjustment is also essential when factors such as ancestry, age, or study site are related to both genotype burden and case status. The score-test construction follows the earlier material on covariate-adjusted score tests: nuisance covariate effects are fitted under the null, and the burden contribution is evaluated after accounting for those covariates.

Definition: Let \(G_{ij}\) be the genotype dosage for individual \(i\) at variant \(j\), coded as the number of copies of the designated allele, and let \(w_j\) be a prespecified weight for variant \(j\). The weighted burden for individual \(i\) is \(B_i=\sum_{j=1}^{m}w_jG_{ij}\). In a logistic burden model, \(\operatorname{logit}\{P(Y_i=1\mid G_i,C_i)\}=C_i^\mathsf{T}\gamma+\beta B_i\), where \(Y_i\) is the binary outcome, \(C_i\) is the vector of covariates (including an intercept), \(\gamma\) contains nuisance coefficients, and \(\beta\) is the burden coefficient tested by \(H_0:\beta=0\).

A positive \(\beta\) means that greater weighted burden is associated with higher odds of being a case, conditional on the covariates; a negative \(\beta\) means lower odds. The odds ratio for a one-unit increase in burden is \(\exp(\beta)\). A burden unit depends on the chosen weights and dosage coding, so it is not automatically a one-allele change at a particular variant.

Model and Score-Test Construction

Under the null hypothesis, fit the logistic regression model without the burden term: \(\operatorname{logit}(\mu_i)=C_i^\mathsf{T}\gamma\), where \(\mu_i=P(Y_i=1\mid C_i)\). Write \(\hat{\mu}_i\) for the fitted null probability and \(W\) for the diagonal matrix with entries \(\hat{\mu}_i(1-\hat{\mu}_i)\). The burden score before nuisance adjustment is \(B^\mathsf{T}(Y-\hat{\mu})\), where \(Y\) is the vector of outcomes.

The fitted nuisance parameters affect the variability of this score. The appropriate variance is therefore not simply the variance of the burden in the sample. It is the information for the burden after projecting out the covariates using the logistic information weights \(W\).

$$ a=(C^\mathsf{T}WC)^{-1}C^\mathsf{T}WB,\qquad R=B-Ca,\qquad U=B^\mathsf{T}(Y-\hat{\mu}),\qquad V=R^\mathsf{T}WR. $$

Here \(R\) is the burden residual after its weighted linear projection onto the covariates. The inverse requires \(C^\mathsf{T}WC\) to be nonsingular. The score can also be written \(U=R^\mathsf{T}(Y-\hat{\mu})\), because the null-model score equations give \(C^\mathsf{T}(Y-\hat{\mu})=0\). This representation makes clear that only the portion of burden not explained by the covariates contributes to the test.

Theorem: Under the null model, suppose observations are independent, the logistic mean model is correctly specified, the null maximum-likelihood fit exists and satisfies the usual regularity conditions, and a central limit theorem applies to the efficient score. In particular, the efficient information must grow (\(V\to\infty\)); a sufficient condition that no observation dominates the bounded Bernoulli score contributions is \(\max_i |R_i|/\sqrt{V}\to0\). Then the standardized score \(U/\sqrt{V}\) is approximately standard normal in large samples. The two-sided score statistic \(Q=U^2/V\) is therefore approximately distributed as a chi-square random variable with one degree of freedom. The approximation is asymptotic; it can be poor when the data are very sparse or the fitted probabilities are extreme.

Why the Covariate Adjustment Works

The score for \(\beta\) at the null is \(U_\beta=B^\mathsf{T}(Y-\hat{\mu})\). In the full model, the information matrix has a burden block, a covariate block, and cross-information between the burden and covariates. Removing the part of the burden score predictable from the nuisance scores gives the efficient score variance \(\,I_{\beta\beta}-I_{\beta\gamma}I_{\gamma\gamma}^{-1}I_{\gamma\beta}\,\). For logistic regression, these information blocks are \(B^\mathsf{T}WB\), \(C^\mathsf{T}WC\), and \(B^\mathsf{T}WC\). Substitution gives \(B^\mathsf{T}WB-B^\mathsf{T}WC(C^\mathsf{T}WC)^{-1}C^\mathsf{T}WB\), which equals \(R^\mathsf{T}WR\).

The null score equations also justify the adjusted score representation: since \(C^\mathsf{T}(Y-\hat{\mu})=0\), \((B-Ca)^\mathsf{T}(Y-\hat{\mu})=B^\mathsf{T}(Y-\hat{\mu})\). Thus nuisance adjustment changes the variance through residualization, while the fitted-null score has the same numerical value whether written with \(B\) or \(R\).

Worked Example: An Intercept-Only Case-Control Test

Suppose eight hypothetical, unrelated participants are sampled, with the first four cases and the last four controls. Their already-defined burden scores are \(0,0,1,1,0,0,1,2\). Assume an intercept-only null model and use a two-sided test.

The null estimate of the case probability in this sample is \(\hat{\mu}=4/8=0.5\). The case burdens sum to \(0+0+1+1=2\), and the total burden is \(0+0+1+1+0+0+1+2=5\). Therefore \(U=\sum_i(Y_i-\hat{\mu})B_i=2-0.5(5)=-0.5\).

The mean burden is \(5/8=0.625\). With an intercept as the only covariate, the residuals are the burdens centered at this mean. Their squared sum is \(4(0-0.625)^2+3(1-0.625)^2+(2-0.625)^2=3.875\). Since \(\hat{\mu}(1-\hat{\mu})=0.25\), the efficient variance is \(V=0.25(3.875)=0.96875\). Thus \(Q=(-0.5)^2/0.96875\approx0.2581\), with one degree of freedom, giving a two-sided asymptotic p-value of approximately \(0.6115\).

This result is not evidence against the null at conventional significance levels. The negative score indicates that the observed burden-case association points in the lower-odds direction under this coding, but its magnitude relative to its standard error is small. The p-value is the probability, under the null and the stated approximation, of a score statistic at least this extreme; it is not the probability that the null hypothesis is true.

Worked Example: Adjusting for a Binary Covariate

Consider eight hypothetical participants in two covariate strata, A and B. Each stratum has two cases and two controls, so the null logistic fit with an intercept and stratum indicator has \(\hat{\mu}_i=0.5\) for everyone. In stratum A, burdens are \(0,1,0,1\), with the first two participants cases. In stratum B, burdens are \(0,0,1,2\), with the first and third participants cases.

In stratum A, case burden is \(0+1=1\), while half the total stratum burden is \(0.5(2)=1\); its score contribution is zero. In stratum B, case burden is \(0+1=1\), while half the total burden is \(0.5(3)=1.5\); its contribution is \(-0.5\). Hence \(U=-0.5\).

The covariate adjustment centers burdens within each stratum because the fitted probabilities are equal. The stratum A centered squared sum is \(4(0.5)^2=1\). The stratum B mean burden is \(3/4=0.75\), giving centered squared sum \((0-0.75)^2+(0-0.75)^2+(1-0.75)^2+(2-0.75)^2=2.75\). Thus \(V=0.25(1+2.75)=0.9375\), and \(Q=0.25/0.9375\approx0.2667\). The one-degree-of-freedom two-sided p-value is approximately \(0.606\). The adjustment tests the within-stratum burden pattern rather than treating differences in mean burden between strata as association evidence.

Weights and Their Consequences

Weights define the burden being tested; they do not merely rescale a universal test. For example, weights may give greater contribution to variants considered more likely to affect function, or may downweight variants according to a prespecified frequency rule. The weighting scheme should be chosen before examining the association results, or the additional selection must be reflected in the inferential procedure. If weights are changed, then \(B\), its covariate projection, \(U\), and \(V\) all change.

A burden test is most aligned with alternatives in which weighted variant effects contribute in a broadly common direction. Opposing effects can cancel in the aggregate. Separate-variant logistic regressions avoid this particular aggregation assumption but may be imprecise for rare alleles and require attention to multiple testing. An unweighted burden is a simpler alternative, but it gives each dosage unit the same contribution. These approaches answer different questions; a significant burden test supports association of the specified weighted summary, not proof that every included variant is associated.

Worked Example: A Prespecified Weighted Burden

For six hypothetical participants, suppose two variant dosages are \(G_1=(0,1,0,0,0,1)\) and \(G_2=(0,0,1,0,1,0)\). Prespecify weights \(w_1=2\) and \(w_2=1\), so \(B_i=2G_{i1}+G_{i2}\), giving burdens \((0,2,1,0,1,2)\). Let the outcomes be \(Y=(1,1,0,0,0,1)\), with three cases. Use an intercept-only null model.

The fitted null probability is \(3/6=0.5\). Total burden is 6 and case burden is \(0+2+2=4\), so \(U=4-0.5(6)=1\). The mean burden is 1, and the centered squared sum is \((0-1)^2+(2-1)^2+(1-1)^2+(0-1)^2+(1-1)^2+(2-1)^2=4\). Therefore \(V=0.25(4)=1\), \(Q=1^2/1=1\), and the one-degree-of-freedom p-value is approximately \(0.3173\).

The test does not provide strong evidence of association for this specified weighted burden. Its unit is one unit of \(2G_1+G_2\), not one allele at either variant in isolation. The changed weight on variant 1 changes the hypothesis being tested as well as the numerical scale of the burden.

Assumptions, Diagnostics, and Common Pitfalls

The calculation relies on a suitable logistic model for the null mean, independent observations under the stated model, and a valid covariate design. In a case-control study sampled according to outcome, logistic slope associations can be estimated from the sampled data when selection depends on outcome and does not further depend on genotype or covariates within outcome groups. The fitted intercept then reflects the sampled case fraction rather than population prevalence. If selection depends additionally on covariates or genotype, the sampling scheme may require a more appropriate model or weighting.

In genomic studies, inspect the burden and covariates for near-collinearity, check that the null model fits and its information matrix is estimable, and verify the genotype and weight coding. A burden almost completely explained by covariates leaves \(V\) near zero and little information for testing. Extreme fitted probabilities, very few cases or controls in covariate strata, and highly sparse carrier patterns can make the large-sample chi-square approximation unreliable. Population structure, relatedness, or other dependence not represented in the model can also invalidate the nominal calibration. The asymptotic score test does not repair these design or model problems.

  • Mistake: using the unadjusted burden variance after fitting covariates. This ignores the covariate-burden overlap. Use \(V=R^\mathsf{T}WR\), with projection based on the fitted null model.
  • Mistake: interpreting a weighted-burden result as evidence for every variant. The test concerns one prespecified weighted aggregate and can be influenced by multiple variants together.
  • Mistake: treating the sample case fraction as population prevalence in an ascertained study. Outcome-based sampling changes the intercept. It does not, by itself, make the sampled case fraction a population risk estimate.
  • Mistake: reading a nonsignificant result as evidence of no genetic effect. A burden test can have limited power when effects oppose one another, weights are poorly aligned with effects, or the effective information is small.
Key takeaway: A logistic burden score test evaluates one prespecified weighted genotype summary in a null logistic model. Its score is the null residual-outcome association, and its variance must account for the covariates through weighted nuisance projection.

Check Your Understanding

Work through these questions, then compare your reasoning with the solutions below.

  1. In a logistic burden model, what does the coefficient \(\beta\) describe, and what is the odds ratio for a one-unit burden increase?
  2. Why is the efficient variance based on a covariate-residualized burden rather than the raw burden?
  3. Suppose a test has \(U=2\) and \(V=4\). Calculate \(Q\) and identify its approximate null distribution.
  4. Why does a case-control sample's fitted intercept generally not estimate population prevalence?
  5. Give one reason a weighted burden test may fail to detect an association even if some variants have effects.

Solutions. (1) \(\beta\) is the conditional change in log odds per one-unit increase in the specified burden; the corresponding odds ratio is \(\exp(\beta)\). (2) Covariates can explain part of the burden, so that part supplies no independent association information; projection removes it and yields \(V=R^\mathsf{T}WR\). (3) \(Q=U^2/V=2^2/4=1\), approximately chi-square with one degree of freedom under the null. (4) Outcome-based sampling sets the case fraction in the analyzed sample, so its intercept reflects that sampling fraction rather than population prevalence. (5) Effects with opposite directions can cancel in the aggregate; alternatively, weights may be poorly aligned with effect sizes or the data may provide little information.

Clinical Biostats