Variable-Threshold Burden Testing

Why Search Across Frequency Cutoffs?

A burden test combines qualifying variants into a score and tests whether that score is associated with a phenotype. Its result can depend on which variants qualify. A stringent minor allele frequency (MAF) cutoff may exclude variants whose effects contribute useful signal; a more inclusive cutoff may add variants with weak or opposing effects. If the appropriate cutoff is not known in advance, testing several cutoffs and reporting the strongest result is tempting—but that search itself must be included in the inference.

A variable-threshold (VT) burden test addresses this issue by calculating a burden statistic at each of several MAF thresholds, selecting the most extreme statistic according to a prespecified rule, and calibrating that selected statistic under the null while repeating the same selection. The adjustment is essential: the best of several noisy results tends to look stronger than any one result chosen in advance.

Definition: A variable-threshold burden test evaluates a burden-based association statistic over a prespecified set of minor-allele-frequency thresholds and tests a summary, usually the maximum statistic, with a null distribution that accounts for searching across those thresholds.

Burden Scores at Nested Thresholds

Suppose a set contains $m$ variants and $G_{ij}$ is the genotype dosage for individual $i$ at variant $j$, coded as 0, 1, or 2 copies of the allele of interest. Let $q_j$ be the minor allele frequency of variant $j$, and let $w_j$ be a prespecified weight. For a threshold $t$, define the individual burden

B_i(t)=\sum_{j:q_j\leq t} w_jG_{ij}.

The thresholds are usually ordered, $t_1<t_2<\cdots<t_K$. Because the included variant sets are nested, increasing the threshold can add variants but does not remove those already included. For a simple quantitative-phenotype illustration, suppose the null model is $Y_i=\mu+\varepsilon_i$, with independent errors of mean zero and variance $\sigma^2$. Define $r_i=Y_i-\bar{x}$, where here $\bar{x}$ is not used for phenotype centering; to avoid confusing the course notation, write the centered phenotype explicitly as $Y_i-\bar{Y}$, with $\bar{Y}=n^{-1}\sum_iY_i$. A score and its null variance in this model are

U(t)=\sum_{i=1}^{n}(Y_i-\bar{Y})B_i(t), \qquad V(t)=\sigma^2\sum_{i=1}^{n}\{B_i(t)-\bar{B}(t)\}^2, \qquad Z(t)=\frac{U(t)}{\sqrt{V(t)}},

where $\bar{B}(t)=n^{-1}\sum_iB_i(t)$. The standardized score $Z(t)$ is approximately standard normal under the null when the score approximation is adequate. In applied settings, the score and variance must match the specified phenotype model; the simple formula above is for an intercept-only quantitative-phenotype model, not a universal formula for all designs.

A directional test for an increasing phenotype selects $T_{\max}=\max_{k=1,\ldots,K}Z(t_k)$. If the direction is not specified, a two-sided rule can instead use $\max_k|Z(t_k)|$. The direction and selection rule should be fixed before inspecting phenotype associations. Choosing a direction after seeing the results is another selection and must also be accounted for.

Worked Example: Building Burdens at Three Cutoffs

Consider four hypothetical individuals, three variants with MAFs 0.004, 0.009, and 0.018, and thresholds 0.005, 0.010, and 0.020. For illustration, take all weights as 1 and phenotype values $Y=(0,0,2,2)$, so $\bar{Y}=1$. The genotype dosages at the three variants are:

Individual	Variant 1	Variant 2	Variant 3	$Y_i-\bar{Y}$
1	0	0	1	-1
2	0	1	0	-1
3	1	0	0	1
4	0	0	0	1

At $t=0.005$, only variant 1 is included, giving burdens $B=(0,0,1,0)$. Thus $U= (-1)(0)+(-1)(0)+(1)(1)+(1)(0)=1$. At $t=0.010$, variants 1 and 2 are included, giving $B=(0,1,1,0)$ and $U=0-1+1+0=0$. At $t=0.020$, all three are included, giving $B=(1,1,1,0)$ and $U=-1-1+1+0=-1$. This shows that adding a variant can change a burden score in either direction; nested inclusion does not guarantee a larger association statistic.

