A New Perspective on Robust $M$-Estimation: Finite Sample Theory and Applications to Dependence-Adjusted Multiple Testing

Wen-Xin Zhou, Koushiki Bose, Jianqing Fan, Han Liu

Introduction

High dimensional data are often automatically collected with low quality. For each feature, the samples drawn from a moderate-tailed distribution may comprise one or two very large outliers in the measurements. When dealing with thousands or tens of thousands of features simultaneously, the chance of including a fair amount of outliers is high. Therefore, the development of robust procedures is arguably even more important for high dimensional problems. In this paper, we develop a finite sample theory for robust $M$ -estimation from a new perspective. Such a finite sample theory is motivated by contemporary statistical problems of simultaneously testing many hypotheses. In these problems, the goal is either to control the false discovery rate (FDR)/false discovery proportion (FDP), or to control the familywise error rate (FWER).

The main contributions of this paper are described and summarized in the following two subsections.

The first contribution of this paper is to develop a new finite sample theory for the Huber estimator. Recall the Huber loss [Huber (1964)]:

For every $\tau>0$ , by definition $\widehat{\bm{\theta}}$ is an $M$ -estimator of

which typically differs from the target parameter

Note that the total error $\|\widehat{\bm{\theta}}-\bm{\theta}^{*}\|$ can be divided into two parts:

where $R_{n}(\tau)$ is a finite sample error bound which characterizes the accuracy of such a linear approximation, and $d$ is the number of covariates that may grow with $n$ . We refer to Theorem 2.1 for a rigorous description of the result (1.4), where we obtain an exponential-type deviation inequality for this Bahadur representation.

Many asymptotic Bahadur-type representations for robust $M$ -estimators have been obtained in the literature; see, for example, Portnoy (1985), Mammen (1989) and He and Shao (1996, 2000), among others. Our result (1.4), however, is nonasymptotic and provides an explicit tail bound for the remainder term $R_{n}(\tau)$ . To obtain such a result, we first derive a sub-Gaussian-type deviation bound for $\widehat{\bm{\theta}}$ , and then conduct a careful analysis on the higher-order remainder term using this bound and techniques from empirical process theory. The expansion (1.4) further yields two classical normal approximation results, the Berry-Esseen inequality and Cramér-type moderate deviation. These results have important applications to large-scale inference [Fan, Hall and Yao (2007), Delaigle, Hall and Jin (2011), Liu and Shao (2014), Chang, Shao and Zhou (2016)]. Consider the statistical problems of simultaneously testing many hypotheses with FDR/FDP control or globally inferring a high dimensional parameter. For multiple testing, the obtained Berry-Esseen bound and Cramér-type moderate deviation result can be used to investigate the robustness and accuracy of the $P$ -values and critical values. For globally testing a high dimensional parameter, the expansion (1.4), combined with the parametric bootstrap, can be used to construct a valid test. In this paper, we only focus on the large-scale multiple testing problem and leave the global testing in high dimensions for future research.

2 FDP control for robust dependent tests

We apply the Bahadur representation (1.4) to construct robust dependence-adjusted test statistics for simultaneous inference. Conventional tasks of large-scale multiple testing, including controlling the FDR/FDP or FWER, have been extensively explored and are now well understood when the test statistics are independent [Benjamini and Hochberg (1995), Storey (2002), Genovese and Wasserman (2004), Lehmann and Romano (2005)]. It is becoming increasingly important to understand and incorporate the dependence information among multiple test statistics. Under the positive regression dependence condition, the FDR control can be conducted in the same manner as that for the independent case [Benjamini and Yekutieli (2001)], which provides a conservative upper bound. For more general dependence, directly applying standard FDR control procedures developed for independent $P$ -values may lead to inaccurate false discovery rate control and include too many spurious discoveries [Efron (2004, 2007), Sun and Cai (2009), Clarke and Hall (2009), Schwartzman and Lin (2011), Fan, Han and Gu (2012)]. In this more challenging situation, various multi-factor models have been used to investigate the dependence structure in high dimensional data; see, for example, Leek and Storey (2008), Friguet, Kloareg and Causeur (2009), Desai and Storey (2012) and Fan, Han and Gu (2012).

The multi-factor model relies on the identification of a linear space of random vectors capturing the dependence structure of the data. Friguet, Kloareg and Causeur (2009) and Desai and Storey (2012) assume that the data are drawn from a strict factor model with independent idiosyncratic errors. They use the expectation-maximization algorithm to estimate the factor loadings and realized factors in the model, and then obtain an estimator for the FDP by subtracting out realized common factors. These methods, however, depend on stringent model assumptions, including the independence of idiosyncratic errors and joint normality of the factor and noise. In contrast, Fan, Han and Gu (2012) and Fan and Han (2017) use a more general approximate factor model that allows dependent noise.

Let $\bm{X}=(X_{1},\ldots,X_{p})^{\intercal}$ be a $p$ -dimensional random vector with mean $\bm{\mu}=(\mu_{1},\ldots,\mu_{p})^{\intercal}$ and covariance matrix $\bm{\Sigma}=(\sigma_{jk})_{1\leq j,k\leq p}$ . We aim to simultaneously test

We are interested in the case where there is strong dependence across the components of $\bm{X}$ . The approximate factor model assumes that the dependence of a high dimensional random vector $\bm{X}$ can be captured by a few factors, that is,

from which we observe independent random samples $(\bm{X}_{1},\bm{f}_{1}),\ldots,(\bm{X}_{n},\bm{f}_{n})$ satisfying

For testing the hypotheses in (1.5) under model (1.6), when the factor $\bm{f}$ is unobserved, a popular and natural approach is based on (marginal) sample averages of $\{\bm{X}_{i}\}$ with a focus on the valid control of FDR [Benjamini and Yekutieli (2001), Genovese and Wasserman (2004), Efron (2007), Sun and Cai (2009), Schwartzman and Lin (2011), Fan and Han (2017)]. As pointed out in Fan, Han and Gu (2012), the power of such an approach is dominated by the factor-adjusted approach that produces alternative rankings of statistical significance from those of the marginal statistics. They focus on a Gaussian model where both $\bm{f}$ and $\bm{u}$ follow multivariate normal distributions, while statistical properties of the corresponding test procedure on FDR/FDP control remain unclear. The normality assumption, however, is really an idealization which provides insights into the key issues underlying the problems. Data subject to heavy-tailed and asymmetric errors are repeatedly observed in many fields of research [Finkenstadt and Rootzén (2003)]. For example, it is known that financial returns typically exhibit heavy tails. The important papers by Mandelbrot (1963) and Fama (1963) provide evidence of power-law behavior in asset prices in the early 1960s. Since then, the non-Gaussian character of the distribution of price changes has been widely observed in various market data. Cont (2001) provides further evidence showing that a Student’s $t$ -distribution with four degrees of freedom displays a tail behavior similar to many asset returns.

For multiple testing with heavy-tailed data, the least squares based test statistics are sensitive to outliers and thus lack robustness. This issue is amplified further by high dimensionality: When the dimension is large, even moderate tails may lead to significant false discoveries. This motivates us to develop new test statistics that are robust to the tails of error distributions. Also, since the multiple testing problem is more complicated with dependent data, theoretical guarantees of the FDP control for the existing dependence-adjusted methods remain unclear.

