Accurate and Ultra-Efficient p-Value Calculation for Higher Criticism Tests

Abstract

In modern data science, higher criticism (HC) method is effective for detecting rare and weak signals. The computation, however, has long been an issue when the number of p-values combined ( $K$ ) and/or the number of repeated HC tests ( $N$ ) are large. Some computing methods have been developed, but they all have significant shortcomings, especially when a stringent significance level is required. In this paper, we propose an accurate and highly efficient computing strategy for four variations of HC. Specifically, we propose an unbiased cross-entropy-based importance sampling method ( ${\text{IS}}_{CE}$ ) to benchmark all existing computing methods, and develop a modified SetTest method (MST) that resolves numerical issues of the existing SetTest approach. We further develop an ultra-fast approach (UFI) combining pre-calculated statistical tables and cubic spline interpolation. Finally, following extensive simulations, we provide a computing strategy integrating MST, UFI and other existing methods with R package "HCp" for virtually any $K$ and small p-values ( $\sim {10}^{-20}$ ). The method is applied to a COVID-19 disease surveillance example for spatio-temporal outbreak detection from case numbers of 804 days in 3,342 counties in the United States. Results confirm viability of the computing strategy for large-scale inferences. Supplementary materials for this article are available online.

Keywords: analytical approximation; asymptotic rare and weak model; higher criticism; importance sampling; p-value computation.

Fig. 1

(A) Left panel: the MST is recommended for $T_{HC}^{R_T}$ of arbitrary truncation with $K\le 1000$, while a hybrid of MST and Li-Siegmund is recommended for targeted p-value $\ge {10}^{-3}$ and p-value $<{10}^{-3}$ respectively when $K$ is above 1000. Right panel: a hybrid strategy combining the MST for $K\le 100$, naive Monte Carlo and the Li-Siegmund for targeted p-value $\ge {10}^{-3}$ and p-value $<{10}^{-3}$ respectively is recommended for $T_{HC}^{R_{TM}}$ of arbitrary truncation under $K>100$. (B) In our “HCp” R package, statistical tables are pre-calculated when $K\le 2000$ and target p-value $>{10}^{-14}$ for the specific HC tests: (i) $T_{HC}^{R_F}$ and (ii) $T_{HC}^{R_T}$ with $k_0=1,k_1=\lfloor K∕2\rfloor$ (left panel), as well as (iii) $T_{HC}^{R_M}$ and (iv) $T_{HC}^{R_{TM}}$ With $k_0=1,k_1=\lfloor K∕2\rfloor$ of target p-value $>{10}^{-12}$ (right panel), where the UFI method is applicable and prefered.

Fig. 2

(A) illustrates the accuracy of the Barnett-Lin, SetTest, and MST method in computing small p-values of $T_{HC}^{R_F}$ with $K=500$ benchmarked by ${\text{IS}}_{CE}$ method (the methods started with * in the legend is our proposed methods). The x-axis is the $\log$ value of $T_{HC}^{R_F}$ statistic, and the y-axis is the corresponding $-{\log }_{10}$ p-value. Similarly, (B) evaluates the accuracy of the SetTest and MST method benchmarked by the ${\text{IS}}_{CE}$ method for case (ii) $T_{HC}^{R_T}$ of $K=500$. (C) and (D) shows the accuracy of the MST benchmarked by the ${\text{IS}}_{CE}$ method for $T_{HC}^{R_M}$ and case (iv) $T_{HC}^{R_{TM}}$ respectively under $K=100$. The sampling size for ${\text{IS}}_{CE}$ is $M={10}^4$.

Fig. 3

(A) and (C) illustrate the performance of the Li-Siegmund, MST, and naive Monte Carlo method with sampling size $M={10}^6$ in estimating large p-values for $T_{HC}^{R_M}$ and case (iv) $T_{HC}^{R_{TM}}$ of $K=100$ respectively. The error bar of the naive Monte Carlo estimate represents the square root of mean square error at the original scale over 50 replications. The y-axis is the p-value at original scale. (B) and (D) illustrate the performance of Li-Siegmund, MST, and ${\text{IS}}_{CE}$ method with sampling size $M={10}^6$ in estimating small p-values for $T_{HC}^{R_M}$ and case (iv) $T_{HC}^{R_{TM}}$ of $K=100$ respectively. The error bar of the ${\text{IS}}_{CE}$ estimate is the square root of mean square error at ${\log }_{10}$ scale over 50 replications and the y-axis represents the $-{\log }_{10}$ p-value.

Fig. 4

(A) compares the time (in second) consumed of computing the p-values of 100 different $T_{HC}^{R_F}$ tests by the Barnett-Lin, SetTest, Li-Siegmund, ${\text{IS}}_{CE}$, MST and UFI method with respect to varying $K$. The formula of the fitted polynomial curve for each method is labeled in the figure. (B) removes the Barnett-Lin method and compares the other five methods. (C) illustrates the consistent p-value estimate by the UFI method with the analytic truth by the MST method for $T_{HC}^{R_F}$ of various $K$. The black straight line $y=x$ is the reference representing that the UFI estimates are exactly the same as the analytic truth.

Fig. 5

The maps show the significance level ($-{\log }_{10}$ p-value) of the COVID-19 outbreak testing by the UFI method for each of the 3,342 counties in the US at three time periods: 03/22/2021-04/20/2021 (top); 11/17/2021-12/16/2021 (middle); 02/15/2022-03/16/2022 (bottom). There are more significant counties outbreaking COVID-19 at the end of 2021 and during the period from mid-February to mid-March in 2022.