Information Theoretic Methods for Variable Selection—A Review

Abstract

We review the principal information theoretic tools and their use for feature selection, with the main emphasis on classification problems with discrete features. Since it is known that empirical versions of conditional mutual information perform poorly for high-dimensional problems, we focus on various ways of constructing its counterparts and the properties and limitations of such methods. We present a unified way of constructing such measures based on truncation, or truncation and weighing, for the Möbius expansion of conditional mutual information. We also discuss the main approaches to feature selection which apply the introduced measures of conditional dependence, together with the ways of assessing the quality of the obtained vector of predictors. This involves discussion of recent results on asymptotic distributions of empirical counterparts of criteria, as well as advances in resampling.

1. Introduction

Conditional independence is one of the main concepts in statistics which plays a fundamental role in such areas as causal inference, dependence analysis, and graphical modelling. In this review, we discuss how measures of dependence considered in information theory are used in classification and regression problems to choose predictors which significantly influence the outcome. This is a vital application of the information theoretic approach to variable selection, also known as feature selection. Let us stress, however, that information-theoretic measures are frequently applied in many other problems in Machine Learning and Statistics, in contexts other than feature selection. Some representative examples include data visualisation (t-SNE, [1]), clustering [2], Independent Component Analysis (ICA [3], Chapter 15), Variational Inference ([4], Chapter 10), and Natural Language Modelling [5], among others.

In the times of the big data challenge, the problem of a choice of a small group of variables from the pool of all potential variables, all of which would accurately describe the changes in response, is gaining in importance. This is needed for better understanding of a studied phenomena, as well as for construction of tractable models for them. Moreover, feature selection is instrumental in avoiding curse of dimensionality problem when a prohibitive amount of data are required for adequate fitting of the model ([6], Chapter 2) and avoiding overfitting (ibid., Chapter 7). This is also important for prediction, as classifiers which avoid using inactive predictors are usually less variable. Feature selection is frequently applied because of cost considerations, especially when measuring some of the potentially useful predictors is costly or inconvenient.

Another competing direction of research with a similar aim is dimensionality reduction, in particular variable extraction, which transforms given variables to obtain a lean group of new predictors. Primary examples of such methods are Principal Components Regression and Partial Least Squares, see, e.g., [6], and more recent approaches based on neural networks [7].

Here, we focus on an important group of selectors, called filters, which, during the selection process, do not take into account the classifier (or regression estimator) which will incorporate selected variables. In contrast, wrappers select features to optimise performance of a specific classifier under consideration (see, e.g., methods which also yield ranking of features, such as [8,9,10]). The property that filters are model-agnostic makes them a universal tool which can be used for any classifier or regression estimator.

In this paper, we try to present critical and, unavoidably, partial assessment of the state of the art for this problem, focusing on results which can be formally proved, and showing motivation, advantages, and limitations of the presented solutions, including the most recent ones. Additionally, we show that the majority of information-based feature selection criteria can be viewed as truncated or truncated and weighted expansions of Möbius decomposition. This sheds new light on similarities and differences between the frequently employed criteria. Moreover, based on the recent results, asymptotic distributions of their plug-in empirical counterparts are discussed. Such results are needed to construct a test of conditional independence which control probability of false alarms. The recent advances in resampling yielding conditionally independent samples which provide alternative method to solve this problem are also discussed. By this, hopefully, a new insight is added to existing reviews (see e.g., [11,12,13,14]).

The paper is organised as follows. We first discuss in Section 2 the main information theoretic objects and their interplay, focusing on Möbius decomposition of conditional mutual information ($CMI$). In Section 3, we introduce and discuss a feature selection problem from the information theory perspective based on measures of dependence introduced in the previous section. The concept of Markov Blanket of a target variable, being the aim of feature selection as the minimal set of predictors containing the whole information about it, is investigated here. In Section 4, we study feature selection criteria related to $CMI$ stressing that most of them are naturally related to Möbius decomposition. Moreover, this section introduces another group of selectors based on variational bounds of information measures. In the following Section 5, we discuss the interplay between various $CMI$-related measures. The next sections cover the approximate distributions of the introduced measures (Section 6) and an alternative method of assessing their distributions using resampling schemes (Section 7). These properties are vital for performance of conditional independence tests, which are building blocks of feature selection algorithms. In Section 8, Markov Blanket discovery algorithms are described. Section 9 discusses feature selection in continuous case, Section 10 covers related problem of interaction detection stressing how these can be incorporated into feature selection approach.

2. Conditional Mutual Information and Related Measures

We start with discussing properties and the role that conditional mutual information ($CMI$) plays in variable selection based on information theoretic concepts. We will focus on a discrete finite case, meaning that the considered variables take a finite number of discrete values. The continuous case is reviewed shortly in Section 9. In the following, Y will denote the class variable whereas X and Z (possibly multivariate) are features which will be used to predict Y. For basic information theoretic concepts we refer to [15,16].

2.1. Mutual Information $MI$

Definition 1

([15], p. 46). The Mutual Information ($MI$) between Y and X is defined as

$$I(Y;X)=\sum _{x_{,}y}P(X=x,Y=y)log\frac{P(X=x,Y=y)}{P(X=x)P(Y=y)}=H\left(Y\right)-H\left(Y\right|X),$$

(1)

where $H\left(Y\right)=-{\sum }_{y}P(Y=y)logP(Y=y)$ and $H\left(Y\right|X)={\sum }_{x}P(X=x)H\left(Y\right|X=x)$ are the entropy and the conditional entropy, respectively.

The second equality in (1) motivates the name Information Gain also used for $MI$. It also yields intuitive meaning for $MI$ as the decrease in variability of Y is measured by its entropy when the information about another variable X is available. We remark that $I(Y;X)$ is frequently denoted as $MI(Y;X)$. Note that the semicolon in $I(Y;X)$ determines between which variables dependence is considered: either between Y and X in the case of $I(Y;X)$ or between Y and $(X,Z)$ in the case of $I(Y;X,Z)$. $I(Y;X)$ evaluates how similar the joint distribution $P_{YX}$ of $(Y,X)$ is to the product $P_Y\otimes P_X$ of their marginal distributions. As $P_Y\otimes P_X$ corresponds to independence of X and Y, $I(Y;X)$ can be considered a measure of strength of dependence between Y and X. If follows from the definition that ([15], Section 2.3)

$$I(Y;X)=KL\left(P_{Y,X}\right||P_Y\otimes P_X),$$

where Kullback–Leibler (KL) divergence between distributions $P_W$ and $P_Z$ is defined as

$$KL(P_W\left|\right|P_Z)=\sum _{w}P(W=w)log\{P(W=w)/P(Z=w)\}.$$

We note that KL divergence is closely related to the Maximum Likelihood (ML) method and popular feature selection method Akaike Information Criterion (AIC) is derived as the bias-corrected ML [17]. We also remark that Kullback–Leibler divergence can be replaced by other pseudo-distance between probability distributions resulting in a different measure of dependence.

It follows from the properties of KL divergence that $I(Y;X)$ is non-negative and is equal to zero if, and only if, $P_{Y,X}=P_Y\otimes P_X$, i.e., when Y and X are independent. Moreover, the definition (1) that Mutual Information is symmetric: $I(Y;X)=I(X;Y)$. It is also easily seen that (see formula (2.450) in [15])

$$I(Y;X)=H\left(Y\right)+H\left(X\right)-H(Y,X).$$

(2)

2.2. Conditional Mutual Information $CMI$

We denote by $P_{X|Z=z}$ conditional distribution of X given $Z=z$ and define $I(Y;X|Z=z)=KL\left(P_{Y,X|Z=z}\right||P_{Y|Z=z}\otimes P_{X|Z=z})$ as the strength of dependence between Y and X given $Z=z$. Thus, $I(Y;X|Z=z)=0$ means that Y and X are independent given that Z equals z. Now we define conditional mutual information.

Definition 2

([15], p. 49). The conditional mutual information ($CMI$) is

$$\begin{array}{ccc}\hfill I(Y;X|Z)& =& E_{Z=z}I(Y;X|Z=z)\hfill \\ & =& \sum _{z}P(Z=z)\sum _{x,y}P(X=x,Y=y|Z=z)log\frac{P(X=x,Y=y|Z=z)}{P(X=x|Z=z\left)P\right(Y=y|Z=z)}.\hfill \end{array}$$

(3)

Thus, the conditional mutual information is the mutual information of Y and X given $Z=z$ averaged over the values of Z. It is measure of dependence between Y and X given the knowledge of X. From the properties of Kullback–Leibler divergence it follows that

$$I(Y;X|Z)=0\iff X\mathrm{and}Y\mathrm{are}\mathrm{conditionally}\mathrm{independent}\mathrm{given}Z.$$

This is a powerful property, not satisfied for other measures for dependence, such as partial correlation coefficient in the case of continuous random variables. The conditional independence of X and Y given Z will be denoted by $X\perp \perp Y|Z$ and abbreviated to CI. We note that since $I(Y;X|Z)$ is defined as a probabilistic average of $I(Y;X|Z=z)$ over $Z=z$, it follows that

$$I(Y;X|Z)=0\iff I(Y;X|Z=z)=0\mathrm{for}\mathrm{any}z\mathrm{in}\mathrm{the}\mathrm{support}\mathrm{of}Z.$$

Thus, formally, testing conditional independence of X and Y given Z is equivalent to testing unconditional dependence of X and Y on every stratum $Z=z$. This, however, does not make the problem an easy task, as sometimes claimed, since performing such tests simultaneously on each strata (global null), even when we have sufficient data to do it, would result in lack of control of the probability of false signals (multiple testing problem).

The above definitions can be naturally extended to the case of random vectors (i.e., when X, Y, and Z are multivariate) by using a multivariate probability mass functions instead of a univariate one. It is also easily seen that the following chain formula holds

$$I(Y;X,Z)=I(Y;X)+I(Y;Z|X).$$

(4)

2.3. Interaction Information $II$

The 3-way interaction information $II(Y;X,Z)$ of, possibly multivariate, $X,Y$, and Z plays also important role in feature selection.

Definition 3.

The 3-way interaction information $II(Y;X,Z)$ is defined as

$$II(Y;X;Z)=I(Y;X|Z)-I(Y;X)=I(Y;X,Z)-I(Y;X)-I(Y;Z).$$

(5)