To illustrate the impact of heavy tailedness, we generate independent and identically distributed (i.i.d.) random variables $\{X_{ij},\,i=1,\ldots,30,\,j=1,\ldots,10000\}$ from a normalized $t$ -distribution with $2.5$ degrees of freedom. In Figure 1, we compare the histogram of the empirical means $\overline{X}_{j}$ with that of the robust mean estimates constructed using (1.3) without covariates, after rescaling both estimators by $\sqrt{30}$ . For a standard normal distribution, we expect $99.73\%$ data points to lie within three standard deviations of the mean or inside $ $. Hence for this experiment, if the distribution of the estimator is indeed approximately normal, we would expect about$ 27 $out of$ 10000 $realizations to lie outside <math xmlns="http://www.w3.org/1998/Math/MathML" display="block"><semantics><mrow><mi mathvariant="normal">.</mi><mi>F</mi><mi>r</mi><mi>o</mi><mi>m</mi><mi>F</mi><mi>i</mi><mi>g</mi><mi>u</mi><mi>r</mi><mi>e</mi><mn>1</mn><mi>w</mi><mi>e</mi><mi>s</mi><mi>e</mi><mi>e</mi><mi>t</mi><mi>h</mi><mi>a</mi><mi>t</mi><mi>t</mi><mi>h</mi><mi>e</mi><mi>r</mi><mi>o</mi><mi>b</mi><mi>u</mi><mi>s</mi><mi>t</mi><mi>p</mi><mi>r</mi><mi>o</mi><mi>c</mi><mi>e</mi><mi>d</mi><mi>u</mi><mi>r</mi><mi>e</mi><mi>g</mi><mi>i</mi><mi>v</mi><mi>e</mi><mi>s</mi><mn>28</mn><mi>p</mi><mi>o</mi><mi>i</mi><mi>n</mi><mi>t</mi><mi>s</mi><mi>t</mi><mi>h</mi><mi>a</mi><mi>t</mi><mi>f</mi><mi>a</mi><mi>l</mi><mi>l</mi><mi>o</mi><mi>u</mi><mi>t</mi><mi>s</mi><mi>i</mi><mi>d</mi><mi>e</mi><mi>t</mi><mi>h</mi><mi>i</mi><mi>s</mi><mi>i</mi><mi>n</mi><mi>t</mi><mi>e</mi><mi>r</mi><mi>v</mi><mi>a</mi><mi>l</mi><mo separator="true">,</mo><mi>w</mi><mi>h</mi><mi>e</mi><mi>r</mi><mi>e</mi><mi>a</mi><mi>s</mi><mi>t</mi><mi>h</mi><mi>e</mi><mi>s</mi><mi>a</mi><mi>m</mi><mi>p</mi><mi>l</mi><mi>e</mi><mi>a</mi><mi>v</mi><mi>e</mi><mi>r</mi><mi>a</mi><mi>g</mi><mi>e</mi><mi>g</mi><mi>i</mi><mi>v</mi><mi>e</mi><mi>s</mi><mi>a</mi><mi>m</mi><mi>u</mi><mi>c</mi><mi>h</mi><mi>l</mi><mi>a</mi><mi>r</mi><mi>g</mi><mi>e</mi><mi>r</mi><mi>n</mi><mi>u</mi><mi>m</mi><mi>b</mi><mi>e</mi><mi>r</mi><mo separator="true">,</mo><mn>79</mn><mo separator="true">,</mo><mi>m</mi><mi>a</mi><mi>n</mi><mi>y</mi><mi>o</mi><mi>f</mi><mi>w</mi><mi>h</mi><mi>i</mi><mi>c</mi><mi>h</mi><mi>w</mi><mi>o</mi><mi>u</mi><mi>l</mi><mi>d</mi><mi>s</mi><mi>u</mi><mi>r</mi><mi>e</mi><mi>l</mi><mi>y</mi><mi>b</mi><mi>e</mi><mi>f</mi><mi>a</mi><mi>l</mi><mi>s</mi><mi>e</mi><mi>l</mi><mi>y</mi><mi>r</mi><mi>e</mi><mi>g</mi><mi>a</mi><mi>r</mi><mi>d</mi><mi>e</mi><mi>d</mi><mi>a</mi><mi>s</mi><mi>s</mi><mi>i</mi><mi>g</mi><mi>n</mi><mi>a</mi><mi>l</mi><mi>s</mi><mi mathvariant="normal">.</mi><mi>W</mi><mi>e</mi><mi>s</mi><mi>e</mi><mi>e</mi><mi>t</mi><mi>h</mi><mi>a</mi><mi>t</mi><mi>i</mi><mi>n</mi><mi>t</mi><mi>h</mi><mi>e</mi><mi>p</mi><mi>r</mi><mi>e</mi><mi>s</mi><mi>e</mi><mi>n</mi><mi>c</mi><mi>e</mi><mi>o</mi><mi>f</mi><mi>h</mi><mi>e</mi><mi>a</mi><mi>v</mi><mi>y</mi><mi>t</mi><mi>a</mi><mi>i</mi><mi>l</mi><mi>s</mi><mi>a</mi><mi>n</mi><mi>d</mi><mi>h</mi><mi>i</mi><mi>g</mi><mi>h</mi><mi>d</mi><mi>i</mi><mi>m</mi><mi>e</mi><mi>n</mi><mi>s</mi><mi>i</mi><mi>o</mi><mi>n</mi><mi>s</mi><mo separator="true">,</mo><mi>t</mi><mi>h</mi><mi>e</mi><mi>r</mi><mi>o</mi><mi>b</mi><mi>u</mi><mi>s</mi><mi>t</mi><mi>m</mi><mi>e</mi><mi>t</mi><mi>h</mi><mi>o</mi><mi>d</mi><mi>l</mi><mi>e</mi><mi>a</mi><mi>d</mi><mi>s</mi><mi>t</mi><mi>o</mi><mi>a</mi><mi>m</mi><mi>o</mi><mi>r</mi><mi>e</mi><mi>a</mi><mi>c</mi><mi>c</mi><mi>u</mi><mi>r</mi><mi>a</mi><mi>t</mi><mi>e</mi><mi>n</mi><mi>o</mi><mi>r</mi><mi>m</mi><mi>a</mi><mi>l</mi><mi>t</mi><mi>a</mi><mi>i</mi><mi>l</mi><mi>a</mi><mi>p</mi><mi>p</mi><mi>r</mi><mi>o</mi><mi>x</mi><mi>i</mi><mi>m</mi><mi>a</mi><mi>t</mi><mi>i</mi><mi>o</mi><mi>n</mi><mi>t</mi><mi>h</mi><mi>a</mi><mi>n</mi><mi>u</mi><mi>s</mi><mi>i</mi><mi>n</mi><mi>g</mi><mi>a</mi><mi>n</mi><mi>o</mi><mi>n</mi><mi>r</mi><mi>o</mi><mi>b</mi><mi>u</mi><mi>s</mi><mi>t</mi><mi>o</mi><mi>n</mi><mi>e</mi><mi mathvariant="normal">.</mi><mi>I</mi><mi>n</mi><mi>f</mi><mi>a</mi><mi>c</mi><mi>t</mi><mo separator="true">,</mo><mi>f</mi><mi>o</mi><mi>r</mi><mi>t</mi><mi>h</mi><mi>e</mi><mi>e</mi><mi>m</mi><mi>p</mi><mi>i</mi><mi>r</mi><mi>i</mi><mi>c</mi><mi>a</mi><mi>l</mi><mi>m</mi><mi>e</mi><mi>a</mi><mi>n</mi><mi>s</mi><mo separator="true">,</mo><mi>m</mi><mi>a</mi><mi>n</mi><mi>y</mi><mi>o</mi><mi>f</mi><mi>t</mi><mi>h</mi><mi>e</mi><mi>m</mi><mi>e</mi><mi>v</mi><mi>e</mi><mi>n</mi><mi>f</mi><mi>a</mi><mi>l</mi><mi>l</mi><mi>o</mi><mi>u</mi><mi>t</mi><mi>s</mi><mi>i</mi><mi>d</mi><mi>e</mi></mrow><annotation encoding="application/x-tex">. From Figure 1 we see that the robust procedure gives 28 points that fall outside this interval, whereas the sample average gives a much larger number, 79, many of which would surely be falsely regarded as signals. We see that in the presence of heavy tails and high dimensions, the robust method leads to a more accurate normal tail approximation than using a nonrobust one. In fact, for the empirical means, many of them even fall outside</annotation></semantics></math>.FromFigure1weseethattherobustproceduregives28pointsthatfalloutsidethisinterval,whereasthesampleaveragegivesamuchlargernumber,79,manyofwhichwouldsurelybefalselyregardedassignals.Weseethatinthepresenceofheavytailsandhighdimensions,therobustmethodleadstoamoreaccuratenormaltailapproximationthanusinganonrobustone.Infact,fortheempiricalmeans,manyofthemevenfalloutside. This inaccuracy in tail approximation for the nonrobust estimator gives rise to false discoveries. In summary, outliers from the test statistics$ \overline{X}_{j}$ can be so large that they are mistakenly regarded as discoveries, whereas the robust approach results in fewer outliers.

In Section 3, we develop robust dependence-adjusted multiple testing procedures with solid theoretical guarantees. We use the approximate factor model (1.6) with an observable factor and relatively heavy-tailed errors to characterize the dependence structure in high dimensional data. Assuming such a model, we construct robust test statistics based on the Huber estimator with a diverging tuning parameter, denoted by $T_{1},\ldots,T_{p}$ , for testing the individual hypotheses. At a prespecified level $0<\alpha<1$ , we apply a family of FDP controlling procedures to the dependence-adjusted $P$ -values $\{P_{j}=2\Phi(-|T_{j}|)\}_{j=1}^{p}$ to decide which null hypotheses are rejected, where $\Phi$ is the standard normal distribution function. To justify the validity of the resulting procedure on FDP control, a delicate analysis of the impact of dependence-adjustment on the distribution of the $P$ -values is required. We show that, under mild moment and regularity conditions, the robust multiple testing procedure controls the FDP at any prespecified level asymptotically. Specifically, applying Storey’s procedure [Storey (2002)] to the above $P$ -values gives a data-driven rejection threshold $\widehat{z}_{{\rm N}}$ such that $H_{0j}$ is rejected whenever $|T_{j}|\geq\widehat{z}_{{\rm N}}$ . Let ${\rm FDP}(z)=V(z)/\max\{1,R(z)\}$ be the FDP at threshold $z\geq 0$ , where $V(z)=\sum_{j=1}^{p}1(|T_{j}|\geq z,\mu_{j}=0)$ and $R(z)=\sum_{j=1}^{p}1(|T_{j}|\geq z)$ are the number of false discoveries and the number of total discoveries, respectively. In the ultra-high dimensional setting that $p$ can be as large as $e^{n^{c}}$ for some $0<c<1$ , we prove that

as $(n,p)\to\infty$ , where $p_{0}=\sum_{j=1}^{p}1(\mu_{j}=0)$ is the number of true null hypotheses. We also illustrate the usefulness of the robust techniques by contrasting the performances of robust and least squares based inference procedures through synthetic numerical experiments.

Key technical tools in proving (1.7) are the Berry-Esseen bound and Cramér-type moderate deviation for the marginal statistic $T_{j}$ . These results are built upon the nonasymptotic Bahadur representation (1.4), and may be of independent interest for other statistical applications. For example, Delaigle, Hall and Jin (2011) explore moderate and large deviations of the $t$ -statistic in a variety of high dimensional settings.

3 Organization of the paper

