Feature selection by Higher Criticism thresholding: optimal phase diagram

Introduction

The modern era of high-throughput data collection creates data in abundance. Some devices – spectrometers and gene chips come to mind – automatically generate measurements on thousands of standard features of each observational unit.

What hasn’t changed in science is the difficulty of obtaining good observational units. High-throughput devices don’t help us to find and enroll qualified patients, or catch exotic butterflies, or observe primate mating behaviors. Hence, in many fields the number of observational units – eg. patients, butterflies, or matings – has not increased, and stays today in the dozens or hundreds. But each of those few observational units can now be subjected to a large battery of automatic feature measurements

Many of those automatically measured features will have little relevance to any given project. This new era of feature glut poses a needle-in-a-haystack problem: we must detect a relatively few valuable features among many useless ones. Unfortunately, the combination of small sample sizes (few observational units) and high dimensions (many feature measurements) makes it hard to tell needles from straw.

Orthodox statistical methods assumed a quite different set of conditions: more observations than features, and all features highly relevant. Modern statistical research is intensively developing new tools and theory to address the new unorthodox setting; such research comprised much of the activity in the recent 6-month Newton Institute program Statistical Theory and Methods for Complex, High-Dimensional Data.

In this article we focus on this new setting, this time addressing the challenges that modern high-dimensional data pose to linear classification schemes. New data analysis tools and new types of mathematical analysis of those tools will be introduced.

Consider a simple model of linear classifier training. We have a set of labelled training samples $(Y_{i},X_{i})$, $i=1,\dots,n$, where each label $Y_{i}$ is $\pm 1$ and each feature vector $X_{i}\in R^{p}$. For simplicity, we assume the training set contains equal numbers of $1$’s and $-1$’s and that the feature vectors $X_{i}\in R^{p}$ obey $X_{i}\sim N(Y_{i}\mu,\Sigma)$, $i=1,\dots,n$, for an unknown mean contrast vector $\mu\in R^{p}$; here $\Sigma$ denotes the feature covariance matrix and $n$ is the training set size. In this simple setting, one ordinarily uses linear classifiers to classify an unlabeled test vector $X$, taking the general form $L(X)=\sum_{j=1}^{p}w(j)X(j)$, for a sequence of ‘feature weights’ $w=(w(j):j=1,\dots,p)$.

Classical theory going back to RA Fisher shows that the optimal classifier has feature weights $w\propto\Sigma^{-1}\mu$; at first glance linear classifier design seems straightforward and settled. However, in many of today’s most active application areas, it is a major challenge to construct linear classifiers which work well.

2 p𝑝p larger than n𝑛n

The culprit can be called the “$p>n$ problem”. A large number $p$ of measurements is automatically made on thousands of standard features, but in a given project, the number of observational units, $n$, might be in the dozens or hundreds. The fact that $p\gg n$ makes it difficult or impossible to estimate the feature covariance $\Sigma$ with any precision.

It is well known that naive application of the formula $w\propto\Sigma^{-1}\mu$ to empirical data in the $p>n$ setting is problematic; at a minimum, because the matrix of empirical feature covariances $\widehat{Cov}_{n,p}(X)$ is not invertible. But even if we use the generalized inverse $\widehat{Cov}_{n,p}(X)^{\dagger}$, the resulting naive classification weights, $\hat{w}_{naive}\propto\widehat{Cov}_{n,p}(X)^{{\dagger}}\widehat{Cov}_{n,p}(Y,X)$, often give very “noisy” classifiers with low accuracy. The modern feature glut thus seriously damages the applicability of ‘textbook’ approaches.

A by-now standard response to this problem is to simply ignore feature covariances, and standardize the features to mean zero and variance one. One in effect pretends that the feature covariance matrix $\Sigma$ is the identity matrix, and uses the formula $w(j)\propto Cov(Y,X(j))$ . Even after this reduction, further challenges remain.

3 When features are rare and weak

In fields such as genomics, with a glut of feature measurements per observational unit, it is expected that few measured features will be useful in any given project; nevertheless, they all get measured, because researchers can’t say in advance which ones will be useful. Moreover, reported misclassification rates are relatively high. Hence there may be numerous useful features, but they are relatively rare and individually quite weak.

Such thinking motivated the following rare/weak feature model (RW Feature Model) in . Under this model,

Useful features are rare: the contrast vector $\mu$ is nonzero in only $k$ out of $p$ elements, where $\epsilon=k/p$ is small, i.e. close to zero. As an example, think of $p=10,000$, $k=100$, so $\epsilon=k/p=.01$. In addition,

Useful features are weak: the nonzero elements of $\mu$ have common amplitude $\mu_{0}$, which is assumed not to be ‘large’. ‘Large’ can be measured using $\tau=\sqrt{n}\mu_{0}$; values of $\tau$ in the range $2$ to $4$ imply that corresponding values of $\mu$ are not large.

Since the elements $X(j)$ of the feature vector where the class contrast $\mu(j)=0$ are entirely uninformative about the value of $Y(j)$, only the $k$ features where $\mu(j)=\mu_{0}$ are useful. The problem is how to identify and benefit from those rare, weak features.

We speak of $\epsilon$ and $\tau$ as the sparsity and strength parameters for the Rare/Weak model, and refer to the $RW(\epsilon,\tau)$ model. Models with a ‘sparsity’ parameter $\epsilon$ are common in estimation settings , but not with the feature strength constraint $\tau$. Also closely related to the RW model is work in multiple testing by Ingster and the authors .

4 Feature selection by thresholding

Feature selection - i.e. working only with an empirically-selected subset of features - is a standard response to feature glut. We are supposing, as announced in Section 1.2, that feature correlations can be ignored and that features are already standardized to variance one. We therefore focus on the vector of feature $Z$-scores, with components $Z(j)=n^{-1/2}\sum_{i}Y_{i}X_{i}(j)$, $j=1,\dots,p$. These are the $Z$-scores of two-sided tests of $H_{0,j}$: $Cov(Y,X(j))=0$. Under our assumptions $Z\sim N(\theta,I_{p})$ with $\theta=\sqrt{n}\mu$ and $\mu$ the feature contrast vector. Features with nonzero $\mu(j)$ typically have significantly nonzero $Z(j)$ while all other features will have $Z(j)$ values largely consistent with the null hypothesis $\mu(j)=0$. In such a setting, selecting features with $Z$-scores above a threshold makes sense.

