Tippett’s method (Tippett, 1931), also known as the Dunn-Šidák correction for multiple testing (Šidák, 1957; Dunn, 1958), follows from the fact that the familywise type I error rate for k independent tests, α s = 1 − ( 1 − α p ) k , will equal a desired nominal rate, α, if we set α _p = 1 − (1 − α)^1/k . Hence, the joint null distribution can be rejected if min(p _i ) ≤ 1 − (1 − α)^1/k . Analogously, we can use p = 1 − 1 − min p 1 , … , p k k (2)

as the p-value for the gene. As opposed to the Bonferroni method, which is slightly conservative even when all tests are independent, the method provides exact control of the type I error rate, but only under independence.

Under the joint null hypothesis, r ∼Binomial (k, α _p ) where r denotes the number of tests that are significant at α _p . Therefore, we can reject the joint null hypothesis if p = ∑ x = r k k x α p x 1 − α p k − x (3)

is equal to or less than the desired type I error rate (Wilkinson, 1951). Intuitively, we can interpret this method as a test of “excess significance” of the SNPs within a gene. For example, the chances of finding r = 10 or more significant SNPs at α _p = 0.05 in a gene with 100 independent SNPs is approximately p = 0.028, which would be significant at α = 0.05.

Assuming that p _i ∼Uniform (0, 1) under the null hypothesis of no association, it is easy to show that −2 ln(p _i ) is chi-square distributed with 2 degrees of freedom. Hence, the combined test statistic X 2 = − 2 ∑ i = 1 k ln p i (4)

follows a chi-squared distribution with 2k degrees of freedom under the joint null hypothesis (Fisher, 1932). The p-value for the gene can therefore be computed with p = 1 − F (X ², 2k), where F (⋅, 2k) denotes the cumulative distribution function of a chi-square distribution with 2k degrees of freedom.

Let Φ(⋅) denote the cumulative distribution function of the standard normal distribution and Φ⁻¹(⋅) its inverse. Since Φ⁻¹ (1 − p _i ) follows a standard normal distribution under H _0i , z = ∑ i = 1 k Φ − 1 ( 1 − p i ) / k ∼ Normal ( 0,1 ) under the joint null hypothesis (Stouffer et al., 1949). The p-value for a gene is then computed with p = 1 − Φ(z).

One potential approach to adjust the previous methods is to quantify the degree of dependence between the tests, estimate the effective number of independent tests based on this information, and incorporate this estimate into the methods described above.

The degree of dependence between the tests is closely related to the strength of the association between the SNPs. The latter can be quantified with various statistics (e.g., D, D′, r, or r ²) expressing the degree of LD between pairs of SNPs (Laird and Lange, 2010). We can use one of these measures to construct a k × k association matrix for all SNPs, sometimes called an “LD map”. The effective number of tests can then be estimated based on this association matrix. A variety of approaches have been described in the literature for this purpose (Cheverud, 2001; Nyholt, 2004; Li and Ji, 2005; Gao et al., 2008; Galwey, 2009). A common feature of all methods is that they start by applying PCA to the association matrix. We use λ _i to denote the ith eigenvalue extracted from the PCA.

The method proposed by Cheverud (2001) and Nyholt (2004) estimates the effective number of tests with k eff CN = 1 + k − 1 1 − Var λ k , (5) where Var(λ) is the variance of the k eigenvalues. On the other hand, Li and Ji (2005) suggested the formula k eff LJ = ∑ i = 1 k f | λ i | , (6) where f x = I x ≥ 1 + x − ⌊ x ⌋ (7) and ⌊⋅⌋ is the floor function. According to the method by Gao et al. (2008), we first sort the eigenvalues in decreasing order, letting λ ₍₁₎ denote the largest and λ _(k) the smallest eigenvalue. Then the effective number of tests is defined as k eff GAO = min x such that ∑ i = 1 x λ i ∑ i = 1 k λ i > C , (8) where C is a user-defined parameter and usually chosen to be 0.995. Finally, Galwey (2009) proposed to estimate the effective number of tests with k eff GAL = ∑ i = 1 k λ i ′ 2 ∑ i = 1 k λ i ′ , (9) where λ i ′ = max ( 0 , λ i ) .