The lay-out of the paper is as follows. In Section 2, we develop a general finite sample theory for Huber’s robust $M$ -estimator from a new perspective where a diverging tuning parameter is involved. In Section 3, we propose a robust dependence-adjusted multiple testing procedure with rigorous theoretical guarantees. Section 4 consists of numerical studies and real data analysis. The simulation study provides empirical evidence that the proposed robust inference procedure improves performance in the presence of asymmetric and heavy-tailed errors, and maintains efficiency under light-tailed situations. A discussion is given in Section 5. Proofs of the theoretical results in Sections 2 and 3 are provided in the supplemental material [Zhou et al. (2017)].

Robust M𝑀M-estimation: A finite sample theory

Consider a heteroscedastic linear regression model $Y=\mu^{*}+\bm{X}^{\intercal}\bm{\beta}^{*}+\sigma(\bm{X})\varepsilon$ , from which we observe independent samples $\{(Y_{i},\bm{X}_{i})\}_{i=1}^{n}$ satisfying

In this section, we study the robust estimator of $\bm{\theta}^{*}$ defined in (1.2). In particular, we show that it admits an exponential-type deviation bound even for heavy-tailed error distributions. Note that, under the heteroscedastic model (2.1), $\bm{\theta}^{*}$ differs from the median effect of $Y$ conditioning on $\bm{X}$ , so that the LAD-based methods are not applicable to estimate $\bm{\theta}^{*}$ . Instead, we focus on Huber’s robust estimator $\widehat{\bm{\theta}}$ given in (1.2) with a diverging tuning parameter $\tau=\tau_{n}$ that balances the approximation error and robustness of the estimator. To begin with, we make the following conditions on the linear model (2.1).

Condition 2.1 allows a family of conditional heteroscedastic models with heavy-tailed error $\varepsilon$ . Specifically, it only requires the second moment of $\nu=\sigma(\bm{X})\varepsilon$ to be finite. Under this condition, our first result, Theorem 2.1, provides an exponential-type deviation bound and a nonasymptotic Bahadur representation for the robust estimator $\widehat{\bm{\theta}}=(\widehat{\mu},\widehat{\bm{\beta}}^{\intercal})^{\intercal}$ defined in (1.2).

Under the linear model (2.1) with Condition 2.1 satisfied, we have for any $w>0$ that, the robust estimator $\widehat{\bm{\theta}}$ in (1.2) with $\tau=\tau_{n}=\tau_{0}\sqrt{n}(d+1+w)^{-1/2}$ and $\tau_{0}\geq\sigma$ satisfies

An important message of Theorem 2.1 is that, even for heavy-tailed errors with only finite second moment, the robust estimator $\widehat{\bm{\theta}}$ with properly chosen $\tau$ has sub-Gaussian tails. See inequality (2.2). To some extent, the tuning parameter $\tau$ plays a similar role as the bandwidth in constructing nonparametric estimators. Furthermore, we show in (2.3) that the remainder of the Bahadur representation for $\widehat{\bm{\theta}}$ exhibits sub-exponential tails. To the best of our knowledge, no nonasymptotic results of this type exist in the literature, and classical asymptotic results can only be used to derive polynomial-type deviation bounds.

Together, Theorems 2.1 and 2.2 lead to a Berry-Esseen type bound for $T:=\sqrt{n}(\widehat{\mu}-\mu^{*})/\sigma$ for properly chosen $\tau$ . In addition, the following theorem gives a Cramér-type moderation deviation result for $T$ , which quantifies the relative error of the normal approximation.

uniformly for $0\leq z=o\{\min(\sqrt{w_{n}},\sqrt{n}w_{n}^{-1})\}$ as $n\to\infty$ , where $G\sim N(0,1)$ ,

and $C>0$ is a constant independent of $n$ . In particular, we have

Motivated by an application to large-scale simultaneous inference considered in Section 3, we only focus on the robust intercept estimator $\widehat{\mu}$ of $\mu^{*}$ in Theorems 2.2 and 2.3. In fact, similar results can be obtained for $\widehat{\bm{\beta}}$ or a specific coordinate of $\widehat{\bm{\beta}}$ based on the Bahadur representation (2.3).

Large-scale multiple testing for heavy-tailed dependent data

In this section, we propose and analyze a robust dependence-adjusted procedure for simultaneously testing the means $\mu_{1},\ldots,\mu_{p}$ in model (1.6), based on independent observations from the population vector $\bm{X}$ which exhibits strong dependence and heavy tails.

Suppose we are given independent random samples $\{(\bm{X}_{i},\bm{f}_{i})\}_{i=1}^{n}$ from model (1.6). We are interested in the simultaneous testing of mean effects (1.5). A naive approach under normality is to directly use the information $X_{ij}\sim N(\mu_{j},\sigma_{jj})$ for the dependent case as was done in the literature, where $\sigma_{jj}=\mbox{var}(X_{j})$ . Such an approach is very natural and popular when the factors are unobservable and focus is on the valid control of FDR, but is inefficient as noted in Fan, Han and Gu (2012). Indeed, if the loading matrix $\mathbf{B}$ is known and the factors are observed (otherwise, replaced by their estimates), for each $j$ , we can construct the marginal test statistic using dependence-adjusted observations $\{X_{ij}-\bm{b}_{j}^{\intercal}\bm{f}_{i}\}_{i=1}^{n}$ from $\mu_{j}+\sigma_{j}(\bm{f})u_{j}$ for testing the $j$ th hypothesis $H_{0j}:\mu_{j}=0$ .

We consider the approximate factor model (1.6) and write

Let $\bm{\Sigma}_{f}$ and $\bm{\Sigma}_{\nu}=(\sigma_{\nu,jk})_{1\leq j,k\leq p}$ denote the covariance matrices of $\bm{f}$ and $\bm{\nu}$ , respectively. Under certain sparsity condition on $\bm{\Sigma}_{\nu}$ (see Section 3.4 for an elaboration), $\nu_{1},\ldots,\nu_{p}$ are weakly dependent random variables with higher signal-to-noise ratios since $\mbox{var}(\nu_{j})=\sigma_{jj}-\|\bm{\Sigma}_{f}^{1/2}\bm{b}_{j}\|^{2}<\sigma_{jj}$ . Therefore, subtracting common factors out makes the resulting FDP control procedure more efficient and powerful. It provides an alternative ranking of the significance of hypothesis from the tests based on marginal statistics.

For each $j=1,\ldots,p$ , we have a linear regression model

A natural approach is to estimate $\mu_{j}$ and $\bm{b}_{j}$ by the method of least squares. However, the least squares method is sensitive to the tails of the error distributions. Also, as noted in Fan, Li and Wang (2017), the LAD-based methods are not applicable in the presence of asymmetric and heteroscedastic errors. Hence, we suggest a robust method that simultaneously estimates $\mu_{j}$ and $\bm{b}_{j}$ by solving

To construct a test statistic for the individual hypothesis $H_{0j}:\mu_{j}=0$ with pivotal limiting distribution, we need to estimate $\sigma_{\nu,jj}=\sigma_{jj}-\mbox{var}(\bm{b}_{j}^{\intercal}\bm{f})$ . For $\mbox{var}(\bm{b}_{j}^{\intercal}\bm{f})$ , a natural and simple estimator is $\widehat{\bm{b}}_{j}^{\intercal}\widehat{\bm{\Sigma}}_{f}\widehat{\bm{b}}_{j}$ , where $\widehat{\bm{\Sigma}}_{f}:=n^{-1}\sum_{i=1}^{n}\bm{f}_{i}\bm{f}_{i}^{\intercal}$ . Let $\widehat{\sigma}_{jj}$ and $\widehat{\sigma}_{\nu,jj}$ be generic estimators of $\sigma_{jj}$ and $\sigma_{\nu,jj}$ , respectively. To simultaneously infer all the hypotheses of interest, we need the following uniform convergence results

For $\sigma_{jj}=\mbox{var}(X_{j})$ , it is known that the sample variance $n^{-1}\sum_{i=1}^{n}(X_{ij}-\overline{X}_{j})^{2}$ performs poorly when $X_{j}$ has heavy tails. Based on the recent developments of robust mean estimation for heavy-tailed data [Catoni (2012), Joly and Lugosi (2016), Fan, Li and Wang (2017)], we consider the following two types of robust variance estimators.

Back to the current problem, we aim to estimate $\sigma_{jj}$ based on independent observations $X_{1j},\ldots,X_{nj}$ . Let $V=V_{n}<n$ be an integer and decompose $n$ as $n=Vm+r$ for some integer $0\leq r<V$ . Let $B_{1},\ldots,B_{V}$ be a partition of $\{1,\ldots,n\}$ defined by

Given robust mean and variance estimators of each type, we construct dependence-adjusted test statistics

2 Dependence-adjusted FDP control procedure

To conduct multiple testing of (1.5) using the test statistics $T_{j}$ ’s, let $z>0$ be the critical value to be determined such that $H_{0j}$ is rejected whenever $|T_{j}|\geq z$ . The main object of interest in this paper is the false discovery proportion

where $V(z)=\sum_{j\in\mathcal{H}_{0}}1(|T_{j}|\geq z)$ is the number of false discoveries, $R(z)=\sum_{j=1}^{p}1(|T_{j}|\geq z)$ and $\mathcal{H}_{0}=\{j:1\leq j\leq p,\mu_{j}=0\}$ represents the set of true null hypotheses. There is substantial interest in controlling the FDP at a prespecified level $0<\alpha<1$ for which the ideal rejection threshold is $z_{{\rm oracle}}=\inf\{z\geq 0:{\rm FDP}(z)\leq\alpha\}$ .

