Sample-Size Planning for Rare-Variant Set Tests

Why Sample-Size Planning Depends on Genetic Architecture

A rare-variant set test can have little power even when some variants have substantial effects. The reason is that power depends not only on how many participants are analyzed, but also on allele frequencies, effect directions, the test's variant weights, and the relationship among variants. A burden test combines variant scores, so effects in opposite directions can cancel. A variance-component test, such as SKAT, aggregates squared scores and can retain signal from effects in either direction. Planning should therefore reflect a plausible genetic architecture rather than assume every test responds equally well.

The calculations below use a quantitative trait to make the relationship between sample size and genetic signal explicit. The same planning logic applies to binary outcomes, but the information per participant must then be obtained under the specified case-control or logistic model; substituting quantitative-trait formulas without adjustment is not justified.

Model and Assumptions

Suppose \(n\) unrelated participants have a quantitative outcome \(Y\), a vector of covariates, and genotypes at \(m\) variants. Genotype \(G_j\) is coded as the allele count \(0,1,\) or \(2\), and \(\beta_j\) is the change in the outcome mean per additional effect allele at variant \(j\). For this tutorial, effects are measured in outcome units and the residual standard deviation is \(\sigma\). The model assumes independent participants, a correctly specified mean model and covariate adjustment, and a residual variance that is adequately represented by \(\sigma^2\). The score-test approximations also rely on enough information for the score vector to be approximately multivariate normal.

Let \(S\) denote the covariance matrix of the genotype score vector after covariate adjustment and scaling by residual variance. Under the null, the score vector \(U\) is approximately \(N(0,S)\); under an alternative with effect vector \(\beta\), it is approximately \(N(S\beta,S)\). If participants contribute similar information independently, write \(S \approx nS_0\), where \(S_0\) is the per-participant information matrix. This linear scaling with \(n\) is the basis of the sample-size calculations.

For a simple analytic calculation, assume variants are in linkage equilibrium, are in Hardy-Weinberg equilibrium, and are not strongly associated with the covariates. If \(p_j\) is the minor-allele frequency, the genotype variance is \(v_j=2p_j(1-p_j)\). With independent variants and no additional covariate adjustment, \(S_0\) is approximately diagonal with entries \(v_j/\sigma^2\). These assumptions are conveniences, not requirements of the tests: with linkage disequilibrium (LD), use the full genotype covariance matrix rather than the diagonal approximation.

Planning principle: First specify the significance level, target power, phenotype model, allele frequencies, effect sizes, and test weights. Then calculate the information per participant under those assumptions. A sample-size estimate is only as credible as these inputs.

Burden-Test Sample Size

A burden test forms a weighted sum of variant scores. Let \(a_j\) be the prespecified burden weight for variant \(j\), and let \(a\) be the vector of these weights. The score is \(B=a^\mathsf{T}U\). Under the alternative, its mean is \(a^\mathsf{T}S\beta\), and its variance is \(a^\mathsf{T}Sa\). Thus its squared, standardized statistic has approximately a noncentral chi-square distribution with one degree of freedom:

T_B=\frac{(a^\mathsf{T}U)^2}{a^\mathsf{T}Sa} \ \dot{\sim}\ \chi^2_1(\lambda_B), \qquad \lambda_B=\frac{(a^\mathsf{T}S\beta)^2}{a^\mathsf{T}Sa}.

Here \(\lambda_B\) is the noncentrality parameter: it measures the expected signal relative to score variability. Since \(S\approx nS_0\), the parameter grows in proportion to sample size. Under independent variants with genotype variances \(v_j\), the expression becomes

\lambda_B \approx n\, \frac{\left(\sum_{j=1}^{m}a_jv_j\beta_j\right)^2} {\sigma^2\sum_{j=1}^{m}a_j^2v_j}.

For a two-sided one-degree-of-freedom test, choose \(\lambda_*\) to give the desired power at the prespecified \(\alpha\), using the noncentral chi-square distribution. Then \(n\approx\lambda_*/(\lambda_B/n)\). A normal approximation sometimes uses \(\lambda_*\approx(z_{1-\alpha/2}+z_{\text{power}})^2\); for stringent significance levels or final planning, numerical inversion of the noncentral chi-square distribution is preferable.

Worked Example: Planning a Burden Test for Aligned Effects