We identify three useful threshold functions: $\eta_{t}^{\star}(z)$, $\star\in\{clip,hard,soft\}$. These are: Clipping – $\eta^{clip}_{t}(z)=\mbox{sgn}(z)\cdot 1_{\{|z|>t\}}$, which ignores the size of the $Z$-score, provided it is large; Hard Thresholding – $\eta^{hard}_{t}(z)=z\cdot 1_{\{|z|>t\}}$, which uses the size of the $Z$-score, provided it is large; and Soft Thresholding – $\eta^{soft}_{t}(z)=\mbox{sgn}(z)(|z|-t)_{+}$, which uses a shrunken $Z$-score, provided it is large.

Let $\star\in\{soft,hard,clip\}$. The threshold feature selection classifier makes its decision based on $L_{t}^{\star}<>0$ where $\hat{L}_{t}^{\star}(X)=\sum_{j=1}^{p}\hat{w}^{\star}_{t}(j)X(j)$, and $\hat{w}^{\star}_{t}(j)=\eta^{\star}_{t}(Z(j)),j=1,\dots,p$.

In words, the classifier sums across features with large training-set $Z$-scores, and a simple function of the $Z$-score generates the feature weight.

Several methods for linear classification in bioinformatics follow this approach: the Shrunken Centroids method is a variant of soft thresholding in this two-class setting; the highly-cited methods in and are variants of hard thresholding. Clipping makes sense in the theoretical setting of the RW model (since then the useful features have all the same strength) and is simpler to analyse than the other nonlinearities.

Thresholding has been popular in estimation for more than a decade ; it is known to be succesful in ‘sparse’ settings where the estimand has many coordinates, of which only a relatively few coordinates are significantly nonzero. However classification is not the same as estimation, and performance characteristics are driven by quite different considerations.

One crucial question remains: how to choose the threshold based on the data? Popular methods for threshold choice include cross-validation ; control of the false discovery rate ; and control of the local false discovery rate .

5 Higher Criticism

In we proposed a method of threshold choice based on recent work in the field of multiple comparisons. We now very briefly mention work in that field and then introduce the threshold choice method.

(HC Testing) . The Higher Criticism objective is

Fix $\alpha_{0}\in(0,1)$ (eg $\alpha_{0}=1/10$). The HC test statistic is $HC^{*}=\max_{1\leq i\leq\alpha_{0}N}HC(i;\pi_{(i)})$.

HC seems insensitive to the selection of $\alpha$, in Rare/Weak situations; here we always use $\alpha_{0}=.10$.

In words, we look for the largest standardized discrepancy for any $\pi_{(i)}$ between the observed behavior and the expected behavior under the null. When this is large, the whole collection of $P$-values is not consistent with the global null hypothesis. The phrase “Higher Criticism” is due to John Tukey, and reflects the shift in emphasis from single test results to the whole collection of tests; see discussion in . Note: there are several variants of HC statistic; we discuss only one variant in this brief note; the main results of still apply to this variant; for full discussion see .

5.2 HC thresholding

Return to the classification setting in previous sections. We have a vector of feature $Z$-scores $(Z(j),j=1,\dots,p)$. To apply HC notions, translate $Z$-scores into two-sided $P$-values, and maximizes the HC objective over index $i$ in the appropriate range. Define the feature $P$-values $\pi_{i}=Prob\{|N(0,1)|>|Z(i)|\}$, $i=1,\dots,p$; and define the increasing rearrangement $\pi_{(i)}$, the HC objective function $HC(i;\pi_{(i)})$, and the increasing rearrangement $|Z|_{(i)}$ correspondingly. Here is our proposal.

(HC Thresholding). Apply the HC procedure to the feature $P$-values. Let the maximum HC objective be achieved at index $\hat{i}$. The Higher Criticism threshold (HCT) is the value $\hat{t}^{HC}=|Z|_{(\hat{i})}$. The HC threshold feature selector selects features with $Z$-scores exceeding $\hat{t}^{HC}$ in magnitude.

Figure 1 illustrates the procedure. Panel (a) shows a sample of $Z$-scores, Panel (b) shows a PP-plot of the corresponding ordered $P$-values versus $i/p$ and Panel (c) shows a standardized PP-plot. The standardized PP-Plot has its largest deviation from zero at $\hat{i}$; and this generates the threshold value.

5.3 Previously-reported results for HCT

Our article reported several findings about behavior of HCT based on numerical and empirical evidence. In the RW model, we can define an ideal threshold, i.e. a threshold based on full knowledge of the RW parameters $\epsilon$ and $\mu$ and chosen to minimize the misclassification rate of the threshold classifier – see Section 2 below. We showed in that:

HCT gives a threshold value which is numerically very close to the ideal threshold.

In the case of very weak feature z-scores, HCT has a False Feature Discovery Rate (FDR) substantially higher than other popular approaches, but a Feature Missed Detection Rate (MDR) substantially lower than those other approaches;

At the same time, HCT has FDR and MDR very closely matching those of the ideal threshold.

In short, HCT has very different operating characteristics from those other thresholding schemes like FDR thresholding and Bonferroni thresholding, but very similar operating characteristics to the ideal threshold.

6 Asymptotic RW model, and the phase diagram

In this paper we further support the findings reported in , this time using an asymptotic analysis. In our analysis the number of observations $n$ and the number of features $p$ tend to infinity in a linked fashion, with $n$ remaining very small compared to $p$. (Empirical results in show our large-$p$ theory is applicable at moderate $n$ and $p$).

More precisely, we consider a sequence of problems with increasingly more features, increasingly more rare useful features, and relatively small numbers of observations compared to the number of features.

The phrase asymptotic RW model refers to the following combined assumptions.

Asymptotic Setting. We consider a sequence of problems, where the number of observations $n$ and the number of features $p$ both tend to $\infty$ along the sequence.

$p$ dramatically larger than $n$. Along this sequence, $n\sim c\cdot\log(p)^{\gamma}$, so there are dramatically more features per observational unit than there are observational units.

Increasing Rarity. The sparsity $\epsilon$ varies with $n$ and $p$ according to $\epsilon=p^{-\beta}$, $0<\beta<1$.

Decreasing Strength. The strength $\tau$ varies with $n$ and $p$ according to $\tau=\sqrt{2r\log(p)}$, $0<r<1$.