The statistical behavior of ${\rm FDP}(z)$ is the center of interest in multiple testing. However, the realization of $V(z)$ for a given experiment is unknown and thus needs to be estimated. When the sample size is large, it is natural to approximate $V(z)$ by its expectation $2p_{0}\Phi(-z)$ , where $p_{0}={\rm Card}(\mathcal{H}_{0})$ . In the high dimensional sparse setting, both $p$ and $p_{0}$ are large and $p_{1}=p-p_{0}=o(p)$ is relatively small. Therefore, we can use $p$ as a slightly conservative surrogate for $p_{0}$ , so that ${\rm FDP}(z)$ can be approximated by

We will prove in Theorem 3.2 that under mild conditions, ${\rm FDP}_{{\rm N}}(z)$ provides a consistent estimate of ${\rm FDP}(z)$ uniformly in $0\leq z\leq\Phi^{-1}(1-m_{p}/(2p))$ for any sequence of positive number $m_{p}\leq 2p$ satisfying $m_{p}\to\infty$ .

In the non-sparse case where $\pi_{0}=p_{0}/p$ is bounded away from 0 and 1 as $p\to\infty$ , ${\rm FDP}_{{\rm N}}$ given in (3.9) tends to overestimate the true FDP. Therefore, we need to estimate the proportion $\pi_{0}$ , which has been studied by Efron et al. (2001), Storey (2002), Genovese and Wasserman (2004), Langaas and Lindqvist (2005) and Meinshausen and Rice (2006), among others. For simplicity, we focus on Storey’s approach. Let $\{P_{j}=2\Phi(-|T_{j}|)\}_{j=1}^{p}$ be the approximate $P$ -values. For a predetermined $\lambda\in[0,1)$ , Storey (2002) suggests the following conservative estimate of $\pi_{0}$ :

The intuition of such an estimator is as follows. Since most of the large $P$ -values correspond to the null and thus are uniformly distributed, for a sufficiently large $\lambda$ , we expect about $(1-\lambda)\pi_{0}$ of the $P$ -values to lie in $(\lambda,1]$ . Hence, the proportion of $P$ -values that exceed $\lambda$ , $p^{-1}\sum_{j=1}^{p}1(P_{j}>\lambda)$ , should be close to $(1-\lambda)\pi_{0}$ . This gives rise to Storey’s procedure.

Incorporating such an estimate of $\pi_{0}$ , we obtain a modified estimate of ${\rm FDP}(z)$ by

In view of (3.9)–(3.11) and the fact $\widehat{\pi}_{0}(0)=1$ , we have ${\rm FDP}_{{\rm N},0}(z)={\rm FDP}_{{\rm N}}(z)$ .

By replacing the unknown quantity ${\rm FDP}(z)$ by ${\rm FDP}_{{\rm N},\lambda}(z)$ given in (3.11) for some $\lambda\in[0,1)$ , we reject $H_{0j}$ whenever $|T_{j}|\geq\widehat{z}_{{\rm N},\lambda}$ , where

By Lemmas 1 and 2 in Storey, Taylor and Siegmund (2004), this procedure is equivalent to a variant of the seminal Benjamini-Hochberg (B-H) procedure [Benjamini and Hochberg (1995)] for selecting $\mathcal{S}=\{j:1\leq j\leq p,P_{j}\leq P_{(k_{p}(\lambda))}\}$ based on the $P$ -values $P_{j}=2\Phi(-|T_{j}|)$ , where $k_{p}(\lambda):=\max\{j:1\leq j\leq p,P_{(j)}\leq\frac{\alpha j}{\widehat{\pi}_{0}(\lambda)p}\}$ and $P_{(1)}\leq\cdots\leq P_{(p)}$ are the ordered $P$ -values. Theorem 3.3 shows that under weak moment conditions, the FDP of this dependence-adjusted procedure with $\lambda=0$ converges to $\alpha$ in the ultra-high dimensional sparse setting.

Note that ${\rm FDP}_{{\rm N},0}(z)$ is the most conservatively biased estimate of ${\rm FDP}(z)$ among all $\lambda\in[0,1)$ using normal calibration. The statistical power of the corresponding procedure can be compromised if $\pi_{0}$ is much smaller than 1. In general, the procedure requires the choice of a tuning parameter $\lambda$ in the estimate $\widehat{\pi}(\lambda)$ , which leads to an inherent bias-variance trade-off. We refer to Section 9 in Storey (2002) and Section 6 in Storey, Taylor and Siegmund (2004) for two data-driven methods for automatically choosing $\lambda$ .

3 Bootstrap calibration

For $j=1,\ldots,p$ , define empirical tail distributions

The bootstrap $P$ -values are thus given by $\{P^{*}_{j}=G^{*}_{j,B}(|\widehat{\mu}_{j}|)\}_{j=1}^{p}$ , to which either the B-H procedure or Storey’s procedure can be applied. For the former, we reject $H_{0j}$ whenever $P^{*}_{j}\leq P^{*}_{(k_{p}^{*})}$ , where $k_{p}^{*}=\max\{j:1\leq j\leq p,P^{*}_{(j)}\leq j\alpha/p\}$ for a predetermined $0<\alpha<1$ and $P^{*}_{(1)}\leq\cdots\leq P^{*}_{(p)}$ are the ordered bootstrap $P$ -values. For the distribution of the bootstrap weights, it is common to choose $W\sim 2{\rm Bernoulli}(0.5)$ , $W\sim{\rm exp}(1)$ or $W\sim N(1,1)$ in practice. Using nonnegative random weights has the advantage that the weighted objective function is guaranteed to be convex.

Weighted bootstrap procedure serves as an alternative method to normal calibration in multiple testing. We refer to Spokoiny and Zhilova (2015) and Zhilova (2016) for the most advanced recent results of weighted bootstrap and a comprehensive literature review. We leave the theoretical guarantee of this procedure for future research.

4 Theoretical properties

First, we impose some conditions on the distribution of $\bm{X}$ and the tuning parameters $\tau$ and $\gamma$ that are used in the robust regression and robust estimation of the second moment.

$(\tau,\gamma)=(\tau_{n},\gamma_{n})$ satisfies $\tau=\tau_{0}\sqrt{n}w_{n}^{-1/2}$ and $\gamma=\gamma_{0}\sqrt{n}w_{n}^{-1/2}$ for some constants $\tau_{0}\geq\max_{1\leq j\leq p}\sigma_{\nu,jj}^{1/2}$ and $\gamma_{0}\geq\max_{1\leq j\leq p}{\rm var}^{1/2}(X_{j}^{2})$ , where the sequence $w_{n}$ is such that $w_{n}\to\infty$ and $w_{n}=o(\sqrt{n})$ .

$\max_{1\leq j<k\leq p}|\rho_{\nu,jk}|\leq\rho$ and

for some $0<\rho<1$ , $\kappa>0$ and $0<r<(1-\rho)/(1+\rho)$ . As $n,p\to\infty$ , $p_{0}/p\to\pi_{0}\in(0,1]$ , $\log p=o(n^{1/5})$ and $w_{n}\asymp n^{1/5}$ , where $w_{n}$ is as in Condition (C2).

${\rm Card}\big{\{}j:1\leq j\leq p,\sigma_{\nu,jj}^{-1/2}|\mu_{j}|\geq\lambda\sqrt{(\log p)/n}\big{\}}\to\infty$ as $n,p\to\infty$ for some $\lambda>2\sqrt{2}$ .

Condition (C3) allows weak dependence among $\nu_{1},\ldots,\nu_{p}$ in the sense that each variable is moderately correlated with $s_{p}$ other variables and weakly correlated with the remaining ones. The technical assumption (C4) imposes a constraint on the number of significant true alternatives, which is slightly stronger than $p_{1}\to\infty$ . According to Proposition 2.1 in Liu and Shao (2014), this condition is nearly optimal for the results on FDP control in the sense that if $p_{1}$ is fixed, the B-H method fails to control the FDP at any level $0<\beta<1$ with overwhelming probability even if the true $P$ -values were known.

where $\mathcal{S}_{j}=\{P_{j}^{{\rm true}}>\alpha/p\}$ , $j=1,\ldots,p$ .

Assume that Conditions (C1) and (C2) hold and $\log p=o\{\min(w_{n},nw_{n}^{-2})\}$ . Then (3.13) holds.

Theorem 3.1 shows that, to ensure the accuracy of the normal distribution calibration, the number of simultaneous tests can be as large as $\exp\{o(n^{1/3})\}$ , when taking $w_{n}\asymp n^{1/3}$ . We are also interested in estimating FDP in the high dimensional sparse setting, that is, $p$ is large, but the number of $\mu_{j}\neq 0$ is relatively small. The following result indicates that ${\rm FDP}_{{\rm N}}(z)$ given in (3.9) provides a consistent estimator of the realized FDP in a uniform sense.

Assume that Conditions (C1)–(C3) hold. Then, for any sequence of positive numbers $m_{p}\leq p$ satisfying $m_{p}\to\infty$ , we have as $(n,p)\to\infty$ ,

Further, Theorem 3.3 shows that the proposed robust dependence-adjusted inference procedure controls the FDP at a given level $\alpha$ asymptotically with $P$ -values estimated from the standard normal distribution.

Assume that Conditions (C1)–(C4) hold. Then, for any prespecified $0<\alpha<1$ ,

as $(n,p)\to\infty$ , where $\widehat{z}_{{\rm N},0}$ is defined in (3.12).

The constraint on $p$ , as a function of $n$ , imposed in Theorems 3.2 and 3.3 can be relaxed in a strict factor model with independent idiosyncratic errors.

$\nu_{1},\ldots,\nu_{p}$ in model (1.6) are independent. As $n,p\to\infty$ , $p_{0}/p\to\pi_{0}\in(0,1]$ , $\log p=o(w_{n})$ and $w_{n}=O(n^{1/3})$ , where $w_{n}$ is as in Condition (C2).