The second equality, stemming from (4), neatly explains the idea of interaction information: we want to evaluate the synergistic effect (positive or negative) of both X and Z influencing Y, disregarding the sum of their individual effects. Thus, from the evaluated strength of overall dependence between Y and a pair $(X,Z)$ measured by $I(Y;X,Z)$, we subtract the strength of individual dependences of Y on X and Y on Z measured by $I(Y;X)$ and $I(Y;Z)$, respectively. Although it is not immediately clear from the definition, $II(Y;X;Z)$ is actually a symmetric function of its arguments. However, it can be positive, negative, or 0. If it is positive, X and Z are said to be complementary wrt to Y or interacting positively in influencing Y. A smaller than 0 value of $II$ indicates redundancy among features or their inhibition.

We define 2-way as $II$ as

$$II(Y;X)=I(Y;X)\mathrm{and}II(Y;X|Z)=I(Y;X\left|Z\right).$$

Definition (5) can be generalised in a recursive way to define k-way interaction information (see [18,19,20] for alternative equivalent definition). For any subset of indices $S=\{s_1,s_2,\dots ,s_{\left|S\right|}\}$, $Z_S:=(Z_{s_1},Z_{s_2},\dots ,Z_{s_{\left|S\right|}})$ will stand for the sub-vector of $Z_i$s with indices in S. Abusing the notion, $Z\in Z_S$ will be mean that $Z\in \{Z_{s_1},Z_{s_2},\dots ,Z_{s_{\left|S\right|}}\}$. Let $S=\{1,\dots ,|S\left|\right\}$ and $k=\left|S\right|$ be the cardinality of S. We define k-way interaction information $II(Z_1;Z_2;\dots ;Z_{\left|S\right|})$ so that the chain formula holds.

Definition 4.

k-way interaction information $II(Z_1;Z_2;\dots ;Z_{\left|S\right|})$ is defined as

$$II(Z_1;Z_2;\dots ;Z_{\left|S\right|})=II(Z_1;\dots ;Z_{\left|S\right|-1}|Z_{\left|S\right|})-II(Z_1;\dots ;Z_{\left|S\right|-1}),$$

(6)

where, consistently with the definition of the conditional mutual information in (3), we define

$$II(Z_1;\dots ;Z_{\left|S\right|-1}|Z_{\left|S\right|})=\sum _{z_{\left|S\right|}}p\left(z_{\left|s\right|}\right)II(Z_1;\dots ;Z_{\left|S\right|-1}|Z_{\left|S\right|}=z_{\left|s\right|}).$$

Using expansions of $MI$ and $CMI$ in terms of entropies, the following formula for k-way $II$ is obtained

$$II(Z_1;\dots ;Z_{\left|S\right|})=-\sum _{i=1}^{\left|S\right|}\sum _{T\subseteq S,\left|T\right|=i}{(-1)}^{\left|S\right|-\left|T\right|}H\left(Z_T\right),$$

(7)

where $Z_T=(Z_{t_1},\dots ,Z_{t_i})$ for $T=\{t_1,\dots ,t_i\}$. In particular, we have

$$II(Z_1;Z_2;Z_3)=H(Z_1,Z_2)+H(Z_2,Z_3)+H(Z_1,Z_3)-H(Z_1,Z_2,Z_3)-H\left(Z_1\right)-H\left(Z_2\right)-H\left(Z_3\right).$$

Note that when, e.g., $Z_3=XOR(Z_1,Z_2)=I\{Z_1\ne Z_2\}$ and $(Z_1,Z_2)$ are independent copies of a random variable having Bernoulli distribution with $p=1/2$, we obtain $II(Z_1;Z_2;Z_3)=log2$, as in this case $I(Z_3,Z_2|Z_1)=H\left(Z_3\right)$. In the following, so-called Möbius expansion (see, e.g., [20]) plays a crucial role.

Theorem 1.

The conditional mutual information satisfies the equation

$$I(X;Y|Z_S)=I(Y;X|Z_1,\dots ,Z_{\left|S\right|})=\sum _{i=0}^{\left|S\right|}\sum _{T\subseteq S,\left|T\right|=i}II(X;Y;Z_{t_1},\dots ;Z_{t_i}),$$

(8)

where $T=\{t_1,\dots ,t_i\}$ and, conversely,

$$II(X;Y;Z_{s_1},\dots ,Z_{s_{\left|S\right|}})=\sum _{i=0}^{\left|S\right|}\sum _{T\subseteq S,\left|T\right|=i}{(-1)}^{\left|S\right|-\left|T\right|}I(X;Y|Z_T).$$

(9)

We stress that the inner sum ranges over all subsets of S. For a description of set functions for which (8) and (9) are equivalent, see [20]. Note that we carefully distinguish between 3-way interaction information $II(X;Y;Z_S)$ with multivariate $Z_S$ as the third component and $II(X;Y;Z_1,\dots ;Z_{\left|S\right|})$ being $\left(\right|S|+2)$-way interaction information between $X,Y$ and all components $Z_1,\dots ,Z_{\left|S\right|}$. Equality (8) can be restated, in view of $II(X;Y)=I(X;Y)$ and (6) as

$$I(X;Y|Z_S)=I(X;Y)+\sum _{i=1}^{\left|S\right|}\sum _{T\subseteq S,\left|T\right|=i}[II(X;Z_{t_1};\dots ;Z_{t_i}|Y)-II(X;Z_{t_1};\dots ;Z_{t_i})].$$

(10)

Finally, note that it follows from (9) that $II(X;Y;Z_1,\dots ;Z_{\left|S\right|})=0$ provided X and Y are conditionally independent given any subvector $Z_T$ of $Z_S$ including $Z_{\varnothing }$. As we shall see in Section 4.1, Möbius expansion (8) is a natural starting point to introduce $CMI$-related criteria for feature selection.

Let us indicate that for any classifier $\widehat{Y}=Y\left(X\right)$ of class variable Y for g classes, its unconditional probability of error $P(Y\ne \widehat{Y})$ can be related to conditional entropy $H\left(Y\right|X)$ by means of Fano’s inequality ([15,21])

$$H\left(Y\right|X)\le 1+P(Y\ne \widehat{Y})log(g-1),$$

which, using $I(X;Y)=H\left(Y\right)-H\left(Y\right|X)$, can be written as

$$P(Y\ne \widehat{Y})\ge \frac{H\left(Y\right)-I(X;Y)-1}{log(g-1)}.$$

This shows that when $I(Y;X)$ decreases, i.e., Y and X become less dependent, the lower bound on probability of error of any classifier increases.

3. Feature Selection

In this section we discuss an objective of feature selection, the concept of Markov Blanket and its properties, as well as a greedy search for active features.

3.1. General Considerations: Characterisations of Markow Blanket

Suppose now that we consider class variable Y and p-dimensional discrete vector $(X_1,\dots ,X_p)=:(Z_1,\dots ,Z_p)$ of all predictors available. Our aim is to describe Y using those features among $X_1,\dots ,X_p$ which jointly influence Y. Note that this is a different task than selecting features which individually affect Y. Namely, we want to find out which features can be discarded without loss of information about Y when the remaining features are taken into account. Thus, our aim is to check which features become redundant in the presence of other features. Moreover, we would like to investigate synergy between active features, that is a possible additional effect due to their simultaneous acting. Joint informativeness is, thus, crucial for feature selection. In this context, we define a minimal set of active predictors

Definition 5.

We define a minimal subset of active predictors as a minimal subset $S^{*}\subset \{1,\dots ,p\}=:F$, such that

$$I(Y;X_{S^{*}})=\underset{T\subseteq F}{max}I(Y;X_T)=I(Y;X_F),$$

(11)

where $X_T$ denotes the subvector of $(X_1,\dots ,X_p)$ with indices in set T. Minimality is meant in the set-theoretic sense i.e., $S^{*}$ is minimal if there is no proper subset $W\subset S$ which satisfies (11).

Note that the second equality in (11) is due to chain formula as it follows that $I(Y;X_T)\le I(Y;X_{T^{\prime }})$ for $T\subseteq T^{\prime }$.

The problem of determining $S^{*}$ is a feature selection problem stated in information-theoretic setting, as $X_{S^{*}}$ may be considered as the minimal set describing adequately the overall dependence of Y on the available vector of predictors. The other possible way of defining a small subset of predictors which contain ’almost’ all information about target Y is, for a given value of hyperparameter $\gamma >0$, to consider the smallest subset $S_{\gamma }^{*}$, such that

$$I(Y;S_{\gamma }^{*})\ge I(Y;X_F)-\gamma .$$

The other variant of the problem is its constrained version when a solution is sought for specific cardinality k of the feature set

$$I(Y;X_{S^{*}})=\underset{T\subseteq F,\left|T\right|=k}{max}I(Y;X_T),$$

(12)

which involves $\left(\genfrac{}{}{0pt}{}{\left|F\right|}{k}\right)$ evaluations of mutual information. We also refer to the related bottleneck problem in which information in X is compressed to a random variable M satisfying a given compression condition making it easy to transmit, which has maximal mutual information with Y, see, e.g., [22].