The symbol $ARW(r,\beta;c,\gamma)$ refers to the model combining these assumptions.

In this model, because $r<1$, useful features are individually too weak to detect and because $0<\beta<1$, useful features are increasingly rare with increasing $p$, while increasing in total number with $p$. It turns out that $c$ and $\gamma$ are incidental, while $r$ and $\beta$ are the driving parameters. Hence we always simply write $ARW(r,\beta)$ below.

There is a large family of choices of $(r,\beta)$ where successful classification is possible, and another large family of choices where it is impossible. To understand this fact, we use the concept of phase space, the two-dimensional domain $0<r,\beta<1$. We show that this domain is partitioned into two regions or ‘phases’. In the “impossible” phase, useful features are so rare and so weak that classification is asymptotically impossible even with the ideal choice of threshold. In the “possible” phase, successfully separating the two groups is indeed possible - if one has access to the ideal threshold. Figure 2 displays this domain and its partition into phases. Because of the partition into two phases, we also call this display the phase diagram. An explicit formula for the graph $r=\rho^{*}(\beta)$ bounding these phases is given in (3.2) below.

The phase diagram provides a convenient platform for comparing different procedures. A threshold choice is optimal if it gives the same partition of phase space as the one obtained with the ideal choice of threshold.

How does HCT compare to the ideal threshold, and what partition in the phase space does HCT yield? For reasons of space, we focus in this paper on the Ideal HC threshold, which is obtained upon replacing the empirical distribution of feature $Z$-scores by it expected value. The Ideal HC threshold is thus the threshold which HCT is ‘trying to estimate’; in the companion paper we give a full analysis showing that the ideal HCT and HCT are close.

The central surprise of our story is that HC behaves surprisingly well: the partition of phase space describing the two regions where ideal thresholding fails and/or succeeds also describes the two regions where Ideal HCT fails and/or succeeds in classifying accurately. The situation is depicted in the table below:

Here by ‘succeeds’, we mean asymptotically zero misclassification rate and by ‘fails’, we mean asymptotically 50% misclassification rate.

In this sense of size of regions of success, HCT is just as good as the ideal threshold. Such statements cannot be made for some other popular thresholding schemes, such as False Discovery threshold selection. As will be shown in even the very popular Cross-Validated choice of Threshold will fail if the training set size is bounded, while HCT will still succeed in the RW model in that case.

The full proof of a broader set of claims – with a considerably more general treatment – will appear elsewhere. To elaborate the whole story on HCT needs three connected papers including , , and the current one. In , we reported numerical results both with simulated data from the RW model and with certain real data often used as standard benchmarks for classifier performance. In , we will develop a more mathematical treatment of many results we cite here and in . The current article, logically second in the triology, develops an analysis of Ideal HCT which is both transparent and which provides the key insights underlying our lengthy arguments in . We also take the time to explain the notions of phase diagram and phase regions. We believe this paper will be helpful to readers who want to understand HCT and its performance, but who would be overwhelmed by the epsilontics of the analysis in .

The paper is organized as follows. Section 2 introduces a functional framework and several ideal quantities. These include the proxy classification error where Fisher’s separation (SEP) plays a key role, the ideal threshold as a proxy for the optimal threshold, and the ideal HCT as a proxy of the HCT. Section 3 introduces the main results on the asymptotic behavior of the HC threshold under the asymptotic RW model, and the focal point is the phase diagram. Section 4 outlines the basic idea behind the main results followed by the proofs. Section 5 discuss the connection between the ideal threshold and the ideal HCT. Section 6 discusses the ideal behavior of Bonferroni threshold feature selection and FDR-controlling feature selection. Section 7 discusses the link between ideal HCT and ordinary HCT, the finite $p$ phase diagram, and other appearances of HC in recent literature.

Sep functional and ideal threshold

Suppose $L$ is a fixed, nonrandom linear classifier, with decision boundary $L<>0$. Will $L$ correctly classify the future realization $(Y,X)$ from simple model of Section 1.3. Then $Y=\pm 1$ equiprobable and $X\sim N(Y\mu,I_{p})$. The misclassification probability can be written

where $\Phi$ denotes the standard normal distribution function and $Sep(L;\mu)$ measures the standardized interclass distance:

The ideal linear classifier $L_{\mu}$ with feature weights $w\propto\mu$, and decision threshold $L_{\mu}<>0$ implements the likelihood ratio test. It also maximizes $Sep$, since for every other linear classifier $L$, $Sep(L;\mu)\leq Sep(L_{\mu};\mu)=2\|\mu\|_{2}$.

Threshold selection rules give random linear classifiers: the classifier weight vector $w$ is a random variable, because it depends on the $Z$-scores of the realized training sample. If $L_{Z}$ denotes a linear classifier constructed based on such a realized vector of $Z$-scores, then the misclassification error can be written as

this is a random variable depending on $Z$ and on $\mu$. Heuristically, because there is a large number of coordinates, some statistical regularity appears, and we anticipate that random quantities can be replaced by expected values. We proceed as if

where the expectation is over the conditional distribution of $Z$ conditioned on $\mu$. Our next step derives from the fact that $\mu$ itself is random, having about $\epsilon\cdot p$ nonzero coordinates, in random locations. Now as $w_{i}=\eta_{t}(Z_{i})$ we write heuristically

Let the threshold $t$ be fixed and chosen independently of the training set. In the $RW(\epsilon,\tau)$ model we use the following expressions for proxy separation

and $W$ denotes a standard normal random variable. By proxy classification error we mean

Normalizations are chosen here so that, in large samples

While ordinarily, we expect averages to “behave like expectations” in large samples, we use the word proxy to remind us that there is a difference (presumably small). Software to compute these proxy expressions has been developed by the authors, and some numerical results using them were reported in .

Of course, the rationale for our interest in these proxy expressions is our heuristic understanding that they accurately describe the exact large-sample behavior of certain threshold selection schemes. This issue is settled in the affirmative, after considerable effort, in .

2 Certainty-equivalent threshold functionals

In general the best threshold to use in a given instance of the RW model depends on both $\epsilon$ and $\tau$. It also depends on the specific realization of $\mu$ and even of $Z$. However, dependence on $\mu$ and $Z$ is simply “noise” that goes away in large samples, while the dependence on $\epsilon$ and $\tau$ remains.

