Biomarker discovery aims to identify biological markers (genes or gene products) that lead to improved diagnostic or prognostic tests, better treatment recommendations, or advances in our understanding of the disease or biological condition of interest (Ilyin et al., 2004; Rifai et al., 2006; Ramachandran et al., 2008). Reliable and reproducible biomarker discovery has proven to be difficult (Diamandis, 2010); while one reason is that preliminary small-sample datasets are inherently difficult to analyze, another factor is that most popular algorithms and methods employed in bioinformatics have inherent limitations that make them unsuitable for many discovery applications.

Consider univariate filter methods like t-tests, which are perhaps the most ubiquitous throughout the bioinformatics literature. Since they focus on only one feature at a time, filter methods cannot take advantage of potential correlations between markers. In particular, they cannot identify pairs (or sets) of markers with tell-tale behaviors only when observed in tandem. Multivariate methods, on the other hand, can account for correlations, but the vast majority are wrapper or embedded feature selection algorithms designed to aid in classification or regression model reduction. Model construction is not a reasonable goal in most small-scale exploratory studies, particularly in biology where the expected number of variables is large and the nature of their interactions can be highly complex and context dependent. Furthermore, feature selection methods designed for model reduction intrinsically penalize redundant features and reward smaller feature sets (Sima and Dougherty, 2008; Awada et al., 2012; Ang et al., 2016; Li et al., 2017). This can be counterproductive in biomarker discovery; so much so that filter methods often far outperform multivariate methods (Sima and Dougherty, 2006, 2008; Foroughi pour and Dalton, 2018b).

Optimal Bayesian feature selection (OBFS) is a supervised multivariate screening method designed to address these problems (Dalton, 2013; Foroughi pour and Dalton, 2019). The OBFS modeling framework, discussed in detail in Section 2, assumes that features can be partitioned into two types of independent blocks: blocks of correlated “good” features that have distinct joint distributions between two classes, and blocks of correlated “bad” features that have identical joint distributions in all classes (Foroughi pour and Dalton, 2018a). Given distributional assumptions on each block, and priors on the distribution parameters, OBFS can be implemented using one of many selection criteria, for example the maximum number correct (MNC) criterion, which maximizes the posterior expected number of features correctly identified as belonging to good versus bad blocks, or the constrained MNC (CMNC) criterion, which additionally constrains the number of selected features. In this way, OBFS aims to optimally identify and rank the set of all features with distinct distributions between the classes, which represents the predicted set of biomarkers, while accounting for correlations.

Optimal Bayesian feature filtering (OBF) is a special case of OBFS that assumes that features are independent (blocks of size one) and the parameters governing each feature are independent (Foroughi pour and Dalton, 2015). Gaussian OBF (which assumes independent Gaussian features with conjugate priors, henceforth referred to as simply OBF), has very low computation cost, robust performance when its modeling assumptions are violated, and particularly excels at identifying strong markers with low correlations (Foroughi pour and Dalton, 2017a, 2018a). Furthermore, it has been shown that OBF is strongly consistent under very mild conditions, including cases where the features are dependent and non-Gaussian (Foroughi pour and Dalton, 2019). In particular, OBF, in the long run, selects precisely the set of features with distinct means or variances between the classes. The consistency theorem in (Foroughi pour and Dalton, 2019) formalizes our intuition that OBF cannot take advantage of correlations; although OBF identifies features that are individually discriminating, we see that it has no capacity to identify features that only discriminate when grouped together, or features that are merely correlated with (or linked through a chain of correlation with) other discriminating features.

In bioinformatics, features that only discriminate when grouped together and features that are correlated with other strong markers are of tremendous interest, since they may represent important actors in the biological condition under study. While multivariate methods have the potential to discover such features, as discussed above, methods focused on model reduction tend to do so unreliably and inconsistently. Rather than involve classification or regression models, OBFS searches for intrinsic differences between features. Like most multivariate methods, one difficulty with OBFS is that it is computationally expensive (except in certain special cases, like OBF). If an exact solution is desired, even Gaussian OBFS (which assumes independent blocks of Gaussian features with conjugate priors, henceforth referred to as simply OBFS) is currently only tractable in small-scale problems with up to approximately ten features. This has lead to the development of a number of promising computationally efficient heuristic algorithms based on OBFS theory, for example 2MNC-Robust (Foroughi pour and Dalton, 2014), POFAC, REMAIN, and SPM (Foroughi pour and Dalton, 2018a).

Although OBFS and heuristic algorithms based on OBFS have demonstrated better performance than other multivariate methods in biomarker discovery, it is currently unknown precisely what features they select in the long run. Our main contribution is a theorem presented in Section 3, analogous to the consistency theorem presented in (Foroughi pour and Dalton, 2019) for OBF, that shows that OBFS is strongly consistent under mild conditions. The consistency proof for OBFS utilizes nine lemmas and is significantly more complex than that of OBF. This theorem identifies precisely what features are selected by OBFS asymptotically, provides conditions that guarantee convergence, justifies the use of OBFS in non-design settings where its internal assumptions are invalid, and characterizes rates of convergence for key posteriors in the framework, including marginal posteriors on different types of features.

The asymptotic behavior of optimal feature selection provides a frame of reference to understand the performance of heuristic algorithms based on OBFS; for instance, we may compare the sets of features selected asymptotically, the conditions required to guarantee convergence, and rates of convergence. Furthermore, while numerous works emphasize the need to account for gene interactions and to detect families of interacting biomarkers [e.g., see Han et al. (2017), Xi et al. (2018) and Fang et al. (2019)], typically the focus is on simple statistics that measure pairwise interactions in the data, rather than on establishing intrinsic characteristics of complete marker families, or on quantifying the performance and properties of selection algorithms designed to identify marker families. Thus, this work is also important because it proposes a formal definition for marker families (correlated blocks of features), and shows that these marker families are identifiable (via OBFS).