Note that since $I(Y;X_S)$ is monotone in the second coordinate ($I(Y;X_S)\ge I(I(Y;X_T)$ for $T\subseteq S$), we have

$$I(Y;X_F)\ge \frac{1}{\left(\genfrac{}{}{0pt}{}{\left|F\right|}{k}\right)}\sum _{T\subseteq F,\left|T\right|=k}I(Y;X_T),$$

thus the RHS can be used as a criterion when looking for $S^{*}$ ([23]).

It is shown in [24] that when the distribution of $(X_1,\dots ,X_p)$ is determined by a subvector $X_{S^{*}}$ corresponding to a certain $S^{*}\subseteq F$ and is parametrised by ${\tau }^{*}$, finding a maximiser of $I(Y;X_T)$ is a first step of maximising the log-likelihood $\mathcal{L}(T,\tau )$ over T and $\tau$ under so called filter assumption stating that optimisations over T and over $\tau$ are independent. The problem of uniqueness of $S^{*}$ satisfying (11) is important and sometimes disregarded as a problem in feature selection.

Remark 1.

Note that although (11) is trivially satisfied for $S^{*}=F$ it does not mean that the minimal set in the sense of inclusion is uniquely defined. The obvious example is constructed by taking any $Y,X$, such that $I(Y;X)>0$ and letting $X_1=X$ and $X_2=f\left(X_1\right)$ where f is any $1--1$ transform which maps a set of values of X onto itself. In this case, $I(Y;X_1)=I(Y;X_2)=I(Y;X_1,X_2)$, thus both subsets $\left\{X_1\right\}$ and $\left\{X_2\right\}$ satisfy (11) and are minimal. Moreover, note that $Y\perp \perp X_2|X_1$. However, uniqueness of $S^{*}$ satisfying (11) holds for the case of continuous X and Y binary under strict positiveness assumptions stating that density $p\left(x\right)$ is positive almost everywhere with respect to Lebesgue measure and $P_X(p(Y=1|X)=1/2)=0$ (see [25], p. 97).

We discuss, now, properties of the set $S^{*}$. The first one is obtained by noticing that in view of the chain formula for $S^{*}$ defined in (11) we have $I(Y;X_{S^{*c}}|X_{S^{*}})=0$, where $T^c$ denotes complement of T in F. This is equivalent to stating that $S^{*}$ is so called Markov Blanket of Y.

Definition 6.

Markov Blanket of target variable Y is the minimal subset $MB\left(Y\right)\subseteq F$, such that

$$X_{MB{\left(Y\right)}^c}\perp \perp Y|MB\left(Y\right),$$

(13)

where $MB{\left(Y\right)}^c=\{1,\dots ,p\}\backslash MB\left(Y\right)$ (Minimal set satisfying (13) is also called Markov boundary).

Thus $MB\left(Y\right)$ shields Y from the rest of predictors in the sense that they become irrelevant once $MB\left(Y\right)$ is known. Again, $MB\left(Y\right)$ does not need to be uniquely defined. In order to discuss properties of $MB\left(Y\right)$ we introduce the concept of strong relevancy.

Definition 7.

$X_i$ is strongly relevant feature provided

$$I(Y;X_i|X_{F\backslash \left\{i\right\}})>0.$$

Note that the last property is equivalent to $Y\cancel{⫫}{X}_i|X_{F\backslash \left\{i\right\}}$ and it can be restated as $P_{Y|X_F}\ne P_{Y|X_{F\backslash \left\{i\right\}}}$ [24]. Intuitively, strongly relevant features should be included in $S^{*}$ as, after exclusion of such features, the dependence of Y on $X_F$ is not adequately described. In the class of General Linear Models, under minimal conditions, features which are not strongly relevant are exactly those for which corresponding regression coefficient equals 0 (see [26], Proposition 2.2). The following facts concerning $S^{*}$ have been formally proved:

Theorem 2.

Assume that Markov Blanket $MB\left(Y\right)$ exists. Then, we have:

(i) 

$S^{*}$ coincides with Markov Blanket $MB\left(Y\right)$;

(ii) 

$S^{*}$ satisfies

$$\forall T\subset F,T\subseteq F\backslash S^{*}\iff Y\perp \perp X_T|X_{S^{*}\backslash T};$$

(14)

(iii) 

Every strongly relevant feature belongs to $S^{*}$.

The first statement is justified above. Part (ii) is due to [27] and indicates that all subsets of F are partitioned into two parts: the first consisting of sets T disjoint with $S^{*}$ for which $Y\perp \perp X_T|X_S^{*}$ holds and the remaining ones which are conditionally dependent with Y given $X_{S^{*}\backslash T}$.

We note that (iii) can be justified by the following simple reasoning: it is enough to show that if $j\notin S^{*}$ then $X_j$ is not strongly relevant. Indeed, we have

$$\begin{array}{ccc}& & I(Y;X_{S^{*}})=I(Y;X_F)=I(Y;X_{F\backslash \left\{j\right\}})+I(Y;X_j|X_{F\backslash \left\{j\right\}})\hfill \\ & =& I(Y;X_{S^{*}})+I(Y;X_{F\backslash \{\left\{j\right\}\cup S^{*}\}}|X_{S^{*}})+I(Y;X_j|X_{F\backslash \left\{j\right\}})\ge I(Y;X_{S^{*}}),\hfill \end{array}$$

where the second equality follows from the chain formula and the third from the fact that $j\notin S^{*}$ and the chain formula again. Thus, from the second equality it follows that $I(Y;X_j|X_{F\backslash \left\{j\right\}})=0$ whence $X_j$ is not strongly relevant. It is not true in general that $MB\left(Y\right)$ consists only of strongly relevant features. Note that in the last example both $X_1$ and $X_2$ substituted for $MB\left(Y\right)$ satisfy (13), but neither of them is strongly relevant as $Y\perp \perp X_1|X_2$ and $Y\perp \perp X_2|X_1$. In the case when X is continuous and Y is binary, this can be proved for strictly positive distributions defined in Remark 1 (cf. Theorem 10 in [28]).

Ref. [27] introduces a concept of m-Markov Blanket $S_m^{*}$ for $m\le p$ for which equivalence in (14) is satisfied for any subset $T\subseteq F\backslash S_m^{*}$, such that $\left|T\right|\le m$. As for $m=p$ the condition $\left|T\right|\le m$ is vacuous, p-Markov Blanket of Y is $MB\left(Y\right)$.

3.2. Greedy Feature Selection

In the view of (i) of Theorem 2, finding $S^{*}$ is equivalent to finding the minimal set satisfying

$$I(Y;X_{{S^{*}}^c}|X_{S^{*}})=0.$$

This is a difficult task when p is large as it would involve checking this condition for all $2^p$ subsets of $\{1,\dots ,p\}$. The problem is frequently replaced by its greedy selection analogue: given S defined as the set of indices of features already selected, one determines

$${\arg \max }_{j\in S^c}[I(X_{S\cup \left\{j\right\}};Y)-I(X_S;Y)]={\arg \max }_{j\in S^c}I(X_j;Y|X_S),$$

(15)

where the second equality follows from (4). The feature having the largest $MI$ with Y is chosen first and the sequential search is stopped when for all $j\in S^c$ we have $I(X_j;Y|X_S)=0$. Note that, in view of chain equality, the maximised criterion equals

$$I(X_j;Y|X_S)=I(X_j;Y)-I(X_j;X_S)+I(X_j;X_S|Y)=I(X_j;Y)+II(X_j;X_S;Y),$$

(16)

where $II(X_j,X_S,Y)$ is three-way interaction information of $X_j,X_S$, and Y. The first two terms of the first equality are called relevance and redundancy of $X_j$ wrt Y and $X_S$, respectively. Note that in the maximisation process of $I(X_j;Y|X_S)$ over $X_j$, redundancy of $X_j$ with respect to $X_S$ may be outweighed by the magnitude of conditional information which $X_j$ contains about $X_S$ within classes: $I(X_j;X_S|Y)$. RHS of (16) involves $II(X_j,X_S,Y)$. In the next section, we show that $II(X_j,X_S,Y)$ is frequently replaced by algebraic expressions involving $II(X_j;X_{t_1};\dots ;X_{t_k};Y)$ where k does not exceed 2. This is motivated by (8) and allows for easier estimation of the expression.

4. Feature Selection Criteria Related to $\mathit{CMI}$

4.1. Criteria Based on Möbius Expansion

As estimation of conditional quantities, such as $I(Y;X|X_S)$ and its effective use for detecting conditional dependence in case when $X_S$ is high-dimensional, requires large sample sizes (see, e.g., [29], Section 3.2), the common approach is to modify expressions, such as RHS of (10) to define analogues of $CMI$. Consider two approaches to achieve this aim. The first is based on truncation of the sum in (10) at a certain order, the second consists of weighing the corresponding terms and truncating them. In both cases, conditioning by random variable $X_S$ is replaced with conditioning by low-dimensional sub-vectors of these vector, which drastically reduces the need for large sample sizes. This makes them easier to estimate and analyse. As a motivational example, consider the situation when $X_S$ consists of $\left|S\right|=10$ binary coordinates, uniformly distributed on $\{0,1\}$ and the sample size $n=1000$. Then, we expect around $1000/2^{10}\approx 1$ observation for each value $X_S=x_S$ of the conditioning variable $X_S$ in $I(Y;X|X_S)$ whereas the respective number is around 500 for a single conditioning variable $X_i$ appearing in the definition of, e.g., $JMI$ in (21) below. The problem that low dimensionality of conditioning set facilitates estimation is recognised (see, e.g., [29]).

More specifically, we will consider below the following criteria based on truncation:

• Mutual Information Maximisation criterion $MIM$;
• Conditional Infomax Feature Selection criterion $CIFE$ of order two and three.

Additionally, the following criteria based on truncation and weighing of Möbius expansion will be considered:

• Generalised Information Criterion $GIC$;
• Joint Mutual Information criterion $JMI$ of order two and three;
• Mutual Information Feature Selection criterion $MIFS$;
• Minimum-Redundancy Maximum-Relevance criterion $mRMR$.

The simplest approach is to consider only the first term of expansion (10) leading to Mutual Information Maximisation criterion $MIM\left(X\right):=I(X;Y)$ (cf. [30]). This completely ignores dependence structure between Y and $X_S$, thus taking into account only feature relevance and disregarding its redundancy. However, it is useful and frequently applied for preliminary screening of predictors which are than subjected to more precise scrutiny. Actually, the name ’filters’ is frequently used for exactly such criteria when $MI$ is replaced by dependence measure between Y and X of user’s choice.

Consideration of the first two summands of expansion (10) leads to Conditional Infomax Feature Selection ($CIFE$, [31], see also [24])

$$CIFE\left(X\right)=CIFE(X,Y|X_S)=I(X;Y)+\sum _{i\in S}II(Y;X;X_i),$$

(17)

which is also called Short Expansion of $CMI$ of order 2 ($SECMI2$) in [29].

We stress that, in the definition of $CIFE$ and in definitions of the following criteria, the argument of the criterion is $X:=X_i$ for $i\in S^c$ which is the variable over which maximisation is performed.

Analogously, taking into account the first three summands in (8) yield [32]

$$CIFE3\left(X\right)=I(X;Y)+\sum _{i\in S}II(Y;X;X_i)+\sum _{i<j,i,j\in S}II(Y;X;X_i;X_j),$$

(18)

which is also called, stressing its relation to $CMI$, $SECMI3$. Thus, for $\left|S\right|\le 2$, $CIFE3$ in greedy selection rule yields the same results as (15).

Incorporating weights into (10) leads to the following definition of Generalised Information Criterion ($GIC$): for any $\beta ,\gamma \in R^{\left|S\right|}$ we define

$$I_{\beta ,\gamma }\left(X\right)=I(X;Y)+\sum _{k=1}^{\left|S\right|}\sum _{T=\{t_1,\dots ,t_k\}\subseteq S}\left[\gamma \left(k\right)II(X,X_{t_1};\dots ;X_{t_k})\right|Y)-\beta \left(k\right)II(X,X_{t_1};\dots ;X_{t_k})\left]\right),$$