The ideal threshold functional $T_{ideal}(\epsilon,\tau)$ maximizes the proxy separation

Heuristically, $T_{ideal}$ represents a near-optimal threshold in all sufficiently large samples; it is what we “ought” to be attempting to use.

Folding. The following concepts and notations will be used in connection with distributions of absolute values of random variables.

The Half Normal distribution function $\Psi(t)=P\{|N(0,1)|\leq t\}$.

The noncentral Half-Normal distribution $\Psi_{\tau}(t)=P\{|N(\tau,1)|\leq t\}$.

Given a distribution function $F$, the folded distribution is $G(t)=F(t)-F(-t)$. The Half Normal is the folded version of the standard Normal, and the noncentral Half Normal is the folded version of a Normal with unit standard deviation and nonzero mean equal to the noncentrality parameter.

Let $F_{\epsilon,\tau}$ denote the 2-point mixture

$G_{\epsilon,\tau}$ denotes the corresponding folded distribution:

We now define an HCT functional representing the target that HC thresholding aims for.

Let $F$ be a distribution function which is not the standard normal $\Phi$. At such a distribution, we define the HCT functional by

here $G$ is the folding of $F$, and $t_{0}$ is a fixed parameter of the HCT method (eg. $t_{0}=\Phi^{-1}(0.1)$). The HC threshold in the $RW(\epsilon,\tau)$ model may be written, in an abuse of notation,

Let $F_{n,p}$ denote the usual empirical distribution of the feature $Z$-scores $Z_{i}$. The HCT of Definition 1.3 can be written as $\hat{t}_{n,p}^{HC}=T_{HC}(F_{n,p})$. Let $F$ denote the expected value of $F_{n,p}$; then $T_{HC}(F)$ will be called the ideal HC threshold. Heuristically, we expect the usual sampling fluctuations and that

with a discrepancy decaying as $p$ and $n$ increase. This issue is carefully considered in the companion paper , which shows that the empirical HC threshold in the ARW model indeed closely matches the ideal HC threshold.

For comparison purposes, we considered two other threshold schemes. First, (ideal) False-Discovery Rate thresholding. For a threshold $t$, and parameters $(p,\epsilon,\tau)$, the expected number of useful features selected is

and the expected number of useless features selected is

Let $TPR(t)=p^{-1}E(TP)(t)$ denote the expected rate of useful features above threshold and $FPR(t)=p^{-1}E(FP)(t)$ denote the expected rate of usless features above threshold. In analogy with our earlier heuristic, we define the proxy False Discovery Rate (FDR)

although for large $p$ the difference will often be small.)

We define the FDRT-$\alpha$ functional by

Heuristically, this is the threshold that FDRT is ‘trying’ to learn from noisy empirical data. We will also need the proxy Local FDR.

Here $FPR^{\prime}$ denotes the derivative of $FPR$, which exists, using smoothness properties of $\Psi_{0}$; similarly for $TPR^{\prime}$ and $\Psi_{\tau}$. Intuitively, $\widetilde{Lfdr}(t)$ denotes the expected fraction of useless features among those features having observed $Z$-scores near level $t$.

Second, we considered Bonferroni-based thresholding.

This threshold level is set at the level that would cause on average one false alarm in a set of $p$ null cases.

In [11, Figures 2-3], we presented numerical calculations of all these functionals and their separation behavior in two cases.

Although our calculations are exact numerical finite-$p$ calculations, we remark that they correspond to sparsity exponents $\beta=1/2$ and $\beta=2/3$, respectively. The figures show the following.

There is a very close numerical approximation of the HCT to the ideal threshold, not just at large $\tau$ but also even at quite small $2<\tau<3$.

FDR and Bonferroni thresholds behave very differently from the ideal and from HC.

The separation behavior of the HCT is nearly ideal. For the constant FDR rules, the separation behavior is close to ideal at some $\tau$ but becomes noticeably sub-ideal at other $\tau$.

The False discovery rate behavior of HCT and Ideal thresholding depends on $\tau$. At small $\tau$, both rules tolerate a high FDR while at large $\tau$, both rules obtain a small FDR.

The Missed detection rate of HCT and Ideal thresholding also depends on $\tau$. At small $\tau$, the missed detection rate is high, but noticeably less than 100%. At large $\tau$, the missed detection rate falls, but remains noticeably above $0\%$. In contrast the MDR for FDR procedures is essentially 100% for small $\tau$ and falls below that of HCT/ideal for large $\tau$.

These numerical examples illustrate the idealized behavior of different procedures. We can think of the HCT functional as the threshold which is being estimated by the actual HCT rule. On an actual dataset sampled from the underlying $F$, the HC threshold will behave differently, primarly due to stochastic fluctuations $F_{n,p}\approx F$. Nevertheless, the close approximation of the HCT threshold to the ideal one is striking and, to us, compelling.

Behavior of ideal threshold, asymptotic RW model

We now study the ideal threshold in the asymptotic RW model of Definition 1.4. That is, we fix parameters $r$,$\beta$ in that model and study the choice of threshold $t$ maximizing class separation.

We first make precise a structural fact about the ideal threshold, first observed informally in .

ROC Curve. The feature detection receiver operating characteristic curve (ROC) is the curve parameterized by $(FPR(t),TPR(t))$. The tangent to this curve at $t$ is

Note that in the $RW(\epsilon,\tau)$ model, $TPR$, $FPR$, $tan$ and $sec$ all depend on $t$, $\epsilon$,$\tau$ and $p$, although we may, as here, indicate only dependence on $t$.

Tangent-Secant Rule. In the $ARW(r,\beta)$ model. we have

Here $\epsilon=p^{-\beta}$, $\tau=\sqrt{2r\log(p)}$ and $n\sim c\log(p)^{\gamma}$, $p\rightarrow\infty$ as in Definition 1.4.

Success Region. The region of asymptotically successful ideal threshold feature selection in the $(\beta,r)$ plane is the interior of the subset where the ideal threshold choice $T_{ideal}(\epsilon,\tau)$ obeys

here we are in the $ARW(r,\beta)$ model of Definition 1.4.

The interesting range involves $(\beta,r)\in^{2}$. The following function is important for our analysis, and has previously appeared in important roles in other (seemingly unrelated) problems; see Section 7.4.