We begin with an overview of the Gaussian modeling framework presented in Foroughi pour and Dalton (2018a), and a derivation of the corresponding optimal selection rule, OBFS. Let F be the set of feature indices, each corresponding to a real-valued feature. We call feature indices we wish to select “true good features,” denoted by ¯G, and we call feature indices we wish to not select “true bad features,” denoted by ¯. In this model, good features assign (w.p. 1 over the prior) distinct distributions between two classes, labeled y = 0 and y = 1, while bad features assign the same distribution across the whole sample. We further assume that ¯G and ¯B can each be partitioned into sub-blocks of interacting features. We call the set of all sub-blocks a “true feature partition,” denoted by ¯P. If ¯G and ¯B are non-empty, we write ¯P = (), where ¯u and ¯v are positive integers, and the set of all ¯’s and ¯’s are disjoint such that ¯¯, ¯¯, and all ¯’s and ¯’s are non-empty. If ¯G is empty, then we still denote ¯P this way, but also define , and define sums and products from 1 to ¯u to be 0 and 1, respectively. We define similar conventions when ¯B is empty.

In the Bayesian model, ¯G, ¯B and ¯P are random sets, and we denote instantiations of these random sets by G, B, and P = (), respectively. We define a prior on the true feature partition, P( ¯P = P), which induces a prior on the true good features, P( ¯G = G) = ).

Given a true feature partition ¯P = P, define for each class y = 0, 1 and each good block index i = 1, . . . , u, and for each bad block index j = 1, . . . v. Let

denote the collection of all ’s and ’s. Assume the ’s and ’s are mutually independent, i.e., we have a prior of the form

For good block , assume = [] for each class y = 0, 1, where is a length |A| column vector, Σis an matrix, and |A| denotes the cardinality of set A. We assume ) = )) is normal-inverse-Wishart:

) = exp5 tr((Σ)(2.2) ) = exp()(Σ))(2.3)

where |X| is the determinant, tr(X) is the trace, and is the transpose of matrix X. For a proper prior, 1, is a symmetric positive-definite matrix, 0, is an 1 vector,

and Γdenotes the multivariate gamma function, where d is a positive integer. Likewise, for bad block assume = [], and that ) is normal-inverse-Wishart with hyperparameters , and , and scaling constants and .

Let be a sample composed of n labeled points with points in class y, where labels are determined by a process independent from ¯P and (for instance, using random sampling or separate sampling). Given ¯, and the labels, assume all sample points are mutually independent, and features in different blocks are also independent. Assume that features in block = 1, . . . , u, are jointly Gaussian with mean and covariance matrix Σunder class y, and that features in block , j = 1, . . . , v, are jointly Gaussian with mean and covariance Σacross the whole sample.

The posterior on ¯P over the set of all valid feature partitions is given by the normalized product of the prior and likelihood function. It can be shown that P( ¯P = ), where ) is normalized to have unit sum, and

For each good block = 1, . . . , u, and class y = 0, 1,

and ˆand ˆΣare the sample mean and unbiased sample covariance of features in A under class y. Similarly, for each bad block , we find and using (2.8) through (2.12) with all subscript y’s removed, where ˆand ˆΣare the sample mean and unbiased sample covariance of features in A across the whole sample.

The posterior on ¯P induces a posterior on ¯G over all subsets ,

as well as posterior probabilities that each individual feature is in ¯G,

Note ) = P(¯) and ) = P( ¯) are distinct.

The objective of OBFS is to identify the set of true good features, ¯G. We will consider two objective criteria: MNC and CMNC. MNC maximizes the expected number of correctly identified features; that is, MNC outputs the set with complement B = F\G that maximizes E[¯¯B|], which is given by = : (Foroughi pour and Dalton, 2014). CMNC maximizes the expected number of correctly identified features under the constraint of selecting a specified number of features, D, and is found by picking the D features with largest ) (Foroughi pour and Dalton, 2017b). Both MNC and CMNC require computing ) for all , which generally requires computing ) for all valid feature partitions P, and is generally intractable unless |F| is small.

Under proper priors, Gaussian OBFS takes the following modeling parameters as input: (1) ) for each valid feature partition, P, (2) 0, , 1, and symmetric positive-definite for all y and all possible good blocks A in valid P, and (3) 0, 1, and symmetric positive-definite for all possible bad blocks A in valid P. If CMNC is used, it also takes D as input. When ¯P = P) is improper, the above derivations are invalid, but ) can still be taken as a definition and computed as long as: (1) ) is proper, (2) 0, 1, and is symmetric positive-definite, (3) 0, 1, and is symmetric positive-definite, and (4) and are no longer given by (2.4) and (2.5), and similarly and are no longer given by analogous equations; instead these are positive input parameters specified by the user.

Gaussian OBF is a special case where ) assumes that all blocks ¯¯and ¯are of size one, and the events are mutually independent. For each , OBF takes as input (1) marginal priors P(¯G), (2) scalars and for all y, (3) scalars and , and (4) ) if improper priors are used. OBF also takes D as input if CMNC is used. Under OBF, it can be shown that ), defined in (2.14), simplifies to ) = h(f)/(1 + h(f)), where

and are defined in (2.9), (2.10) and (2.11), respectively, and , and are defined similarly. Rather than evaluate ) for all feature partitions, OBF under both MNC and CMNC reduces to simple filtering, where features are scored by the h(f) given in (2.15).

Let be an arbitrary sampling distribution on an infinite sample, . Each sample point in consists of a binary label, y = 0, 1, and a set of features corresponding to a set of feature indices, F. For each n = 1, 2, . . ., let denote a sample consisting of the first n points in , let denote the number of points in class y, and define .

The goal of feature selection is to identify a specific subset of features (say those with different distributions between two classes), which we denote by ¯G. Proving strong consistency for a feature selection algorithm thus amounts to showing that