All of the methods described above have the following desirable properties. When applied to an identity matrix (i.e., when there is no association between any pair of SNPs), then k _eff = k, so that the effective number of tests is equal to the number of SNPs. An exception to this property can occur with k eff GAO . Depending on the value of C and the number of tests, it can happen that the effective number of tests is then estimated to be less than k (i.e., when k (1 − C) > 1 then k eff G A O < k ). On the other hand, when all of the SNPs are perfectly associated (i.e., the correlation matrix is equal to a k × k matrix of 1’s), then k _eff = 1. In essence, the same test is then repeated k times, yielding identical results, so that effectively only a single test has been carried out. However, the methods differ in how association matrices that fall in between these two extremes are handled, yielding varying estimates of the effective number of tests between 1 and k.

Once k _eff has been estimated with one of these approaches, it can be used to adjust each of the methods for combining p-values described earlier. For the Bonferroni and Tippett’s methods, we substitute k _eff for k so that p = min 1 , min p 1 , … , p k × k eff . (10) and p = 1 − 1 − min p 1 , … , p k k eff (11) are then the p-values for the gene. For the binomial test, we first define r ̃ = ⌊ r × k eff k ⌋ as the adjusted (i.e., effective) number of significant SNPs within the gene. Then the p-value for the gene is computed with p = ∑ x = r ̃ k eff k eff x α p x 1 − α p k eff − x . (12) Use of the floor function for computing r ̃ may be conservative, but we consider this preferable over rounding and the risk of a too liberal test. Fisher’s method can be adjusted by replacing the degrees of freedom of the chi-square distribution with 2k _eff and adjusting the test statistic with X ̃ 2 = k eff k × X 2 . Hence, the p-value for the gene is then computed with p = 1 − F ( X ̃ 2 , 2 k eff ) . Finally, for Stouffer’s method, we let z ̃ = k eff k × z denote the adjusted test statistic and hence p = 1 − Φ ( z ̃ ) is then the p-value for the gene.

Another approach to account for dependence is to make use of permutation testing (Johnson et al., 2010; Moskvina et al., 2011). The idea is to empirically derive the null distribution of the test statistic of interest by reshuffling the data in such a way that relevant features of the data structure are preserved except for the actual association being tested. For example, when testing for the association between each SNP and case-control status, reshuffling the status variable breaks any existing associations, but keeps the LD structure of the SNPs intact. Hence, any dependence among the p-values to be combined using one of the methods described earlier is automatically incorporated into the null distribution. The p-value for a gene is then computed from the percentile of the actually observed test statistic under the empirical null distribution. Note that in the present case, the test statistic of interest is actually a p-value itself, so letting p _j denote the combined p-value based on the jth permutation of the data (with j = 1, … , s) and p _obs the observed combined p-value, the p-value for a gene is then given by p = ∑ j = 1 s I ( p j ≤ p obs ) / s .

Permuting the data in the manner described above requires access to the raw data, so that the phenotype variable can be reshuffled and the test of association can be conducted for each SNP. In addition, repeatedly computing the test of association for each SNP within a gene can be computationally demanding. We can reduce the computational burden and eliminate the dependence upon the raw data by directly generating p-values based on an association matrix that reflects the degree of LD among the SNPs (Liu et al., 2010) which may be obtained from a reference population and not necessarily the given data.

In particular, let R denote the LD association matrix constructed from the correlations among the SNPs. We can then quickly generate a large number (s) of samples from a multivariate normal distribution with a true mean vector equal to zeros and covariance matrix R. Let Z denote the s × k matrix of these values and P = 2 (1 − Φ(|Z|)) the matrix of two-sided p-values obtained by applying Φ(⋅) element-wise. For each row in P, we then apply one of the methods for combining p-values, yielding p _j . The p-value for a gene is then again computed as described above.

The last set of methods we will consider are modifications of Fisher’s and Stouffer’s method so that dependence among the tests is directly taken into consideration.

Figure 1 shows boxplots of the type I error rates of all methods observed on all 30,910 genes applied to the original HapMap dataset (i.e., based on the non-bootstrapped data). Individual genes are indicated as points when their rate was more than 1.5 times the interquartile range above the third or below the first quartile. The mean rejection rates (and SDs) are also indicated in the figure.