As it turns out, it marks the boundary between success and failure for threshold feature selection.

Existence of Phases. The success region is precisely $r>\rho^{*}(\beta)$, $0<\beta<1$. In the interior of the complementary region $r<\rho^{*}(\beta)$, $1/2<\beta<1$, even the ideal threshold cannot send the proxy separation to infinity with increasing $(n,p)$.

Regions I,II, III. The Success Region can be split into three regions, referred to here and below as Regions I-III. The interiors of the regions are as follows:

$\beta-1/2<r\leq\beta/3$ and $1/2<\beta<3/4$; $r>\rho^{*}(\beta)$.

$\beta/3<r\leq\beta$ and $1/2<\beta<1$; $r>\rho^{*}(\beta)$.

$\beta<r<1$ and $1/2<\beta<1$; $r>\rho^{*}(\beta)$.

In the asymptotic RW model, the optimal threshold must behave asymptotically like $\sqrt{2q\log(p)}$ for a certain $q=q(r,\beta)$. Surprisingly we need not have $q=r$.

Formula for Ideal Threshold. Under the Asymptotic RW model $ARW(r,\beta)$, with $r>\rho^{*}(\beta)$, the ideal threshold has the form $T_{ideal}(\epsilon,\tau)\sim\sqrt{2q^{*}\log(p)}$ where

Note in particular that in Regions I and II, $q^{*}>r$, and hence $T_{ideal}(\epsilon,\tau)>\tau$. Although the features truly have strength $\tau$, the threshold is best set higher than $\tau$.

We now turn to FDR properties. The Tangent-Secant rule implies immediately

Hence any result about FDR is tied to one about local FDR, and vice versa.

Under the Asymptotic RW model $ARW(r,\beta)$, at the ideal threshold $T_{ideal}(\epsilon,\tau)$ proxy FDR obeys

as $p\rightarrow\infty$, and the proxy local FDR obeys

Several aspects of the above solution are of interest.

Threshold Elevation. The threshold $\sqrt{2q^{*}\log(p)}$ is significantly higher than $\sqrt{2r\log(p)}$ in Regions I and II. Instead of looking for features at the amplitude they can be expected to have, we look for them at much higher amplitudes.

Fractional Harvesting. Outside of Region III, we are selecting only a small fraction of the truly useful features.

False Discovery Rate. Outside Region III, we actually have a very large false discovery rate, which is very close to $1$ in Region I. Surprisingly even though most of the selected features are useless, we still correctly classify!

Training versus Test performance. The quantity $\sqrt{q}$ can be interpreted as a ratio: $\sqrt{q^{*}/r}=\min\{\frac{\beta+r}{2r},2\}=$ strength of useful features in training / strength of those features in test. From (3.3) we learn that, in Region I, the selected useful features perform about half as well in training as we might expect from their performance in test.

Behavior of ideal clipping threshold

We now sketch some of the arguments involved in the full proof of the theorems stated above. In the RW model, it makes particular sense to use the clipping threshold function $\eta_{t}^{clip}$, since all nonzeros are known to have the same amplitude. The ideal clipping threshold is also very easy to analyze heuristically. But it turns out that all the statements in Theorems 1-3 are equally valid for all three types of threshold functions, so we prefer to explain the derivations using clipping.

In the RW model, we can express the components of the proxy separation very simply when using the clipping threshold:

where $W$ denotes an $N(0,1)$ random variable. Recall the definitions of useful selections TP and useless selections FP; we also must count Inverted Detections, for the case where the $\mu_{i}>0$ but $\eta_{t}^{clip}(Z_{i})<0$. Put

with again $\Phi$ the standard normal distribution, and define the inverted detecton rate by $IDR=p^{-1}E(ID)$. Then

We arrive at an identity for $\widetilde{Sep}$ in the case of clipping:

We now explain Theorem 1, the Tangent-Secant rule. Consider the alternate proxy

i.e. drop the term $IDR$. It turns out that for the alternate proxy, the tangent secant rule and resulting FDR-Lfdr balance equation are exact identites.

Let $\epsilon>0$ and $\tau>0$. The threshold $t_{alt}$ maximizing $\overline{Sep}(t)$ as a function of $t$ satisfies the Tangent-Secant rule as an exact identity; at this threshold we have

Proof. Now, $A$ and $B$ are both smooth functions of $t$, so at the $t$ optimizing $AB^{-1/2}$ we have

By inspection $B^{\prime}(t)<0$ for every $t>0$. Hence,

The Tangent-Secant Rule follows. We now remark that

The full proof of Theorem 1, which we omit, simply shows that the discrepancy caused by $\overline{Sep}\neq\widetilde{Sep}$ has an asymptotically negligible effect; the two objectives have very similar maximizers.

2 Analysis in the asymptotic RW model

We now invoke the $ARW(r,\beta)$ model: $\epsilon=p^{-\beta}$, $\tau=\sqrt{2r\log(p)}$, $(n,p)\rightarrow\infty$, $n\sim c\log(p)^{\gamma}$, $p\rightarrow\infty$. Let $t_{p}(q)=\sqrt{2q\log(p)}$. The classical Mills’ ratio can be written in terms of the Normal survival function as:

We also need a notation for poly-log terms.

Any occurrence of the symbol $\mbox{\sc PL}(p)$ denotes a term which is $O(\log(p)^{\zeta})$ and $\Omega(\log(p)^{-\zeta})$ as $p\rightarrow\infty$ for some $\zeta>0$. Different occurrences of this symbol may stand for different such terms.

In particular, we may well have $T_{1}(p)=\mbox{\sc PL}(p)$, $T_{2}(p)=\mbox{\sc PL}(p)$, as $p\rightarrow\infty$, and yet $\frac{T_{1}(p)}{T_{2}(p)}\not\rightarrow 1$ as $p\rightarrow\infty$. However, certainly $\frac{T_{1}(p)}{T_{2}(p)}=\mbox{\sc PL}(p)$, $p\rightarrow\infty$.

The following Lemma exposes the main phenomena driving Theorems 1-4. It follows by simple algebra, and several uses of Mills’ Ratio (4.3) in the convenient form $\bar{\Psi}_{0}(t_{p}(q))=PL(p)\cdot p^{-\beta}$.

In the asymptotic RW model $ARW(r,\beta)$, we have:

Quasi power-law for useful feature discoveries:

where the useful feature discovery exponent $\delta$ obeys

Quasi power-law for useless feature discoveries:

As an immediate corollary, under $ARW(r,\beta)$, we have:

On the right side of this display, the poly-log term is relatively unimportant. The driving effect is the power-law behavior of the fraction. The following Lemma contains the core idea behind the appearance of $\rho^{*}$ in Theorems 1 and 2, and the distinction between Regions I and Region II,III.

Let $(\beta,r)\in(\frac{1}{2},1)^{2}$. Let $\gamma(q;r,\beta)$ denote the rate at which

tends to $\infty$ as $p\rightarrow\infty$, for fixed $q$, $r$, and $\beta$. Then $\gamma>0$ if and only if $r>\rho^{*}(\beta)$. A choice of $q$ maximizing this ratio is given by (3.3).

Proof Sketch. By inspection of $\delta$, it is enough to consider $q\leq 1$. The ratio grows to infinity with $p$ like $p^{\gamma}$, where

provided there exists $q\in$ obeying $\gamma(q;r,\beta)>0$.

Let $q^{*}(r,\beta)$ denote the value of $q$ maximizing the rate of separation:

in the event of ties, we take the smallest value of $q$. Let’s define $\rho^{*}(\beta)$ without recourse to the earlier formula (3.2) but instead by the functional role claimed for it by this lemma:

We will derive the earlier formula (3.2) from this. Now $\gamma=\min(\gamma_{1},\gamma_{2})$ where

In dealing with this maximin, two special choices of $q$ will recur below.

$q_{1}$: Viewed as a function of $q$, $\gamma_{1}$ is maximized on the interval $[r,1]$ (use calculus!) at $q_{1}(r,\beta)\equiv 4r$, and is monotone on either side of the maximum.

$q_{2}$: On the other hand, $\gamma_{2}$ is monotone decreasing as a function of $q$ on $[r,1]$. Hence the maximizing value of $q$ in (4.9) over the set of $q$-values where $\gamma_{2}$ achieves the minimum will occur at the minimal value of $q$ achieving the minimum, i.e. at the solution to:

(4.10) is satisfied uniquely on $[r,1]$ by $q_{2}(r,\beta)=(\beta+r)^{2}/4r$.

The behavior of $\min(\gamma_{1},\gamma_{2})$ varies by cases; see Table 1. To see how the table was derived, note that

Consider the first and second rows. For $q>r$, $\delta=1-q-\beta-r+2\sqrt{rq}$. Hence $\delta<1-q$ on $[r,1]$ iff $-\beta-r+2\sqrt{rq}<0$ iff $q<q_{2}$. Consider the third and fourth rows. For $q<r$, $\delta=1-\beta$. Hence $\delta<1-q$ on $[0,r]$ iff $\beta>q$.

Derivatives $\dot{\gamma}_{i}=\frac{\partial}{\partial q}\gamma_{i}$, $i=1,2$, are laid out in Table 2.

Table 3 presents results of formally combining the two previous tables. There are four different cases, depending on the ordering of $q_{1}$,$q_{2}$,$\beta$ and $r$. In only one case does the above information leave $q^{*}$ undefined. (We note that this is a purely formal calculation; Lemma 4.4, Display (4.13) below shows that rows 2 and 4 never occur.)

To see how Table 3 is derived, consider the first row. Using the derivative table above, we see that $\min(\gamma_{1},\gamma_{2})$ is increasing on $[0,\beta]$, constant on $[\beta,r]$ and decreasing on $[r,1]$. Hence the maximin value is achieved at any $q\in[\beta,r]$. For row 2, $\min(\gamma_{1},\gamma_{2})$ is increasing on $[0,\beta]$, constant on on $[\beta,r]$ increasing on $[r,\min(q_{1},q_{2})]$ and decreasing on $[\min(q_{1},q_{2}),1]$. For row 3, $\min(\gamma_{1},\gamma_{2})$ is constant on $[0,r]$, increasing on $[r,\min(q_{1},q_{2})]$ and monotone decreasing on $[\min(q_{1},q_{2}),1]$. For row 4, $\min(\gamma_{1},\gamma_{2})$ is increasing on $[0,r]$ and decreasing on $[r,1]$.

We are trying to find circumstances where $\gamma\leq 0$. In the above table, we remarked that the hypotheses of rows 2 and 4 can never occur. We can see that in row 1, $\gamma_{1}(\beta;r,\beta)>0$ for $\beta\in$, $r>\beta$. This leaves only row 3 where we might have $\gamma\leq 0$; in that case either $q^{*}=q_{1}$ or $q^{*}=q_{2}$. Writing out explicitly

Hence $\gamma_{1}(q_{1};r,\beta)=0$ along $r=\beta-1/2$, and $\gamma_{1}(q_{1};r,\beta)<0$ for $r<\beta-1/2$.

Consider $1/2<\beta<3/4$. In this range, Lemma 4.4 shows $r<q_{1}(\beta-1/2,\beta)<q_{2}(\beta-1/2,\beta)<1$. Hence $q^{*}(r,\beta)=q_{1}(r,\beta)$ for $r\leq\beta-1/2$, $\beta\in(1/2,3/4)$. We conclude that

We have $\gamma_{1}(q_{2};r,\beta)=0$ along $r=(1-\sqrt{1-\beta})^{2}$, and $\gamma_{1}(q_{2};r,\beta)<0$ for $0<r<(1-\sqrt{1-\beta})^{2}$.

Consider $3/4<\beta<1$. In this range, Lemma 4.4 shows shows $r<q_{2}((1-\sqrt{1-\beta})^{2},\beta)<q_{1}((1-\sqrt{1-\beta})^{2},\beta)<1$. We conclude that

Together (4.8), (4.11), and (4.12) establish that the formula (3.2) has the properties implied by the Lemma.

To complete the Lemma, we need to validate formula (3.3). This follows from rows 1 and 3 of Table 3 and Lemma 4.4. $\Box$

Let $q_{1}(r,\beta)=4r$ and $q_{2}(r,\beta)=(\beta+r)^{2}/4r$ just as in the previous lemma. For $0<\beta<1$ we have