Assume that Conditions (C1), (C2), (C4) and (C5) hold. Then, for any prespecified $0<\alpha<1$ , $(p_{0}/p)^{-1}{\rm FDP}(\widehat{z}_{{\rm N},0})\to\alpha$ in probability as $(n,p)\to\infty$ , where $\widehat{z}_{{\rm N},0}$ is defined in (3.12).

Theorems 3.2–3.4 provide theoretical guarantees on the FDP control for the B-H procedure with dependence-adjusted $P$ -values $P_{j}=2\Phi(-|T_{j}|)$ , $j=1,\ldots,p$ . A similar approach can be defined by using the median-of-means approach, namely, replacing $T_{j}$ ’s with $S_{j}$ ’s in the definition of ${\rm FDP}_{{\rm N}}(z)$ in (3.9), which is equivalent to the B-H procedure with $P$ -values $Q_{j}=2\Phi(-|S_{j}|)$ , $j=1,\ldots,p$ . Under similar conditions, the theoretical results on the FDP control remain valid.

Let ${\rm FDP}(z)$ and $\widehat{z}_{{\rm N},0}$ be defined in (3.8) and (3.12) with $T_{j}$ ’s replaced by $S_{j}$ ’s, and let $V=V_{n}$ in (3.5) satisfy $V\asymp w_{n}$ for $w_{n}$ as in Condition (C2). Moreover, let $\tau=\tau_{n}$ be as in Condition (C2).

Under Conditions (C1), (C3) and (C4), (3.15) holds for any prespecified $0<\alpha<1$ .

Under Conditions (C1), (C4) and (C5), (3.15) holds for any prespecified $0<\alpha<1$ .

Numerical study

For each $1\leq j\leq p$ , we apply the above algorithm with $\tau=\tau_{j}:=c\,\widehat{\sigma}_{j}\sqrt{n/\log(np)}$ to obtain $(\widehat{\mu}_{j},\widehat{\bm{b}}_{j}^{\intercal})^{\intercal}$ , where $\widehat{\sigma}_{j}^{2}$ denotes the sample variance of the fitted residuals using OLS and $c>0$ is a control parameter. We take $c=2$ in all the simulations reported below. In practice, we can use a cross-validation procedure to pick $c$ from only a few candidates, say $\{0.5,1,2\}$ .

2 Simulations via a synthetic factor model

In this section, we perform Monte Carlo simulations to illustrate the performance of the robust test statistic under approximate factor models with general errors. Consider the Fama-French three factor model:

where $\bm{u}_{i}=(u_{i1},\ldots,u_{ip})^{\intercal}$ are i.i.d. copies of $\bm{u}=(u_{1},\ldots,u_{p})^{\intercal}$ . We simulate $\{\bm{b}_{j}\}_{j=1}^{p}$ and $\{\bm{f}_{i}\}_{i=1}^{n}$ independently from $N_{3}(\bm{\mu}_{B},\bm{\Sigma}_{B})$ and $N_{3}(\mathbf{0},\bm{\Sigma}_{f})$ , respectively. To make the model more realistic, parameters are calibrated from the daily returns of S&P 500’s top 100 constituents (chosen by market cap), for the period July 1st, 2008 to June 29th, 2012.

To generate dependent errors, we set $\bm{\Sigma}_{u}=\operatorname{cov}(\bm{u})$ to be a block diagonal matrix where each block is four-by-four correlation matrix with equal off-diagonal entries generated from Uniform $[0,0.5]$ . The hypothesis testing is carried out under the alternative: $\mu_{j}=\mu$ for $1\leq j\leq\pi_{1}p$ and $\mu_{j}=0$ otherwise. In the simulations reported here, the ambient dimension $p=2000$ , the proportion of true alternatives $\pi_{1}=0.25$ and the sample size $n$ takes values in $\{80,120\}$ . For simplicity, we set $\lambda=0.5$ in our procedure and use the Matlab package mafdr to compute the estimate $\widehat{\pi}_{0}(\lambda)$ of $\pi_{0}=1-\pi_{1}$ . For each test, the empirical false discovery rate (FDR) is calculated based on $500$ replications with FDR level $\alpha$ taking values in $\{5\%,10\%,20\%\}$ . The errors $\{\bm{u}_{i}\}_{i=1}^{n}$ are generated independently from the following distributions:

Model 1. $\bm{u}\sim N(\mathbf{0},\bm{\Sigma}_{u})$ : Centered normal random errors with covariance matrix $\bm{\Sigma}_{u}$ ;

Model 2. $\bm{u}\sim(1/\sqrt{5})\,t_{2.5}(\mathbf{0},\bm{\Sigma}_{u})$ : Symmetric and heavy-tailed errors following a multivariate $t$ -distribution with degrees of freedom 2.5 and covariance matrix $\bm{\Sigma}_{u}$ ;

For weakly dependent errors following the normal distribution and $t$ -distribution, Table 1 shows that the RD-A procedure consistently outperforms the OD-A method, in the sense that the RD-A method provides a much better control of the FDR at the expense of slight compromises of the FNR and TPR. Should the FDR being controlled at the same level, the robust method will be more powerful. In Table 2, when the errors are both asymmetric and heavy-tailed, the RD-A procedure has the biggest advantage in that it significantly outperforms the OD-A on controlling the FDR at all levels while maintaining low FNR and high TPR. Together, these results show that the RD-A procedure is indeed robust to outliers and does not lose efficiency when the errors are symmetric and light-tailed. In terms of controlling FDR, Models 3 and 4 present more challenges than Models 1 and 2 due to being both heavy-tailed and asymmetric. In Table 2 we see that although both the RD-A and OD-A methods achieve near-perfect power, the empirical FDR is higher than the desired level across all settings and much more higher for OD-A. Hence we compare the FDR of the RD-A and OD-A methods for various sample sizes in Figure 2. We see that the empirical FDR decreases with increase in sample size, while consistently outperforming the OD-A procedure. The difference between the two methods is greater for lower sample sizes, reinforcing the usefulness of our method for high dimensional heavy-tailed data with moderate sample sizes.

The naive procedure suffers from a significant loss in FNR and TPR. The reasons are twofold: (a) the naive procedure ignores the actual dependency structure among the variables; (b) the signal-to-noise ratio of $H_{0j}:\mu_{j}=0$ for the naive procedure is $\sigma_{jj}^{-1/2}|\mu_{j}|$ , which can be much smaller than $\sigma_{\nu,jj}^{-1/2}|\mu_{j}|$ for the dependence-adjusted procedure.

3 Stock market data

In this section, we apply our proposed robust dependence-adjusted multiple testing procedure to monthly stock market data. Consider Carhart’s four-factor model [Carhart (1997)] on S&P 500 index, where the excess return of a stock has the following representation:

for $j=1,\ldots,p$ and $t=1,\ldots,T$ . Here $r_{jt}$ is the excess returns of stock $j$ at month $t$ , $r_{ft}$ is the risk free interest rate at month $t$ , and MKT, HML, SMB and UMD represent the market, value, size, and momentum factors respectively. We are interested in the value $\mu_{j}$ , which represents the alpha of stock $j$ . A stock can be said to have excess returns if its alpha is positive, or in other words, the stock exhibits returns higher than those that can be accounted for by the four factors. If the alpha is negative, the stock is consistently underperforming, given the level of risk it undertakes. Detecting nonzero alpha is important since it is directly related to the efficient equity market hypothesis. When the market is inefficient, we can conduct multiple hypothesis testing to identify those stocks in the market that have statistically significant alphas. When the returns of mutual fund data are used, the test is related to test whether the fund manager has skills or not [Barras, Scaillet and Wermers (2010)]. All the data in this section was obtained from Kenneth French’s website and the COMPUSTAT and CRSP databases.

We obtain monthly data for $393$ S&P 500 constituents over the time period from January 2005 to December 2013, after removing those stocks that have missing values or have discontinuous inclusion in the index. The stock returns exhibit severely heavy-tails, as illustrated by the histogram of the excess kurtosis of the data in Figure 3. Among the $393$ series, $112$ have distributions whose tails are fatter than the $t$ -distribution with $5$ degrees of freedom.

The regression in (4.2) is carried out over rolling windows: for each month, we evaluate the model using data from the preceding three years. For each rolling window, we simultaneously test the hypotheses $H_{0j}:\mu_{j}=0$ versus $H_{1j}:\mu_{j}\neq 0$ for $j=1,\ldots,p$ , using the proposed robust dependence-adjusted procedure. We see that out of a portfolio of size $393$ , only a few stocks exhibit statistically significant nonzero alphas at the FDR threshold of $5\%$ , $10\%$ and $20\%$ . Table 3 summarizes the results for the number of selected stocks and estimated alphas with the FDR controlled at 5%. For the robust dependence-adjusted procedure, we see that no stock is selected on average, with maximum $2$ stocks selected over the entire time period. In particular, our method does not select any stocks from the third quarter of 2008 to the third quarter of 2010, coinciding with the financial crisis during which the market volatility is much higher. Moreover, the estimated alphas for the selected stocks are much higher than those not selected by the robust multiple testing procedure. This is represented as $|\widehat{\mu}_{j}|$ in Table 3. The naive method, which directly performs multiple $t$ -tests ignoring the common factors, appears to be unstable with the number of stocks selected being extremely variable. Additionally, a tremendously large number of stocks are selected in a few time periods, pointing towards false discoveries. In summary, the robust dependence-adjusted multiple testing procedure is particularly suited for the problem of finding a few stocks with nonzero alphas, which is explained by the focus on a balanced panel of highly traded stocks with large capitalizations, namely, the constituents of the S&P 500.