with probability 1 (w.p. 1) over the infinite sampling distribution, where n is the sample size, is a sample of size n, and ˆ) is the output of the selection algorithm. Here, ˆ) and ¯G are sets, and we define lim) = ¯G to mean that ˆ) = ¯G for all but a finite number of n.

OBFS and OBF fix the sample and model the set of features we wish to select and the sampling distribution as random. To study consistency, we reverse this: now the sampling distribution is fixed and the sample is random. We begin by reviewing a few definitions and a special case of the strong consistency theorem for OBF.

Definition 1 (Foroughi pour and Dalton, 2019). ¯G is an independent unambiguous set of good features if the following hold, where and are the mean and variance of feature f under class y, respectively:

1. For each ¯and exist for all y such that either or . 2. For each ¯and exist for all y such that = and = .

Definition 2 (Foroughi pour and Dalton, 2019). An infinite sample, , is a balanced sample if lim inf0, lim sup1, and, conditioned on the labels, sample points are independent with points belonging to the same class identically distributed.

Theorem 1 (Foroughi pour and Dalton, 2019). Let be a fixed infinite sample and let G be a fixed subset of F. If lim) = 1, then limˆ) = G and limˆ) = G if D = |G|, where ˆ) and ˆ) are the MNC and CMNC feature sets under ), respectively.

By Theorems 1 and 2, Gaussian OBF under MNC and Gaussian OBF under CMNC with “correct” D (the user knows in advance how many features to select) are strongly consistent if the conditions of Theorem 2 hold. Condition (i) uses Definition 1, and characterizes the features that OBF aims to select. Condition (iii) uses Definition 2, and characterizes requirements of the sample. Conditions (i)–(iii) impose very mild assumptions on the data generating process, which is not required to satisfy the Gaussian OBF modeling assumptions in Condition (iv).

Let us now turn to OBFS. Denote the true mean and covariance of features in feature set A and class y by and Σ, respectively (the features need not be Gaussian). OBFS aims to select the smallest set of features, ¯G, with different means or covariances between the classes, where ¯B = F\ ¯G has the same means and covariances between the classes and is uncorrelated with ¯G. This is formalized in the following definition.

Definition 3. ¯G is an unambiguous set of good features if the following hold:

4. Each feature in ¯G is uncorrelated with each feature in ¯B in both classes. 5. Conditions 1-4 do not all hold for any strict subset ¯G.

Assuming all first and second order moments exist, an unambiguous set of good features always exists and is unique. To prove uniqueness, let ¯and ¯be arbitrary unambiguous sets of good features, and let ¯= ¯¯. By condition 4, Σ, has a block diagonal structure for each y = 0, 1 corresponding to the sets of features ¯, ¯, ¯, and F\( ¯). Thus, condition 4 holds for ¯. By condition 3, = and Σ= Σfor all of these blocks except ¯. Thus, condition 3 holds for ¯. If , then (since the means for each class are equal for ¯. Alternatively, if Σ= Σ, then (since ¯and ¯are uncorrelated for each class, and the covariances between each class are equal for ¯¯) Σ= Σ. In either case, condition 2 holds for ¯. Since ¯¯and ¯¯, by condition 5 we must have ¯= ¯= ¯. We denote the unique unambiguous set of good features by ¯G, and its complement by ¯B, throughout.

Similar to the Bayesian model, define a feature partition to be an ordered set of the form P = (), where the set of ’s and ’s partition F. We call each a “good block” and each a “bad block”, keeping in mind that these terms are always relative to a specified arbitrary partition. Feature partitions with no good blocks or no bad blocks are permitted, with appropriate conventions for unions, sums, and products over the (empty set of) corresponding blocks. The following definitions formalize a non-Bayesian analog of the true feature partition.

Definition 4. Let and be arbitrary feature partitions. is a mesh of if every block in is a subset of a block in is a refinement of if every good block in is a subset of a good block in and every bad block in is a subset of a bad block in is a strict refinement of if it is a refinement and .

Definition 5. ¯P = ({ ¯¯¯¯unambiguous feature partition

1. ¯¯is an unambiguous set of good features. 2. Each feature in any arbitrary block is uncorrelated with each feature in any other block in both classes. 3. Conditions 1 and 2 do not hold for any feature partition P that is a strict refine-ment of ¯P.

An unambiguous feature partition always exists and is unique, and we denote it by ¯P throughout. By definition, the unambiguous feature partition ¯P induces the unambiguous set of good features ¯G.

Our main result is given in the following theorem (Theorem 3), which provides sufficient conditions for the (almost sure) convergence of ) to a point mass at ¯P, thereby guaranteeing the (almost sure) convergence of ) to a point mass at ¯G. By Theorem 1, Gaussian OBFS under MNC and Gaussian OBFS under CMNC with “correct” D are strongly consistent if the conditions of Theorem 3 hold, which are very mild. Condition (i) is based on Definition 3 and essentially guarantees that ¯G really represents the type of features OBFS aims to select (those with different means or covariances between the classes). Conditions (i) and (ii) ensure certain moments exist, and Condition (iii) ensures that all covariances are full rank. There are no other distributional requirements. Condition (iv) is based on Definition 2 and characterizes the sampling strategy; the proportion of points observed in any class must not converge to zero, and sample points should be independent. Condition (v) places constraints on the inputs to the OBFS algorithm; we must assign a non-zero probability to the true feature partition.

Finally, from the proof of Theorem 3 we have that, under the conditions stated in the theorem, w.p. 1 there exist 0 < r < 1 and N > 0 such that

for all n > N and all = ¯P that label any good features as bad or that put features that are correlated in separate blocks. Also, w.p. 1 there exist s, N > 0 such that

for all n > N and all = ¯P. Fix ¯G. By (3.2), for each feature partition P such that , w.p. 1 there exist 0 1 and 0 such that ( ¯for all . Therefore, w.p. 1,

for all n > N, where max1 and N > max. Thus, w.p. 1,

for all n > N. The marginal posterior of good features converges to 1 at least exponentially w.p. 1. Now fix ¯B. By (3.3), w.p. 1 there exist s, N > 0 such that

for all n > N. In other words, the marginal posterior of bad features converges to 0 at least polynomially w.p. 1. This characterizes rates of convergence of the posterior on feature partitions, and the marginal posterior probabilities on individual features.

Throughout, “f and g are asymptotically equivalent” and “as ” mean that lim) = 1, v(i) denotes the ith element of vector v, M(i, j) denotes the ith row, jth column element of matrix M, 0denotes the all-zero matrix, and the sample mean and unbiased sample covariance of features in block A and class y are denoted by ˆand ˆΣ, respectively.

(i) ¯P is the unambiguous feature partition. (ii) For all , the fourth order moment across both classes exists and is finite. (iii) Σand Σare full rank. (iv) is a balanced sample w.p. 1. (v) ) assumes a fixed Gaussian OBFS model for all n with ( ¯P) > 0.

