Robust High-dimensional Tuning Free Multiple Testing^††thanks: Supported by NSF grants DMS-2210833, DMS-2053832, DMS-2052926 and ONR grant N00014-22-1-2340

Large-scale, high-dimensional data with rich structures have been widely collected in almost all scientific disciplines and humanities, thanks to the advancements of modern technologies. Massive developments have been made in statistics over the past two decades on extracting valuable information from these high dimensional data; see Bühlmann and Van De Geer (2011); Hastie et al. (2009, 2015); Wainwright (2019); Fan et al. (2020a); Chen et al. (2021) for a detailed account and references therein.

Despite convenience for theoretical analysis, the sub-Gaussian tails condition is not realistic in many applications that involve high-dimensional variables. For instance, it is well known that heavy-tailed distributions is a stylized feature for financial returns and macroeconomic variables (Cont, 2001; Stock and Watson, 2002; Fan and Yao, 2017; Fan et al., 2021). Therefore, tools designed for sub-Gaussian data can lead to erroneous scientific conclusions. Asking thousands of gene expressions to have all sub-Gaussian tails is a mathematical dream, not a reality that data scientists face. For example, comparing gene expression profiles between various cell sub-populations, especially after treatments and therapies, is an essential statistical task (Nagalakshmi et al., 2008; Shendure and Ji, 2008; Wang et al., 2009; Li et al., 2012; Li and Tibshirani, 2013; Gupta et al., 2014; Finotello and Di Camillo, 2015). However, it is unrealistic to hope all thousands of gene expressions have sub-Gaussian distributions: outliers in non-sub-Gaussian distributions can have a significant impact on nonrobust procedures and lead to many false positives and negatives (Gupta et al., 2014; Wang et al., 2015). The situation also arises for inferences using functional magnetic resonance imaging (fMRI) data since the data do not conform to the assumed Gaussian distribution (Eklund et al., 2016). These practical challenges demand for developing efficient and reliable robust inference methods.

Recently, robust statistical methods have gained popularity as a mean of resolving outliers and heavy-tailed noises. Many preceding arts have taken a significant stride toward effective statistical estimation under heavy-tailed distributions. For instance, aiming at dealing with the heavy-tailed noise contamination, the Huber regression is proposed (Huber, 1973), and subsequent publications along these lines include Yohai and Maronna (1979); Mammen (1989); He and Shao (1996, 2000), where the asymptotic properties of the Huber estimator have been thoroughly investigated. From a non-asymptotic perspective, the Huber-type estimator was recently revisited by Sun et al. (2020), in which the authors propose an adaptive Huber regression method and establish its non-asymptotic deviation bounds by only requiring finite $(1+\delta)$-th moment of the noise with any $\delta>0$. Moreover, using a similar idea for making the correspondent $M$-estimator insensitive to extreme values, Catoni (2012) developed a novel approach through minimizing a robust empirical loss. It is demonstrated that the estimator has exponential concentration around the true mean and enjoys the same statistical rate as the sample average for sub-Gaussian distributions when the population only has a bounded second moment. Brownlees et al. (2015) further investigates empirical risk minimization based on the robust estimator proposed in Catoni (2012). Additionally, the so-called median of means strategy, which can be traced back to Nemirovsky and Yudin (1983), is another successful method for handling heavy-tailed distributions. By only requiring bounded second moment, it achieves the sub-Gaussian type of concentration around the population mean parameter. Minsker (2015) and Hsu and Sabato (2016) further generalize this idea to multivariate cases. Moreover, there also exists a series of works that focus on solving the issue caused by heavy-tailed noises using quantile-based robust estimation; see Arcones (1995); Koenker and Hallock (2001); Belloni and Chernozhukov (2011); Fan et al. (2014); Zheng et al. (2015) for more details. Furthermore, recently, in Fan et al. (2021); Yang et al. (2017); Fan et al. (2022c), under heavy-tailed contamination, they proposed a novel principle by simply truncating or shrinking the response variables appropriately to achieve sub-Gaussian rates, and they only require bounded second moment of the measurements. Additionally, the aforementioned methodologies can also be applied to a wide range of problems, such as matrix sensing, matrix completion, robust PCA, factor analysis, and neural networks. For interested readers, we refer to Minsker (2015); Hsu and Sabato (2016); Fan et al. (2017); Loh (2017); Minsker (2018); Goldstein et al. (2018); Wang et al. (2020); Fan et al. (2022b); Wang and Fan (2022); Fan et al. (2022a) for more details.