4 Gene expression data

In this section, we apply the proposed procedure to the analysis of a neuroblastoma data set reported in Oberthuer et al. (2006) to identify differentially expressed genes between the group of patients who had 3-year event-free survival after the diagnosis of neuroblastoma and the group of patients who did not. This data set consists of 251 patients of the German Neuroblastoma Trials NB90-NB2004, diagnosed between 1989 and 2004. The complete data set, obtained via the MicroArray Quality Control phase-II (MAQC-II) project [Shi et al. (2012)], includes gene expression over 10,707 probe sites. There are 246 subjects with 3-year event-free survival information available (56 positive and 190 negative). See Oberthuer et al. (2006) for more details about the data sets.

In the first stage, we use standard principal component analysis on the two samples to obtain the factors, based on which we construct dependence-adjusted $P$ -values to conduct multiple testing in the second step. Note that the test statistic given in (3.7) can be directly generalized to the two-sample case: Given two groups of $p$ -dimensional ( $p=10,707$ ) observations with sizes $n_{1}=56$ and $n_{2}=190$ , we compute robust mean and variance estimators $(\widehat{\mu}_{1j},\widehat{\mu}_{2j})$ and $(\widehat{\sigma}_{1\nu,jj},\widehat{\sigma}_{2\nu,jj})$ for $j=1,\ldots,p$ . Define two-sample test statistics $T_{j}=(\widehat{\mu}_{1j}-\widehat{\mu}_{2j})/(\widehat{\sigma}_{1\nu,jj}/n_{1}+\widehat{\sigma}_{2\nu,jj}/n_{2})^{1/2}$ so that the corresponding $P$ -values are $\{2\Phi(-|T_{j}|)\}_{j=1}^{p}$ .

The number of factors is estimated by the eigenvalue ratio estimator proposed in Ahn and Horenstein (2013), which was also used in the context of factor-adjusted multiple testing in Fan and Han (2017). The estimator is defined as $\widehat{K}=\text{argmax}_{1<k<k_{\max}}(\widehat{\lambda}_{k}/\widehat{\lambda}_{k+1})$ , where $\widehat{\lambda}_{j}$ is the $j$ th eigenvalue of the sample covariance matrix and $k_{\max}$ is the maximum possible number of factors. Following this procedure, we use $K=2$ to model the latent structure in the data.

Next, we conduct multiple testing using the proposed robust dependence-adjusted procedure and the naive procedure based on two-sample $t$ -tests. At FDR level 1%, we detect $3779$ genes and the naive procedure detects $3236$ genes; while at FDR level 5%, we discover $5223$ genes and the naive procedure discovers $4685$ genes. In general, taking the latent structure into account causes a visible increase in the number of genes that are declared statistically significant regardless of the prechosen FDR level, reflecting the improved power of our method. This phenomenon is in accord with that in Desai and Storey (2012). These results may serve as an exploratory step for more refined analyses regarding those significant genes.

Summary and discussion

This paper consists of two main parts with each one being of independent interest. In the first part, we study the conventional robust $M$ -estimation [Huber (1973)] from a new perspective by allowing the robustification parameter $\tau$ to diverge with the sample size to balance the bias and robustness of the estimator. Our main theoretical contribution (Theorem 2.1) is a nonasymptotic Bahadur representation of the proposed robust estimator along with a sub-Gaussian-type deviation bound if the error variable has a finite second moment. As by-products, we prove the Berry-Esseen inequality and a Cramér-type moderate deviation theorem for the estimator. These probabilistic results are particularly useful in investigating robustness and accuracy of the $P$ -values in multiple testing, among other high dimensional statistical inference problems [Fan, Hall and Yao (2007), Delaigle, Hall and Jin (2011), Chang, Shao and Zhou (2016)].

In the second part, we focus on large-scale multiple testing for dependent and heavy-tailed data. To characterize the dependence, we employ a multi-factor model similar to that used in Desai and Storey (2012), Fan, Han and Gu (2012) and Fan and Han (2017) but with an observable factor. To achieve robustness, we propose a Huber loss based approach to construct test statistics for testing the individual hypotheses. Under mild conditions, our procedure asymptotically controls the overall false discovery proportion at the nominal level. Thorough numerical results on both simulated and real world datasets are also provided to back up our theory. It is shown that the newly proposed robust dependence-adjusted method performs well numerically in terms of both the size and power. It significantly outperforms the multiple $t$ -tests under strong dependence, and is applicable even when the true error distribution deviates wildly from the normal distribution. A more interesting and challenging problem is when the dependence structure is characterized by latent factors. In this case, robust estimators of the unobservable factors along with the loadings are required. Large-scale simultaneous inference for latent factor models with heavy-tailed errors is our ongoing work. We leave the details of the results elsewhere in the future.

References

Appendix A Proofs of the results in Section 2

In this section, we present the proofs to the theoretical results from Section 2. Throughout, we use $C,C_{1},C_{2},\ldots$ and $c,c_{1},c_{2},\ldots$ to denote positive constants independent of $n$ and $p$ , which may take different values at each occurrence. First, we collect two useful propositions in Section A.1.

Proposition A.1 reveals that the approximation error vanishes as $\tau$ diverges.

Under Condition 2.1, it holds as long as $\tau\geq 8K_{1}^{2}\,\sigma$ that

where $\widetilde{\bm{\theta}}_{\tau}=\lambda\bm{\theta}^{*}+(1-\lambda)\bm{\theta}^{*}_{\tau}$ for some $0\leq\lambda\leq 1$ and $\nu_{i}=\sigma(\bm{X}_{i})\varepsilon_{i}$ . Moreover, write

provided that $\tau\geq 8K_{1}^{2}\,\sigma$ . This, together with (A.2) and (A.3), proves (A.1). ∎

The following lemma is borrowed from Fan et al. (2015). It provides a localized analysis mechanism, which is a crucial element in the proof of Theorem 2.1.

A.2 Proof of Theorem 2.1

Throughout, let $c_{0},c_{1},c_{2},\ldots$ be positive constants depending only on $\tau_{0}$ , $K_{0}$ and $\sigma^{2}$ . Moreover, we write

where the last step follows from the first order condition that $\nabla L(\widehat{\bm{\theta}})=0$ . By the mean value theorem, $\nabla L(\widehat{\bm{\theta}}_{\eta})-\nabla L(\bm{\theta}^{*})=\nabla^{2}L(\widetilde{}\bm{\theta}_{\eta})(\widehat{\bm{\theta}}_{\eta}-\bm{\theta}^{*})$ , where $\widetilde{}\bm{\theta}_{\eta}$ is a convex combination of $\bm{\theta}^{*}$ and $\widehat{\bm{\theta}}_{\eta}$ and thus satisfies $\|\mathbf{S}^{1/2}(\widetilde{}\bm{\theta}_{\eta}-\bm{\theta}^{*})\|\leq r$ . Together with (A.6), this indicates that

To bound the left-hand side of (A.7) from below, note that

where $\mathbf{A}=(\mathbf{a}_{1},\ldots,\mathbf{a}_{n})^{\intercal}$ with $\mathbf{a}_{i}=\bm{\zeta}_{i}1(|\nu_{i}|+\|\bm{\zeta}_{i}\|r\leq\tau)$ . For any $w>0$ , applying Remark 5.40 in Vershynin (2012) to $\mathbf{A}$ yields that, with probability greater than $1-2e^{-w}$ ,

It then follows that with probability greater than $1-2e^{-w}$ ,

Putting the above calculations together, we conclude that with probability at least $1-2e^{-w}$ ,

as long as $n\geq 16c_{0}^{2}(d+1+w)$ and $\tau\geq 8K_{1}^{2}\{\sigma+(d+1)^{1/2}r\}$ .

Next we bound the quadratic form $\|\mathbf{S}^{-1/2}\nabla L(\bm{\theta}^{*})\|$ . Define

where $\nu_{0}\geq 1$ is a constant depending on $\tau_{0}$ , $K_{0}$ and $\|\bm{\Delta}\|$ . Hence, condition (1.18) in the supplement of Spokoiny (2012) holds with $V_{0}=\mathbf{I}_{d+1}$ and $\mathbf{g}=\sqrt{2(d+1+w)}$ . Moreover, put $D_{0}=\mathbf{I}_{d+1}$ so that $D_{0}^{-1}V_{0}^{2}D_{0}^{-1}=\mathbf{I}_{d+1}$ . Applying Corollary 1.13 there implies that for any $\sqrt{2(d+1)}/18<x\leq x_{c}$ ,

where $x_{c}=(1-0.5\log 3)(d+1)+1.5w\geq 0.45(d+1)+1.5w$ . In particular, taking $x=(d+1)/3+w$ in the preceding inequality we obtain that, with probability greater than $1-5e^{-w}$ ,

Together, the last two displays imply that, with probability at least $1-5e^{-w}$ ,

Taking $r_{1}=4.1r_{0}$ , then it follows from (A.7)–(A.9) that with probability at least $1-7e^{-w}$ , $\|\mathbf{S}^{1/2}(\widehat{\bm{\theta}}_{\eta}-\bm{\theta}^{*})\|\leq 4r_{0}<r_{1}$ whenever $n\geq c_{1}(d+w)^{3/2}$ . By the definition of $\widehat{\bm{\theta}}_{\eta}$ in the beginning of the proof, we must have $\eta=1$ and thus (2.2) follows.

This verifies Condition $(\mathcal{L}_{0})$ by taking