(19)

where $\beta =\left(\beta \right(1),\dots ,\beta (k\left)\right)$ and $\gamma =\left(\gamma \right(1),\dots ,\gamma (k\left)\right).$ Usually $\beta \left(l\right)=\gamma \left(l\right)=0$ for $l\ge l_0$, where $l_0$ is a predefined small integer. Letting $\beta \left(l\right)=\gamma \left(l\right)=0$ for $l\ge 2$ one obtains criteria introduced in [24] parametrised by $\beta$ and $\gamma$, where, abusing the notion slightly, $\beta :=\beta \left(1\right)$ and $\gamma :=\gamma \left(1\right)$:

$$J_{\beta ,\gamma }\left(X\right)=I(X;Y)+\gamma \sum _{i=1}^{\left|S\right|}I(X;X_i|Y)-\beta \sum _{i=1}^{\left|S\right|}I(X;X_i).$$

(20)

For $\beta =\gamma ={\left|S\right|}^{-1}$Joint Mutual Information criterion ($JMI$) is obtained [33]:

$$JMI\left(X\right)=I(Y;X)+\frac{1}{\left|S\right|}\sum _{i=1}^{\left|S\right|}(I(X;X_i|Y)-I(X,Z_i))=I(Y;X)+\frac{1}{\left|S\right|}\sum _{k=1}^{\left|S\right|}II(Y;X;X_i).$$

(21)

Note that, in comparison to $SECMI2$ in (17), interaction terms are down-weighted by a factor ${\left|S\right|}^{-1}$.

Remembering that $II$ is symmetric and writing $II(Y;X,Z_i)=I(Y;X|Z_i)-I(Y;X)$ one arrives at an useful form of $JMI$

$$JMI\left(X\right)=\frac{1}{\left|S\right|}\sum _{i=1}^{\left|S\right|}I(Y;X|X_i),$$

(22)

which, in particular, shows that $JMI=0$ is equivalent to $Y\perp \perp X|X_i$ for any $i\in S$.

Remark 2.

Note that $JMI\left(X\right)$ is always non-negative, but it is not the case for $CIFE\left(X\right)$ (and neither for $CIFE3\left(X\right))$. Indeed, taking arbitrary X such that $I(Y,X_1)>0$ and letting $X_2=\cdots =X_p=X_1$ it is easy to check that w have for for $p>2$:

$$CIFE(X_1,Y|X_2,\dots ,X_p)=I(X_1;Y)+\sum _{j=2}^p(I(X_j;Y|X_1)-I(X_j;Y))=-(p-2)I(X_1;Y)<0$$

as $I(X_j;Y|X_1)=0$ (cf. also [34]).

Letting $\gamma =0$ and $\beta \in [0,1]$ one obtains Mutual Information Feature Selection ($MIFS$) criterion [35]

$$MIF{S}_{\beta }\left(X\right)=I(X;Y)-\beta \sum _{k=1}^{\left|S\right|}I(X;X_i),$$

with a special case $\beta =1/\left|S\right|$ called Minimum-Redundancy Maximum-Relevance criterion (mRMR) introduced in [36]. A normalised version of the $MIFS$ criterion was considered in [37].

In [38], 3-way $JMI$ has been introduced by starting from the equality

$$\sum _{\{i,j\}\subseteq S}I(Y;X,X_i,X_j)=\sum _{\{i,j\}\subseteq S}I(Y;X_i,X_j)+\sum _{\{i,j\}\subseteq S}I(Y;X|X_i,X_j).$$