Lemmas 4.3 and 4.4 show that (3.3) gives us one choice of threshold maximizing the rate at which $\widetilde{Sep}$ tends to infinity; is it the only choice? Except in the case, $r>\beta$ and $q_{2}<r$, this is indeed the only choice. As Table 3 shows, in case $r>\beta$ and $q_{2}<r$, any $q\in[\beta,r]$ optimizes the of separation. It turns out that in that case, our formula $q^{*}$ not only maximizes the rate of separation, it correctly describes the leading-order asymptotic behavior of $T_{Ideal}$. The key point is the Tangent-Secant Formula, which picks out from among all $q\in[\beta,r]$ uniquely $q_{2}$. This is shown by the next two lemmas, which thereby complete the proof of Theorem 3.

Set $\gamma_{0}(q;r,\beta)=-\beta-r+2\sqrt{rq}$, for $q\in(0,1)$. In the $ARW(r,\beta)$ model, consider the threshold $t_{q}(p)$. Suppose that $q>r$. Then

Proof. Simple manipulations with Mills’ ratio, this time not simply grouping polylog terms together with the symbol $PL$, give that for the threshold under the ARW model, if $q>0$,

Display (4.16) follows. If $q\neq r$ we have the exact identities:

$q_{2}(r,\beta)$ is the unique solution of $\gamma_{0}(q;r,\beta)=0$. Suppose $r>\beta$ and $\beta<q_{2}<r$.

Proof. By inspection for $q<q_{2}$, $\gamma_{0}<0$ while for $q>q_{2}$, $\gamma_{0}>0$. So from (4.17), $Lfdr(t_{q}(p))$ tends to or $1$ depending on $q<q_{2}$ or $q>q_{2}$. Fix $\eta>0$. If $T_{Ideal}<t_{q_{2}-\eta}$ infinitely often as $p\rightarrow\infty$, then would be a cluster point of $Lfdr(T_{Ideal})$. Similarly, if $T_{Ideal}>t_{q_{2}+\eta}$ infinitely often as $p\rightarrow\infty$, then $1$ would be a cluster point of $Lfdr(T_{Ideal})$. From $q>\beta$ and (4.16) we know that $FDR(T_{Ideal})\rightarrow 0$. ¿From the Tangent-Secant rule we know that $Lfdr(T_{Ideal})\rightarrow 1/2$. (4.18) follows. $\Box$

3 FDR/Lfdr properties of ideal threshold

We now turn to Theorem 4. Important observation:

The asymptotic $T_{Ideal}=t_{q^{*}}(p)\cdot(1+o(1))$ is simply too crude to determine the FDR and Lfdr properties of $T_{Ideal}$; it is necessary to consider the second-order effects implicit in the $(1+o(1))$ term. For this, the Tangent-Secant formula is essential.

Indeed (4.17) shows that the only possibilities for limiting local FDR of a threshold of the exact form $t_{q}(p)$ are $0,1/2,1$. The actual local FDR of $T_{Ideal}$ spans a continuum from $$, due to the fact that small perturbations of$t_{q}(p)(1+o(1))$implicit in the$o(1)$can cause a change in the local FDR. To understand this, for$q\neq r$and$s\in(0,\infty)$, put

Choosing $s$ appropriately, we can therefore obtain a perturbed $q$ which perturbs the Lfdr and FDR. In fact there is a unique choice of $s$ needed to ensure the Tangent-Secant formula.

For given values $F_{1}$, $F^{\prime}_{1}$, $T_{1}\neq 0$, $T^{\prime}_{1}\neq 0$, put

This choice of $s$ obeys the Tangent-Secant rule:

To use this recall (4.15)-(4.16)-(4.17). These formulas give expressions for $T_{1}/F_{1}$ and $T^{\prime}_{1}/F^{\prime}_{1}$. Plugging in $q=q^{*}$ we get

where $s^{*}$ is obtained by setting $T_{1}=(\sqrt{q^{*}}-\sqrt{r})^{-1}$, $F_{1}=2/\sqrt{q}$, $T^{\prime}_{1}=1$ and $F^{\prime}_{1}=2$. Moreover, if $\beta/3<r<\beta$,

Connection of HC objective with S​e​p~~𝑆𝑒𝑝\widetilde{Sep}

Let $F=F_{\epsilon,\tau}$ be the two-point mixture of Definition 2.3 and $G=G_{\epsilon,\tau}$ the corresponding folded distribution. Then in the asymptotic RW model we have, for $t=t_{q}(p)$:

Step (5) follows from $G_{\epsilon,\tau}(t_{q}(p))\rightarrow 1$ as $p\rightarrow\infty$. Step (5.2) follows from $\epsilon FPR(t_{q}(p))=o(TPR(t_{q}(p)))$ as $p\rightarrow\infty$. Step (5.3) follows from $TPR(t_{q}(p))=o(TPR(t_{q}(p)))$ as $p\rightarrow\infty$.

This derivation can be made rigorously correct when $q=q^{*}(r,\beta)$, where $q^{*}$ is as announced in Theorem 2. With extra work, not shown here, one obtains:

The statements made for the ideal threshold in Theorems 1-3 are equally valid for the ideal HC threshold.

Suboptimality of phase diagram for other methods

Setting thresholding by control of the False Discovery rate is a popular approach. Another approach, more conservative classical and probably even more popular, is the Bonferroni approach, which controls the expected total number of false features selected. Theorem 2 and 3 implicitly show the suboptimality of these approaches. The implications include an attractive relationship to the regions of Definition 3.3.

Under the Asymptotic RW model $ARW(r,\beta)$, the regions of the $(r,\beta)$ phase space for successful classification are:

Remark. The successful region of Ideal FDRT and Bonferroni is smaller than that of HCT. However, even for regions when FDRT or Bonferroni would be successful, HCT stil yields more accurate classifications in terms of the convergence rate. This is especially important for finite-p performance. In short, Bonferroni and FDRT fail to adapt the difficulty level of the classification problem measured by $(\epsilon,\tau)$, one picks a fixed threshold, another picks a fixed false discovery rate.

Proof Sketch. We continue to use the notations $q^{*}(r,\beta)$, $q_{i}(r,\beta),i=1,2$, $\gamma(r,\beta)$ $\gamma_{i}(r,\beta)$, $i=1,2$ from the proof of Lemma 4.2. The Bonferroni threshold $-\Phi^{-1}(1/p)=t_{1}(p)(1+o(1))$ as $p\rightarrow\infty$. In this proof sketch, we analyze $t_{1}(p)$ as if it were the Bonferroni threshold. Applying (4.7), and the definition of $\gamma$ we have