where $\widetilde{\nu}_{0}>0$ is a constant depending only on $K_{0}$ . This verifies Condition ( $ED_{2}$ ) by taking $\omega=n^{-1/2}$ , $\nu_{0}=\widetilde{\nu}_{0}$ and $\mathbf{g}(r)=c_{2}\sqrt{n}$ for all $r>0$ . Then, applying Proposition 3.1 in Spokoiny (2013) with $D_{0}$ replaced by $\mathbf{S}^{1/2}$ yields that, as long as $n\geq c_{2}^{-2}\{4(d+1)+2w\}$ ,

with probability greater than $1-e^{-w}$ , where $\delta(r)$ is as in (A.10). This, together with (2.2) and the fact $\nabla L(\widehat{\bm{\theta}})=0$ , proves (2.3) by taking $r=4r_{0}$ . The proof of Theorem 2.1 is then complete. ∎

A.3 Proof of Theorem 2.2

where $c_{1}>0$ is an absolute constant. By Proposition A.2, $\sigma_{\tau}^{2}\geq\sigma^{2}-\sigma^{4}\tau^{-2}-2(\kappa-2)^{-1}v_{\kappa}\tau^{2-\kappa}$ . Hence, the right-hand side of (A.11) can be further bounded by $2\sqrt{2}\,c_{1}\sigma^{-3}v_{3}\,n^{-1/2}$ , provided that $\tau\geq 2\sigma\vee\{8(\kappa-2)^{-1}\sigma^{-2}v_{\kappa}\}^{1/(\kappa-2)}$ .

Next we prove that the distributions of Gaussian random variables $G_{0}$ and $G_{\tau}$ are also close. Again, using Proposition A.2 we deduce that $|\sigma^{-2}\sigma_{\tau}^{2}-1|\leq\sigma^{-2}\{2(\kappa-2)^{-1}v_{\kappa}\tau^{2-\kappa}+\sigma^{4}\tau^{-2}\}\leq 1/2$ for sufficiently large $\tau$ as above. Then, applying Lemma A.7 in the supplement of Spokoiny and Zhilova (2015) yields

A.4 Proof of Theorem 2.3

Keeping the notations in the proof of Theorem 2.2, we define $T_{0}=\sigma_{\tau}^{-1}W_{0n}$ . First we prove that $T$ and $T_{0}$ are sufficiently close with overwhelmingly high probability. Note that

By Proposition A.2, we have $|\sigma^{2}-\sigma_{\tau}^{2}|\leq 2v_{3}\tau^{-1}$ and $|m_{\tau}|\leq v_{3}\tau^{-2}$ . This, together with Theorem 2.2 yields that

The conclusion (2.4) for $0\leq z\leq 1$ thus follows immediately.

It suffices to prove (2.4) for $z\in[1,o\{\min(\sqrt{w_{n}},\sqrt{n}w_{n}^{-1})\})$ . By (A.13), we have for every $z\geq 1$ and for all sufficiently large $n$ ,

For $|T_{0}|$ , applying Lemma 3.1 in the supplement of Liu and Shao (2014) with $d=1$ , $B_{n}=n$ , $c_{n}\asymp w_{n}^{-1/2}$ , $b_{n}=n^{-1}$ , $d_{n}=n^{-3/10}$ , $t_{n}=(C_{3,1}^{-1/2}\vee 4)(\sqrt{\log n}+B_{n}^{-1/2}x)$ and $x=B_{n}^{1/2}t$ , we deduce that for all sufficiently large $n$ ,

uniformly for $0\leq t\leq c\min(\sqrt{w_{n}},n^{1/6})$ . As a direct consequence, we have

uniformly for $0\leq t\leq c\min(\sqrt{w_{n}},n^{1/6})$ , where $|C_{n,t}|\leq C\{(\sqrt{\log n}+t)^{3}n^{-1/2}+(1+t)n^{-3/10}\}$ . Moreover, note that for any $t>0$ , $t(1+t^{2})^{-1}e^{-t^{2}/2}\leq\sqrt{2\pi}\,\{1-\Phi(t)\}\leq t^{-1}e^{-t^{2}/2}$ . Therefore, for every $z\geq 1$ and all sufficiently large $n$ such that $z-\delta_{n}>0$ ,

Combining (A.14), (A.15) and (A.16) we deduce that the convergence in (2.4) holds uniformly for $1\leq z\leq o\{\min(\sqrt{w_{n}},\sqrt{n}w_{n}^{-1})\}$ , which completes the proof of (2.4). ∎

Appendix B Proofs of the results in Section 3.4

We present here the proofs to the main theorems in Section 3.4, starting with a few essential technical results stated as propositions and proved in Section C below. Throughout, we use $C$ and $c$ to denote positive constants independent of $n$ and $p$ , which may take different values at each occurrence.

Proposition B.1 is a direct consequence of Theorem 2.3. The proof is thus omitted.

with probability at least $1-Cpe^{-w_{n}}$ for all sufficiently large $n$ .

The next two propositions give an uniform law of large numbers for $p_{0}^{-1}\sum_{j\in\mathcal{H}_{0}}1(|T_{j}|\geq z)$ under dependence and independence, respectively, where $\mathcal{H}_{0}=\{j:1\leq j\leq p,\mu_{j}=0\}$ and $p_{0}={\rm Card}(\mathcal{H}_{0})$ .

Assume Conditions (C1)–(C3) hold. Then, for any sequence of positive numbers $m_{p}\leq p$ satisfying $m_{p}\to\infty$ , we have as $(n,p)\to\infty$ ,

Assume Conditions (C1), (C2) and (C5) hold. Then (B.3) remains valid for any sequence of positive numbers $m_{p}\leq p$ satisfying $m_{p}\to\infty$ .

Proofs of Propositions B.2–B.4 are provided in Section C.

B.2 Proof of Theorem 3.1

First, using Propositions B.1, B.2 and the inequality $t(1+t^{2})^{-1}e^{-t^{2}/2}\leq\int_{t}^{\infty}e^{-u^{2}/2}\,du\leq t^{-1}e^{-t^{2}/2}$ for $t>0$ , we deduce that as $n\to\infty$ ,

Let $z_{n}>0$ satisfy $1-\Phi(z_{n})=\alpha/(3p)$ . Then it is easy to see that $z_{n}=\{1+o(1)\}\sqrt{2\log p}$ . This, together with (B.4) yields $F_{j,n}(-z_{n})+1-F_{j,n}(z_{n})=\{1+o(1)\}\{\Phi(-z_{n})+1-\Phi(z_{n})\}=\{1+o(1)\}2\alpha/(3p)$ . On the event $\mathcal{S}_{j}$ , we have $P_{j}^{{\rm true}}>\alpha/p\geq F_{j,n}(-z_{n})+1-F_{j,n}(z_{n})$ and hence $|T_{j}|\leq z_{n}$ for all sufficiently large $n$ . This, together with (B.4) proves (3.13). ∎

B.3 Proof of Theorem 3.2

the conclusion (3.14) follows immediately from Proposition B.3. ∎

B.4 Proof of Theorem 3.3

Recall the definition of $\widehat{z}_{{\rm N},0}$ in (3.12). As in Liu and Shao (2014), using the monotonicity of the indicator function and the continuity of $\Phi$ , it can be shown that

To proceed, we need to derive concentration inequalities for $\widehat{\mu}_{j}$ ’s and $\widehat{\sigma}_{\nu,jj}$ ’s, which follow from Propositions B.1 and B.2. For any $\varepsilon>0$ , define the event

Since $w_{n}\asymp n^{1/5}$ and $\log p=o(n^{1/5})$ , it follows immediately from (B.1) that

Together, (B.5), (B.8) and Proposition B.3 with $m_{p}=\alpha c_{p}$ imply that with probability converging to 1 as $(n,p)\to\infty$ ,

B.5 Proof of Theorem 3.4

Theorem 3.4 is proved similarly to Theorem 3.3, and so is not derived in detail here. The only difference in the argument is to use Proposition B.4 instead of Proposition B.3. ∎

B.6 Proof of Theorem 3.5

The proof follows an argument similar to that in the proof of Theorem 3.3. It suffices to show that Propositions B.1, B.2 and B.3 remain valid for test statistics $S_{j}$ ’s and variance estimators $\widetilde{\sigma}_{\nu,jj}$ ’s. Note that the only difference between $T_{j}$ and $S_{j}$ is on the variance estimation. The former uses $\widehat{\sigma}_{\nu,jj}$ as an estimator of $\sigma_{\nu,jj}=\sigma_{jj}-\bm{b}_{j}^{\intercal}\bm{\Sigma}_{f}\bm{b}_{j}$ , while the latter uses $\widetilde{\sigma}_{\nu,jj}$ defined in (3.6).

Appendix C Proofs of Propositions B.2–B.4

The proof is based on combining exponential bounds for $\widehat{\theta}_{j}$ ’s and $\widehat{\bm{b}}_{j}^{\intercal}\widehat{\bm{\Sigma}}_{f}\widehat{\bm{b}}_{j}$ ’s. For $\widehat{\theta}_{j}$ ’s, using Theorem 5 in Fan, Li and Wang (2017) we deduce that, with $\gamma=\gamma_{n}$ as in Condition (C2),

with probability greater than $1-2pe^{-w_{n}}$ whenever $n\geq 8w_{n}$ . Next, note that for each $j$ ,