While many effective solutions have been developed to address the problem of heavy-tailedness, these solutions still have some potential shortcomings. In specific, the developments call for the second moments to be bounded and are primarily based on shrinkage data, Huber-type of loss, median of means (Huber, 1973; Nemirovsky and Yudin, 1983; Catoni, 2012; Fan et al., 2021). More critically, Huberization, Winsorization and sample splitting introduce additional tuning parameters and these tuning paprameters should be variable-dependent that makes large-scale applications difficult and damages the fidelity of empirical results . Although, the quantile estimators (Koenker and Hallock, 2001) such as the median and the Hodge-Lehmann (HL) estimator (Hodges and Lehmann, 1963) can eliminate the restriction on moment conditions and tuning parameters selection, the empirical median is often less efficient and requires stronger distribution assumptions. In addition, the existing literature on the HL estimator focuses mainly on low-dimensional asymptotic analysis and can not be applied to large-scale inferences.

In this paper, we revisit the celebrated HL estimator (Hodges and Lehmann, 1963) and conduct non-asymptotic and large-scale theoretical studies for both one-sample and two-sample problems. For one-sample location estimation in the univariate case, we let $X_{1},\ldots,X_{n}\in\mathbb{R}$ be independent and identically distributed (i.i.d.) random variables with

where $\theta$ represents the location parameter of interest and $\xi_{1},\ldots,\xi_{n}$ are i.i.d. random variables drawn from some unknown distribution. In this scenario, the HL estimator of $\theta$ is defined by

By assuming the pseudomedian of $\xi_{1}$ to be zero, we derive non-asymptotic Bahadur representation of $\widehat{\theta}.$ To the best of our knowledge, this is the first study of its kind on the non-asymptotic expansion of the HL estimator. From there, we also establish the Berry-Essen bound and moderate deviation for $\widehat{\theta}$ in the widest range. Furthermore, as there are multiple unknowns in the asymptotic distribution of $\widehat{\theta}$, including the density function of $\xi_{1}$ and the unknown location parameter $\theta$, we then propose a weighted bootstrap approach to construct confidence intervals for $\theta$ based on data. These results and methods are essential for the large-scale inference.

In addition to the study of one-sample problem, two-sample location shift problems arise frequently in many scientific studies, including choosing genes that are expressed differently in normal and injured spinal cord, determining the effects of treatment between treated and control groups, finding change points, etc. To this end, we let $Y_{1},\ldots,Y_{m}\in\mathbb{R}$ be another independent sample of i.i.d. random variables satisfying

The primary goal is to conduct statistical inference for $\Theta=\theta-\theta^{\circ}$. In the sequel, following Hodges and Lehmann (1963), the two-sample HL estimator for $\Theta$ is given by

Instead of assuming the noises are generated from the same distribution (Hodges and Lehmann, 1963), we only require median$(\xi_{1}-\varepsilon_{1})=0$, which is more general and allows random noises to have different distributions. In a similar vein, we establish the non-asymptotic expansion of $\widehat{\Theta}$, investigate its asymptotic distributions, and calculate the confidence interval via bootstrap techniques. Again, the techniques and results developed can be applied to large-scale multiple testing problems.

There is a rich literature on large-scale multiple testing problems for location parameters (Benjamini and Hochberg, 1995; Storey, 2002, 2003; Genovese and Wasserman, 2004; Ferreira and Zwinderman, 2006; Chi, 2007; Blanchard and Roquain, 2009). However, most of these works assume the noise distributions sub-Gaussian. Moving away for sub-Gaussian assumptions, Fan et al. (2019) propose estimating mean vector via minimizing the Huber type loss and perform the false discovery proportion (FDP) control. However, leveraging Huber type estimators necessitates moment limits and introduces tuning parameters, making it hard to be applied in large-scale inference, as the tuning parameters should ideally be variable-dependent. Additionally, while the HL estimator enjoys tuning-free and moment-free qualities in the univariate setting, its behavior in high dimensions is largely unknown. To this end, in this research, we further expand the HL estimator to high-dimensional regimes prompted by the lack of tuning free large-scale multiple testing problems for heavy-tailed distributions.

