Combining dependent p-values to assess the global null hypothesis has long been a fundamental challenge in statistical inference. A common scenario arises when integrating the results of various methods on the same dataset to enhance signal detection power (Wu et al., 2016;Rosenbaum, 2012). When individual p-values have arbitrary dependence, the Bonferroni test is the most common approach with a theoretical guarantee. However, it is often criticized for being overly conservative in practical applications.
Specifically, consider n individual p-values P 1 , . . . , P n . To test the global null hypothesis, i.e., all n null hypotheses are true, the Bonferroni test calculates the combined p-value as n × min (P 1 , . . . , P n ). Due to the scaling factor n, the Bonferroni combined p-value may exceed any of the individual p-values, leading to a loss of power during the combination process.
Recently, a novel approach gaining traction involves the combination of p-values through transformations based on heavy-tailed distributions (Liu et al., 2019;Wilson, 2019a). Let X i be defined as Q F (1 -P i ), where F(•) represents the cumulative distribution function of a heavytailed distribution and Q F is its quantile function. The core idea is to compute the combined p-value based on the tail distribution of S n = n i=1 X i , which under the global null is robust to dependence among the heavy-tailed variables X 1 , . . . , X n . The Cauchy combination test, which sets F as the standard Cauchy distribution, was first introduced in Liu et al. ( 2019) for genomewide association studies (GWAS) and has since been applied in genetic and genomic research, including spatial transcriptomics (Sun et al., 2020), ChIP-seq data (Qin et al., 2020), and singlecell genomics (Cai et al., 2022). Another popular method, the harmonic mean p-value (Wilson, 2019a), employs the Pareto distribution with shape parameter γ = 1 as F.
Despite the growing popularity of these heavy-tailed combination tests in practical applications, there has been limited theoretical investigation and empirical evaluation of these methods. Existing studies (Liu & Xie, 2020;Fang et al., 2023) have provided asymptotic validity of these tests as the significance level α → 0 for pairwise bivariate normal test statistics. These results closely related to earlier findings on sums of regularly varying tail variables, showing that pr (S n > x) and n {1 -F(x)} are asymptotically equivalent as x → +∞, provided that the variables X 1 , . . . , X n are pairwise quasi-asymptotically independent (Chen & Yuen, 2009). Intuitively, for heavytail distributed X 1 , . . . , X n , their maximum typically dominates the sum, making the latter less sensitive to dependence among X 1 , . . . , X n . Yet this same intuition raises doubts about the true benefits of these tests compared to the Bonferroni test. Additionally, the assumption of quasiasymptotic independence, while covering any bivariate normal variables that are not perfectly correlated, remains more stringent than allowing arbitrary dependence. For example, bivariate t-distributed variables, which are frequently used as test statistics, are not quasi-asymptotically independent. This raises questions about the robustness of these tests when faced with unknown dependence structures.
This paper addresses these concerns through theoretical and empirical analyses. Many applications employ heavy-tailed combination tests to aggregate results from different methods or studies, often in settings where the number of base hypotheses, n, is moderate rather than excessively large. Accordingly, we focus on scenarios where n is fixed and analyze the asymptotic regime as the significance level α → 0. Our theoretical investigation shows that when test statistics are quasi-asymptotically independent, particularly when they follow a bivariate normal distribution with imperfect correlation, the rejection regions of heavy-tailed combination tests are asymptotically equivalent to those of the Bonferroni test as α approaches zero. This suggests that in the same asymptotic regime where combination tests have proven to be valid, they offer no real power advantage over Bonferroni's approach. However, when the assumption of asymptotic independence is violated, such as when test statistics follow a multivariate t distribution, our empirical results indicate that combination tests still appear to be asymptotically valid when the tail index γ 1, despite the lack of a theoretical guarantee. More strikingly, they exhibit significantly greater power than the Bonferroni test, highlighting their potential advantages in settings where p-values are strongly dependent, a scenario that often arises when aggregating results from different methods applied to the same dataset. Furthermore, through simulations and real-world case studies, we observe that the empirical validity and power of these tests are affected by both the heaviness and support of the heavy-tail distribution.
Consider n test statistics T 1 , . . . ,T n , where each T i is for a base null hypothesis H 0,i . For each base hypothesis, we construct a one-sided or two-sided base p-value P i based on the distribution of T i under H 0,i . We are interested in testing the global null hypothesis
The test statistics T 1 , . . . ,T n may exhibit unknown dependence structures among each other.
For the heavy-tailed combination tests, we apply a transformation of the p-values into quantiles of heavy-tailed distributions. Specifically, let F denote the cumulative distribution function (CDF) of the heavy-tailed distribution and Q F represent its quantile function, defined as
We define the individual transformed test statistics as
A combination test can then be constructed based on the sum
or more generally, any weighted sum S n, ì ω = n i=1 ω i X i with non-random positive weights ω i s.
We begin by reviewing existing theoretical results on the tail properties of S n . If X 1 , . . . , X n belong to the sub-exponential family, a major class of heavy-tailed distributions, it is well-known that the tail probability of S n = X 1 + • • • + X n is asymptotically equivalent to the sum of individual tail probabilities under the assumption that the X i s are mutually independent. That is,
where F = 1 -F denotes the tail probability (Embrechts et al., 2013). When the independence assumption fails, previous works (Chen & Yuen, 2009;Asmussen et al., 2011;Albrecher et al., 2006;Kortschak & Albrecher, 2009;Geluk & Ng, 2006;Tang, 2008) have shown that (1) still holds for different subclasses of sub-exponential distributions under certain assumptions of the dependence structure.
Here, we restate several key results that form the foundation of the theoretical properties of the heavy-tailed combination tests, which will be detailed in Section 2.3. For any variable X, we denote X + = max (X, 0) and X -= max (-X, 0). To begin, we introduce the concepts of quasi-asymptotic independence and the consistently-varying subclass C of sub-exponential distributions, following Chen & Yuen (2009).
Definition 1 (Quasi-asymptotic independence). Two non-negative random variables X 1 and X 2 with cumulative distribution functions F 1 and F 2 , are quasi-asymptotically independent if
More generally, two real-valued random variables, X 1 and X 2 , are quasi-asymptotically indepen-
When X 1 and X 2 have the same marginal distribution, (2) can be rewritten as pr(X
Theorem 3.1 in Chen & Yuen (2009) established the asymptotic tail probability of S n for distributions within C , provided that quasi-asymptotic independence holds.
Theorem 1 (Theorem 3.1 of Chen & Yuen (2009)). Let X 1 , . . . , X n be n pairwise quasiasymptotically independent real-valued random variables with distributions F 1 , . . . ,
The asymptotic equivalence (3) can hold for broader subclasses of heavy-tailed distributions beyond C under stronger dependence assumptions. For instance, Geluk & Tang (2009) provided the necessary dependence structure requirements for this equivalence to hold for dominatedvarying tailed and long-tailed random variables. Additionally, Asmussen et al. (2011) verified this for log-normal distributions when coupled with a Gaussian copula. However, Botev & L'Ecuyer (2017) showed that convergence in (3) can be extremely slow for log-normal distributions, requiring the tail probability to be as small as 10 -233 to achieve reasonable approximations.
Moreover, researchers have observed asymptotic equivalence between the tail probability of the S n and that of max(X 1 , . . . , X n ).
Corollary 1. With the same setting as in Theorem 1, the tail probability of the sum and the maximum has the following relationship
Remark 1. We provide a proof of Corollary 1 in Supplementary S3.2, which essentially restates earlier results (Geluk & Ng, 2006;Tang, 2008;Ko & Tang, 2008), to facilitate understanding for interested readers.
Table 1 presents a list of common distributions in C . All of these distributions also belong to a smaller subclass, the regularly varying tailed distributions R, defined as follows:
Table 1: Regularly varying tailed distributions and their tail indices. Φ is the cumulative distribution function of a standard normal distribution. Γ is the gamma function. J(s, x) = ∫ ∞ x t s-1 e -t dt is the incomplete gamma function and
dt is the regularized incomplete eta function, Ft (c) is the survival function at c of the corresponding t distribution with the same degree of freedom γ Distributions: Survival Function Tail index Support Cauchy:
Following Cline (1983), the parameter γ is referred to as the tail index, characterizing the tail heaviness (Teugels et al., 1987) of a distribution. Distributions with a smaller γ exhibit heavier tails. For example, for the Student's t distribution, γ is the same as the degree of freedom, with the Cauchy distribution being a special case with γ = 1.
In Table 1, all distributions, except for the Student's t distributions that includes the Cauchy distribution, have a lower bound in their support. In contrast, the Student's t distributions have symmetric densities around the origin, and their supports cover the entire real line. As a consequence, when p i approaches 1, the transformed test statistics X i can become substantially negative, which may affect both the power and type-I error control in the associated combination tests. To address this issue, we introduce the left-truncated Student's t distribution in Table 1, defined as a conditional Student's t distribution with a left-bounded support interval of [c, +∞). Specifically, we define F t,γ (x) = pr(X x) with X following a Student's t distribution with degree of freedom γ. The cumulative distribution function of the left-truncated t distribution is
With this definition, the left-truncated t distribution remains a regularly varying tailed distribution with the same tail index γ, as proved in Proposition S1. In our experiments, we vary the truncation level c by setting c as the 1p 0 quantile of the t distribution with the same tail index γ, and we refer to p 0 as the truncation threshold. This approach allows us to explore the effects of different levels of truncation on the performance of combination tests in practice.
The asymptotic validity of heavy-tailed transformation-based combination tests can be established based on Theorem 1. In particular, Liu & Xie (2020) demonstrated the asymptotic validity L. Gui et al.
of the Cauchy combination test. Extending this work, Fang et al. (2023) expanded these results to cover regularly varying distributions under additional constraints. However, both results are only limited to two-sided p-values, which are always positively dependent. In this section, we present a unified theory for the asymptotic validity of the heavy-tailed combination tests that accommodates both one-sided and two-sided p-values.
We first define combination tests applying the sum S n , directly inspired by Theorem 1.
Definition 4 (Combination test). Let F be the cumulative distribution function of a distribution in R -γ . The combination test approximates the tail probability pr(S n > x) by n F(x). Specifically, the combined p-value is defined as n F(S n ), and the corresponding decision function at the significance level α is
In addition to the sum S n , the widely accepted Cauchy and harmonic combination tests, as introduced by Liu et al. (2019) and Wilson (2019a), utilize the average M n and directly approximate the tail probability pr(M n > x) using F(x). Indeed, any regularly varying tailed distribution with tail index γ = 1 can be used to define a similar average-based combination test:
Definition 5 (Average-based combination test). Let F be the cumulative distribution function of a distribution in R -1 . The average-based combination test approximates the tail probability pr(M n > x) by F(x). Specifically, the combined p-value is defined as F(M n ) and the corresponding decision function at the significance level α is
More generally, one can define a weighted combination test, which includes both the tests defined in Definitions 4 and 5 as special cases. As noted in Liu & Xie (2020) and Fang et al. (2023), the weighted test can incorporate prior information on the importance of each base hypothesis to enhance power.
Definition 6 (Weighted combination test). Let F be the cumulative distribution function of a distribution in R -γ and let ì ω = (ω 1 , . . . , ω n ) ∈ R n + be a non-random weight vector associated with each hypothesis. Define the weighted sum as S n, ì ω = n i=1 ω i X i and let κ = n i=1 ω γ i where ω γ i is the γth power of ω i . Then the weighted combination test approximates the tail probability pr(S n, ì ω > x) by κ F(x). Specifically, the combined p-value is defined as κ F(S n, ì ω ) and the corresponding decision function at the significance level α is
Remark 2. The sum-based and average-based combination test in Definitions 4 and 5 are special cases of the weighted combination tests with uniform weights ω i = 1 or ω i = 1/n. Although the weighted combination test is not scale-free regarding the weights, empirical simulations suggest that the weight scaling has minimal practical impact.
The asymptotic validity of the combination tests in Liu & Xie (2020) and Fang et al. (2023) relies on pairwise bivariate normality of the test statistics {T i } n i=1 , ensuring pairwise quasiasymptotic independence as required by Theorem 1. Under the same assumption, we can establish the asymptotic validity for the combination tests defined in Definitions 4 to 6. Additionally, the asymptotic result is uniform in the nuisance parameters, particularly pairwise correlation ρ i j s, if we impose mild constraints on them.
Theorem 2. Assume that the test statistics {T i } n i=1 are pairwise normal with correlations ρ i j ∈ [-ρ 0 , ρ 0 ] (ρ 0 > 0) and are marginally following standard normal distributions under the global null. Then, the type-I error of the tests defined in Definitions 4 to 6 using two-sided p-values
where φ F comb is the test's decision function defined in (4) to (6). For the combination tests with onesided p-values
i=1 , the relationship (7) still holds with an additional assumption that the cumulative distribution function F(•) satisfies that F(x) F(-x) for sufficiently large x.
Remark 3. Our analysis considers fixed n. Prior work (Liu & Xie, 2020;Long et al., 2023) established the asymptotic validity of the Cauchy combination test as n → ∞, assuming n grows at a slower rate than the decay of α → 0. In addition, Vovk & Wang (2020) introduced an adjusted rejection threshold for the harmonic mean p-value to ensure validity as n → ∞, even under arbitrary dependence among the p-values.
The asymptotic validity of combination tests hinges on proving the pairwise asymptotic independence of the transformed statistics {X i } n i=1 . Theorem 2 provides a stronger asymptotic validity than previous studies (Liu & Xie, 2020;Fang et al., 2023) as uniform convergence is guaranteed over the set of correlation matrices. It also imposes minimal distributional requirements on F and further addresses one-sided p-values. Unlike two-sided p-values, which are always non-negatively correlated under bivariate normality as stated in Proposition S2, one-sided p-values can exhibit negative correlations. To establish the test's asymptotic validity, an additional constraint is required that F(x) F(-x) for sufficiently large x. This condition, met by all distributions in Table 1, ensures that the left tail is either absent or lighter than the right tail.
Theorem 2 requires no (i, j) pair has perfect correlaion. When ρ i j = ±1, though the transformed statistics X i and X j are no longer quasi-asymptotically independent, a weaker form of asymptotic validity still holds when the tail index γ 1, as stated below.
Corollary 2. Under assumptions of Theorem 2 while allowing ρ i j = ±1 for any (i, j) pairs, if additionally the tail index γ 1, then the combination tests defined in Definitions 4 to 6 using two-sided p-values are still asymptotically valid satisfying
where Σ = (ρ i j ) n×n is the correlation matrix of test statistics T i s, and
As a special case of Corollary 2, when ρ i j ≡ 1 for all pairs of test statistics, it holds that Corollary 3. Under conditions of Corollary 2, if ρ i j ≡ 1 for all (i, j) pairs, then the combination tests defined in Definitions 4 to 6 using either one-sided or two-sided p-values satisfies
In particular, when all weights are 1, the limit is n γ-1 .
L. Gui et al.
In this subsection, we explore the relationship between the heavy-tailed combination tests and the Bonferroni test. We begin by defining the weighted Bonferroni test, a generalization of the standard Bonferroni test that incorporates pre-chosen weights.
Definition 7 (Weighted Bonferroni test). Let P 1 , . . . , P n be the p-values and ì ω = (ω 1 , . . . , ω n ) ∈ R n + be a non-random weight vector satisfying n i=1 ω i = 1, then the weighted Bonferroni test at the significance level α has the decision function
It has been shown that the weighted Bonferroni test controls type-I error under any dependence structure (Genovese et al., 2006). The standard Bonferroni test is a special case where ω i = 1/n.
Given that the set
, the decision function of the Bonferroni test can be rewritten as
Thus, Corollary 1 implies that the type-I error of the Bonferroni test and the standard combination tests are asymptotically the same. Given this, we investigate whether the combination tests are indeed asymptotically equivalent to the Bonferroni test. We find out that the rejection regions of the weighted combination tests converge to those of a weighted Bonferroni test as α → 0.
Theorem 3. Assume that the test statistics {T i } n i=1 are pairwise normal with correlations ρ i j ∈ [-ρ 0 , ρ 0 ] and have a common marginal variance 1. Means of marginal normals are all finite. Then for two-sided p-values, when α → 0, any weighted heavy-tailed combination test defined in Definition 6 is asymptotically equivalent to a weighted Bonferroni test. Namely,
For one-sided p-values, the conclusion retains when further assuming that the cumulative distribution function F(•) satisfies that F(x) F(-x) for sufficiently large x.
Theorem 3 establishes the asymptotic equivalence between the combination tests and the Bonferroni test under any hypothesis configuration, provided that the test statistics are pairwise normal and not perfectly correlated. As the significance level α approaches zero, the rejection regions of both the combination tests and the Bonferroni test shrink, and the differences between these rejection regions diminish at a higher order. This equivalence does require that the test statistics are not perfectly correlated, so that they are quasi-asymptotically independent.
To provide an intuitive understanding of Theorem 3, Fig. 1 compares the rejection regions of various tests in the test statistics space for two-sided p-values with n = 2. The key takeaway is that the heavy-tailed nature of the transformation distribution yields nearly square rejection regions, which closely resemble those of the Bonferroni test as α decreases. In contrast, for combination tests relying on light-tailed distributions, such as Fisher's combination method, different rejection region shapes persist regardless of how small α becomes. Thus, in the asymptotic regime where these heavy-tailed combination tests are proven valid and when the individual test statistics are not perfectly correlated, there is no power gain over the Bonferroni test. 3. Empirical evaluations of the heavy-tailed combination tests under asymptotic independence 3.1
The theoretical results in Section 2 provide valuable insight into the heavy-tailed combination tests. However, it is unclear to what extent these asymptotic results align with their practical performance at finite significance levels. We aim to conduct an empirical evaluation of the tests' validity, focusing on commonly used finite significance levels.
For a comprehensive study, we vary the significance level α, number of hypotheses, tail heaviness and support of the distribution, and the level of dependence among the p-values. Specifically, we generate test statistics as z-values sampled from a multivariate normal distribution with mean ì µ = ì 0 n and covariance matrix Σ ρ . The covariance matrix Σ ρ ∈ R n×n has 1s on the diagonal and a common value ρ off the diagonal, representing varying degrees of dependence. We assess performance at three values of ρ, 0, 0.5, and 0.99, in line with no, moderate, and strong dependence. We calculate two-sided p-values from the z-values and conduct the combination tests based on different heavy-tailed distributions from four distribution families, the Student's t, Fréchet, Pareto, and inverse Gamma distributions. Each family has have a tunable tail index γ quantifying the tail heaviness, with a larger γ corresponding to a lighter tail. We vary this γ from 0.7 to 1.5 by 0.01 for all four distribution families. We also include the Bonferroni test and the Cauchy combination test as baselines. For significance levels, we adopt α = 0.05 and 5 × 10 -4 to account for different testing scenarios. The standard 0.05 is commonly used for a single global null hypothesis, while 5 × 10 -4 reflects the stricter threshold needed in genetic applications, where multiple testing adjustments lower the effective significance level for individual p-values. For the number of hypotheses, we consider n = 5 and 100. Each scenario is replicated 10 6 times to calculate the empirical type-I errors of the tests.
Figure 2 and S1 present the results for n = 5 and 100. When α = 0.05 and γ = 1, only the Cauchy combination test can strictly control empirical type-I error under independence, and no method achieves strict control when correlation ρ i j = 0.5. Smaller α improves error control and leads to a flatter curve across γ, consistent with the theoretical limit. Regarding the impact of tail heaviness on validity, differences between various distribution families diminish as α decreases, making the empirical validity of the tests primarily dependent on the tail index γ. Fig. 2: The type-I error of the combination test when n = 5 with different distributions: Cauchy (star point), inverse Gamma (blue), Fréchet (green), Pareto (purple), student t (red), left-truncated t with truncation threshold p 0 = 0.9 (dark orange), left-truncated t with truncation threshold p 0 = 0.7 (orange), left-truncated t with truncation treshold p 0 = 0.5 (light orange). The vertical axis represents the empirical type-I error, and the horizontal axis stands for the tail index γ.
error control is approximately achieved when γ 1. Distribution support also plays a role in type-I error control. Tests based on t distributions, which allow negative transformed statistics, outperform those using distributions with only positive support at α = 0.05. To examine this further, left-truncated t-distributions with different truncation thresholds p 0 = 0.5, 0.7 and 0.9 are adopted. As shown in Fig. 2 and S1, their empirical type-I errors fall between those of the original t distributions and other distribution families. This suggests that a wider support to the left of the real line tends to reduce the type-I error of the combination tests.
Additionally, we have investigated the type-I error control of the combination tests when the base p-values are negatively correlated by generating one-sided p-values. Results are shown in Table S1. We observe that when the p-values are negatively correlated, the Cauchy combination test can be even more conservative than the Bonferroni test due to its unbounded support. This undesired conservativeness can be mitigated by using a left-truncated t-distribution with a moderate truncation threshold. For more details, see Supplementary Section S1.
Theoretically, we have shown that the combination tests are asymptotically equivalent to the Bonferroni test for pairwise normal test statistics. Empirically, we aim to compare their power at finite significance levels and determine how small α needs to be for the asymptotic results to appear. Specifically, we evaluate significance levels α = 0.05 and 5 × 10 -4 while also approximating the asymptotic setting by letting α approach 0.
We start with assessing the power of the combination tests and the Bonferroni test at finite αs. Specifically, we define power as pr H 1,global (reject global null). We adopt the same simulation settings as in Section 3.1, generating one-sided p-values to obtain both positive and negative correlated p-values. We introduce both sparse and dense signals in the mean vector ì µ and consider three different numbers of hypotheses n = 5, 20, 100. The dense signals are generated as ì µ = ì µ n = (µ, µ, . . . , µ) ∈ R n . For sparse signals, we employ ì and ì µ = ( ì 0 95 , ì µ 5 ) ∈ R 100 as signal vectors. The parameter µ ranges from 0 to 6, ensuring that all testing methods can reach a power of 1, in increments of 0.5. For the covariance matrix Σ ρ , we select ρ = 0, 0.5, 0.9, 0.99 and also consider the negative correlation ρ = -0.2, to ensure the covariance matrix is positive definite, for n = 5. Each scenario is replicated 10 6 times to calculate the empirical power of the tests.
Figure 3 displays the maximum power difference between the combination tests using the Cauchy, truncated t 1 , Pareto, Fréchet, and Levy distributions, compared to the Bonferroni test when allowing µ to increase until all methods reach a power of 1. The truncation threshold for the t 1 distribution is set at p 0 = 0.9. The Cauchy, truncated t 1 , Fréchet, and Pareto distributions share a tail index γ = 1, whereas the Levy distribution has a tail index of 0.5, resulting in a smaller power difference compared to the Bonferroni test.
Our findings reveal that combination tests can achieve higher power at finite significance levels, particularly in situations where signals are dense. This remains the case for the Cauchy combination test even when p-values are negatively correlated, a setting in which it tends to be overly conservative. This likely stems from the nature of the combination test, which synthesizes L. Gui et al. signals from multiple sources rather than relying on a single dominant signal. These results suggest that the onset of asymptotic equivalence may occur at much smaller values of α compared to that for asymptotic validity, especially when signals are dense.
To further investigate the asymptotic equivalence between the combination tests and the Bonferroni test, we examine how the size of their non-overlapping rejection regions evolves as α approaches 0. Using the same settings as earlier in this section with n = 5 and ρ = 0.5, we fix the signal level µ = 2 to ensure the power difference between the two tests is not negligible.
We consider three mean vectors: ì µ = ì 0 5 (global null),( ì 0 4 , 2) (sparse signal), ì 2 5 (dense signal), allowing us to compare their performance under different scenarios. As shown in Fig. 4, the difference, quantified by the probability ratio between the overlapping rejection region and individual rejection regions, converges to zero as α decreases, being consistent with the asymptotic equivalence established in Theorem 3.
Since the Bonferroni test is known to suffer under strong dependence, we also compare the combination tests against the adjusted Bonferroni method, minP. Specifically, we calibrate the cutoff for min(p 1 , . . . , p n ) using Monte Carlo sampling from the true data-generating model to ensure the actual type I error matches the nominal level α (Table S4). We replicate the simulation settings from Fig. 3, replacing Bonferroni with minP as the baseline. As shown in Fig. S2, combination tests outperform minP when signals are dense and test statistics are weakly correlated, consistent with findings in Liu & Xie (2020). However, minP relies on knowledge of the dependence structure among p-values, limiting its practicality in many applications and making it computationally intensive.
Although heavy-tailed combination tests are typically employed when p-values have unknown dependence, they do not guarantee control of the type-I error under arbitrary dependence structures, even asymptotically. One key assumption for ensuring asymptotic type-I error control in Section 2.3 is the requirement of quasi-asymptotic independence, which can be restrictive in practice. For instance, when the sample size is small, test statistics are likely to follow a t-distribution rather than a normal distribution. Additionally, even when the sample size is large, it can still be challenging to ensure that two dependent test statistics are pairwise normal.
The strength of asymptotic dependence between any two variables (X 1 , X 2 ) with the same marginal distribution F can be quantified by the upper tail dependence coefficient (Joe, 1997)
As discussed earlier, if X 1 and X 2 are bivariate normal and are not perfectly correlated, they are quasi-asymptotic independent, and hence λ = 0. However, many dependent variables do not satisfy quasi-asymptotic independence. For instance, for bivariate t-distributed variables (T 1 ,T 2 ) with degree of freedom ν, variances 1 and correlation ρ, their tail dependent coefficient (Demarta & McNeil, 2005) is
, where t ν (•) is the cumulative distribution function of the t distribution. As a result, T 1 and T 2 are never quasi-asymptotically independent, even when ρ = 0, due to shared covariance estimation.
To understand the sensitivity of the combination tests to violations of quasi-asymptotic independence, we generate test statistics (T 1 , . . . ,T n ) from a multivariate t distribution t ν ( ì 0 n , Σ ρ ), where Σ ρ is defined in Section 3. We choose an extreme degree of freedom ν = 2 and set the correlation ρ to 0, 0.5, 0.9, and 0.99, resulting in tail dependence indices ranging from 0.18 to 0.91. All base p-values are one-sided and derived from the test statistics.
Table 2 and S3 compare the empirical type-I errors of different combination tests at the significance level α = 0.05 and 5 × 10 -4 , and for n = 5 and n = 100. Surprisingly, the results indicate that type-I errors remain well-controlled regardless of the tail dependence coefficient, demonstrating the robustness of the combination tests to violations of the pairwise normal assumption for the test statistics.
Furthermore, Table 2 and S3 suggest that the Bonferroni test tends to be exceedingly conservative when the dependence coefficient λ > 0, especially when both n and λ are large. In contrast, the combination tests based on heavy-tailed distributions with γ = 1 consistently maintain a type-I error rate close to the specified significance level. Thus, we hypothesize that when test statistics are quasi-asymptotically dependent, the combination tests with a tail index γ 1 are still asymptotically valid when α → 0, but they will not be asymptotically equivalent to the Bonferroni test. While the Bonferroni test can exhibit excessive conservatism, the combination tests with γ = 1 display neither conservatism nor inflation in their type-I error rates. For example, as discussed in Corollary 2, in situations where test statistics are perfectly correlated with ρ = 1, resulting in a tail dependence coefficient of λ = 1, the combination tests with γ = 1 maintain an asymptotic type-I error of α, whereas the true type-I error of the Bonferroni test is only α/n.
We further investigate the power gain of the combination test over the Bonferroni test when test statistics follow a multivariate t-distribution. Compared to the power comparison in Section 3.2, we replace the distribution of the test statistics from a multivariate normal distribution to a multivariate t distribution with ν = 2, while keeping all other settings the same. Figure 5 displays the maximum power gain of each combination test over the Bonferroni test as the power of both tests grows from 0 to 1 as signal strength increases. Compared to the subtle power improvement we observed for multivariate normally distributed test statistics in Fig. 3, the maximum power difference for multivariate t-distributed test statistics can be as large as 1 even when signals are sparse. The power difference does not diminish even when the significance level decreases from 0.05 to 5 × 10 -4 .
Table 2: Type-I error control of the combination tests when test statistics follow a multivariate t-distribution when n = 5. Values in parentheses are the corresponding standard errors. For the Fréchet and Pareto distributions, the tail index γ = 1. For truncated t 1 , the truncation threshold p 0 = 0.9
ρ Cauchy Pareto Truncated t 1 Fréchet Levy Bonferroni Fisher 0 0.18 2.90E-02 5.30E-02 4.73E-02 5.17E-02 3.89E-02 3.56E-02 6.26E-02 (1.68E-04) (2.24E-04) (2.21E-04) (1.93E-04) (2.12E-04) (1.85E-04) (2.42E-04) 0.5 0.39 4.48E-02 5.24E-02 5.01E-02 5.13E-02 3.19E-02 2.65E-02 1.16E-01 (2.07E-04) (2.23E-04) (2.18E-04) (2.21E-02) (1.76E-04) (1.61E-04) (3.20E-04) 0.9 0.72 5.00E-02 5.09E-02 5.06E-02 4.99E-02 2.50E-02 1.67E-02 1.51E-01 (2.18E-04) (2.20E-04) (2.19E-04) (2.18E-04) (1.56E-04) (1.28E-04) (3.58E-04) 0.99 0.91 5.02E-02 5.03E-02 5.02E-02 4.92E-02 2.27E-02 1.19E-02 1.59E-01 (2.18E-04) (2.18E-04) (2.18E-04) (2.16E-04) (1.49E-04) (1.09E-04) (3.66E-04) (b) α = 5 × 10 -4 ρ λ 2,ρ Cauchy Pareto Truncated t 1 Fréchet Levy Bonferroni Fisher 0 0.18 2.48E-04 4.57E-04 4.57E-04 4.57E-04 3.49E-04 3.18E-04 2.17E-02 (1.57E-05) (2.14E-05) (2.14E-05) (2.14E-05) (1.88E-05) (1.78E-05) (1.46E-04) 0.5 0.39 3.94E-04 4.65E-04 4.65E-04 4.65E-04 3.08E-04 2.67E-04 2.63E-02 (1.98E-05) (2.16E-05) (2.16E-05) (2.16E-05) (1.75E-05) (1.63E-05) (1.60E-04) 0.9 0.72 5.20E-04 5.28E-04 5.28E-04 5.28E-04 2.37E-04 1.65E-04 3.82E-02 (2.28E-05) (2.30E-05) (2.30E-05) (2.30E-05) (1.54E-05) (1.28E-05) (1.92E-04) 0.99 0.91 5.24E-04 5.24E-04 5.24E-04 5.24E-04 2.22E-04 1.16E-04 4.25E-02 (2.29E-05) (2.29E-05) (2.29E-05) (2.29E-05) (1.49E-05) (1.08E-05) (2.02E-04)
These findings indicate a potential power advantage of the combination tests over the Bonferroni test, even in the asymptotic regime where α → 0, when test statistics are pairwise asymptotically dependent. Our empirical results indicate that, unlike in the case of asymptotic independence, combination tests can remain asymptotically valid while achieving a nontrivial power improvement over the Bonferroni test under asymptotic dependence. This highlights the potential of asymptotic dependence as a valuable framework for advancing both the theoretical and practical understanding of combination tests.
Circadian rhythms, which are oscillations of behavior, physiology, and metabolism, are observed in almost all living organisms (Pittendrigh, 1960). Recent advances in omics technologies, such as microarray and next-generation sequencing, provide powerful platforms for identifying circadian genes that encode molecular clocks crucial for health and diseases (Rijo-Ferreira & Takahashi, 2019). In this case study, we focus on a gene expression dataset obtained from mouse liver samples, collected every hour across 48 different circadian time points, denoted as CT points, ranging from CT18 to CT65, under complete darkness conditions (Hughes et al., 2009). At each time point, the expression levels of approximately 13,000 mouse genes were profiled by microarray. The objective of this case study is to identify genes that exhibit significant oscillatory behavior by aggregating results across all measured time points. One of the most widely used methods is JTK_CYCLE (Hughes et al., 2010). JTK_CYCLE determines whether a gene exhibits significant cyclic behavior by performing a Kendall's tau test. It compares the observed gene expression measurements across 48 time points to expected patterns with specific phases and periods using a rank-based correlation test. This process involves testing 216 combinations of phase and period, resulting in 216 correlated base p-values for each gene. By default, JTK_CYCLE combines these p-values using the Bonferroni test, though this approach has been shown to lack power in benchmarking studies (Mei et al., 2021).
In place of the Bonferroni test, we use the heavy-tailed combination tests to aggregate the 216 correlated p-values for each gene. For comparison, we also include Fisher's method. To assess the performance of different tests, we utilize a set of the 60 positive control, i.e., cyclic genes, and 61 negative control, i.e., non-cyclic genes, from Wu et al. (2014) as ground truth. Figure 6 displays the box plots of the combined p-values for the positive and negative controls. Compared to the Bonferroni method, the combined p-values from heavy-tailed combination tests have higher detection power of the true signals, while avoiding false positives in negative controls compared to Fisher's method. Truncated" refers to using the t 1 distribution with truncation threshold p 0 = 0.9. For Fréchet and Pareto distributions, the tail index is set to γ = 1.
In the second real data analysis, similar to Liu et al. (2019), we combine correlated p-values to identify genes that are significantly associated with diseases in genome-wide association studies, referred to as GWAS for brevity. A gene of interest may contain multiple single-nucleotide polymorphisms, referred to as SNPs, each tested individually against the trait, e.g., disease status, using a simple regression framework, resulting in SNP-level p-values. Then, p-values from the SNPs within the same gene region are further combined via a gene-level test. SNPs that are close to each other on the genome are highly correlated due to linkage disequilibrium, leading to highly correlated SNP-level p-values for the same gene. Several methods have been developed for gene-level association testing, such as EPIC (Wang et al., 2022) andMAGMA (de Leeuw et al., 2015), which account for SNP-SNP correlations within the same gene. However, these methods can be computationally intensive. For example, deriving gene-level test statistics in these methods often requires inverting large covariance matrices.
In this analysis, we apply heavy-tailed combination tests to test for each gene's association with schizophrenia, referred to as SCZ (Ripke et al., 2013). To adjust for multiple testing errors, we apply the Benjamini-Hochberg procedure (Benjamini & Hochberg, 1995) on the gene-level combined p-values to control the false discovery rate, referred to as FDR for simplicity. Figure 7 shows the number of overlapping genes rejected by each method compared when FDR is controlled at 0.05 and 0.2. As illustrated, the number of genes detected by the combination tests is comparable to or even higher than those identified by Epic and Magma. Notably, the combination tests are highly computationally efficient, completing analyses almost instantly compared to domainspecific methods that require modeling the correlation structure.
Compared to the Bonferroni test, the combination tests identify 25% more significant genes, even at a low nominal FDR level of α = 0.05. Figure S4 summarizes the number of SNPs for each gene, showing that most genes have fewer than 100 SNPs, suggesting that the significant power gain is not due to combining an excessively large number of p-values, which could lead to inflation of type-I errors. Figure S5 displays that even when focusing solely on genes with 50 or "Truncated" refers to the truncated t 1 distribution with truncation threshold p 0 = 0.9. For Fréchet and Pareto distributions, the tail index is set to γ = 1.
fewer SNPs, the combination tests still identify substantially more genes than the Bonferroni test.
Compared to the simulation results, the substantial power gain in this real data analysis likely results from the violations of quasi-asymptotic independence of the SNP-level p-values. To evaluate whether the additional genes detected by the heavy-tailed combination tests are biologically meaningful, we analyze the set of 939 genes detected at the FDR level α = 0.2 by the Cauchy, truncated t 1 , Fréchet, or Pareto combination tests but not by the Bonferroni test. We conduct a gene-set enrichment analysis using DAVID (Sherman et al., 2022). Results are shown in Fig. S6. The top two significantly enriched gene ontology terms are "regulation of ion transmembrane transport" and "chemical synaptic transmission", both of which have been reported and confirmed by independent studies (Favalli et al., 2012;Liu et al., 2022). These findings underscore the enhanced statistical power of transformation tests compared to the Bonferroni test in practical genetic applications.
In this section, we examine the extensions and limitations of our results and discuss related literature. The asymptotic validity of the heavy-tailed combination tests can be generalized to cases where p-values are only valid, i.e., they satisfy pr(p α) α, as long as the p-values are pairwise quasi-asymptotically independent. Though the transformed test statistics X i derived from these valid p-values may lack a regularly varying tailed distribution, the combination tests should maintain control over type-I errors. Intuitively, this is because we can always construct uniformly distributed variables that are stochastically smaller than these valid p-values.
In the context of multiple testing, the combination tests can be applied within a closed testing procedure to identify individual non-null hypotheses. In Supplementary Section S2, we provide a shortcut algorithm for applying closed testing with combination tests. Goeman et al. (2019a) L. Gui et al. demonstrated that, as n → +∞, the closed testing procedure using harmonic mean p-values is significantly more powerful than the one based on Bonferroni corrections. However, for finite n and when the family-wise error rate approaches zero, the equivalence between combination tests and the Bonferroni test may extend to their respective closed testing procedures.
To balance validity and power, we recommend using a truncated t 1 distribution with truncation threshold p 0 = 0.9, based on the empirical results. This definition differs slightly from the truncated Cauchy distribution proposed by Fang et al. (2023), which assigns a point mass at the truncation threshold rather than rescaling the distribution. Notably, the half-Cauchy distribution in Long et al. (2023) is a special case of our definition with p 0 = 0.5. While our focus is on establishing the asymptotic validity of combination tests using the truncated t 1 distribution under an unknown dependence structure, both Fang et al. (2023) and Long et al. (2023) have also provided adjustments that ensure exact validity when p-values are independent.
While our results establish the asymptotic validity of the heavy-tailed combination tests under quasi-asymptotically independent test statistics, the combination tests can exhibit noticeable inflation in type-I error rates under arbitrary dependence and finite α. Exact control over type-I errors may be achieved with additional adjustments. For the harmonic mean pvalues, Vovk & Wang (2020) demonstrated that it is valid under arbitrary dependence when scaled by a factor a n = (y n + n) 2 /(ny n + n) where y n is the unique solution to the equation y 2 n = n {(y n + 1) log(y n + 1)y n }. This factor asymptotically approaches log n when n increases and the test can be further improved through randomization techniques (Gasparin et al., 2024).
Other studies, such as Wilson (2019a) and subsequent works (Held, 2019;Wilson, 2019b;Goeman et al., 2019b) have provided empirically calibrated thresholds for harmonic mean p-value. Additionally, Chen et al. ( 2024) establish an adjustment of the harmonic mean p-value to guarantee its validity when the individual p-values follow a Clayton copula.
Numerous alternative methods for combining dependent p-values exist, each with distinct tradeoffs. Some approaches, such as those by Goeman et al. (2004) and Edelmann et al. (2020), model specific dependence structures, which can be powerful but require strong model assumptions and can be computationally intensive. Other methods, like those by Hommel (1983) and Vovk & Wang (2020), guarantee type-I error control under arbitrary dependence. However, as discussed in earlier studies (Fang et al., 2023;Chen et al., 2023), combination methods with proven validity guarantees under arbitrary dependence may have limited power in practical applications.
The R package facilitating the implementation of heavy-tailed combination tests is accessible at https://github.com/gl-ybnbxb/heavytailcombtest. The code to reproduce figures and tables is available at https://github.com/gl-ybnbxb/ combination-test-reproduce-code. Time-series circadian gene expression data of mouse liver is downloaded from the Gene Expression Omnibus (GEO) database with accession number GSE11923. GWAS summary statistics of schizophrenia (SCZ) is downloaded from the Psychiatric Genomics Consortium at https://pgc.unc.edu/for-researchers/download-results/.
This work was supported by the National Science Foundation under grants DMS-2113646 and DMS-2238656, and by the National Institute of Health under grant R35 GM138342. We thank Nancy R. Zhang and Jian Ding for their helpful discussions. We acknowledge the University of