Consider three hypothetical independent variants with minor-allele frequencies \(0.002\), \(0.004\), and \(0.006\). Suppose all have the same effect, \(\beta_j=0.10193\) outcome standard deviations per allele, the residual standard deviation is one outcome standard deviation, and the burden weights are all one. Assume the model and covariance conditions above hold.

The genotype variances are \(v_1=2(0.002)(0.998)=0.003992\), \(v_2=2(0.004)(0.996)=0.007968\), and \(v_3=2(0.006)(0.994)=0.011928\). Their sum is \(0.023888\). Because the effects have the same direction, the per-participant burden noncentrality is

\frac{\lambda_B}{n} =\frac{(0.023888\times0.10193)^2}{0.023888} =(0.10193)^2(0.023888) \approx 0.0002482.

Suppose the design requires \(\alpha=5.6\times10^{-7}\) and 80% power. Inverting the one-degree-of-freedom noncentral chi-square power calculation gives a target noncentrality of approximately \(\lambda_*=34.19\). Therefore,

n\approx\frac{34.19}{0.0002482} =137{,}752 \approx 1.378\times10^5.

The estimate is about 137,800 analyzable participants under the stated assumptions. It is not a universal sample size for these allele frequencies: different effects, weights, residual variance, significance level, or LD would change the result. Enrollment should also account for expected exclusions or missing data, which are not included in this calculation.

Variance-Component Tests and Mixed Directions

A variance-component test uses a quadratic form in the variant scores. Write \(W\) for a diagonal matrix of nonnegative SKAT weights; the statistic is \(Q=U^\mathsf{T}WU\). Its null distribution is generally a mixture of chi-square distributions, with mixture weights determined by the eigenvalues of \(S^{1/2}WS^{1/2}\). Under an alternative, the distribution is also affected by the shifted mean \(S\beta\). Consequently, a single one-degree-of-freedom noncentrality parameter is generally not enough to describe SKAT power.

For independent variants, \(Q\) can be expressed as a weighted sum of independent squared normal variables. For unit SKAT weights and a common residual variance, the noncentrality in the standardized \(j\)th component is \(n v_j\beta_j^2/\sigma^2\). Unlike a burden score, the squared contributions do not cancel just because effects have opposite signs. For general LD, unequal weights, or a stringent significance threshold, calculate power from the appropriate mixture distribution or by simulation under the planned model.

Worked Example: Opposite Effects Favor a Variance-Component Test

Suppose two hypothetical independent variants each have minor-allele frequency \(0.005\), residual standard deviation one, and unit SKAT weights. Their effects are \(+0.20\) and \(-0.20\) outcome standard deviations per allele. The genotype variance at each variant is \(v=2(0.005)(0.995)=0.00995\).

For an unweighted burden test, the two terms in the expected score are \(v(0.20)\) and \(v(-0.20)\), so their sum is zero. Its noncentrality is therefore zero under this symmetric model, even though both variants affect the outcome. For SKAT, the total standardized component noncentrality is

\lambda_{\text{SKAT}} =n\frac{v(0.20)^2+v(-0.20)^2}{\sigma^2} =n\{0.00995(0.04)+0.00995(0.04)\} =0.000796n.

For example, at \(n=10{,}000\), this gives \(\lambda_{\text{SKAT}}=7.96\), while the burden noncentrality remains zero in this idealized calculation. The SKAT statistic has a mixture distribution, so 7.96 alone does not give its power; the null eigenvalues and the chosen \(\alpha\) are also needed. This calculation does show why a burden-based sample-size estimate would be misleading for this architecture.

Worked Example: How Burden Weights Change the Signal

Consider two hypothetical independent variants, with \(p_1=0.01\), \(p_2=0.002\), effects \(\beta_1=0.10\), \(\beta_2=-0.30\), and \(\sigma=1\). Their genotype variances are \(v_1=0.0198\) and \(v_2=0.003992\). At \(n=10{,}000\), compare an unweighted burden, a burden with \(a_j=1/\sqrt{v_j}\), and a unit-weight SKAT signal.

For the unweighted burden, the per-participant noncentrality is

\frac{(0.0198(0.10)+0.003992(-0.30))^2}{0.0198+0.003992} =\frac{(0.0007824)^2}{0.023792} \approx 0.00002573.

Thus its noncentrality at 10,000 participants is approximately \(0.2573\). For the frequency-weighted burden, the numerator before squaring is \(\sqrt{v_1}(0.10)+\sqrt{v_2}(-0.30)\approx0.14071(0.10)-0.06318(0.30)=-0.004883\), and the denominator is \(a_1^2v_1+a_2^2v_2=2\). The per-participant noncentrality is approximately \((-0.004883)^2/2=0.00001192\), or \(0.1192\) at \(n=10{,}000\). These weights do not improve power for this particular pattern of opposing effects.

