On the Consistency of Optimal Bayesian Feature Selection in the Presence of Correlations

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 $\bar{G}$, and we call feature indices we wish to not select “true bad features,” denoted by $\bar{B}=F\backslash\bar{G}$. 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 $\bar{G}$ and $\bar{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 $\bar{P}$. If $\bar{G}$ and $\bar{B}$ are non-empty, we write $\bar{P}=(\{{\bar{G}_{1}},\ldots,{\bar{G}_{\bar{u}}}\},\{{\bar{B}_{1}},\ldots,{\bar{B}_{\bar{v}}}\})$, where $\bar{u}$ and $\bar{v}$ are positive integers, and the set of all $\bar{G}_{i}$’s and $\bar{B}_{j}$’s are disjoint such that $\bar{G}=\cup_{i=1}^{\bar{u}}\bar{G}_{i}$, $\bar{B}=\cup_{j=1}^{\bar{v}}\bar{B}_{j}$, and all $\bar{G}_{i}$’s and $\bar{B}_{j}$’s are non-empty. If $\bar{G}$ is empty, then we still denote $\bar{P}$ this way, but also define $\cup G_{i}=\varnothing$, and define sums and products from $1$ to $\bar{u}$ to be $0$ and $1$, respectively. We define similar conventions when $\bar{B}$ is empty.

In the Bayesian model, $\bar{G}$, $\bar{B}$ and $\bar{P}$ are random sets, and we denote instantiations of these random sets by $G$, $B$, and $P=(\{{G_{1}},\ldots,{G_{u}}\},\{{B_{1}},\ldots,{B_{v}}\})$, respectively. We define a prior on the true feature partition, $\pi(P)\equiv\text{P}(\bar{P}=P)$, which induces a prior on the true good features, $\pi(G)\equiv\text{P}(\bar{G}=G)=\sum_{P:\cup G_{i}=G}\pi(P)$.

Given a true feature partition $\bar{P}=P$, define $\theta^{G_{i}}_{y}$ for each class $y=0,1$ and each good block index $i=1,\ldots,u$, and $\theta^{B_{j}}$ for each bad block index $j=1,\ldots v$. Let $\theta^{P}$ denote the collection of all $\theta^{G_{i}}_{y}$’s and $\theta^{B_{j}}$’s. Assume the $\theta^{G_{i}}_{y}$’s and $\theta^{B_{j}}$’s are mutually independent, i.e., we have a prior of the form

For good block $A=G_{i}$, assume $\theta^{A}_{y}=[\mu^{A}_{y},\Sigma^{A}_{y}]$ for each class $y=0,1$, where $\mu^{A}_{y}$ is a length $|A|$ column vector, $\Sigma^{A}_{y}$ is an $|A|\times|A|$ matrix, and $|A|$ denotes the cardinality of set $A$. We assume $\pi(\theta^{A}_{y})=\pi(\mu^{A}_{y}|\Sigma^{A}_{y})\pi(\Sigma^{A}_{y})$ is normal-inverse-Wishart:

where $|X|$ is the determinant, $\operatorname{tr}(X)$ is the trace, and $X^{T}$ is the transpose of matrix $X$. For a proper prior, $\kappa^{A}_{y}>|A|-1$, $S^{A}_{y}$ is a symmetric positive-definite $|A|\times|A|$ matrix, $\nu^{A}_{y}>0$, $m^{A}_{y}$ is an $|{A}|\times 1$ vector,

and $\Gamma_{d}$ denotes the multivariate gamma function, where $d$ is a positive integer. Likewise, for bad block $A=B_{j}$ assume $\theta^{A}=[\mu^{A},\Sigma^{A}]$, and that $\pi(\theta^{A})$ is normal-inverse-Wishart with hyperparameters $\kappa^{A}$, $S^{A}$, $\nu^{A}$, and $m^{A}$, and scaling constants $K^{A}$ and $L^{A}$.

Let $S_{n}$ be a sample composed of $n$ labeled points with $n_{y}$ points in class $y$, where labels are determined by a process independent from $\bar{P}$ and $\theta^{\bar{P}}$ (for instance, using random sampling or separate sampling). Given $\bar{P}=P$, $\theta^{P}$, and the labels, assume all sample points are mutually independent, and features in different blocks are also independent. Assume that features in block $G_{i}$, $i=1,\ldots,u$, are jointly Gaussian with mean $\mu^{G_{i}}_{y}$ and covariance matrix $\Sigma^{G_{i}}_{y}$ under class $y$, and that features in block $B_{j}$, $j=1,\ldots,v$, are jointly Gaussian with mean $\mu^{B_{j}}$ and covariance $\Sigma^{B_{j}}$ across the whole sample.

The posterior on $\bar{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 $\pi^{*}(P)\equiv\text{P}(\bar{P}=P|S_{n})\propto\pi(P)q(P)a(P)$, where $\pi^{*}(P)$ is normalized to have unit sum, and

For each good block $A=G_{i}$, $i=1,\ldots,u$, and class $y=0,1$,

and $\hat{\mu}^{A}_{y}$ and $\hat{\Sigma}^{A}_{y}$ are the sample mean and unbiased sample covariance of features in $A$ under class $y$. Similarly, for each bad block $A=B_{j}$, we find $\nu^{A*}$, $\kappa^{A*}$, $Q^{A}$, $S^{A*}$ and $C^{A}$ using (2.8) through (2.12) with all subscript $y$’s removed, where $\hat{\mu}^{A}$ and $\hat{\Sigma}^{A}$ are the sample mean and unbiased sample covariance of features in $A$ across the whole sample.

The posterior on $\bar{P}$ induces a posterior on $\bar{G}$ over all subsets $G\subseteq F$,

as well as posterior probabilities that each individual feature $f\in F$ is in $\bar{G}$,

Note $\pi^{*}(f)=\text{P}(f\in\bar{G}|S_{n})$ and $\pi^{*}(\{f\})=\text{P}(\bar{G}=\{f\}|S_{n})$ are distinct.

The objective of OBFS is to identify the set of true good features, $\bar{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 $G\subseteq F$ with complement $B=F\backslash G$ that maximizes $\text{E}[|G\cap\bar{G}|+|B\cap\bar{B}|]$, which is given by $G^{\textrm{MNC}}=\{f\in F:\pi^{*}(f)>0.5\}$ (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 $G^{\textrm{CMNC}}$ is found by picking the $D$ features with largest $\pi^{*}(f)$ (Foroughi pour and Dalton, 2017b). Both MNC and CMNC require computing $\pi^{*}(f)$ for all $f\in F$, which generally requires computing $\pi(P)$ 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) $\pi(P)$ for each valid feature partition, $P$, (2) $\nu^{A}_{y}>0$, $m^{A}_{y}$, $\kappa^{A}_{y}>|A|-1$, and symmetric positive-definite $S^{A}_{y}$ for all $y$ and all possible good blocks $A$ in valid $P$, and (3) $\nu^{A}>0$, $m^{A}$, $\kappa^{A}>|A|-1$, and symmetric positive-definite $S^{A}$ for all possible bad blocks $A$ in valid $P$. If CMNC is used, it also takes $D$ as input. When $\pi(\theta^{P}|\bar{P}=P)$ is improper, the above derivations are invalid, but $\pi^{*}(P)\propto\pi(P)q(P)a(P)$ can still be taken as a definition and computed as long as: (1) $\pi(P)$ is proper, (2) $\nu^{A*}_{y}>0$, $\kappa^{A*}_{y}>|A|-1$, and $S^{A*}_{y}$ is symmetric positive-definite, (3) $\nu^{A*}>0$, $\kappa^{A*}>|A|-1$, and $S^{A*}$ is symmetric positive-definite, and (4) $K^{A}_{y}$ and $L^{A}_{y}$ are no longer given by (2.4) and (2.5), and similarly $K^{A}$ and $L^{A}$ are no longer given by analogous equations; instead these are positive input parameters specified by the user.

Gaussian OBF is a special case where $\pi(P)$ assumes that all blocks $\bar{G}_{1},\ldots,\bar{G}_{\bar{u}}$ and $\bar{B}_{1},\ldots,\bar{B}_{\bar{v}}$ are of size one, and the events $\{f\in\bar{G}\}$ are mutually independent. For each $f\in F$, OBF takes as input (1) marginal priors $\pi(f)\equiv\text{P}(f\in\bar{G})$, (2) scalars $\nu_{y}^{\{f\}}$, $m_{y}^{\{f\}}$, $\kappa_{y}^{\{f\}}$ and $S_{y}^{\{f\}}$ for all $y$, (3) scalars $\nu^{\{f\}}$, $m^{\{f\}}$, $\kappa^{\{f\}}$ and $S^{\{f\}}$, and (4) $L^{f}\equiv K_{0}^{\{f\}}L_{0}^{\{f\}}K_{1}^{\{f\}}L_{1}^{\{f\}}/(K^{\{f\}}L^{\{f\}})$ if improper priors are used. OBF also takes $D$ as input if CMNC is used. Under OBF, it can be shown that $\pi^{*}(f)$, defined in (2.14), simplifies to $\pi^{*}(f)=h(f)/(1+h(f))$, where

$\kappa_{y}^{\{f\}*}$, $Q_{y}^{\{f\}}$ and $S_{y}^{\{f\}*}$ are defined in (2.9), (2.10) and (2.11), respectively, and $\kappa^{\{f\}*}$, $Q^{\{f\}}$ and $S^{\{f\}*}$ are defined similarly. Rather than evaluate $\pi^{*}(P)$ 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 $\mathcal{F}_{\infty}$ be an arbitrary sampling distribution on an infinite sample, $S_{\infty}$. Each sample point in $S_{\infty}$ 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,\ldots$, let $S_{n}$ denote a sample consisting of the first $n$ points in $S_{\infty}$, let $n_{y}$ denote the number of points in class $y$, and define $\rho_{n}=n_{0}/n$.

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 $\bar{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, $S_{n}$ is a sample of size $n$, and $\hat{G}(S_{n})$ is the output of the selection algorithm. Here, $\hat{G}(S_{n})$ and $\bar{G}$ are sets, and we define $\lim_{n\to\infty}\hat{G}(S_{n})=\bar{G}$ to mean that $\hat{G}(S_{n})=\bar{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.

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 $\mu^{A}_{y}$ and $\Sigma^{A}_{y}$, respectively (the features need not be Gaussian). OBFS aims to select the smallest set of features, $\bar{G}$, with different means or covariances between the classes, where $\bar{B}=F\backslash\bar{G}$ has the same means and covariances between the classes and is uncorrelated with $\bar{G}$. This is formalized in the following definition.

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

Similar to the Bayesian model, define a feature partition to be an ordered set of the form $P=(\{G_{1},\ldots,G_{u}\},\{B_{1},\ldots,B_{v}\})$, where the set of $G_{i}$’s and $B_{j}$’s partition $F$. We call each $G_{i}$ a “good block” and each $B_{j}$ 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.

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

Our main result is given in the following theorem (Theorem 3), which provides sufficient conditions for the (almost sure) convergence of $\pi^{*}(P)$ to a point mass at $\bar{P}$, thereby guaranteeing the (almost sure) convergence of $\pi^{*}(G)$ to a point mass at $\bar{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 $\bar{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\neq\bar{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\neq\bar{P}$. Fix $f\in\bar{G}$. By (3.2), for each feature partition $P$ such that $f\notin\cup G_{i}$, w.p. $1$ there exist $0<r_{P}<1$ and $N_{P}>0$ such that $\pi^{*}(P)/\pi^{*}(\bar{P})<r_{P}^{n}$ for all $n>N_{P}$. Therefore, w.p. $1$,

for all $n>N$, where $\max_{P:f\notin\cup G_{i}}r_{P}<r<1$ and $N>\max_{P:f\notin\cup G_{i}}N_{P}$. 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 $f\in\bar{B}$. By (3.3), w.p. $1$ there exist $s,N>0$ such that

for all $n>N$. Thus, w.p. $1$,

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 “$f\sim g$ as $n\to\infty$” mean that $\lim_{n\to\infty}f(n)/g(n)=1$, $v(i)$ denotes the $i$th element of vector $v$, $M(i,j)$ denotes the $i$th row, $j$th column element of matrix $M$, $0_{a,b}$ denotes the all-zero $a\times b$ matrix, and the sample mean and unbiased sample covariance of features in block $A$ and class $y$ are denoted by $\hat{\mu}_{y}^{A}$ and $\hat{\Sigma}_{y}^{A}$, respectively.