In specific, we assume $\bm{X}_{i}=\bm{\theta}+\bm{\xi}_{i}$, $i\in[n]$ and $\bm{Y}_{j}=\bm{\theta}^{\circ}+\bm{\varepsilon}_{j},j\in[m]$, where $\bm{\theta}=(\theta_{1},\ldots,\theta_{p})^{\top}$ and $\bm{\theta}^{\circ}=(\theta_{1}^{\circ},\ldots,\theta_{p}^{\circ})^{\top}$ are p-dimensional vectors of unknown parameters and random noises $\bm{\xi}_{i},i\in[n],\bm{\varepsilon}_{j},j\in[m]$. Let $\bm{\Theta}=\bm{\theta}-\bm{\theta}^{\circ}$. For both one- and two-sample problems, we propose a carefully constructed Gaussian multiplier bootstrap to test global null hypotheses

by extending the HL estimator to high-dimensional regimes. When the null hypothesis above is rejected, we then perform multiple testing, allowing weakly dependent measurements, and efficiently control the FDP. Compared with existing literature (Liu and Shao, 2014; Fan et al., 2019), our procedures do not involve any tuning parameters and moment conditions for testing global null and large-scale multiple testing. These theoretical finds are further supported by exhaustive numerical studies.

The main contributions of the paper can be summarized as follows:

• •

  The existing studies on the HL estimator mainly focus on its asymptotic behavior, which is too weak for high-dimensional applications. In practice, however, it is crucial to understand the HL estimator’s performance under finite sample, especially in high-dimensional and large-scale experiments. For this purpose, we first derive the non-asymptotic expansions of the HL estimators for both one-sample and two-sample problems.

• •

  With the non-asymptotic expansions of the HL estimators, for both one- and two-sample problems, we derive its Berry-Essen type bounds and Cramér type moderate deviations, with the widest range. To deal with unknown components in the distribution, we further develop the weighted bootstrap to build data-driven confidence intervals. In addition, we also furnish the non-asymptotic analysis of the bootstrap estimator.

• •

  Existing work on large-scale testing with heavy-tail errors typically involves additional tuning parameters and the moment conditions. In order to address these issues, we generalize the HL estimator to large-scale studies and propose tuning-free and moment-free high-dimensional testing procedures. Additionally, we develop bootstrap methods for calculating critical values for large-scale applications. We show that the resulting false discovery proportion is well controlled.

This section is devoted to studying the non-asymptotic expansions of the Hodges-Lehmann estimator and conducting statistical estimation and inference for population location shift parameters. For both one-sample and two-sample problems, the theoretical properties, which are needed for large-scale inferences, are presented in the following sections.

In this section, we propose a weighted bootstrap method to construct confidence intervals for $\theta$ and $\Theta$, rather than directly estimating those involved unknown terms in asymptotic variances using the brute force methods. The reason is that the direct estimation approach always necessitates the additional selection of tuning parameters and imposes moment conditions. Additionally, the bootstrap calibration performs admirably with finite samples, particularly when the sample size is modest. Therefore, in the sections that follow, we outline the bootstrap procedures for both the one- and two-sample problems.

With the advancement of technology, large-scale, high-dimensional data have been extensively collected over the past two decades in a variety of fields such as medicine, biology, genetics, earth science, and finance. In large-scale regimes, there are inevitably heavy-tailed noises and it is crucial to develop robust statistical inference procedures. Yet, existing research that infers location shifts via Huber-type estimates calls for variable-dependent tuning parameters and moment limitations (Fan et al., 2019; Sun et al., 2020). This type of technique is hard to apply efficiently and faithfully to large-scale inferences due to the choices of tuning values. Moment constraints also exclude a large number of heavy-tailed distributions. To remedy the issues, this section focuses on extending the HL estimation to high dimensions and developing tuning-free and moment-free high-dimensional multiple testing procedures.