Note that by dropping the first sum on RHS which does not depend on X (as the introduced criteria will be used to choose X) and scaling the second term by $2/\left(\right|S\left|\right(\left|S\right|-1)$ one obtains in view of (8)

$$\begin{array}{ccc}\hfill JMI3\left(X\right)& =& \frac{2}{S\left|\right(\left|S\right|-1)}\sum _{\{i,j\}\subseteq S}I(Y;X|X_i,X_j)\hfill \\ & =& I(Y;X)+\frac{2}{\left|S\right|}\sum _{i\in S}II(Y;X;X_i)+\frac{2}{\left|S\right|\left(\right|S|-1)}\sum _{i<j}II(Y;X;X_i;X_j).\hfill \end{array}$$

(23)

This is generalised information criterion (20) with $\beta \left(1\right)=\gamma \left(1\right)=2/\left|S\right|$ and $\beta \left(1\right)=\gamma \left(1\right)=2/\left(\right|S\left|\right(\left|S\right|-1)$ and all other coefficients equal to 0.

The important criterion Conditional Mutual Information Maximisation $CMIM$ using a different approach consisting of considering non-linear function of $I(Y;X|Z_i),i=1,\dots ,|S|$ has been proposed in [39]

$$CMIM\left(X\right)=\underset{j\in S}{min}I(Y;X|X_j)=I(Y;X)-\underset{j\in S}{max}[I(X;X_j)-I(X;X_j|Y)),$$

(24)

which should be maximised over $X:=X_i$ for $i\in S^c$. Thus, we look for the best surrogate of $X_i$ among already chosen variables and the candidate having the worst the best surrogate of X is chosen. Generalisation of the rule based on $CMIM$ is considered in [40].

4.2. Variational Approach

The other promising approach to approximate $MI$ and $CMI$ is based on construction of variational lower bounds of these quantities. The bounds obtained are then used as selection criteria. We discuss two approaches which are similar in nature. The first one is based on Donsker–Varadhan inequality which states that [41]

$$I(Y;X)\ge \underset{f\in \mathcal{F}}{sup}[E_{P_{X,Y}}f(X,Y)-log{E}_{P_X\otimes P_Y}{e}^{f(X,Y)}],$$

(25)

and inequality becomes equality when family $\mathcal{F}$ contains the logarithm of the ratio of densities of $P_{X,Y}$ and $P_X\otimes P_Y$, namely $h(x,y)=log(f_{XY}(x,y)/f_X\left(x\right)f_Y\left(y\right))$. Indeed, note that by plugging $h(x,y)+C$ into the expression on the RHS of (25) we obtain equality. Estimation of $I(X;Y)$ based on the Donsker–Varadhan formula has been proposed in [42] where neural network is applied for $\mathcal{F}$ resulting in MINE method. The approach has been further pursued in [43,44]. Other lower bounds, such as Nguyen–Wainwright–Jordan bound [45] can also be used. Note that in order to evaluate the expected value under independence in (25) one uses permutations of the original sample which consists in permuting X values while keeping the values of Y at their original places (see Section 7). For feature selection, one can either directly evaluate $I(Y;X_{S\cup \left\{j\right\}})$ or $I(Y;X_j|X_S)$. For the second approach resampling methods which would ensure conditional independence for generated data are needed. Note that the efficiency of such methods depends on flexibility of function class $\mathcal{F}$ over which the bound in (25) is optimised. The empirical evidence suggest that relatively simple neural nets are sufficiently flexible to yield satisfactory estimators of $I(Y;X)$. They are also preferable, as large networks will increase the variability of the estimators.

The second bound, similar in flavour to (25), is due to [46]. It is noticed there that for arbitrary density $q(x,y)$ we have

$$\begin{array}{ccc}\hfill I(Y;X)& =& H\left(Y\right)-H\left(Y\right|X)=H\left(Y\right)+E_{P_X}{E}_{P_{Y|X}}logp\left(Y\right|X)\hfill \\ & \ge & H\left(Y\right)+E_{P_X}{E}_{P_{Y|X}}logq\left(Y\right|X)=E_{P_{X,Y}}log\frac{q\left(X\right|Y)}{p\left(Y\right)},\hfill \end{array}$$

where the inequality is due to $D\left(p\right(Y\left|x\right)\left|\right|q\left(Y\right|x\left)\right)\ge 0$. Assuming that $q\left(y\right)=p\left(y\right)$ in view of Bayes theorem it is now sufficient to specify marginal density $q\left(x\right|y)$ to obtain a lower bound on $I(Y;X)$. In [46] $q\left(x\right|y)$ are considered which either satisfy naive Bayes assumption or more general assumption that $q\left(x_t\right|y,x_1,\dots ,x_{t-1})$ is geometric mean of $q\left(x_t\right|y,x_i)$ for $i=1,\dots ,t-1$. Other proposals, using specific tree-representation for $q\left(x\right|y)$ are possible (see, e.g., [47]).

5. Interplay between $\mathit{CMI}$ and $\mathit{CMI}$-Related Criteria

The following result states conditions under which $CMI$ and one of the introduced criteria coincide. As before, X denotes a feature from $X_{F\backslash S}$ considered as a candidate to be added to $\{X_s,s\in S\}$.

Theorem 3.

(i) 

Assume that all features are conditionally independent given class (naive Bayes assumption): and features in S are conditionally independent given any not chosen feature (i.e., belonging to $S^c$). Then, for $X\in X_{F\backslash S}$, $I(Y;X|X_S)$ differs from $MIF{S}_1\left(X\right)=I(Y;X)-{\sum }_{i\in S}I(X;X_i)$ by a factor which does not depend on X.

(ii) 

Assume that all k-way interaction informations $II(Y;X;X_{t_1};\dots ;X_{t_l})=0$ for $k>3$ ($l>1$). Then, $I(Y;X|X_S)=CIFE(Y,X|X_S)$ for any $X\in X_{F\backslash S}$.

(iii) 

Ref. [24] If $X_i,i\in S$ are conditionally independent given X, for $X\in X_{F\backslash S}$ and additionally they are conditionally independent given X and Y, then $I(Y;X|X_S)$ differs from $CIFE(Y,X|X_S)$ by a factor which does not depend on X.

(iv) 

Ref. [48] Assume that for any $i\in S$ we have $X\perp \perp X_{S\backslash \left\{i\right\}}|X_i$ and additionally $X\perp \perp X_{S\backslash \left\{i\right\}}|X_i,Y$. Then, $I(Y;X|X_S)=JMI(Y;X|X_S)$.

For completeness, the proof is included in the Appendix A. Note that the conditions imposed by (i) are very stringent. There are no known probabilistic vectors satisfying $I(Y;X;X_{t_1},\dots ;X_{t_l})=0$ for $l>1$ apart from simple situations, such as $X_S\perp \perp (X,Y)$, which obviously does not hold when $X_S$ is the set of chosen predictors.

We remark that under conditions of (i) $MIF{S}_1\left(X\right)=CIFE\left(X\right)$, as due to the naive Bayes assumption ${\sum }_{i\in S}I(X_i;X|Y)=0$. Additionally, it follows from (21) that $JMI\left(X\right)=mRMR\left(X\right)$ provided that the naive Bayes condition holds, thus in order to have $JMI\left(X\right)=CMI\left(X\right)=mRMR\left(X\right)$ under rather strong assumptions of (iv), we have to assume additionally naive Bayes condition. This indicates that the last equality can hold only under restrictive conditions and is not true in general (see [24], Section 4.1).

Assumptions in (iii) and (iv) are not compatible as they lead to different forms of criteria equivalent to $CMI$ criterion. It is argued in [46] that the only plausible graph representation of probability structure for which conditions of (iii) is graph with edges $E=\left\{(Y,X)\right\}\cup {\left\{(X,X_i)\right\}}_{i\in S}\}$. In this case, due to data-processing inequality, we have $I(Y;X)\ge I(Y;X_S)$. This means, however, that X should have been chosen before any features from S.

Additional results of similar flavours to Theorem 3 are discussed in [38,49].

It follows from preceding discussion that criteria introduced in Section 4.1 are formally introduced by truncation (or truncation and weighing) of terms in Möbius expansion. Their analytical properties concerning how well they approximate $CMI$ have yet to be established.

6. Asymptotic Distributions of Information-Theoretic Measures

Sequential feature selection for predicting target variable Y typically involves checking whether a new candidate feature X is a valuable addition to the set of already chosen features $X_S$. This is usually based on testing whether null hypothesis $H_0:X\perp \perp Y|X_S$ holds and its rejection is interpreted as an indication that X carries an additional information about Y to that provided by $X_S$. To this end, $\widehat{I}(Y;X|X_S)$ or its modified versions are used. The usual strategy following statistical testing approach is to derive asymptotic distribution of $\widehat{I}(Y;X|X_S)$ or its modifications described in Section 4 under the null. This distribution is used as a benchmark for which value of the statistic for the sample under consideration is compared to obtain the asymptotic p-value of the test. In the following, we describe the asymptotic distribution of $\widehat{CMI}$ and $\widehat{CMI}$-related criteria under $H_0$. A competing approach to approximate the distribution of the considered statistic under $H_0$ based on resampling is described in Section 7.

6.1. Asymptotic Distribution of $\widehat{CMI}$

We assume that $X,Y,X_S$ take $I,J,K$ possible values, respectively. As $X_S$ is a $\left|S\right|$-dimensional vector, K is the number of all possible combinations of values of its coordinates. We let $p(x,y,x_S)=P(X=x,Y=y,X_S=x_S)$ and consider the case when all probabilities $p(x,y,x_S)$ are positive. It is assumed throughout that the estimation of $I(Y;X|X_S)$ and related quantities is based on n independent identically distributed (iid) samples from the distribution of $(X,Y,X_S)$. Construction of estimators of $CMI$ relies on plugging-in frequencies in place of unknown probabilities, e.g.,

$$\widehat{I}(Y;X|X_S)=\sum _{x,y,x_S}\widehat{p}(x,y,x_S)log\frac{\widehat{p}(x,y|x_S)}{\widehat{p}\left(x\right|x_S)\widehat{p}\left(y\right|x_S)}=\sum _{x,y,x_S}\widehat{p}(x,y,x_S)log\frac{\widehat{p}(x,y,x_S)\widehat{p}\left(x_S\right)}{\widehat{p}(x,x_S)\widehat{p}(y,x_S)},$$

(26)

where $\widehat{p}(x,y,x_S)=n(x,y,x_S)/n$ and $n(x,y,x_S)$ is a number of samples equal to $(x,y,x_S)$. We will consider only frequencies as estimators of discrete probabilities. Other estimators exist, e.g., regularised versions of sample frequencies, for which regularisation reflects a level of departure from conditional independence, see [38,50]. The following known result is frequently used in dependence analysis (compare [51], see also [52]).

Theorem 4.

(i) 

Assume that $I(Y;X|X_S)\ne 0$. Then, we have

$$n^{1/2}(\widehat{I}(Y;X|X_S)-I(Y;X|X_S))\stackrel{d}{\to }N(0,{\sigma }_{\widehat{CMI}}^2),$$

(27)

where

$${\sigma }_{\widehat{CMI}}^2=\sum _{x,y,x_S}p(x,y,x_S){log}^2\frac{p(x,y,x_S)p\left(x_S\right)}{p(x,x_S)p(y,x_S)}-I^2(X,Y|X_S)=\mathrm{Var}\left(log\frac{p(X,Y,X_S)p\left(X_S\right)}{p(X,X_S)p(Y,X_S)}\right)$$

and ${\sigma }_{\widehat{CMI}}^2>0$.

(ii) 

Assume that $I(Y;X|X_S)=0$. Then,

$$2n\widehat{I}(Y;X|X_S)\stackrel{d}{\to }{\chi }_d^2,$$

(28)

where $d=(I-1)(J-1)K$.

The frequently applied test of conditional independence $X\perp \perp Y|X_S$ is based on the above useful fact (ii) that under independence asymptotic distribution $\widehat{I}(Y;X|X_S)$ does not depend on the distribution of $(X,Y,X_S)$ and is chi-square with the known number of degrees of freedom. On the other hand, when the conditional independence does not hold, the limiting distribution of $\widehat{I}(Y;X|Z)-I(Y;X|X_S)$ is normal with the variance depending on the underlying probability distribution $P_{XYZ}$. We stress that speeds of convergence of $\widehat{I}(Y;X|X_S)$ to $I(Y;X|X_S)$ are different in both cases: they equal $n^{-1}$ in the first case and $n^{-1/2}$ in the second.

The test based on $CMI$ is a popular tool in dependence analysis, in particular for Markov Blanket discovery (see Section 8). It has different names among which $G^2$ test is the most popular (see, e.g., [53]). Additionally, in the literature $X^2$ denotes the second order approximation of $\widehat{CMI}$, which turns out to be the conditional chi-square test and has the same asymptotic distribution as $\widehat{I}(Y;X|X_S)$.

6.2. Asymptotic Distribution of Modified Criteria

Asymptotic distribution of the modified criteria related to $\widehat{CMI}$ can be also derived. Let $\widehat{J^{\beta ,\gamma }}(X,Y|Z)$ be a plug in-version of $J^{\beta ,\gamma }$ defined in (20). Moreover, let $p=\left(p(x,y,x_s)\right)$ be a vector of probabilities of dimension $M=I\times J\times K$, $\widehat{p}$ corresponding vector of fractions and $f:{[0,1]}^M\to R$ be a function which represents $J^{\beta ,\gamma }$ as a function of p, i.e., $J^{\beta ,\gamma }=f\left(p\right)$. For example for $JMI$ criterion (see (22)) the corresponding function is defined as

$$f_{JMI}\left(p\right):=\sum _{x,y,x_S}p(x,y,x_S)\frac{1}{\left|S\right|}\sum _{s\in S}ln\frac{p(x,y,x_s)p\left(x_s\right)}{p(x,x_s)p(y,x_s)},$$

where $x_s$ denotes coordinate of $x_S$. We state here the result on asymptotic distribution of $\widehat{J^{\beta ,\gamma }}(X,Y|X_S)$ proved in [29]. Let ${\mathsf{\Sigma }}_{x,y,x_S}^{x{,}^{\prime }{y}^{\prime },x_S^{\prime }}$ denote an element of matrix $\mathsf{\Sigma }$ with row index $x,y,x_S$ and column index $x^{\prime },y^{\prime },x_S^{\prime }$.

Theorem 5.

(i) 

We have

$$n^{1/2}(\widehat{J^{\beta ,\gamma }}(X,Y|X_S)-J^{\beta ,\gamma }(X,Y|X_S))\stackrel{d}{\to }N(0,{\sigma }_{\widehat{J}}^2),$$

(29)

where ${\sigma }_{\widehat{J}}^2=Df{\left(p\right)}^T\mathsf{\Sigma }Df\left(p\right)=Var\left(Df{\left(p\right)}^T\widehat{p}\right)$ and $\mathsf{\Sigma }=n{\mathsf{\Sigma }}_{\widehat{p}}$ is a matrix consisting of elements ${\mathsf{\Sigma }}_{x,y,x_S}^{x^{\prime },y^{\prime },x_S^{\prime }}=p(x^{\prime },y^{\prime },x_S^{\prime })(I(x=x^{\prime },y=y^{\prime },x_s=x_S^{\prime })-p(x,y,x_S))/n$.

(ii) 

If ${\sigma }_{\widehat{J}}^2=0$ then

$$2n(\widehat{J^{\beta ,\gamma }}(X,Y|X_S)-J^{\beta ,\gamma }(X,Y|X_S))\stackrel{d}{\to }{V}^{T}HV,$$

(30)

where V follows $N(0,\mathsf{\Sigma })$ distribution, and $H=D^{2}f\left(p\right)$ is a Hessian of f.

In particular, we have the following result for $\widehat{JMI}$ [54]:

Corollary 1.

Let Y be binary.

(i) 

If ${\sigma }_{\widehat{JMI}}^2\ne 0$ then

$$n^{1/2}(\widehat{JMI}-JMI)\stackrel{d}{\to }N(0,{\sigma }_{\widehat{JMI}}^2),$$

where

$$\begin{array}{cc}\hfill {\sigma }_{\widehat{JMI}}^2& =\sum _{x,y,x_S}p(x,y,x_S){\left(\frac{1}{\left|S\right|}\sum _{s\in S}ln\frac{p(x,y,x_s)p\left(x_s\right)}{p(x,x_s)p(y,x_s)}\right)}^2\hfill \\ \hfill & -{\left(JMI\right)}^2\hfill \end{array}$$

(31)

(ii) 

We ${\sigma }_{\widehat{JMI}}^2=0\iff JMI=0$. In this case,

$$2n\widehat{JMI}\stackrel{d}{\to }{V}^{T}HV=\sum _{i=1}^M{\lambda }_i{Z}_i^2,$$

where V and H are defined in Theorem 2, $Z_i$ are iid $N(0,1)$ and ${\lambda }_i$ are eigenvalues of the matrix $W=\mathsf{\Sigma }H$ which has the following elements

$$W_{x,y,x_S}^{x^{\prime },y^{\prime },x_S^{\prime }}=\frac{1}{\left|S\right|}\sum _{s=1}^{\left|S\right|}\left[\frac{I(x_s=x_s^{\prime })}{p\left(x_s\right)}-\frac{I(x=x^{\prime },x_s=x_s^{\prime })}{p(x,x_s)}-\frac{I(y=y^{\prime },x_s=x_s^{\prime })}{p(y,x_s)}+\frac{I(x=x^{\prime },y=y^{\prime },x_s=x_s^{\prime })}{p(x,y,x_s)}\right].$$

It follows from Theorem 1 that if for any $i\in S$ we have that $I(Y;X|X_i)=0$, i.e., Y and X are conditionally independent given $X_i$, then asymptotic distribution is that of quadratic form specified in (ii), otherwise the distribution is normal.

The main advantage of using modified criteria instead of $CMI$ is that their estimation does not require as large samples as for $CMI$ itself. Note that for modifications of order k, conditioning involves k-dimensional strata, whereas for $CMI$ p-dimensional strata are considered. However, modified criteria considered in Section 4.1 suffer from the fact that under hypothesis of conditional independence $X\perp \perp Y|Z$ the asymptotic distribution of empirical criterion is not uniquely determined and has to be estimated from the sample. As in the case of $\widehat{JMI}$, the asymptotic distribution can be either normal (when $X\cancel{⫫}Y|X_j$ for at least one $j\in S$ or coincides with distribution of the quadratic form. Which type of asymptotic distribution is valid can be decided using resampling schemes shortly discussed in Section 7. The chosen distribution can used as a benchmark to test the conditional independence hypothesis (see, e.g., [29,55]). Alternatively, testing of $H_0:I(Y;X|X_S)=0$ can be replaced by testing ${\tilde{H}}_0:I(Y;X|X_i)=0$ for any $i\in S$ at each stage of forward feature selection procedure and the benchmark distribution can be obtained by approximating eigenvalues of M specified in (ii) (see [54]) or using scaled and shifted $\chi$ squared distribution $\alpha +\beta {\chi }_d^2$, where $\alpha ,\beta ,d$ are estimated from the sample [56]. Note that, although $H_0$ and ${\tilde{H}}_0$ are not equivalent, in the case of faithful distributions (See Section 9.2 for definition of faithfulness) ${\tilde{H}}_0$ implies $H_0$, as conditional independence is inherited by conditioning supersets in such a case.

7. Resampling Schemes

When using the conditional independence test, which the test statistic is used for, what is very common, the exact or even the asymptotic distribution under conditional independence is not known, the usual practice is to use conditional randomisation (CR) and resampling schemes to approximate this distribution and use the resulting approximation as the benchmark distribution, as described in the beginning of Section 6. Analogously as before, comparison of the value of the statistic based on the observed sample with the benchmark distribution is used to calculate the CR or resampling p-value. It can be also used as alternative way of assessing the distribution of the statistic even when its approximate distribution is known. We briefly discuss these procedures, first for a discrete $X_S$, then for a continuous case.

• Conditional Randomisation (CR).In the case of the $CR$ approach we assume that the probability mass function $p\left(x\right|X_S=x_S)$ or density of X given $X_S=x_S$ is known. Then, given a sample $(\mathbf{X},\mathbf{Y},{\mathbf{X}}_S):={(X_i,Y_i,X_{S,i})}_{i=1}^n$ one generates a CR sample $({\mathbf{X}}^{*},\mathbf{Y},{\mathbf{X}}_S)={(X_i^{*},Y_i,X_{S,i})}_{i=1}^n$, when $X_i^{*}$ are independently sampled from p.m.f. or density $p\left(x\right|X_S=x_{S,i})$ for $i=1,\dots ,n$. Let $T(\mathbf{X},\mathbf{Y},{\mathbf{X}}_S)$ be a test statistic which we would like to use for testing CI and consider M generated $CR$ samples. We define the permutation p-value as

  $$p_{CR}=\frac{\#\{k=1,\dots ,M:T({\mathbf{X}}_k^{*},\mathbf{Y},{\mathbf{X}}_S)\ge T(\mathbf{X},\mathbf{Y},{\mathbf{X}}_S)\}+1}{M+1}.$$

  (32)
  It follows that $p_{CR}$ is a valid p-value, i.e., conditionally on $((\mathbf{Y},{\mathbf{X}}_S)$, i.e., its distribution is uniform on the set $\{1/M,2/M,\dots ,M/(M+1\left)\right\}$ under CI and, thus,

  $$P(p_{CR}\le \alpha |\mathbf{Y},\mathbf{Z})\le \alpha .$$
  This follows from the following simple fact. Suppose that CI holds and ${\mathbf{X}}^{*}$ is such that $P_{{\mathbf{X}}^{*}|{\mathbf{X}}_S}=P_{\mathbf{X}|{\mathbf{X}}_S}$. Then, given $\mathbf{Y},{\mathbf{X}}_S$, when CI holds, we have the following equality in distribution

  $$T(\mathbf{X},\mathbf{Y},{\mathbf{X}}_S)\stackrel{d}{=}T({\mathbf{X}}^{*},\mathbf{Y},{\mathbf{X}}_S);$$
  (see, e.g., [57]). For recent improvements of the CR method, random permutations are appropriately chosen and considered instead of choosing them at random, see [58]. In the case when $p\left(x\right|X_S=x_S)$ is unknown, CR sampling is replaced by Conditional Permutations or Bootstrap X method.
• Conditional permutation ($CP$). Conditional permutation method is similar to $CR$ method and differs in that for value $x_S$ taken by elements of ${\mathbf{X}}_S$ we consider the strata of the sample corresponding to this value, namely

  $$P_i=\{j:(X_j,Y_j,X_{S,j}):X_{S,j}=x_{S,i}\}.$$
  The CP sample is obtained from the original sample by replacing $(X_j,Y_j,X_{S,j})$ for $j\in P_i$ by $(X_{{\pi }^i\left(j\right)},Y_j,X_{S,j})$, where ${\pi }^i$ is a randomly chosen permutation of $P_i$. Thus, on every strata $X_S=x_s$ we randomly permute values of corresponding X independently of values of Y. Once M of CP samples are obtained independently in this fashion, we calculate p-value $p_{PC}$ based on them analogously as before.
• Conditional Bootstrap X (CB.X). Instead of permuting values of X on each strata, we draw a bootstrap sample from X observations, that is why we sample them with replacement as many times as is the size of the strata. The remaining steps are as previously described.

We discuss now the resampling for continuous predictors. Feature selection for this case is covered shortly in Section 9. We remark here that the aim, methodology of solutions, and criteria considered are very similar, the main differences consist of technical problems of adequate estimation of information-theoretic measures in the continuous case.

In the case of $X_S$ being continuous, CR methods work as stated, however the situation is more complicated for $CP$ and $CB.X$ method as those depend on transforming the strata of the sample, which degenerate to single observation points in this case. One possible solution is to sample from some estimate $\widehat{p}\left(x\right|X_S=x_S)$, e.g., kernel estimate, however this requires a large number of observations on each strata. The following proposal of constructing pseudo sample distributions of which is close to $P_{X,Y,X_S}$ under CI has been suggested by [59] and consists in using observations having close values to $X_S=x_S$. More specifically, consider observation $(x_i,y_i,x_{S,i})$ and find $k^{th}$ nearest observation to $x_{S,i}$ in $X_S$ space, say, $x_{S,j}$ with corresponding triple $(x_j,y_j,x_{S,j})$. Then, observation of $(x_i,y_i,x_{S,i})$ is replaced in the resampling sample by $(x_j,y_i,x_{S,i})$.

Another technique for data augmentation, applicable in both discrete and continuous cases, is knock-off construction which we now describe. Consider the regression problem involving vector $X=(X_1,\dots ,X_p)$ of features and response Y. Vector $\tilde{X}=({\tilde{X}}_1,\dots ,{\tilde{X}}_p)$ is the knock-off vector for X in this regression problem, if $Y\perp \perp \tilde{X}$ and, moreover, for any $S\subseteq \{1,\dots ,p\}$ we have the following equality in distribution [57]:

$${(X_1,\dots ,X_p,{\tilde{X}}_1,\dots ,{\tilde{X}}_p)}_{swap\left(S\right)}\stackrel{d}{=}(X_1,\dots ,X_p,{\tilde{X}}_1,\dots ,{\tilde{X}}_p)$$

where $swap\left(S\right)$ means swapping entries $X_j$ and ${\tilde{X}}_j$ for any $j\in S$. Thus, in the above sense ${\tilde{X}}_1,\dots ,{\tilde{X}}_p$ are interchangeable with $X_1,\dots ,X_p$ and independent of Y at the same time.

The construction of knock-offs is complicated in general, feasible only in special cases as of now, such as the Gaussian case, and requires knowledge of distribution of $P_X$. Thus, even more information is needed than in the case of CR approach. However, the gains are considerable, as for natural classes of statistics measuring influence of predictors on the response one can compare the performance of $X_1,\dots ,X_p$ with that of their knock-offs and on this basis decide which ones are not strongly relevant. Intuitively, when the performance of $X_j$ measured, e.g., by an absolute value of estimated regression coefficient, is comparable to that of its knock-off, then such a variable should be discarded. This approach yields a bound on FDR of strongly relevant features (Theorem 3 in [26]). Thus, in cases discussed in Section 3, when the set of strongly relevant features is exactly equal $MB\left(Y\right)$, this approach yields a bound on FDR of its recovery. We stress that, in this way, one can analyse properties of the whole procedure of MB determination, and not only that of an individual steps in the algorithm. An additional advantage is that, unlike for the resampling schemes above, only one sample of knock-offs needs to be generated.

8. Markov Blanket Discovery Algorithms

We discuss now Markov Blanket discovery algorithms. Their aim is to solve a feature selection problem posed in Section 3 (see Equation (11)) applying $CMI$ or $CMI$ approximations introduced in Section 4 for which theoretical guarantees are discussed in Section 5. Such algorithms are used as building blocks for conditional independence tests, based either on asymptotic distributions of test statistics discussed in Section 6 or approximations of their distributions based on resampling covered in Section 7.

We discuss three representative examples of Markov Blanket discovery algorithms. Other examples include [60] and algorithms using assumed faithful Directed Acyclic Graph (DAG) representation of the underlying probability structure (for DAG fathfulness of probability distribution based on d-separation (see e.g., [4]), such as HITON-MB [61], MMMB [62], IPC-MB [63], and STMB [64].

• The GS (Grow and Shrink) algorithm [65]. It consists of two phases. Specific ordering of variables in F is considered, e.g., variables may be ordered according to value of $I(Y;X)$ and, then, the variable having the largest value of the mutual information is the first variable chosen. Denoting by S the current set of chosen variables, we pick the first variable (in the considered ordering) which depends on Y given $X_S$ ($I(Y;X|X_S)>0$) and we add it to S, then repeat the step. When there is no longer such a variable among candidates, the first phase is terminated. We call the resulting chosen set $S^{*}$. In the second phase, we remove from $S^{*}$, again using the considered ordering, any variable $X_j$ which is not strongly relevant with respect to the current $S^{*}$, i.e., such that $X_j\perp \perp Y|X_{S^{*}\backslash \left\{j\right\}}$ ($I(Y;X_j|X_{S^{*}\backslash \left\{j\right\}})=0)$ and we let $S^{*}:=S^{*}\backslash \left\{j\right\}$.
• The IAMB (Incremental Association Markov Blanket) [66]. The algorithm is similar to GS with one important difference in the first step. Namely, it disregards initial ordering and S is augmented by the most plausible candidate, that is the variable realising ${max}_{i\in S^c}I(X_i,Y|X_S)$, provided it is not conditionally independent from Y given $X_S.$
• $G{S}^{\left(m\right)}$ ([27], where $m\le p$. $G{S}^{\left(m\right)}$ differs from GS in the growing phase only. Namely, instead of an individual variables with indices in $F\backslash S$ being considered as possible candidates, all subsets T of size not exceeding m are taken into account and the check is performed whether there are conditionally dependent on Y given current S. If this holds $S:=S\cup T$.

It is proved in [27] that $G{S}^{\left(m\right)}$ algorithm yields the m-Markov Blanket (see Section 3 for definition of m-Markov Blanket) of Y. Thus, for $m=p$ we have that the output of $G{S}^{\left(p\right)}$ is $MB\left(Y\right)$. However, since, in the growing phase, all subsets of $F\backslash S$ of size not exceeding m have to be checked which is computationally intensive, in practice k subsets for k large enough are checked at every step.

Note that for such results we assume that condition $Y\perp \perp X_{T\backslash S}|X_S$ or $Y\perp \perp X|X_S$ can be verified. In practice, it is impossible to check conditional independence of X and Y given the set of variables already chosen without error and this has to be replaced with the appropriate test discussed in Section 6. Obviously, as such a test has to be performed at each step, one does not control probability of including false discoveries. False discoveries are dealt with in the second phase, but no formal results exist concerning how likely recovering $MB\left(Y\right)$ is in the case of practical algorithm. However, see [67] for the results concerning the PC algorithm in the Gaussian case and knock-off construction discussed in Section 7.

Another possibility is to omit phase two and stop early phase one, applying stopping rules devised in multiple testing problems to control Family-Wise Error Rate or False Discovery Rate (see, e.g., [68] where a stopping rule using the Holm procedure has been applied for CIFE criterion). These methods work promisingly in practice, however, due to the fact that the set S augmented at each step is data-dependent, also in this case formal results on the consistency of such procedures are not available.

There is no definitive study of empirical performance of the presented criteria and/or selection procedures; note that any criterion considered can be used for any feature selection method described above, thus creating a large number of filter methods. Moreover, those can be evaluated using their classification performance, as well as for synthetic datasets, according to their ability to choose active predictors. Thus, we discuss only the reports on the simplest greedy procedures consisting on the application of discussed criteria on synthetic datasets with a fixed ahead number of chosen features. Authors of [24] discuss, among other things, results for datasets coming from NIPS Feature Selection Data Challange, Gisette, and Madelon (procedures stopped at 200 variables in the first case and 20 in the second). Two interesting points arise from the study: strong performance of $JMI$, which was the second best and co-winner in those cases, with respect to balanced accuracy (BA) and the number of feature chosen, with a strong performance of $CMIM$ (winner in the case of Gisette and together with $CIFE$ co-winner in the case of Madelon). The overall strong performance of $CIFE$ was also confirmed in [68,69]. The second important observation is the failure of performance of $CMI$, which due to scarcity of data, failed to detect conditionally dependent variables very early. This is due to the fact that for a large conditioning set the test becomes very conservative (see also [70] for discussion of this property).

9. Case of Continuous Distribution ${\mathit{P}}_{\mathit{X},\mathit{Y},\mathit{Z}}$

In the following, we shortly discuss conditional independence testing for continuous distributions. The main motivation to include a continuous case in this review, devoted mainly to discrete case, is to underline the strong similarities between these two cases. In particular, we discuss below that all information-theoretic tools defined for the discrete case have their analogues for continuous distributions. Moreover, we note that the selection methods presented till now do not depend on continuity of the underlying measure, once a specific test for conditional independence $X\perp \perp Y|Z$ which takes this into account is used. Thus, the differences between two cases are mostly due to the fact that estimation of information-theoretic measures are much more difficult for continuous distributions. Moreover, their asymptotic behaviour under CI is not known, thus making construction of corresponding tests difficult. On the other hand, for Gaussian case significant simplifications exists due to the existence of closed formula for $CMI$ discussed in Section 9.2.

9.1. General Considerations

Returning to notation of Section 2, we discuss conditional independent testing in a continuous case. In this context, it is enlightening to mention the result in [71], where it is shown that in the continuous case conditional independence testing is a hard problem in the following sense. For any test on a natural family of $P_{X,Y,Z}$ defined in [71] satisfying $X\perp \perp Y|Z$, which uniformly achieves type I error asymptotically not larger than a prescribed significance level $\alpha$ its power for any alternative will be also no larger than $\alpha$. Thus, such a test is useless for discriminating between CI and Conditional Dependence.

The indices defined in Section 2 have their natural analogues in the continuous case; namely probabilities $P(X=x,Y=y,Z=z)$ are replaced by density functions $p(x,y,z)$ and summations by integrals (for details see [15]). What is important for development here is that Conditional Independence $X\perp \perp Y|Z$ in the continuous case is also equivalent to conditional Mutual Information $I(X;Y|Z)$ being 0, as in the discrete case. Thus, the problem boils down to an accurate estimation of $I(X;Y|Z)$ and corresponding testing of CI using constructed estimator as a test statistic. This, however, is also a much harder task than in a discrete case as plug-in estimators of entropy, and conditional and unconditional information when, e.g., kernel estimators of densities (see, e.g., [6], Chapter 6) are plugged-in are unstable and work satisfactorily only when very large samples are available. Much better results are obtained using the variational approach described in Section 4.2. Additionally, some estimators based on the nearest neighbour idea, such as the Kozachenko–Leonenko estimator of entropy [72,73] and its refinements, such as the Kraskov et al. estimator [74] work better than straightforward plug-in estimator. Then, estimators of $MI$ are constructed using (1). Additionally, the simple approach based on discretising underlying variables is frequently used. There are two problems related to this approach. The first one is a loss of information due to this operation. The second, related one, is a difficulty of choosing an appropriate bin size, especially in many dimensions. Too large bin size may result in hiding interesting characteristics of predictors’ distributions and, consequently, a loss of predictive power.

We mention two other approaches for the continuous case. One which is frequently used is kernel based approach which relies on the following fact [75]:

$$X\perp \perp Y|Z\iff Ef(X,Z)h(Y,Z)=0\mathrm{for}\mathrm{any}f\in {\mathcal{F}}_{X|Z},h\in {\mathcal{F}}_{Y|Z},$$

where ${\mathcal{F}}_{X|Z}=\{h(Y,Z):h(Y,Z)=h_0\left(Y\right)-E\left(h_0\left(Y\right)\right|Z),h_0\in L_Y^2\}$. Thus, although conditional independence is not equivalent to conditional covariance being zero, this is true when transformations of vectors $(X,Z)$ and $(Y,Z)$ are allowed. In view of this result it turns out that its possible to find rich enough function spaces, such that the condition that the generalized conditional covariance of transformed X and Y being equal to 0 is equivalent to conditional independence. Moreover, it is the function space that one can consider to appropriately define Reproducing Kernel Hilbert Spaces (RKHS). This, conceptually, is much more involved than checking $I(X;Y|Z)$ being 0, save for the technical problems of testing this condition in the continuous case.

We also mention in this context the second approach, called distillation method [76] which relies on ideas similar to the construction of partial residual plot and resampling.

9.2. Gaussian Case

We review the case when $(Y,X)$ is Gaussian, where $X=(X_1,\dots ,X_p)$. This case offers significant simplifications, as the quantities on which feature selection is based can be explicitly calculated. Namely, direct calculation yields [15]:

$$H\left(X\right)=\frac{1}{2}log\left|{\mathsf{\Sigma }}_X\right|+\frac{p}{2}log\left(2\pi e\right),$$

if $X\sim N({\mu }_X,{\mathsf{\Sigma }}_X)$. If $(X,Y)\sim N({\mu }_{X,Y},{\mathsf{\Sigma }}_{X,Y})$ it follows from (2) that

$$I(X;Y)=\frac{1}{2}log\frac{|{\mathsf{\Sigma }}_{X,Y}|}{|{\mathsf{\Sigma }}_X|\times |{\mathsf{\Sigma }}_Y|},$$

where $\left|\mathsf{\Sigma }\right|$ stands for determinant of matrix $\mathsf{\Sigma }$. In particular for bivariate normal case when $\rho (X,Y)=\rho$ we have $I(Y;X_1)=-(log(1-{\rho }^2))/2$. Whence, in the case of the conditional distribution of (Y, X), given $X_S$, where $X_S\subseteq \{2,\dots ,p\}$ this yields

$$I(Y;X_1|X_S)=-\frac{1}{2}log(1-{\rho }_{par}^2),$$

where ${\rho }_{par}$ is partial correlation coefficient between Y and $X_1$ given $X_S$. Thus, in Gaussian case information-based dependence indices can be expressed in terms of either correlation or partial correlation coefficients.

In order to present the PC algorithm, we first introduce the concept of the faithfulness of the distribution which plays an important role in its construction (see, e.g., [67]). Let $X_0:=Y$ and consider undirected graph $G=(V,E)$, where vertices $V=\{0,1,\dots ,p\}$ correspond to variables $(Y,X_1,\dots ,X_p)$ and $E\subset V\times V$, such that $(j,k)\in E\iff (k,j)\in E$ is the set of edges.

Definition 8.

Distribution $P_{Y,X}$ is faithful to graph G when for any triple $A,B,C\subset V$ we have

$$A \mathit{and} B \mathit{are} \mathit{separated} \mathit{by} C\iff X_A\perp \perp X_B|X_C,$$

(33)

when separation by C means that every path joining elements of A and B has to pass through C. (Implication from LHS to RHS in (33) is called global Markov property.)

In the case when $(Y,X_1,\dots ,X_p)$ is Gaussian and there exists faithful representation of its distribution, PC algorithm [77] reconstructs its dependence structure with probabilistic guarantees. The main tool is clever usage of one of the consequences of faithfulness, namely

$$X_1\perp \perp X_2|X_{C_1}\Rightarrow X_1\perp \perp X_2|X_{C_2}\mathrm{for}\mathrm{any}{C}_2\supseteq C_1.$$

(34)

The algorithm starts from the fully connected graph and removes edges between any vertices which are marginally uncorrelated (and, thus, in view of (34) and Gaussianity, conditionally independent given any subset of variables). Then, it proceeds to remove edges between vertices which correspond to variables which are conditionally independent given a subset of their neighbours of size l, where $l=1,2,\dots$ is increased stepwise. The conditional independence is tested using an empirical partial correlation coefficient ${\widehat{\rho }}_{par}$ and the property that its Fisher transform $log[(1+{\widehat{\rho }}_{par})/(1-{\widehat{\rho }}_{par})]$ is approximately $N(0,1)$ normal under CI (details in [67], Section 13.7).

The case when Y is discrete and $P_{X|Y=y_i}$ are normal is more complicated, as in this scenario distribution of X is that of normal mixture for which no explicit formulae for its entropy exist (see, however, [34] for special cases of mixtures of summands having the same covariance matrix). Authors of [78] use approximations of $MI$ (see (18) in [78]) to address this problem.

To summarise, we note that the section discusses possible approaches to test conditional independence and discover Markov Blanket of target Y when the underlying distribution is continuous. The first task can be performed by using some available estimators of $CMI$, discussed above, such as Kraskov et al. estimator [74] in conjunction with benchmark distribution based on resampling. The other possibility is a kernel method which reduces the problem to checking whether generalised conditional covariance is 0. For the parametric case when $(X,Y)$ is Gaussian, the testing problem can be reduced to testing whether the partial correlation coefficient is 0. Moreover, in this case under faithfulness assumption on the distribution Markov Blanket can be recovered with theoretical guarantees.

10. Interaction Detection

Detection of existing interactions between features in influencing the outcome Y is an important task of data analysis; primary examples being gene–gene interactions in Genome-Wide Association Studies (GWAS) and gene–environment interactions. In this section, we show how information-theoretic measures introduced before can be applied to solve this problem. Moreover, we indicate that interaction detection problems may be viewed as a special case of feature selection. Various indices have been proposed to measure the strength of interactions, one of the most popular tools being interaction information discussed in Section 2 and Section 4.1. There it was considered as a tool to define new selection criteria by truncating and possibly weighing summands in the Möbius expansion of $CMI$. Here, $II$ is adopted as a main tool in detection of interactions. Its most popular competitor is the value of Likelihood Ratio Test statistics for two nested logistic models: an additive model with no interactions and a saturated model taking all possible interactions into account. However, it has been shown in [68] that for logistic model when predictors are independent and at least one of interaction terms is non-zero then $II\ne 0$ but not vice versa. This shows that there are situations when interactions are present and are detected by $II$ but will go undetected using logistic model approach. $II$ is applied as the main tool to find interactions in AMBIENCE package [79] and BOOST package uses so called Kirkwood approximation discussed below, which is closely related to $II$ [80]. Other competitors to LRT and $II$ include Multifactor Dimensionality Reduction (MDR) [81] and Restricted Partitioning Method (RPM) [82].

Now, we will discuss two additional properties of Interaction Information which are useful in interaction analysis, focusing on its three-way variant applied to predict synergistic effect of two variables $X_1$ and $X_2$, say, in determining the target Y.

From the first equality in (5) it follows that when both variables are jointly independent from Y, i.e., $(X_1,X_2)\perp \perp Y$ then $II(X;X_2;Y)=0$ as $I(X_1;X_2|Y)=I(X_1;X_2)$. Although the converse is not true, it can be shown in adversarially constructed examples only and, thus, testing $H_0:II=0$ is usually replaced by (more stringent) hypothesis ${\tilde{H}}_0:(X_1,X_2)\perp \perp Y$.

The other property is insightful representation

$$II(X_1;X_2;Y)=KL(P_{X_1,X_2,Y}\left|\right|{\tilde{P}}_K),$$

where ${\tilde{P}}_K$ is Kirkwood Superposition Approximation with masses assigned to points equal

$${\tilde{p}}_K(x_1,x_2,y)=\frac{p(x_1,x_2)p(x_1,y)p(x_2,y)}{p\left(x_1\right)p\left(x_2\right)p\left(y\right)}$$

is positive but necessarily summing up to 1, mass function supported on values $(x_1,x_2,y)$ of $(X_1,X_2,Y)$. It can be shown that for $\eta$ denoting the summary mass of ${\tilde{P}}_K$ it follows that when $\eta \le 1$ and $II=0$, then $P_{X_1,X_2,Y}$ is equal to its Kirkwood approximation [68].

Let us discuss two commonly used methods of testing that $II(X;X_2;Y)=0$. The most popular one is based on Han’s approximation [20] derived under a stringent assumption of overall independence, which results in considering as the benchmark distribution chi-squared distribution with $(I-1)(J-1)(K-1)$ degrees of freedom, when $I,J,K$ are numbers of values of $X_1,X_2$, and Y, respectively. This, however, may lead to a large number of false signals when the pertaining test is employed as the overall independence is only a very special situation when $II$ vanishes. It is shown in [69] that when $(X_1,X_2)$ are independent of Y, similarly to other measures of conditional dependence discussed in Section 4.1, the approximate distribution is weighted chi-square distribution. Its complete description has been given for $X_1$ and $X_2$ having three values and Y being binary, which includes typical GWAS situation of two loci with two co-dominant alleles. The properties of the corresponding test have not been investigated yet.

The other method uses permuting Y values of the considered sample and in this way obtaining M random samples satisfying ${\tilde{H}}_0:(X_1,X_2)\perp \perp Y$. Then, permutation p-value can be calculated as ${(M+1)}^{-1}(1+{\sum }_{i=1}^{M}I({\widehat{II}}_i\ge \widehat{II}))$ (compare (32)). One can also use chi-square distribution with the number of degrees of freedom as the benchmark distribution, analogously to [70]. The advantage of the latter approach is that the number of permuted samples M may be much smaller than in the case of permutation test.

The main challenge to detect significant interactions among potential predictors $F=\{1,\dots ,p\}$ is usually computational burden of the procedure. Indeed, any method discussed should in principle include interaction term for any pair of predictors $X_i,X_j,i,j\in F$. This requires enlarging set F by $\left|F\right|\times \left(\right|F|-1)/2$ hot-encoded interactions, i.e., fitted model will contain $\left|F\right|\times \left(\right|F|+1)/2$ nominal predictors. This is practically infeasible, e.g., for 100 Single Nucleotide Polymorphisms (SNP) being our predictors this would yield number of needed tests of order ${10}^6$ at each step of selection procedure. The frequently adapted approach is to screen the set of variables first and then apply more sophisticated procedure to the chosen variables only.

Consider two other possible solutions. The first one uses the premise that significant interactions may arise only among variables which themselves are significant. Thus any of the method described in Section 8 can be used for predictors in F and then two methods described above are applicable for all interactions of the chosen features using Bonferroni adjustment for multiple testing or Holm method [83]. The problematic aspect of this approach is exclusion of possible interactions between two features when only one of them has non-negligible main effect.

The other group of methods, for which BOOST proposed in [80] is a representative example, screens interactions by the modified version of LRT test which consists of replacing likelihood for a fitted additive model by likelihood for normalised Kirkwood approximation ${\tilde{P}}_K/\eta$ which can be very quickly computed. The pairs for which the value of LRT statistic modified in this way exceeds certain threshold. Then, an exact LRT test is performed for all remaining pairs. The weak side of this approach is a choice of a threshold $\tau$ in the first stage which needs to be chosen by a rule-of-thumb as properties of modified LRT test remain unknown.

11. Conclusions

This paper reviews main ideas concerning feature selection using information theoretic tools which revolve around detecting and measuring the strength of conditional independence. As estimation of $CMI,$ which is a natural and powerful tool for this endeavour, encounters a problem when dimensionality of the task is large; a natural way is to look for its good substitutes. Several ways in which such substitutes are constructed are discussed and the current knowledge about their properties is presented. It is argued that selection criteria based on truncation of Möbius decomposition make quite far-reaching compromises on the dependence structures of feature sets, whereas properties of approximations based on variational approaches are yet to be established. Additionally, major Markov Blanket discovery algorithms are constructed under assumptions that conditional independence or its absence can be established without error, when, in reality, the intrinsic features of any practical CI test applied at each stage is its fallibility. Therefore, further efforts, both theoretical and experimental, are needed to understand the advantages, drawbacks, and interrelations between up-to-date developments. The paper presents some recently established tools which can be used for this purpose, such as asymptotic distributions of $CMI$-related criteria discussed in Section 6. They are used to construct tests of conditional independence with approximate control of probability of false signals. The problem of existence of interactions can be similarly approached using results in [69].

Another challenging task is to extend general approaches discussed here, such as construction of feature selection criteria by truncation of Möbius expansion or variational approach to multilabel classification when Y is a multivariate binary vector. Several criteria, such as $CIFE$ or $JMI$, have been generalised to this case already (cf. [84,85], respectively, see also [86] for an approach based on Han’s inequality and [87] for the review), but the general approach, e.g., in the spirit of [24], is still missing. Another important challenge here seems construction of feature selection criteria which will efficiently take into account dependence between coordinates of the response.