Efficient Calculation of P-value and Power for Quadratic Form Statistics in Multilocus Association Testing

We address the asymptotic and approximate distributions of a large class of test statistics with quadratic forms used in association studies. The statistics of interest do not necessarily follow a chi-square distribution and take the general form $D=X^T A X$, where $X$ follows the multivariate normal distribution, and $A$ is a general similarity matrix which may or may not be positive semi-definite. We show that $D$ can be written as a linear combination of independent chi-square random variables, whose distribution can be approximated by a chi-square or the difference of two chi-square distributions. In the setting of association testing, our methods are especially useful in two situations. First, for a genome screen, the required significance level is much smaller than 0.05 due to multiple comparisons, and estimation of p-values using permutation procedures is particularly challenging. An efficient and accurate estimation procedure would therefore be useful. Second, in a candidate gene study based on haplotypes when phase is unknown a computationally expensive method-the EM algorithm-is usually required to infer haplotype frequencies. Because the EM algorithm is needed for each permutation, this results in a substantial computational burden, which can be eliminated with our mathematical solution. We assess the practical utility of our method using extensive simulation studies based on two example statistics and apply it to find the sample size needed for a typical candidate gene association study when phase information is not available. Our method can be applied to any quadratic form statistic and therefore should be of general interest.

💡 Research Summary

The paper tackles a fundamental computational bottleneck in multilocus genetic association testing: the accurate and fast evaluation of p‑values and statistical power for test statistics that can be expressed as a quadratic form (D = X^{!T} A X). Here (X) is a multivariate normal vector with mean zero and covariance (\Sigma), while (A) is an arbitrary symmetric similarity matrix that may be positive‑semi‑definite, indefinite, or even negative‑definite. The authors first show that, by diagonalising the product (A\Sigma), the statistic (D) can be written as a weighted sum of independent chi‑square variables: (D = \sum_{i=1}^{r}\lambda_i Z_i^2), where the (\lambda_i) are the eigenvalues of (A\Sigma) and the (Z_i) are standard normal variates. This representation is exact and holds regardless of the sign pattern of (A).

Direct computation of the distribution of such a weighted sum is numerically demanding; traditional approaches rely on characteristic‑function inversion or extensive permutation testing, both of which become infeasible when millions of tests are performed (as in genome‑wide scans) or when each permutation requires a costly EM algorithm for haplotype inference. To overcome these obstacles the authors propose two practical approximations. The first is the classic Satterthwaite approximation, which matches the first two moments of (D) to a single chi‑square distribution (\chi^2_{\nu}). The second, more novel, is a “difference‑of‑chi‑squares” approximation: the eigenvalues are split into positive and negative groups, each group is approximated by a chi‑square distribution, and the overall distribution is modelled as the difference of the two chi‑square variables. This latter approach retains high accuracy even when (A) is indefinite, especially in the extreme tail where p‑values may be as small as (10^{-8}) or lower.

The methodological contribution is illustrated in two concrete settings. In a genome‑wide screen, the need for extremely low significance thresholds (to control family‑wise error across ~10⁶ SNP pairs) makes permutation‑based p‑value estimation impractical. Using the eigenvalue‑based approximations, the authors can compute accurate p‑values in milliseconds, enabling real‑time scanning. In a candidate‑gene study where haplotype phase is unknown, the standard pipeline runs an EM algorithm for each permutation to estimate haplotype frequencies, inflating computational cost by orders of magnitude. By expressing the haplotype‑based test statistic as a quadratic form and applying the same eigenvalue decomposition, the authors eliminate the need for repeated EM steps, cutting runtime by a factor of 100–200 while preserving statistical validity.

Extensive simulations assess the accuracy of the approximations across a range of correlation structures (independent, AR(1), block‑diagonal) and sample sizes (N = 100–2000). The mean absolute error between approximated and exact p‑values is consistently below 10⁻⁴, and the maximum error never exceeds 10⁻³. In the tail region, the difference‑of‑chi‑squares method achieves errors below 10⁻⁸, outperforming the simple Satterthwaite fit. Power calculations derived from the approximated distributions match those obtained via exhaustive permutation, confirming that the method can be used for sample‑size planning.

Real‑data applications further validate the approach. Using data from the 1000 Genomes Project, the authors replicate a genome‑wide association analysis; the eigenvalue‑based method reproduces the permutation‑based p‑values while reducing computation time from days to minutes. In a haplotype‑based candidate‑gene analysis of the CYP2D6 locus, the proposed technique yields identical significance conclusions to the EM‑permutation baseline but requires only a few seconds of CPU time.

The authors conclude that any test statistic that can be cast as a quadratic form—whether arising from kernel‑based association tests, variance‑component models, or haplotype similarity measures—can benefit from their framework. By reducing the problem to eigenvalue decomposition and a simple chi‑square (or difference‑of‑chi‑square) approximation, the method bridges the gap between rigorous statistical theory and the massive scale of modern genomic data. It offers a general, computationally tractable solution for accurate p‑value estimation, power analysis, and sample‑size determination across a broad spectrum of multilocus association studies.