Fix = ¯P. If ) = 0 then (3.8) holds trivially. In the remainder of this proof, assume = 0. It suffices to show:

We prove this by constructing a sequence of partitions from P to ¯P in six steps. In Step 1, blocks that are labeled bad in P are split into subsets of bad blocks in ¯P and subsets of the unambiguous set of good features ¯G. In Steps 2 and 3, blocks that are labeled bad in the previous partition but that are actually subsets of ¯G are labeled good (and in Step 3 they are also merged with other good blocks). In Step 4, blocks that are labeled good in the previous partition are split into subsets of (good and bad) blocks in ¯P. In Step 5, blocks that are labeled good in the previous partition but that are actually subsets of bad blocks in ¯P are labeled bad. In Step 6, we merge blocks in the previous partition until we arrive at ¯P. An example of this sequence of partitions is illustrated in Fig. 1. Finally, Step 7 uses the sequence of partitions constructed in Steps 1-6 to show (3.9).

Step 1. Let 2, be a sequence of partitions where (1) we let = P, (2) for each k = 11, we let be identical to except with one bad block partitioned into two bad blocks such that one of the smaller blocks is the intersection of the original block with a bad block in ¯P, and (3) we iterate until no more blocks can be split. The order we split blocks does not matter; at the end we always obtain the partition that is identical to P except with each bad block split into smaller blocks that are each either a subset of a bad block in ¯P or a subset of the unambiguous set of good features ¯G. Suppose . Note is a strict refinement of P. By Lemma 1, w.p. 1 there exist 0 (which may depend on the sample) such that for all . Furthermore,

Figure 1: An example of the sequence of partitions constructed in the consistency proof for the OBFS algorithm, starting from the partition P (top) and ending at ¯P (bottom). Green represents good blocks in the current partition (thus ¯G is the set of green blocks at the bottom), and red represents bad blocks in the current partition (thus ¯B is the set of red blocks at the bottom).