For a scale illustration, set $\sigma^2=1$. At the first threshold, $\bar{B}=0.25$, so $V=0.25^2+0.25^2+0.75^2+0.25^2=0.75$, and $Z=1/\sqrt{0.75}=1.155$. At the second threshold, $\bar{B}=0.5$, $V=1$, and $Z=0$. At the third, $\bar{B}=0.75$, $V=0.75$, and $Z=-1/\sqrt{0.75}=-1.155$. These four observations are only a calculation demonstration; they are far too few to justify a normal approximation or a reliable association conclusion.

Why Threshold Selection Requires Calibration

Even if each $Z(t_k)$ has a standard normal null distribution, the maximum over thresholds does not. The thresholds are nested, so their statistics are correlated, often strongly. The null distribution of the maximum depends on that joint correlation, not only on the number of thresholds. Treating the selected statistic as if its threshold had been prespecified generally gives a p-value that is too small.

Null calibration: If $(Z(t_1),\ldots,Z(t_K))$ has an approximate multivariate normal null distribution with mean zero and correlation matrix $R$, then for a one-sided maximum $M=\max_k Z(t_k)$, the selection-adjusted p-value at observed value $m$ is $$ P(M\geq m\mid H_0)=1-P\{Z(t_1)<m,\ldots,Z(t_K)<m\mid H_0\}. $$ For a two-sided maximum, replace the event inside the probability by $\{|Z(t_k)|<m\text{ for every }k\}$.

This is a joint-tail probability. If the statistics were independent, the one-sided probability would be $1-\Phi(m)^K$, where $\Phi$ is the standard normal cumulative distribution function. In a VT test they are usually not independent, so this expression is generally not the correct calibration. A multivariate normal calculation can use an estimated or model-based $R$; permutation provides another route when its assumptions are appropriate.

Permutation Calibration

Under a valid permutation scheme, the phenotype labels (or suitable null residuals) are rearranged relative to the genotypes in a way that preserves the null distribution. For each rearrangement, calculate every threshold statistic and retain the maximum using exactly the same direction, thresholds, and selection rule as for the observed data. This produces a null distribution for the selected statistic rather than for an individual cutoff.

Simple phenotype permutation is justified when observations are exchangeable under the null—for example, in an appropriate independent, identically distributed setting with no structure that makes some labels non-interchangeable. If covariates, relatedness, population structure, or other design features are present, unrestricted permutation may not preserve the null model. The permutation method must respect those features; a generic shuffling instruction does not make an invalid permutation valid.

With $B$ random permutations and $b$ permuted maxima at least as large as the observed maximum, a commonly used Monte Carlo estimate is $(b+1)/(B+1)$. The added one prevents a zero reported p-value and gives a finite-simulation estimate. It does not make a small number of permutations precise, so the number of permutations should be large enough for the tail probability being estimated.

Worked Example: Adjusting a Selected Threshold by Permutation

Suppose a hypothetical analysis tests three prespecified MAF thresholds using a one-sided rule. The observed standardized statistics are $1.20$, $1.65$, and $2.121$, so the selected maximum is $m=2.121$. If this last statistic had been the only prespecified test, its standard normal upper-tail probability would be $P(Z\geq2.121)\approx0.01696$, for $Z\sim N(0,1)$.

Now suppose 10,000 valid phenotype permutations are performed. For each permutation, all three statistics are recomputed and their maximum is recorded. Assume 542 permuted maxima are at least 2.121. The selection-adjusted Monte Carlo p-value is

\hat p_{\mathrm{VT}}=\frac{b+1}{B+1} =\frac{542+1}{10{,}000+1} =\frac{543}{10{,}001} \approx0.05429.

At $\alpha=0.05$, this result is not statistically significant by the permutation-calibrated VT test. The nominal value 0.01696 answers a different question: it treats the observed threshold as though it had been selected in advance. The adjusted value accounts for searching across the three thresholds.

Selection Adjustment and Its Alternatives

Permutation is not the only way to account for threshold selection. A multivariate normal calibration uses the joint null distribution of the threshold scores, including their correlations. Bonferroni correction instead bounds the probability of at least one false rejection by testing each of $K$ thresholds at $\alpha/K$, or by reporting the adjusted value $\min(1,Kp_{\min})$. Bonferroni does not require independent tests, but can be conservative when nested-threshold statistics are highly correlated. A single prespecified cutoff avoids the threshold search, but may be less robust if that cutoff is poorly matched to the set’s allele-frequency spectrum.