None of the unadjusted methods could achieve on average a nominal rejection rate of α = 0.05. As expected, the Bonferroni method tended to be conservative, as was Tippett’s method, which produced very similar results to the former method throughout the entire simulation study and which will therefore not be further discussed (note that rates above 0.05 for the Bonferroni method–which occurred for about 1.5% of the genes–reflect simulation error, since we know that the method guarantees that the type I error rate is equal to or less than α regardless of the degree of dependence). On the other hand, the remaining methods were generally liberal, at times dramatically so, with the binomial test at least providing an average rejection rate closest to the nominal level.

The adjustments for addressing the dependence did bring the rejection rates closer to the nominal level with varying degrees of success. In fact, when adjusted with the Li & Ji method, the Bonferroni method had a nominal average rejection rate, although this came at the cost of increased variability in the type I error rates, and the occurrence of rates well above the nominal level for particular genes. For the binomial test, the average rates fluctuated around the nominal level, being slightly conservative with the Li & Ji and Galwey adjustments and slightly liberal with the Nyholt and Gao adjustments. In contrast, none of the PCA-based adjustments could bring the average type I error rates of the Fisher and Stouffer methods sufficiently close to α = 0.05.

Using empirically-derived null distributions produced rejection rates that were on average reasonably close to the nominal level, especially for Fisher’s and Stouffer’s methods. Moreover, the type I error rates of individual genes had much lower variability than the rates obtained with the PCA-based adjustments. This was also true for the Bonferroni method and the binomial test, but these methods were slightly conservative on average when adjusted in this manner.

The (two-sided) generalization of Fisher’s method to dependent tests (i.e., Brown’s method) yielded a nominal rejection rate on average. Furthermore, the variability (i.e., SD) of the rates for individual genes was lowest compared to all other methods. Quite importantly (mis)application of the one-sided version of Brown’s method (since the p-values for the SNPs were computed from two-sided tests) resulted in worse performance (further references to Brown’s method will therefore pertain to the two-sided version unless otherwise stated). On the other hand, the generalization of Stouffer’s method to dependent tests (i.e., Strube’s method) performed reasonably well, although its type I error rate was on average slightly inflated.

Supplementary Figures S1–S3 show the type I error rates of the methods based on bootstrap samples of sizes 102, 500, and 1,000, respectively. A comparison of Figure 1; Supplementary Figure S1 (both with n = 102) shows that the performance of the methods was similar regardless of whether they were applied to the original data or to bootstrap samples of the same size. The only exception to this was Stouffer’s method, which became slightly more conservative for the bootstrapped data. Also, the patterns in the type I error rates of the methods were not fundamentally altered when applied to larger sample sizes. Using empirical distributions in combination with Fisher’s and Stouffer’s methods and Brown’s method generally resulted in adequate control of the type I error rate on average and comparatively low variability in the rates for individual genes.

To examine whether the performance of the methods was affected by certain characteristics of the genes, we examined their type I error rates as a function of the (log transformed) number of SNPs in the genes, the average correlation in the LD maps, the degree of variability (i.e., standard deviation) of the correlations, the square-root of the mean squared correlations (SRMSC), the average minor allele frequencies (MAF) of the SNPs in the genes, and the standard deviation (SD) of the MAFs. The SRMSC was of particular interest as it distinguishes genes whose SNPs are independent (SRMSC equal to 0) from genes with SNPs in strong LD regardless of the directionality of the association (SRMSC close to 1). We used locally estimated scatterplot smoothing to visualize the relationship between these characteristics and the rejection rates for each method.

Figure 2 shows that the type I error rate of many methods was affected by the number of SNPs in the genes. In particular, the Bonferroni method became increasingly conservative as the number of SNPs increased. Interestingly, this dependence on k was essentially removed when using the adjustment of Li & Ji and, to a slightly lesser extent, the adjustment of Gao. In contrast, the binomial test and Fisher’s method became increasingly liberal as a function of k, whereas Stouffer’s method displayed non-monotonic behaviour. The PCA-based adjustments helped to reduce the inflation in the type I error rates of these methods, but could not eliminate the dependence on k. Furthermore, all methods adjusted based on empirical distributions became increasingly conservative as the number of SNPs increased. Finally, Brown’s method yielded essentially nominal rates regardless of k, except for very large genes, where the method became slightly conservative.