For unit-weight SKAT, the component signal is \(v_1\beta_1^2+v_2\beta_2^2=0.0198(0.01)+0.003992(0.09)=0.00055728\) per participant, or \(5.573\) across 10,000 participants. This is not directly comparable to a one-degree-of-freedom burden noncentrality as a power value; SKAT's mixture null distribution must also be used. The contrast illustrates why weights and effect architecture must be considered together.

Calculation, Diagnostics, and Alternatives

A practical planning calculation proceeds from a specified genetic architecture to a test-specific power calculation. For burden tests, the closed-form noncentrality above is often convenient. For SKAT and other quadratic tests, numerical evaluation of the mixture distribution is usually more appropriate. Simulation is a useful alternative when the model includes LD, covariates, relatedness, non-normal outcomes, or other features not represented by a simple covariance formula.

Specify the design.
Choose the phenotype model, analysis sample size, covariates, significance level, desired power, and test weights before calculating.

Specify the genetic architecture.
Use plausible allele frequencies and effect sizes, including effect directions and any assumptions about which variants are causal.

Calculate information.
Use the burden noncentrality for a burden test; for a quadratic test, calculate its null mixture and alternative distribution, accounting for LD when present.

Check sensitivity and feasibility.
Repeat the calculation across plausible frequencies, effects, and weights, then distinguish analyzable sample size from required recruitment.

A useful diagnostic is to inspect the expected score contribution \(a_jv_j\beta_j\) for each variant. Contributions with opposing signs flag possible cancellation in a burden test. For a quadratic test, inspect the covariance eigenvalues and the variant contributions to the signal; near-zero eigenvalues can indicate little independent information, while very few observed minor-allele carriers make asymptotic approximations questionable. These checks complement, rather than replace, null calibration diagnostics for the planned test.

Analytical formulas are useful for understanding scaling and comparing design assumptions, but they are not guarantees of realized power. Rare alleles may be absent or represented by very few carriers in a finite sample. Estimates of allele frequencies and effects may be uncertain, and weights selected from annotations may not match the true effect pattern. In a case-control design, imbalance between cases and controls changes information; power must be calculated under the specified logistic model and sampling ratio. Use simulations that reproduce the intended analysis when these features are material.

Common Planning Mistakes

Assuming every variant has the same effect direction. That assumption can inflate burden power when real effects oppose one another. Calculate scenarios with plausible direction patterns.
Interpreting a SKAT signal as a burden noncentrality. A quadratic test generally has a mixture null distribution. Use its eigenvalues and alternative distribution to calculate power.
Using allele frequencies but ignoring LD. The diagonal genotype-variance approximation fails when variants are correlated. Use the full covariance matrix or simulate genotypes with realistic LD.
Reporting recruited \(n\) as the analyzable sample size. Missingness, quality control, and exclusions can reduce the number contributing to the test. Inflate recruitment separately using defensible retention assumptions.

Key takeaway: Sample size is determined by test-specific information under a specified genetic architecture. Burden calculations depend on the weighted, signed sum of effects; variance-component calculations retain squared component signals but require their mixture distribution for power.

Check Your Understanding

Try each question, then compare your reasoning with the solutions below.

For an independent-variant burden test, what happens to the expected burden signal if two equally frequent variants have equal and opposite effects and equal weights?
With \(\lambda_B/n=0.0005\) and a target noncentrality of 10, what is the approximate required sample size?
Why is the sum of component noncentralities not by itself a complete SKAT power calculation?
Name two assumptions behind using \(v_j=2p_j(1-p_j)\) as the genotype variance in the analytic formulas.

Solutions. (1) The expected burden contributions cancel: \(a v\beta+a v(-\beta)=0\), so the burden noncentrality is zero in this idealized case. (2) \(n\approx10/0.0005=20{,}000\) analyzable participants. (3) SKAT has a mixture null distribution determined by the covariance and weights; power depends on its eigenvalues and the alternative, not only on a summed signal. (4) The formula assumes Hardy-Weinberg equilibrium and uses the minor-allele frequency for allele-count genotype coding; the diagonal approximation additionally assumes no LD. Covariate adjustment may also change the effective covariance.