For each k = 11, by Lemma 3, w.p. 1 there exist 0 (which may depend on the sample) such that )(log for all . By (3.10), w.p. 1 there exist 0 (namely from above, + , and N = max()) such that ((log for all n > N. Finally, since for 0 dominates (log , ((0 w.p. 1.

Step 2. Construct a sequence of partitions starting from where, in each step, one bad block in the current partition that is (a) contained in ¯G, and (b) has either different means or different covariances between the classes, is labeled good in the next partition, and iterate until no more blocks can be re-labeled. The order we re-label blocks does not matter; we always obtain the same final partition, which we call . Suppose . By Lemma 1, w.p. 1 there exists such that for all n large enough. By Lemma 4, w.p. 1 there exists 0 1 such that for all n large enough. Since for 0 1 dominates 0 w.p. 1.

Step 3. Some bad blocks in may have equal means and covariances between the classes, but may still be contained in ¯G because they are correlated with, or connected through a chain of correlation with, features in ¯G that have either different means or different covariances between the classes. Here we construct a sequence of partitions to correct the label of these features. The construction goes as follows: starting from , in each step merge one bad block in the current partition that is contained in ¯G (which must have the same means and covariances between classes after Step 2) with one good block in the current partition that is correlated in at least one class (because it is correlated it must also be in ¯G), label the resulting block as good in the next partition, and iterate until no more blocks can be merged.

First only blocks directly correlated with good blocks in the original partition can be moved, then only blocks directly correlated with good blocks in the original partition or the block just moved can be moved, etc. All bad blocks in that are actually in ¯G will eventually be moved because: (1) all features that have either different means or different variances have already been correctly labeled good in Step 2, and (2) the definition of ¯P guarantees that every feature in ¯G is connected to at least one feature in ¯G with either different means or different variances via at least one chain of pairwisecorrelated features that are all in the same good block in ¯P. Thus, each bad block in the final partition, call it , is guaranteed to be a subset of a bad block in ¯P.

Suppose . By Lemma 1, w.p. 1 there exists such that for all n large enough. By Lemma 5, w.p. 1 there exists 0 1 such that for all n large enough. Since for 0 1 dominates , 0 w.p. 1.

Step 4. Construct a sequence of partitions starting from where, in each step, one good block in the current partition is split into two good blocks such that one of the smaller blocks is the intersection of the original block with a (good or bad) block in ¯P, and iterate until no more blocks can be split. The order we split blocks does not matter; in the final partition, call it , each bad block is a subset of a bad block in ¯P and each good block is a subset of a good or bad block in ¯P. Suppose . Note that is a strict refinement of . By Lemma 1, w.p. 1 there exists 0 such that for all n large enough. By Lemma 6, w.p. 1 there exists 0 such that (log for all n large enough. Finally, since for 0 dominates (log 0 w.p. 1.

Step 5. Construct a sequence of partitions starting from where, in each step, one good block in the current partition that is contained in a bad block in ¯P is labeled bad in the next partition, and iterate until no more blocks can be re-labeled. The order we re-label blocks does not matter; in the final partition, call it , each block is a subset of a block in ¯P of the same type. Suppose . Note that is a mesh of , every good block in is a good block in , and every bad block in is a block (good or bad) in . By Lemma 1, w.p. 1 there exists 0 such that for all n large enough. By Lemma 7, w.p. 1 there exists 0 such that (log for all n large enough. Since for 0 dominates (log 0 w.p. 1.

Step 6. Construct a sequence of partitions from to ¯P where each partition is a strict refinement of ¯P and, in each step, two blocks of the same type (good or bad) in the current partition are merged into a larger block of the same type in the next partition that is contained in a block of the same type in ¯P. The order we merge blocks does not matter, as long as the pair of blocks merged in each step are correlated in at least one class. Suppose = ¯P. By Lemma 1, w.p. 1 there exists such that ( ¯for all n large enough. By Lemmas 8 and 9, w.p. 1 there exists 0 1 such that ( ¯for all n large enough. Since for 0 1 dominates , ( ¯( ¯0 w.p. 1.

Step 7. We have that

q( ¯P)a( ¯P) = ) ) ) ) ) ( ¯P)a( ¯P) .

(3.11) Since = ¯P, at least one of Steps 1-6 applies. For the steps that apply, the corresponding ratio in (3.11) converges to 0 w.p. 1. For the steps that do not apply, the corresponding ratio equals 1. Thus, (3.9) holds.

Lemma 1. Let P = () and = () be arbitrary feature partitions, and let be a fixed sample such that lim inf0 and lim sup1. Then there exists and N > 0 such that for all n > N. Further, this holds for s > 0 if: (i) is a strict refinement of P, or (ii) is a mesh of P where , every good block in is a subset of a good block in P, and every bad block in is a subset of a good or bad block in P.

Proof. For any feature partition P, we may write

where is a positive constant that does not depend on n. We first focus on each block individually. and are asymptotically equivalent to , and and are asymptotically equivalent to n. Further, for all positive integers d,

and by Stirling’s formula,

as . Thus, for all ,

as , where 0 does not depend on n. Applying (3.18) and (3.19) in (3.12), we have that ((2as , where 0 does not depend on + 0+ 0, and we treat and as functions of n. Applying this in the ratio, we have:

as , where , and + . We always have lim sup25, and by our assumed bounds on the limit inferior and limit superior of we also have lim inf0. Thus, for all ,

Further,

= 14

If is a strict refinement of P, then 0. Also, if is a mesh of P where , then every good block in is a subset of a good block in P and every bad block in is a subset of a good or bad block in P, thus 0. If 0, then (3.21) holds for all 0 .

Lemma 2. Let be any non-empty feature set such that Σand Σare full rank. Suppose for all the fourth order moment across both classes exists and is finite. Let be a balanced sample. Then, w.p. 1 there exist K, N > 0 (which may depend on the sample) such that

A. Foroughi pour and L. A. Dalton 15 ˆΣ(3.24)

Proof. Fix A. Since is a balanced sample, as , and by the strong law of large numbers ˆconverges to and ˆΣconverges to Σw.p. 1 for both y = 0, 1. In the rest of this proof, we only consider the event where ˆand ˆΣconverge.

Since is a balanced sample, there exist 0 1 and 0 such that for all and y = 0, 1. By the law of the iterated logarithm (Hartman and Wintner, 1941), w.p. 1 there exist 0 and such that

for all and i = 1, . . . , |A| (Foroughi pour and Dalton, 2019). This proves (3.23).

We can decompose the sample covariance for class y as follows:

for all and i, j = 1, . . . , |A|. ((i))((j)) is a random variable with mean Σ(i, j) and finite variance ). Suppose 0. By the law of the iterated logarithm, w.p. 1

for all i, j = 1, . . . , |A|. Thus, w.p. 1 there exists 0 and such that 1

log log n n (3.36)

for all and i, j = 1, . . . , |A|. A similar inequality holds when ) = 0. Combining (3.33), (3.34) and (3.36), w.p. 1 there exists 0 such that

for all and i, j = 1, . . . , |A|. Thus (3.25) holds. Since |ˆΣis a polynomial function of the ˆΣ(i, j), where each term has degree |A| and coefficient 1, w.p. 1 there exists 0 such that

for all . Dividing both sides by proves (3.26).

Note that,

Further (3.23) implies that w.p. 1 ˆ((i) is bounded for all . Similarly, (3.25) implies that w.p. 1 ˆΣ(i, j) is bounded for all . Thus, from (3.39), w.p. 1 there exists 0 and such that

for all and i, j = 1, . . . |A|. Again noting that the determinant is a polynomial, w.p. 1 there exists 0 such that

for all . By (3.38), |ˆΣis bounded for all w.p. 1. Thus, (3.24) holds. Applying the triangle inequality on ((i, j) and applying (3.25) and (3.40), we have that (3.27) also holds.

The sample covariance across all samples in both classes can be decomposed as

Define

Again since w.p. 1 ˆ(i) and ˆΣ(i, j) are bounded for all , by the triangle inequality w.p. 1 there exists 0 such that

for all and i, j = 1, . . . , |A|. Further, by the triangle inequality

for all and i, j = 1, . . . , |A| w.p. 1, where

= ˆ, and . ) is constant for all i = 1, . . . , |A|. By (3.23), w.p. 1 there exists such that log log n/n for all = 0, 1 and i = 1, . . . , |A|. Thus, w.p. 1 there exists 0 such that log log n/n for all and i, j = 1, . . . , |A|. Applying this fact and (3.25) to (3.44), w.p. 1 there exists 0 and such that

for all and i, j = 1, . . . , |A|. Combining this fact with (3.45), w.p. 1 there exists 0 such that

for all and i, j = 1, . . . , |A|. Again since the determinant is a polynomial, w.p. 1 there exists such that

for all . Observe from (3.31) that Σ(i, j) is lower bounded by

and upper bounded by

Since is a polynomial function of the Σ(i, j), where each term has degree |A| and coefficient must also be upper and lower bounded when n is large enough. (3.29) follows by dividing both sides of (3.49) by , and applying a bound on on the right hand side.

Observe that

Since ˆ= + (1 , (3.23) implies that w.p. 1 ˆ) is bounded for all . Similarly, (3.48) and the fact that Σ(i, j) is bounded for n large enough

implies that w.p. 1 there exists such that ˆΣ) is bounded for all . Thus, from (3.52), w.p. 1 there exists 0 and such that

for all and i, j = 1, . . . |A|. Again since the determinant is a polynomial, w.p. 1 there exists 0 such that

for all . By (3.49) and the fact that is bounded for n large enough, w.p. 1 there exists such that |ˆΣis bounded for all . Thus, (3.28) holds. Applying the triangle inequality on (i, j) and applying (3.48) and (3.53), we have that (3.30) also holds.

Corollary 1. Let be any non-empty feature set such that Σand Σare full rank. Suppose for all the fourth order moment across both classes exists and is finite. Let be a balanced sample. Then, w.p. 1 the following all hold for all y = 0, 1 and i, j = 1, . . . , |A|, where Σis given in (3.31).

Proof. (1) follows from (3.27). (2) follows by combining (3.24) and (3.26). (5) follows from (3.30). (6) follows by combining (3.28) and (3.29). (3) and (4) are special cases of (5) and (6), where in this case Σ= Σis constant. Recall that in Lemma 2 we showed that Σ(i, j) and are lower and upper bounded when n is large enough. By (3.30), ) must also be upper and lower bounded when n is large enough, so (7) holds. Finally, since by item 6 and are asymptotically equivalent, must also be upper and lower bounded when n is large enough, so (8) holds.

Lemma 3. Let and be any disjoint feature sets such that , ΣΣ= Σ, features in and are uncorrelated in both classes, and Σand Σare full rank, where . Suppose for all the fourth order moment across both classes exists and is finite. Let be a balanced sample. Then, w.p. 1 there exist c, N > 0 such that

for all n > N. The left hand side of (3.55) is ) when P and are identical feature partitions except A is a bad block in P, and and are bad blocks in .

where

By Corollary 1, the first term in (3.56) is upper bounded by a positive finite constant for n large enough w.p. 1. Without loss of generality, assume features in A are ordered such that

for some matrix V . Observe that

where (. Since , Σ= Σ, and features in and are uncorrelated in both classes,

where Σand Σare defined in (3.31). By Lemma 2, w.p. 1 there exist 0 such that for all = 1and j = 1,

Since W is comprised of only quadratic terms in V and ) w.p. 1 (Corollary 1), w.p. 1 there exist 0 such that for all = 1and j = 1,

Further, since is a polynomial function of the elements of W, where each term has a degree between 1 and and a coefficient that is a polynomial function of the elements of of at most degree 1, and since ) is upper and lower bounded for n large enough w.p. 1 for all i, j = 1, . . . , |A| (Corollary 1), w.p. 1 there exists 0 such that for all ,

Therefore,

for w.p. 1, where the first line follows from (3.57) and (3.59) and the second line from (3.63). Further, by Corollary 1, w.p. 1 there exist 0 and such that for all , thus

for all . The second line holds for large enough that x = () log log n/n is between 0 and 0.5 for all , so that (1 1 + 2x. The last line follows from the fact that (1 + for all x, t > 0. From (3.56), the theorem holds with and , where is such that (log exceeds a bound on the first term in (3.56).

Lemma 4. Let G be any feature set such that either or Σ= Σ, and Σand Σare full rank. Let be a balanced sample. Then, w.p. 1 there exist 0 1 and N > 0 such that

for all n > N. The left hand side of (3.66) is ) when P and are identical feature partitions except G is a bad block in P and a good block in .

where

By Corollary 1, the first term in (3.67) is upper bounded by a positive finite constant for n large enough, and

as , both w.p. 1. Suppose Σ= Σand . Then,

The last line follows from the matrix determinant lemma (Harville, 1998). Since is balanced and Σis positive-definite, there exists 0 < T < 1 such that 0 for all n large enough. The theorem holds for any T < r < 1. If Σ= Σ, then

Fan (1950) showed that for any symmetric positive-definite X and Y and 0 1,

Although not mentioned by Fan, if 0 < p < 1 then equality holds if and only if X = Y (by the weighted arithmetic-geometric-mean inequality). Suppose and fix 0 < 5. Then there exists 0 such that + (1 for all ). Applying this fact here, since is balanced there exists 0 < T < 1 such that 0 . The theorem holds for any T < r < 1.

Lemma 5. Let and be any disjoint feature sets such that , ΣΣ= Σ, there exists at least one feature in and one feature in that are correlated in at least one class, and Σand Σare full rank, where . Let be a balanced sample. Then, w.p. 1 there exist 0 1 and N > 0 such that

for all n > N. The left hand side of (3.73) is ) when P and are identical feature partitions except is a good block in is a bad block in P, and G is a good block in .

Proof. We have that

where

By Corollary 1, the first term in (3.74) converges to a positive finite constant w.p. 1, and

where is a matrix of correlations between features in and . Since

we have that

We always have (Σ), and thus 0 1. Since there exists at least one feature in and one feature in that are correlated in at least one class, in this class is non-zero, (Σ), and 1. Since is balanced, there exist 0 1 such that and for n large enough. Thus,

for n large enough, where . Note , where 0 1. The theorem holds for any R < r < 1.

Lemma 6. Let and be any disjoint feature sets such that features in and are uncorrelated in both classes, and Σand Σare full rank, where . Suppose for all the fourth order moment across both classes exists and is finite. Let be a balanced sample. Then, w.p. 1 there exist c, N > 0 such that

for all n > N. The left hand side of (3.81) is ) when P and are identical feature partitions except A is a good block in P, and and are good blocks in .

Proof. We have that

where

By Corollary 1, the first term in (3.82) converges to a positive finite constant w.p. 1. Fix . Without loss of generality, assume features in A are ordered such that

for some matrix . Observe that

where (). Since features in and are uncorrelated in both classes,

By Lemma 2, w.p. 1 there exist 0 such that for all = 0, 1, i = 1and j = 1we have

Since is comprised of only quadratic terms in and ((i, j) w.p. 1 (Corollary 1), w.p. 1 there exist 0 such that for all = 0, 1, i = 1and j = 1

Further, since is a polynomial function of the elements of , where each term has a degree between 1 and and a coefficient that is a polynomial function of the elements of of at most degree 1, and since ((i, j) w.p. 1 (Corollary 1), w.p. 1 there exists 0 such that for all and y = 0, 1,

Therefore,

for w.p. 1, where the first line follows from (3.83) and (3.85) and the second line from (3.89). Further, w.p. 1 there exists such that for all and y = 0, 1,

The first line holds as long as is large enough that 2 for all and y = 0, 1 (this is possible since both converge). The second line holds as long as is large enough that x = (2) log log n/n is between 0 and 0.5 for all and y = 0, 1, so that (11+2x. The third line holds as long as is large enough that for all and y = 0, 1 and some 0 1 (this is possible since is balanced). Finally, the last line follows from the fact that (1 + for all x, t > 0. From (3.82), the theorem holds with +2and , where is such that (log exceeds a bound on the first term in (3.82).

Lemma 7. Let B be any feature set such that = , Σ= Σ, and Σis full rank. Suppose for all the fourth order moment across both classes exists and is finite. Let be a balanced sample. Then, w.p. 1 there exist c, N > 0 such that

for all n > N. The left hand side of (3.92) is ) when P and are identical feature partitions except B is a good block in P and a bad block in .

Proof. We have that

where

Let be such that is between 0 and 0.5 for all , so that (1 1 + 2x. Then for all ,

w.p. 1, where the last line follows from the fact that (1 + for all x, t > 0. From (3.42),

and for y = 0, 1,

for all = 0, 1 and i, j = 1, . . . , |B|. Since is balanced, there exist 0 < 1 and such that for all and y = 0, 1. Since ˆΣas w.p. 1, there exists B > 0 such that (i, j)| < B for all n, y, i and j w.p. 1. By Lemma 2, w.p. 1 there exists 0 and such that log log n/n for all and y = 0, 1. By the triangle inequality,

for all and i, j = 1, . . . , |B| w.p. 1. In particular, there exists 0 such that

Observe that |ˆΣis a polynomial function of the elements of Σ, (1 and E, where each term has a degree of |B| and a coefficient of 1. In particular,

where is the Levi-Civita symbol, equal to +1 if () is an even permutation of (11 if its an odd permutation, and 0 otherwise, = Σ, = (1 , and

From (3.101) and (3.103), w.p. 1 there exists such that

+ log log n n (3.106)

for all , where

) : exactly one k equals 2, the rest equal 1},

) : exactly one k equals 3, the rest equal 1},

for all , where

= Σ= , and we have used the facts that 1) = m(1, . . . , 1) and ) = ) when (. Similarly, w.p. 1 there exists such that

for all . Combining (3.106), (3.107) and (3.109), w.p. 1 there exists such that

for all . Thus, w.p. 1 there exists such that

for all , where is chosen based on the fact that |ˆΣw.p. 1 (an application of Corollary 1 with hyperparameters = 0, = |A| and = 0 in place of the hyperparameters used by the selection rule), and thus |ˆΣmust be bounded for all n. In addition, w.p. 1 there exists 0 and such that

for , where in the second line we have applied the triangle inequality and the fact that |ˆΣw.p. 1, so is chosen such that for all , and the last line follows from Lemma 2. By Lemma S3 in Foroughi pour and Dalton (2019), there exists (0, 1) such that for all (0, 1) and (1 1 + r),

where in the last inequality we use the fact that 25 for all [0, 1]. Thus, w.p. 1 there exists such that

for all , where the last line follows from the fact that (1 + for all x, t > 0. From (3.93), the theorem holds with + 0and , where is such that (log exceeds times a bound on the first term in (3.93).

Lemma 8. Let and be any disjoint feature sets such that there exists at least one feature in and one feature in that are correlated in at least one class, and Σand Σare full rank, where . Let be a balanced sample. Then, w.p. 1 there exist 0 1 and N > 0 such that

for all n > N. The left hand side of (3.116) is ) when P and are strict refinements of the unambiguous feature partition that are identical except and are good blocks in P and G is a good block in .

where

By Corollary 1, the first term in (3.117) converges to a positive finite constant w.p. 1, and (3.76) holds. The rest of the proof proceeds exactly as in Lemma 5.

Lemma 9. Let and be any disjoint feature sets such that = , Σ= Σ, there exists at least one feature in and one feature in that are correlated in at least one class, and Σis full rank, where . Let be a balanced sample. Then, w.p. 1 there exist 0 1 and N > 0 such that

for all n > N. The left hand side of (3.119) is ) when P and are strict refinements of the unambiguous feature partition that are identical except and are bad blocks in P and B is a bad block in .

where

By Corollary 1, the first term in (3.120) converges to a positive finite constant w.p. 1, and

as , both w.p. 1, where Σ= Σand Σ= Σ. Without loss of generality, assume features in B are ordered such that

where V is a matrix of correlations between features in and . Since

we have that

Since there exists at least one feature in and one feature in that are correlated, V is non-zero, (Σ, and 0 1. The theorem holds for any R < r < 1.

The consistency theory presented herein is important because it provides a richer understanding the type of features selected by OBFS. Furthermore, we have characterized rates of convergence for the posterior on feature partitions, and the marginal posterior probabilities on individual features.

Although here we focus on identifying ¯G using the posterior ), the OBFS framework can be used for optimal Bayesian partition selection, which aims to identify ¯P using the full posterior ). Partition selection may be of interest, for instance, if one wishes to identify communities of closely interacting genes. Since Theorem 3 proves that ) converges to a point mass at ¯P, this theorem has direct implications on the consistency of optimal Bayesian partition selection as well.

The conditions in Theorem 3 are sufficient but not necessary. For example, it may be possible to relax Condition (iii) of Theorem 3. It is also possible for ) to converge to a point mass at ¯G, but for ) to not converge to a point mass at ¯P. For example, if non-zero correlations are present in the data generation process, then the OBF variant of OBFS sets ( ¯P ) = 0, which means that Condition (v) of Theorem 3 does not hold. However, if the independent unambiguous set of good features given in Definition 1 and the unambiguous set of good features given in Definition 3 happen to be equal (the latter always contains the former), then OBF can be strongly consistent relative to the unambiguous set of good features.

OBFS searches for the unambiguous set of good features, which expands upon the independent unambiguous set of good features targeted by OBF. In particular, the unambiguous set of good features includes features that are only strongly discriminating when grouped together, features that are individually weak but strongly correlated with strong features, and features that are linked to discriminating features only through a chain of correlation. The unambiguous feature set is complete in the sense that any features that are not included must not have any first or second order distributional differences between the classes and must be uncorrelated with all features that are included, and minimal in the sense that it is the smallest feature set with this property. The unambiguous feature set is thus of great interest and importance in bioinformatics and many other application areas, although we are not aware of any selection algorithms besides OBFS that aim to identify this set.

This work is supported by the National Science Foundation (CCF-1422631 and CCF-1453563).

Ang, J. C., Mirzal, A., Haron, H., and Hamed, H. N. A. (2016). “Supervised, unsuper- vised, and semi-supervised feature selection: A review on gene selection.” IEEE/ACM Transactions on Computational Biology and Bioinformatics, 13(5): 971–989. 2

Awada, W., Khoshgoftaar, T. M., Dittman, D., Wald, R., and Napolitano, A. (2012). “A review of the stability of feature selection techniques for bioinformatics data.” In Proceedings of the 2012 IEEE 13th International Conference on Information Reuse and Integration (IRI), 356–363. 2

Dalton, L. A. (2013). “Optimal Bayesian feature selection.” 65–68. IEEE. 2

Diamandis, E. P. (2010). “Cancer biomarkers: Can we turn recent failures into success?” Journal of the National Cancer Institute, 102(19): 1462–1467. 1

Fan, K. (1950). “On a theorem of Weyl concerning eigenvalues of linear transformations II.” Proceedings of the National Academy of Sciences, 36(1): 31–35. 22

Fang, G., Wang, W., Paunic, V., Heydari, H., Costanzo, M., Liu, X., Liu, X., Vander- Sluis, B., Oately, B., Steinbach, M., et al. (2019). “Discovering genetic interactions bridging pathways in genome-wide association studies.” Nature communications, 10(1): 1–18. 3

Foroughi pour, A. and Dalton, L. A. (2014). “Optimal Bayesian feature selection on high dimensional gene expression data.” In Proceedings of the 2014 IEEE Global Conference on Signal and Information Processing (GlobalSIP), 1402–1405. 3, 5

— (2015). “Optimal Bayesian feature filtering.” In Proceedings of the 6th ACM Conference on Bioinformatics, Computational Biology, and Health Informatics (ACMBCB), 651–652. 2

— (2017a). “Robust feature selection for block covariance Bayesian models.” In Proceedings of the 2017 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2696–2700. 2

— (2017b). “Robust feature selection for block covariance Bayesian models.” 5

— (2018a). “Heuristic algorithms for feature selection under Bayesian models with block-diagonal covariance structure.” BMC Bioinformatics, 19(Suppl 3): 70. 2, 3

— (2018b). “Optimal Bayesian filtering for biomarker discovery: Performance and ro- bustness.” IEEE/ACM Transactions on Computational Biology and Bioinformatics, early access. 2

— (2019). “Theory of Optimal Bayesian Feature Filtering.” Bayesian Analysis, accepted. 2, 3, 7, 15, 29

Han, K., Jeng, E. E., Hess, G. T., Morgens, D. W., Li, A., and Bassik, M. C. (2017). “Synergistic drug combinations for cancer identified in a CRISPR screen for pairwise genetic interactions.” Nature biotechnology, 35(5): 463. 3

Hartman, P. and Wintner, A. (1941). “On the law of the iterated logarithm.” American Journal of Mathematics, 63(1): 169–176. 15

Harville, D. A. (1998). “Matrix algebra from a statistician’s perspective.” Technometrics, 40(2): 164–164. 22

Ilyin, S. E., Belkowski, S. M., and Plata-Salam´an, C. R. (2004). “Biomarker discovery and validation: Technologies and integrative approaches.” Trends in Biotechnology, 22(8): 411–416. 1

Li, J., Cheng, K., Wang, S., Morstatter, F., Trevino, R. P., Tang, J., and Liu, H. (2017). “Feature selection: A data perspective.” ACM Computing Surveys, 50(6): 94. 2

Ramachandran, N., Srivastava, S., and LaBaer, J. (2008). “Applications of protein microarrays for biomarker discovery.” Proteomics - Clinical Applications, 2(10–11): 1444–1459. 1

Rifai, N., Gillette, M. A., and Carr, S. A. (2006). “Protein biomarker discovery and validation: The long and uncertain path to clinical utility.” Nature Biotechnology, 24(8): 971–983. 1

Sima, C. and Dougherty, E. R. (2006). “What should be expected from feature selection in small-sample settings.” Bioinformatics, 22(19): 2430–2436. 2

— (2008). “The peaking phenomenon in the presence of feature-selection.” Pattern Recognition Letters, 29(11): 1667–1674. 2

Xi, J., Li, A., and Wang, M. (2018). “A novel unsupervised learning model for detecting driver genes from pan-cancer data through matrix tri-factorization framework with pairwise similarities constraints.” Neurocomputing, 296: 64–73. 3