Worked Example: Comparing Bonferroni With the VT Result

Use the three-threshold example above. The smallest nominal p-value is 0.01696, and $K=3$. The Bonferroni-adjusted value is

p_{\mathrm{Bonf}}=\min(1,3\times0.01696)=0.05088.

At the 0.05 level, Bonferroni also does not reject the null, although its adjusted value (0.05088) differs from the permutation estimate (0.05429). There is no requirement that the two adjustments agree exactly: Bonferroni uses a general upper bound, while permutation estimates the maximum’s null distribution under the stated exchangeability conditions. The unadjusted 0.01696 is smaller than both because it ignores selection.

Assumptions, Diagnostics, and Failure Modes

A trustworthy VT analysis requires more than computing several scores. The variant set and genotype dosages must be defined consistently; MAF estimates and thresholds should be determined without using the phenotype association results; and weights, thresholds, test direction, and maximum rule should be specified before selecting a result. The score statistic must be calibrated under the null model appropriate to the phenotype and study design. If a normal approximation is used, its accuracy must be adequate, especially when rare alleles yield sparse counts or a small number of carriers.

Useful checks include confirming that thresholds are ordered and that each included variant set is nested; verifying that the same samples and variant quality rules are used across thresholds; and inspecting the threshold-specific scores and their correlation. A sudden statistic change when one variant enters may indicate a genuine contribution, but it can also reveal a genotype, allele-frequency, or coding problem. For permutation analyses, verify that the permutation scheme preserves the design and that increasing the number of permutations stabilizes the estimated tail probability.

Using the smallest nominal p-value. This ignores that the selected threshold was chosen because it looked strongest. Calibrate the maximum or use a valid multiplicity adjustment.
Permuting only the winning threshold. This fails to reproduce the search. Each permutation must repeat the whole threshold scan and selection.
Assuming thresholds are independent. Nested burdens share variants and observations. Use their joint dependence or a valid method that automatically preserves it.
Choosing the direction after seeing the data. A one-sided rule is justified only when the direction was specified in advance. Otherwise use a prespecified two-sided maximum or account for the direction search too.
Interpreting significance as proof of a causal variant. A significant VT burden result is evidence against the tested null for the set-level score procedure; it does not identify which variant is responsible or establish causality.

Interpretation and Limitations

The VT procedure can reduce dependence on a single arbitrary MAF cutoff, but it does not remove all modeling choices. It still depends on the variant set, genotype coding, weights, phenotype model, direction rule, and calibration. A burden statistic also combines variant contributions: opposing effects can cancel, and a selected threshold does not by itself establish that every included variant has the same effect. These are limits of the burden summary, not problems fixed merely by scanning more cutoffs.

Report the tested thresholds, the selection rule, the selected statistic, and the selection-adjusted p-value. State how the null was calibrated and the assumptions supporting that choice. If reporting the winning threshold, describe it as the threshold selected by the procedure, not as an independently confirmed optimum. The inferential claim concerns the prespecified search procedure as a whole.

Key takeaway: A variable-threshold burden test searches across prespecified MAF cutoffs and selects a score, but the p-value must be calibrated for that search. Repeating the entire selection rule under the null—through a joint-distribution calculation or a valid permutation scheme—is what makes the reported evidence interpretable.

Check Your Understanding

Answer each question, then compare with the solutions below.

Why is a nominal p-value at the winning MAF threshold generally not a valid p-value for the VT search?
A directional VT test has four thresholds. In 2,000 valid permutations, 79 maxima equal or exceed the observed maximum. Calculate the plus-one Monte Carlo p-value.
What must be repeated for each permutation in a VT analysis?
Why is the independence formula $1-\Phi(m)^K$ usually unsuitable for nested thresholds?
Name one condition that can make unrestricted phenotype permutation invalid.

Solutions. (1) The winning threshold was selected because its statistic was extreme, so the nominal calculation omits the selection and understates the tail probability. (2) $(79+1)/(2{,}000+1)=80/2{,}001\approx0.03998$. (3) Recalculate statistics at every prespecified threshold and retain the maximum using the same direction and rule as in the observed analysis. (4) Nested scores share variants and are correlated, whereas that formula assumes independent standard normal statistics. (5) Covariates or population structure, for example, can make observations non-exchangeable; the permutation scheme must preserve the relevant design.