Figure 3 shows the type I error rates as a function of the SRMSC values. As expected, the figure points out the increasingly conservative behavior of the Bonferroni method as the SNPs within the genes become more dependent, while Fisher’s and Stouffer’s methods then become liberal. The conservative behavior of the binomial test under independence is also expected (due to the discrete nature of the binomial distribution, the type I error rate of the test will not exceed α = 0.05, but will often fall well below it). More surprising is the fact that the type I error rate of the method was essentially nominal for genes with very strong LD. To understand this phenomenon, consider a gene with k SNPs in perfect LD. In that case, all (two-sided) p-values are identical and hence either none or all k SNPs are significant. Since the latter will happen (under the joint null) with probability α, the test will exhibit nominal performance under this extreme scenario. With respect to the adjustments, the various PCA-based approaches again had the effect of counteracting the conservativeness of the Bonferroni method, while leading to a reduction in the type I error rates of the other methods. For all methods, adjusting with the use of empirical distributions was most successful when LD is strong, while Brown’s method, although slightly conservative under independence, and performed well over the range of SRMSC values. Strube’s method performed similarly, but with some slight inflation for larger SRMSC values.

Supplementary Figures S4, S5 show the type I error rates of the methods as a function of the average correlation and the SD of the correlations of the SNPs within the genes. Here, we again find that changes in these characteristics have essentially no impact on the performance of Brown’s method, as well as when using empirical distributions in combination with Fisher’s and Strube’s methods. Interestingly, the one-sided version of Brown’s method performed similarly to the two-sided one when the average LD was larger than ≈ 0.3 .

Supplementary Figures S6, S7 display the type I error rates of the methods as a function of the average MAFs and their SD within the genes. The performance of Brown’s method and Fisher’s method with the empirical distribution adjustment was again not affected substantially by these factors except that these methods were slightly conservative when the average MAF was below 0.1 within the genes. Stouffer’s method adjusted based on empirical distributions and its generalization to dependence was also relatively robust to both factors.

Figure 4 illustrates the power of the methods (averaged over genes) for the 54 conditions where the joint null hypothesis was false. Each panel corresponds to one of three sample size conditions (i.e., 102, 500, and 1,000), while the x-axis indicates the condition, starting with the two “single SNP” conditions (with effect sizes of 0.2 vs. 0.5) followed by the 16 “multiple SNPs” conditions (with either 5% or 20% of SNPs selected, an effect sizes of 0.2 or 0.5, a non-distributed or distributed effect, and either non-compact or compact SNP positions). We only show the power rates for the (unadjusted) Bonferroni method and those method and adjustment combinations that could control the type I error rate on average (i.e., the Bonferroni method with the Li & Ji adjustment, Brown’s method, and use of empirical distributions in combination with Fisher’s and Stouffer’s methods).

As expected, power increased with the sample size, with the effect size, when a higher percentage of SNPs was associated with the phenotype, and when the effect was not distributed across these SNPs. Whether the selected SNPs fell into a compact region of the gene or were distributed throughout had comparatively little influence on the results. The (unadjusted) Bonferroni method typically had lower power compared to the other methods especially for lower sample sizes. Also, Stouffer’s method adjusted with the use of empirical distributions tended to have slightly lower power compared to the other adjusted methods. While differences between the other methods were often negligible, a slight power advantage could be observed for the Bonferroni method adjusted with the Li & Ji correction when a single SNP or a low percentage of them contained a strong signal. Otherwise, Brown’s method tended to have a slight power advantage.

To compare the computational efficiency of the methods, Supplementary Figure S8 presents the average computation times (based on 500 iterations) of the unadjusted methods along with the Li & Ji, empirical (using 10⁴ samples), and dependence adjustments (the latter only for the Fisher and Stouffer methods) for a set of genes with {10, 25, 50, 100, 199, 254, 492, 897, and 1150} SNPs. While the unadjusted methods show no noteworthy increase in computational times as a function of the gene size, the results show that the use of the adjustments does come at the cost of increased computational times, less so for the Li & Ji adjustment and more so when using the Brown and Strube methods. Finally, as expected, the empirical methods demand the highest amount of computer time (although this can be mitigated to some extent; see Cinar and Viechtbauer (2022) for details). The computation times for the different test types, however, did not appear to vary substantially.