Applying Theorem 5.39 in Vershynin (2012) to i.i.d. random vectors $\bm{f}_{0i}=\bm{\Sigma}_{f}^{-1/2}\bm{f}_{i}$ yields that for every $t>0$ , $\|\bm{\Sigma}_{f}^{-1/2}\widehat{\bm{\Sigma}}_{f}\bm{\Sigma}_{f}^{-1/2}-\mathbf{I}_{K}\|\leq\max(\delta,\delta^{2})$ with probability at least $1-2e^{-t}$ , where $\delta=C(K+t)^{1/2}n^{-1/2}$ . Taking $t=w_{n}$ , we have

with probability greater than $1-2e^{-w_{n}}$ for all sufficiently large $n$ .

Together, (C.1)–(C.3) and Theorem 2.1 prove (B.2). ∎

C.2 Proof of Proposition B.3

The proof is based on a discretization technique used to prove Theorem 2.1 in Liu and Shao (2014) and the following results on joint Gaussian approximations.

Under the Conditions (C1)–(C3), we have for any $0<\rho\leq 1$ and $0<\delta<1$ ,

uniformly for $z\in[0,o(\sqrt{w_{n}}))$ and all all $j,k\in\mathcal{H}_{0}$ satisfying $j\neq k$ and $|\rho_{\nu,jk}|\leq\rho$ . In addition, we have for any $A>0$ ,

uniformly for $0\leq z\leq A\sqrt{\log p}$ and all $j,k\in\mathcal{H}_{0}$ satisfying $j\neq k$ and $|\rho_{\nu,jk}|\leq(\log p)^{-2-\kappa}$ , where $G\sim N(0,1)$ and $|C_{n,z}|\leq C\{(\log p)^{1/2}n^{-3/10}+(\log p)^{-1-\kappa/2}\}$ .

$b_{p}=\sqrt{2\log(p/m_{p})}$ , and note that, as $p\to\infty$ ,

Therefore, $\Psi(b_{p})\leq\Psi(z_{0})=t_{0}\leq\Psi(a_{p})$ for all sufficiently large $p$ , which proves the claim.

In view of the last two displays, it is enough to prove that

in probability. For any $\varepsilon>0$ , applying Boole’s and Markov’s inequalities we deduce that

where $\mathcal{I}_{j}=\{k\in\mathcal{H}_{0}:|\rho_{\nu,jk}|\leq(\log p)^{-2-\kappa}\}$ and

On the other hand, for $j\in\mathcal{H}_{0}$ and $k\in\mathcal{I}_{j}^{{\rm c}}\setminus\{j\}$ , using (C.4) and the inequality $\Psi(z)\geq(2/\pi)^{1/2}(1+z^{2})^{-1}ze^{-z^{2}/2}$ for $z\geq 0$ , we deduce that

as $(n,p)\to\infty$ , where $s_{p}$ , $r$ and $\rho$ are defined in Conditions (C3). This proves (C.6), and hence completes the proof of Proposition B.3. ∎

C.3 Proof of Proposition B.4

for all sufficiently large $n$ , where $\delta_{n}=cn^{-1/2}w_{n}$ . Moreover, note that

uniformly in $0\leq z\leq\Phi^{-1}(1-m_{p}/(2p))$ . In view of (C.10) and (C.11), and by the argument leading to (C.6), we only need to prove that as $(n,p)\to\infty$ ,

where $0=z_{d_{p}}<\cdots<z_{1}<z_{0}=\{1+o(1)\}\sqrt{2\log(p/m_{p})}$ are as in the proof of Proposition B.3.

This proves (C.12), and thus completes the proof of (B.3). ∎

C.4 Proof of Lemma C.1

for all sufficiently large $n$ , where $\delta_{n}=cn^{-1/2}w_{n}$ .

For $z\geq 1$ and $j\neq k\in\mathcal{H}_{0}$ , it follows from (C.13) that for all sufficiently large $n$ ,

Also, it follows from Lemma A.2 that $\max_{j\in\mathcal{H}_{0}}\max(|m_{j,\tau}|,|\sigma_{\nu,jj}-\sigma_{j,\tau}^{2}|)\leq Cn^{-1}w_{n}$ . Substituting this into (C.15) gives

Then, applying Bernstein’s inequality to $T_{0j}+T_{0k}$ we deduce that

where $T_{0j}^{-}=-T_{0j}$ and $T_{0k}^{-}=-T_{0k}$ . Applying Theorem 1.1 in Bentkus (2003) to $\mathbf{T}_{0}=n^{-1/2}\sum_{i=1}^{n}\bm{\xi}_{ni}$ and $(T_{0j}^{-},T_{0k}^{-})^{\intercal}=-n^{-1/2}\sum_{i=1}^{n}\bm{\xi}_{ni}$ , we deduce that

where $\mathbf{W}=(W_{1},W_{2})^{\intercal}\sim N(\mathbf{0},\mathbf{A})$ . Further, it can similarly shown that

By (C.16) and the assumption that $w_{n}\asymp n^{-1/5}$ and $\log p=o(n^{1/5})$ , $W_{1}$ and $W_{2}$ are weakly correlated with $\rho=\operatorname{cov}(W_{1},W_{2})\leq C(\log p)^{-2-\kappa}$ . Therefore, for every $0\leq t\leq 1$ and all sufficiently large $n$ , we deduce that

where $\bm{G}=(G_{1},G_{2})^{\intercal}\sim N(\mathbf{0},\mathbf{I}_{2})$ . Consequently, the conclusion (C.5) for $0\leq z\leq 1$ follows from (C.21), (C.22) and (C.23).

Now it remains to consider the case of $z\geq 1$ . By (C.13),

Hence, it follows from Lemma 3.1 in the supplement of Liu and Shao (2014) by taking $d=2$ , $B_{n}=n$ , $c_{n}\asymp w_{n}^{-1/2}$ , $b_{n}\asymp(\log p)^{-2-\kappa}$ , $d_{n}=(\log p)^{-3/2-\kappa/2}$ , $t_{n}=(C_{3,2}^{-1/2}\vee 4)(\sqrt{\log p}+B_{n}^{-1/2}x)$ and $x=B_{n}^{1/2}t$ that

for all $0\leq t\leq c\min\{\sqrt{w_{n}},(\log p)^{(3+\kappa)/2}\}$ with $n$ sufficiently large.

Together, (C.24), (C.25) and (A.16) prove (C.5) for $1\leq z\leq A\sqrt{\log p}$ . ∎

Appendix D Additional simulation results

In this section, we present numerical results comparing the performance of the RD-AN procedure and the OD-A procedure. Again, we consider the three factor model $X_{ij}=\mu_{j}+\bm{b}_{j}^{\intercal}\bm{f}_{i}+u_{ij}$ for $i=1,\ldots,n$ , where $\bm{u}_{i}=(u_{i1},\ldots,u_{ip})^{\intercal}$ are i.i.d. copies of $\bm{u}=(u_{1},\ldots,u_{p})^{\intercal}$ . We simulate $\{\bm{f}_{i}\}_{i=1}^{n}$ from $N_{3}(\mathbf{0},\bm{\Sigma}_{f})$ as in the main text; independently, we generate the loadings $\{\bm{b}_{j}=(b_{j1},b_{j2},b_{j3})^{\intercal}\}_{j=1}^{p}$ according to $b_{j1},b_{j3}\overset{{\rm i.i.d.}}{\sim}\mbox{Uniform}(0.5,1.5)$ and $b_{j2}\overset{{\rm i.i.d.}}{\sim}\mbox{Uniform}(-2,-1)$ . The errors $\{\bm{u}_{i}\}_{i=1}^{n}$ are generated independently from the following distributions:

Model 5. $\bm{u}\sim(1/\sqrt{2})\,t_{4}(\mathbf{0},\bm{\Sigma}_{u})$ ;

To use normal calibration, we need to estimate the variance $\sigma_{jj}=\mbox{var}(X_{j})$ . We compute the median-of-means estimator $\widetilde{\sigma}_{jj}(V)$ with a universal parameter $V=\lceil 0.5\log(pn)\rceil$ for all $j$ . Since we need to subtract the common variance estimator according to (3.6), the final estimate of $\mbox{var}(u_{j})$ may sometimes be too close to zero. To slightly improve numerical performance, we proceed as follows: (i) For each $j$ , compute $\widetilde{\sigma}_{jj}(v)$ for $v=1,\ldots,V$ ; (ii) remove those $\widetilde{\sigma}_{jj}(v)$ ’s that are smaller than $\widehat{\bm{b}}_{j}^{\intercal}\widehat{\bm{\Sigma}}_{f}\widehat{\bm{b}}_{j}$ ; (iii) taking the 0.75 quantile of the remaining $\widetilde{\sigma}_{jj}(v)$ ’s as the final estimate of $\sigma_{jj}$ . The numerical results indicate that this modified procedure is numerically stable.

In the simulations reported here, we take $p=2000$ , $n=80,120$ , $\mu_{j}=\mu$ for $1\leq j\leq\pi_{1}p$ and $\mu_{j}=0$ otherwise, where $\mu=\sqrt{2(\log p)/n}$ and $\pi_{1}=0.25$ . For simplicity, we set $\lambda=0.5$ in our procedure and use the Matlab package mafdr to compute the estimate $\widehat{\pi}_{0}(\lambda)$ of $\pi_{0}=1-\pi_{1}$ . All the results are based on 500 simulation rounds. From Tables 4 and 5 we see that, in the presence of heavy-tailed errors, the RD-A procedure provides consistently better controls of the FDR at the expense of slight compromises of the FNR and TPR. Figure 4 compares the performance of the three methods, for various sample sizes and signal strengths, in the cases of Model 2 and Model 6. Wee see that the RD-A consistently outperforms the two other methods, across all sample sizes and even for low signal strengths when all the methods exhibit higher errors.