while $\gamma_{1}(1;r,\beta)=\delta(1;r,\beta)/2$ and $\gamma_{2}(1;r,\beta)=\delta(1;r,\beta)$. Hence $\gamma(1;r,\beta)>0$ iff $\delta(1;r,\beta)>0$. But

which is positive iff $r>1-\sqrt{1-\beta}$.

Precise FDR control at level $\alpha$ requires that $TP/FP=(1-\alpha)/\alpha$. Suppose that $r<\beta$. Equation (4.15) relates the TP/FP ratio to $\gamma_{0}$, $q$ and $r$. Clearly, the only way to keep $TP/FP$ bounded away from 0 and infinity is to choose $q$ so that $\gamma_{0}(q;r,\beta)=0$. We note that $q_{2}(r,\beta)$ exactly solves this problem:

It follows that in the $ARW(r,\beta)$ model, the FDRT functional obeys $T_{FDR,\alpha}(\epsilon,\tau)=t_{q_{2}}(p)(1+o(1))$, $p\rightarrow\infty$. In this proof sketch we analyze $t_{q_{2}}(p)$ as if it were exactly the FDR threshold.

while $\gamma_{1}(q_{2};r,\beta)=1-\frac{(\beta+r)^{2}}{4r}$ and $\gamma_{2}(q_{2};r,\beta)=\frac{1}{2}-\frac{(\beta+r)^{2}}{8r}$. Hence $\gamma(q_{2};r,\beta)>0$ is positive iff $r>1-\sqrt{1-\beta}$.

The last paragraph assumed $r<\beta$. On inuitive grounds, the region $r>\beta$ offers even better perfomance, so it must lie entirely above the phase transition. We omit details. $\Box$

Table 4 compare the exponents in SEP for different methods. Also see Figure 3 for a comparison of the exponents for $\beta=1/2$ and $\beta=5/8$.

Discussion and conclusions

We consider here only asymptotic, ideal behavior. Conceptually, the ideal threshold envisions a situation with an oracle, who, knowing $\epsilon$ and $\tau$ and $n$ and $p$, chooses the very best threshold possible under those given parameters. In this paper we have analyzed the behavior of this threshold within a certain asymptotic framework. However, clearly, no empirical procedure can duplicate the performance of the ideal threshold.

We have seen that ideal HC thresholding comes close. This ideal HC threshold does not involve optimal exploitation of knowledge of $\epsilon$ and $\tau$, but merely the availability of the underlying distribution of feature $Z$-scores $F_{\epsilon,\tau}$. In this paper,we analyzed the behavior of $t^{HC}=T_{HC}(F_{\epsilon,\tau})$.

This is an ideal procedure because we never would know $F_{\epsilon,\tau}$; instead we would have the empirical CDF $F_{n,p}$ defined by

The (non-ideal) HCT that we defined in Section 1.5.2 is then simply

Because $F_{\epsilon,\tau}(z)=E(F_{n,p})(z)$, we are conditioned to expect that $t^{HC}\approx\hat{t}^{HC}_{n,p}$; indeed there are generations of experience for other functionals $T$ showing that we typically have $T(F_{n,p})\approx T(E(F_{n,p}))$ for large $n$,$p$ for those functionals. Proving this for $T=T_{HC}$ is more challenging than one might anticipate; the problem is that $T_{HC}$ is not continuous at $F=\Phi$, and yet $F_{\epsilon(p),\tau(p)}\rightarrow\Phi$ as $p\rightarrow\infty$. After considerable effort, we justify the approximation of HCT by ideal HCT in . Hence the analysis presented here only partially proves that HCT gives near-optimal threshold feature selection; it explains the connection at the level of ideal quantities but not at the level of fluctuations in random samples.

2 Other asymptotic settings

In the analysis here we consider only the case that $n\sim c\log(p)^{\gamma}$. The phase diagram will be slightly different in case $n$ is bounded and does not go to infinity with $p$, and will be again slightly different in case $n\sim cp^{\gamma}$. The full details are presented in .

3 Phase diagram for finite sample sizes

While the main focus of our paper has been asymptotic analysis, we mention that the phase diagram also reflects the finite-sample behavior. In Figure 4, we consider $p=3\times 10^{3}\times N$, $N=(1,10,100)$. For such $p$, we take $n=\log(p)/2$ and display the boundary of the set of $(\beta,r)$ where ideal HCT yields a classification error between $10\%$ and $40\%$. The figure illustrates that as $p$ grows, both the upper bound and the lower bound migrate towards the common limit curve $r=\rho^{*}(\beta)$.

4 Other work on HC

HC was originally proposed for use in a detection problem which has nothing to do with threshold feature selection: testing an intersection null hypothesis $\mu(j)=0\;\forall j$ . The literature has developed since then. Papers extends the optimality of Higher Criticism in detection to correlated setting. Wellner and his collaborators investigated the Higher Criticism in the context of Goodness-of-fit, see for example . Hall and his collaborators have investigated Higher Criticism for robustness, see for example . HC has been applied to data analysis in astronomy () and computational biology .

HC has been used as a principle to synthesize new procedures: Cai, Jin, Low in (see also Meinshausen and Rice ) use Higher Criticism to motivate estimators for the proportion $\epsilon$ of non-null effects.

5 HC in classification

Higher Criticism was previously applied to high-dimensional classification in Hall, Pittelkow, Ghosh (2008) , but there is a key conceptual difference from the present paper. Our approach uses HC in classifier design – it selects features in designing a linear classifier; but the actual classification decisions are made by the linear classifier when presented with specific test feature vectors. The approach in uses HC to directly make decisions from feature vector data. Let’s call these two different strategies HC-based feature selection (used here) and HC-based decision (used in .)

In the ARW model, for a given classifier performance level, the HC-based decision strategy requires much stronger feature contrasts to reach that level than the HC-based feature selection strategy. In this paper, we have shown that HC-feature-selection requires useful features to have contrasts exceeding

HC-decision requires useful features to have contrasts exceeding

Therefore, if the number of training samples $n>1$, HC-feature selection has an advantage. Imagine that we have $n=36$ training samples; we will find that HC feature selection can asymptotically detect features roughly 1/6 the strength of using HC directly for decision.

References