---
abstract: 'We present an efficient feature selection method that can find all multiplicative combinations of *continuous* features that are statistically significantly associated with the class variable, while rigorously correcting for multiple testing. The key to overcome the combinatorial explosion in the number of candidates is to derive a lower bound on the $p$-value for each feature combination, which enables us to massively prune combinations that can never be significant and gain more statistical power. While this problem has been addressed for binary features in the past, we here present the first solution for continuous features. In our experiments, our novel approach detects true feature combinations with higher precision and recall than competing methods that require a prior binarization of the data.'
author:
- |
    Mahito Sugiyama\
    *National Institute of Informatics; JST PRESTO*\
    `mahito@nii.ac.jp`
- |
    Karsten M. Borgwardt\
    *D-BSSE, ETH Z[ü]{}rich*\
    `karsten.borgwardt@bsse.ethz.ch`
title: '**Finding Significant Combinations of Continuous Features**'
---

=1

Introduction {#sec:intro}
============

A big challenge in high-dimensional data analysis is the search for features that are statistically significantly associated with the class variable, while accounting for the inherent multiple tesing problem. This problem is relevant in a broad range of applications including natural language processing, statistical genetics, and healthcare. To date, this problem of *feature selection*  [@Guyon03] has been extensively studied in machine learning, including the recent advances in *selective inference* [@Taylor15], a technique that can assess the statistical significance of features selected by linear models such as the Lasso [@Lee2016].

However, current approaches have a crucial limitation: They can only find *single features or linear combinations of features*, but it is still an open problem to find [*patterns*]{}, that is, *combinations of features with multiplicative effect*. A relevant line of research towards this goal is *significant pattern mining* [@Llinares2015KDD; @Papaxanthos16; @Terada13], which tries to find statistically associated feature combinations while controlling the *family-wise error rate* (*FWER*), that is, the probability to detect one or more false positive patterns.

However, all existing methods for significant pattern mining only apply to combinations of *binary or discrete* features, and none of methods can handle real-valued data, although such data is common in many applications. If we binarize data beforehand to use significant pattern mining approaches, a binarization-based method cannot distinguish correlated and uncorrelated features (see Figure \[figure:example\] for an example). *Subgroup discovery* [@Atzmueller15; @Herrera11; @Novak09] also has the same goal of finding associated feature combinations, but the existing methods are also designed for discrete data, which means that binarization is required [@Grosskreutz09] for real-valued data and the above problem still exists.

To date, there is no method that can find all combinations of *continuous* features that are significantly associated with an output variable and that accounts for the inherent multiple testing problem. To solve this problem, one has to address the following three challenges:

1.  \[prob:pval\] How to assess the significance for a combination of continuous features?

2.  \[prob:mcp\] How to perform multiple testing correction? In particular, how to control the FWER (family-wise error rate), the probability to detect one or more false positive combinations?

3.  \[prob:enum\] How to manage the combinatorial explosion of the candidate space, where the number of possible combinations is $2^n$ for $n$ features?

The second problem and the third problem are related to each other: If one can reduce the number of combinations by pruning unnecessary ones, one can gain statistical power and reduce false negative combinations while at the same time controlling the FWER.

![Although there is a clear correlation between Features 1 and 2 in the left panel and none in the right panel, after median-based binarization of Feature 1 and Feature 2, the estimated probability of the occurrence of the feature combination (that is the number of points for which Feature $1 =$ Feature $2 = 1$, or the number of points in the red box) will be exactly the same in both examples.[]{data-label="figure:example"}](example.pdf){width=".75\linewidth"}

> *Our goal in this paper is to present the first method that can find all combinations of continuous features that are statistically significantly associated with the class variable, while controlling the FWER.*

Our approach is to use the *rank order statistics* to directly estimate the probability of occurrence of each feature combination from continuous data, which is known as *copula support* [@Tatti13], and apply a *likelihood ratio test* to assess the significance of association between feature combinations and the class label, which solves problem \[prob:pval\]. We present a lower bound on $p$-values of association, which enables us to prune unnecessary combinations that can never be significant through the technique of *testability* proposed by @Tarone90, and can solve both problems \[prob:mcp\] and \[prob:enum\].

The remainder of this paper is organized as follows: We introduce our method in Section \[sec:method\]. We introduce a likelihood ratio test as a statistical association test for combinations in Section \[subsec:test\], analyze multiple testing correction using the testability in Section \[subsec:mcp\], and present an algorithm that enumerates all significant feature combinations in Section \[subsec:alg\]. We experimentally validate our method in Section \[sec:exp\] and summarize our findings in Section \[sec:conclusion\].

The Proposed Method {#sec:method}
===================

We consider a supervised dataset $D = \{\,({\boldsymbol{v}}_1, y_1), ({\boldsymbol{v}}_2, y_2), \dots, ({\boldsymbol{v}}_N, y_N)\,\}$, where each data point is a pair of an $n$-dimensional vector ${\boldsymbol{v}}_i = (v_{i}^1, v_{i}^2, \dots, v_{i}^n) \in {\mathbb{R}}^n$ and a binary class label $y_i \in \{0, 1\}$[^1]. We denote the set of features by $V = \{1, 2, \dots, n\}$ and the power set of features by $2^V$. For each feature $j \in \{1, 2, \dots, n\}$, we write ${\boldsymbol{v}}^j = (v_1^j, v_2^j, \dots, v_N^j)$, which is the $N$-dimensional vector composed of the $j$th feature of the dataset $D$. We use the standard notation $A {\perp\!\!\!\perp}B$ for the statistical independence between $A$ and $B$.

Our goal is to find *every feature combination that is associated with class labels*. More precisely, let $X_{J}$ be a binary variable induced by a combination $J \in 2^V$, where $X_{J} = 1$ means the occurrence of $J$ and $X_{J} = 0$ otherwise, and $Y$ be an output binary variable of which class labels are realizations. Our objective is to find *all* feature combinations $J \in 2^V$ such that the null hypothesis $X_J {\perp\!\!\!\perp}Y$ is rejected by a statistical association test while rigorously controlling the FWER, the probability of detecting one or more false positive associations, under a predetermined significance level $\alpha$.

To estimate the probability distribution of $X_J$ for each combination $J$, we use the prominent result given by @Tatti13 in the context of frequent pattern mining from continuous data. Given a feature combination $J \in 2^V$, first we convert a dataset $D$ with respect to features $J$ into a single $N$-dimensional vector ${\boldsymbol{x}}_{J} = (x_{J1}, x_{J2}, \dots, x_{JN}) \in [0, 1]^N$ such that $$\begin{aligned}
 \label{eq:rankvector}
 x_{Ji} = \prod_{j \in J} {\pi}(v_i^j),\end{aligned}$$ where ${\pi}({\boldsymbol{v}}^j) = ({\pi}(v_1^j), {\pi}(v_2^j), \dots, {\pi}(v_N^j))$ is a *normalized ranking* with respect to a feature $j$ and defined as ${\pi}(v_i^j) = (k - 1) / (N - 1)$ if $v_i^j$ is the $k$th smallest value among $N$ values $v_1^j, v_2^j, \dots, v_N^j$. @Tatti13 showed that the probability $\Pr(X_J = 1)$ is estimated by the *copula support* defined as $$\begin{aligned}
 \label{eq:copula}
 \eta(J) = \frac{1}{N} \sum_{i = 1}^N x_{Ji}.\end{aligned}$$ The copula support has the *monotonicity* with respect to the inclusion relationship of combinations: $$\begin{aligned}
 \eta(J) \ge \eta(J \cup \{j\}),\quad \forall J \in 2^V,\ j \in V \setminus J.\end{aligned}$$ Note that, although @Tatti13 considered a statistical test for the copula support, it cannot be used in our setting due to the following two reasons: (1) his setting is unsupervised while ours is supervised; (2) multiple testing correction was not considered for the test.

Statistical Testing {#subsec:test}
-------------------

We propose to use a *likelihood ratio test* [@Fisher22], which is a generalized $\chi^2$-test and often called G-test [@Woolf57] for testing the independence of each feature combinations. Although Fisher’s exact test has been used as the standard statistical test in recent studies (e.g. [@Llinares2015KDD; @SugiyamaSDM15; @Terada13]), it can be applied to only discrete test statistics and cannot be used in our setting.

[.49]{}

              $X_J = 1$         $X_J = 0$        Total
  --------- -------------- -------------------- -------
   $Y = 1$   $\eta(J)r_1$   $r_1 - \eta(J)r_1$   $r_1$
   $Y = 0$   $\eta(J)r_0$   $r_0 - \eta(J)r_0$   $r_0$
    Total     $\eta(J)$       $1 - \eta(J)$       $1$

  : Observed distribution.[]{data-label="table:contingency:observed"}

[.49]{}

              $X_J = 1$        $X_J = 0$       Total
  --------- ------------- ------------------- -------
   $Y = 1$   $\eta_1(J)$   $r_1 - \eta_1(J)$   $r_1$
   $Y = 0$   $\eta_0(J)$   $r_0 - \eta_0(J)$   $r_0$
    Total     $\eta(J)$      $1 - \eta(J)$      $1$

  : Observed distribution.[]{data-label="table:contingency:observed"}

Suppose that $\Pr(Y = l) = r_l$ for each class $l \in \{0, 1\}$. From two binary variables $X_J$ and $Y$, we obtain a $2 \times 2$ contingency table, where each cell denotes the joint probability $\Pr(X_j = l, Y = l')$ with $l, l' \in \{0, 1\}$ and can be described as a four-dimensional probability vector ${\boldsymbol{p}}$ such that $$\begin{aligned}
 {\boldsymbol{p}} = \big(\,\Pr(X_j = 1, Y = 1), \Pr(X_j = 1, Y = 0), \Pr(X_j = 0, Y = 1), \Pr(X_j = 0, Y = 0)\,\big).\end{aligned}$$ Let ${\boldsymbol{p}}_{\mathrm{E}}$ be the expected probability vector under the null hypothesis $X_J {\perp\!\!\!\perp}Y$ and ${\boldsymbol{p}}_{\mathrm{O}}$ be the vector obtained from $N$ observations. The difference between two distributions ${\boldsymbol{p}}_{\mathrm{O}}$ and ${\boldsymbol{p}}_{\mathrm{E}}$ can be measured by the Kullback–Leibler (KL) divergence ${D_{\mathrm{KL}}}({\boldsymbol{p}}_{\mathrm{O}}, {\boldsymbol{p}}_{\mathrm{E}}) = \sum_i p_{\mathrm{O}i} \log (p_{\mathrm{O}i} / p_{\mathrm{E}i})$, and the independence $X_J {\perp\!\!\!\perp}Y$ can be translated into the condition ${D_{\mathrm{KL}}}({\boldsymbol{p}}_{\mathrm{O}}, {\boldsymbol{p}}_{\mathrm{E}}) = 0$. Hence the null hypothesis $H_0$ and the alternative hypothesis $H_1$ are $$\begin{aligned}
 H_0: {D_{\mathrm{KL}}}({\boldsymbol{p}}_{\mathrm{O}}, {\boldsymbol{p}}_{\mathrm{E}}) = 0,\quad\quad
 H_1: {D_{\mathrm{KL}}}({\boldsymbol{p}}_{\mathrm{O}}, {\boldsymbol{p}}_{\mathrm{E}}) \not= 0.\end{aligned}$$ This is known as G-test, a special case of likelihood ratio test, and the test statistic $\lambda$ is given as $\lambda = 2N{D_{\mathrm{KL}}}({\boldsymbol{p}}_{\mathrm{O}}, {\boldsymbol{p}}_{\mathrm{E}})$, which follows the $\chi^2$-distribution with the degree of freedom $1$.

In our case, for each combination $J \in 2^V$, the expected probability vector ${\boldsymbol{p}}_{\mathrm{E}}$ is given as $$\begin{aligned}
 {\boldsymbol{p}}_{\mathrm{E}} = \Big(\,\eta(J)r_1,\, \eta(J)r_0,\, r_1 - \eta(J)r_1,\, r_0 - \eta(J)r_0\,\Big),\end{aligned}$$ where $\eta(J)$ is the copula support and $r_l = (1/N) \sum_{i = 1}^N \mathbf{1}_{y_i = l}$, the estimated probability for each class $l \in \{0, 1\}$. In contrast, the observed probability vector ${\boldsymbol{p}}_{\mathrm{O}}$ is given as $$\begin{aligned}
 {\boldsymbol{p}}_{\mathrm{O}} = \Big(\,\eta_1(J),\, \eta_0(J),\, r_1 - \eta_1(J),\, r_0 - \eta_0(J)\,\Big),\end{aligned}$$ where $\eta_1(J) = {\boldsymbol{x}}_J{\boldsymbol{y}}$ and $\eta_1(J) = {\boldsymbol{x}}_J(1 - {\boldsymbol{y}})$ with ${\boldsymbol{x}}_J$ defined in Equation . These two distributions are shown in Table \[table:contingency\]. While is it known that the statistic $\lambda$ does not exactly follow the $\chi^2$-distribution if one of components of ${\boldsymbol{p}}_{\mathrm{E}}$ is too small, the situation does not usually occur since $\eta(J)$ is not too large as $\eta(J) \le 0.5$ by definition and not too small as significant combinations are always included in the largest $k$ combinations with respect to $\eta$ for some $k$, which will be shown in Section \[subsec:alg\].

In the following, we write the set of all possible probability vectors by ${\boldsymbol{{\mathcal{P}}}} = \{\,{\boldsymbol{p}} \mid p_i \ge 0, \sum p_i = 1\,\}$, and its subset given marginals $a$ and $b$ as ${\boldsymbol{{\mathcal{P}}}}(a, b) = \{\,{\boldsymbol{p}} \in {\boldsymbol{{\mathcal{P}}}} \mid p_1 + p_2 = a,\, p_1 + p_3 = b\,\}$.

Multiple testing correction {#subsec:mcp}
---------------------------

In statistical testing, we need to perform *multiple testing correction* to control the FWER (family-wise error rate), otherwise we will find massive false positive combinations as we have $2^n$ hypotheses, all combinations of features. The most popular multiple testing correction method is *Bonferroni correction* [@Bonferroni36]. In the method, the predetermined significance level $\alpha$ (e.g. $\alpha = 0.05$) is corrected as $\delta_{{\mathrm{Bon}}} = \alpha / m$, where $m$ is called a *correction factor* and is the number of hypotheses $m = 2^n$ in our case, and each hypothesis $J$ is declared as significant only if ${\mathrm{pval}}(J) < \delta_{{\mathrm{Bon}}} = \alpha / 2^n$. Then the resulting FWER $< \alpha$ is guaranteed. However, it is well known that the method is too conservative. In particular in our case, the correction factor $2^n$ is too massive due to the combinatorial effect and it is almost impossible to find significant combinations, which will generate massive *false negatives*.

Here we use the *testability* of hypotheses introduced by @Tarone90 and widely used in significant pattern mining [@Llinares2015KDD; @Papaxanthos16; @SugiyamaSDM15; @Terada13], which allows as to prune unnecessary feature combinations that can never be significant while controlling the FWER at the same level. For each feature combination $J \in 2^V$, we use the *lower bound of the $p$-value* of association between $X_J$ and $Y$ given marginals $\eta(J)$ and $r_1$, denoted by $\psi(J)$, which is defined as the $p$-value for the upper bound of the KL divergence in our statistical test. Suppose that $J_1, J_2, \dots, J_{2^n}$ be the sequence of all combinations such that the relationship $$\begin{aligned}
 \psi(J_1) \le \psi(J_2) \le \psi(J_3) \le \dots \le \psi(J_{2^n})\end{aligned}$$ is satisfied. Let $m$ be the threshold such that $$\begin{aligned}
 \label{eq:testability}
 m \cdot \psi(J_m) < \alpha \quad\text{and}\quad (m + 1) \cdot \psi(J_{m + 1}) \ge \alpha.\end{aligned}$$ @Tarone90 showed that the FWER is controlled under $\alpha$ with the correction factor $m$. The set ${\mathcal{T}}= \{J_1, J_2, \dots, J_m\}$ is called *testable* combinations, and each $J \in {\mathcal{T}}$ is significant if ${\mathrm{pval}}(J) < \delta_{{\mathrm{Tar}}} = \alpha / m$. On the other hand, combinations $J_{m + 1}, \dots, J_{2^n}$ are called *untestable* as they can never be significant. Since $m \ll 2^n$ usually holds, we can expect to obtain higher statistical power in Tarone’s method compared to Bonferroni method. Moreover, if we can efficiently enumerate such testable combinations without seeing untestable ones, we can overcome combinatorial explosion.

Our key theoretical contribution is that the tight upper bound of the KL divergence to compute $\psi(J)$ is analytically obtained in the following.

\[theorem:KLmax\] For $a, b \in [0, 1]$ with $a \le b \le 1/2$ and a probability vector ${\boldsymbol{p}}_{\mathrm{E}} = (ab, a(1 - b), (1 - a)b, (1 - a)(1 - b))$, we have $$\begin{aligned}
  \label{eq:KLmax}
  {D_{\mathrm{KL}}}({\boldsymbol{p}}, {\boldsymbol{p}}_{\mathrm{E}}) < a\log\frac{1}{b} + (b - a)\log\frac{b - a}{(1 - a)b} + (1 - b)\log\frac{1}{(1 - a)}
 \end{aligned}$$ for all ${\boldsymbol{p}} \in {\boldsymbol{{\mathcal{P}}}}(a, b)$ and this is tight.

Let $f(x) = {D_{\mathrm{KL}}}({\boldsymbol{p}}, {\boldsymbol{p}}_{\mathrm{E}})$ with ${\boldsymbol{p}} = (a, a - x, b - x, (1 - b) - (a - x))$ and $a' = 1 - a$, $b' = 1 - b$. We have $$\begin{aligned}
 f(x) = x\log\frac{x}{ab} + (a - x)\log\frac{a - x}{ab'} + (b - x)\log\frac{b - x}{a'b} + (b' - a + x)\log\frac{b' - a + x}{a'b'}.\end{aligned}$$ The second derivative of $f(x)$ is given as $$\begin{aligned}
 \frac{\partial^2 f}{\partial x^2} = \frac{1}{b' - x} + \frac{1}{b - a + 1} + \frac{1}{a - x} + \frac{1}{x},\end{aligned}$$ which is always positive from the constraint $0 < x < \min\{a, b\} = a$. Thus $f(x)$ is maximized when $x$ goes to $0$ or $a$. The limit of $f(x)$ is obtained as $$\begin{aligned}
 \lim_{x \to c} f(x) =
 \left\{
 \begin{array}{*2{>{\displaystyle}l}}
  a\log\frac{1}{b'} + b\log\frac{1}{a'} + (b' - a)\log\frac{b' - a}{a'b'} & \!\!\text{if } c = 0, \vspace*{5pt}\\
  \eta\log\frac{1}{b} + (b - a)\log\frac{b - a}{a'b} + (1 - b)\log\frac{1}{a'} & \!\!\text{if } c = a.
 \end{array}
 \right.\end{aligned}$$ The difference $\delta = \lim_{x \to a} f(x) - \lim_{x \to 0} f(x)$ is obtained as $$\begin{aligned}
 \delta = - b\log b + b'\log b' + (b - a)\log(b - a) - (b' - a)\log(b' - a).\end{aligned}$$ The partial derivative of $\delta$ with respect to $\eta$ is $$\begin{aligned}
\frac{\partial \delta}{\partial \eta} = \log(1 - a - b) - \log(b - a),\end{aligned}$$ which is always positive as $(1 - a - b) - (b - a) = 1 - 2b \ge 0$ with the condition $b \le 1/2$. Hence the difference $\delta$ takes the minimum value at $a = 0$ and we obtain $\delta \ge - b\log b + b'\log b' + b'\log b' - b'\log b' = 0$. Therefore $f(x)$ is the tight upper bound when $x \to a$.

![The upper bound $B(a, b)$ of the KL divergence (left-hand side) and the corresponding *p*-value (right-hand side) with $N = 100$ with respect to changes in $a$ when $b = 0.3$.[]{data-label="figure:klmax_pmin"}](klmax_pmin.pdf){width=".7\linewidth"}

Let $B(a, b)$ be the bound shown in Equation . We analyze the behavior of the bound $B(a, b)$ as a function of $a$ with fixed $b$. Assume that $b \le 1/2$. When $a < b$, we have $$\begin{aligned}
 \frac{\partial B(a, b)}{\partial a} = \log\frac{1}{b} - \log\frac{b - a}{b(1 - a)} = \log\frac{1 - a}{b - a} > 0,\end{aligned}$$ hence it is monotonically increases as $a$ increases. When $b < a < 1/2$, we have $$\begin{aligned}
 \frac{\partial B(a, b)}{\partial a} = \log\frac{a - b}{a(1 - b)} - \log\frac{1}{1 - b} = \log\frac{a - b}{a} < 0,\end{aligned}$$ thereby it monotonically decreases as $a$ increases. We illustrate the bound $B(a, b)$ with $b = 0.3$ and the corresponding minimum achievable $p$-value with the sample size $N = 100$ in Figure \[figure:klmax\_pmin\].

Algorithm {#subsec:alg}
---------

We present a depth-first search algorithm that finds $J_1, J_2, \dots$, $J_m$ such that $\psi(J_1) \le \psi(J_2) \le \dots \le \psi(J_m)$, where $m$ satisfies the condition . Let $r_1 \le 1/2 \le r_0$. Since the lower bound of the $p$-value $\psi(J)$ takes the minimum value when $\eta(J) = r_1$ (see Figure \[figure:klmax\_pmin\]) and is monotonically decreasing as $\eta(J)$ decreases, for any $\eta(J) \le r_1$, $$\begin{aligned}
 \label{eq:testable_cand}
 {\mathcal{C}}= \Set{J \in 2^V \mid \eta(J) \ge \sigma} \supseteq {\mathcal{T}}= \Set{J \in 2^V | \psi(J) \le \psi_{\sigma}}\end{aligned}$$ is always guaranteed, where $\sigma$ is some threshold and $\psi_{\sigma}$ is the lower bound of the $p$-value for $\sigma$. Thus if the condition  is satisfied for some $m \le |{\mathcal{T}}|$, the $m$ smallest combinations in ${\mathcal{T}}$ are the testable feature combinations.

Moreover, since $\eta(J)$ has the monotonicity with respect to the inclusion relationship, that is, $\eta(J) \ge \eta(J \cup \{j\})$ for all $j \in V \setminus J$, finding the set ${\mathcal{C}}$ is simply achieved by DFS (depth-first search): Starting from the smallest combination $\emptyset$ with assuming $\eta(\emptyset) = 1$, if $\eta(J) \ge \sigma$ for a combination $J$, we update $J$ by adding $j \in V \setminus J$ and recursively check $J$. Otherwise if $\eta(J) < \sigma$, we prune all combinations $K \supseteq J$ as $\eta(K) < \sigma$ holds.

Our algorithm dynamically updates $\sigma$ to enumerate testable combinations while pruning massive unnecessary combinations. We summarize our search strategy in Algorithm \[alg:enum\]. First we set $\sigma = 0$, which means that all combinations are testable. Whenever we find a new combination $J$, we update ${\mathcal{C}}$ and ${\mathcal{T}}$ and check the condition . If $|{\mathcal{T}}|\psi_{\sigma} > \alpha$ holds, $\sigma$ is too low and ${\mathcal{T}}$ will become too large with the current $\sigma$, hence we update $\sigma$ to $\min_{J \in {\mathcal{C}}} \eta(J)$ and remove the combination ${\mathop{\mathrm{argmin}}}_{J \in {\mathcal{C}}} \eta(J)$. Finally, when the algorithm stops, it is clear that the set ${\mathcal{T}}$ coincides with the set of testable combinations, hence each combination $J \in {\mathcal{T}}$ is significant if the $p$-value of $J$ is smaller than $\delta_{{\mathrm{Tar}}} = \alpha / |{\mathcal{T}}|$.

Experiments {#sec:exp}
===========

We examine the efficiency and the effectiveness of our method using synthetic and real-world datasets. We used Amazon Linux AMI release 2016.09 and ran all experiments on a single core of 2.3 GHz Intel Xeon CPU E7-8880 v3 and 2.0 TB of memory. All methods were implemented in `C`/`C++` and compiled with [gcc]{} 4.8.3.

There exists no method that can enumerate significant feature combinations from continuous data with multiple testing correction. Thus we compare our method to *significant itemset mining* method with prior *binarization* of a given dataset since significant itemset mining offers to enumerate all significant feature combinations while controlling the FWER from binary data. We employ median-based binarization as a preprocessing. For each feature $j \in V$, we pick up the median of ${\boldsymbol{v}}^j = (v_1^j, v_2^j, \dots, v_N^j)$, denoted by ${\mathop{\mathrm{med}}}(j)$, and binarize each value $v_i^j$ as a pair $(v_i^{\le {\mathop{\mathrm{med}}}(j)}, v_i^{> {\mathop{\mathrm{med}}}(j)})$, where $v_i^{\le {\mathop{\mathrm{med}}}(j)} = \mathtt{1}$ if $v_i^{j} \le {\mathop{\mathrm{med}}}(j)$ and $\mathtt{0}$ otherwise, and $v_i^{> {\mathop{\mathrm{med}}}(j)} = \mathtt{1}$ if $v_i^{j} > {\mathop{\mathrm{med}}}(j)$ and $\mathtt{0}$ otherwise. Thus a given dataset is converted to the binarized dataset with $2n$ features. We use the state-of-the-art significant itemset mining algorithm LAMP ver.2[^2] [@Minato14], which incorporates the fastest frequent pattern mining algorithm LCM [@Uno04] to enumerate testable feature combinations from binary data.

We also tried other binarization approaches, *interordinal scaling* used in numerical pattern mining @Kaytoue11 and *interval binarization* in subgroup discovery [@Grosskreutz09] as a preprocessing of significant itemset mining. However, these methods generate too many binarized dense features, resulting in the lack of scalability in the itemset mining step. They can finish their computation for a tiny dataset with $N = 50$ and $n = 5$ in approximately 24 hours, but did not finish after 48 hours for $N = 100$ and $n = 10$, and they exceed the memory limit (2.0 TB) for larger datasets.

![Results on synthetic data with the minor class ratio $r_1 = 0.5$ (**a**) and $r_1 = 0.1$ (**b**). In each figure, the number of features is $n = 20$ in the left column and the sample size is $N = 200$ in the right column. Missing values in the plots of precision and recall mean that no significant combination is detected. Both x- and y-axes are in logarithmic scale. Our proposal is shown in blue circles, the binarization approach in green crosses.[]{data-label="figure:synth"}](synth_fig.pdf){width="\linewidth"}

To examine the quality of detected combinations, we compute precision and recall. Let $A \subseteq V$ be the features that are truely associated with the class labels (ground truth). In the binarization-based method, $A$ is the set of binarized features such that each of them is generated from each feature in the original $A$. Suppose that ${\mathcal{S}}$ is the set of significant combinations. For each combination $J \in {\mathcal{S}}$, $J$ is true positive if $J \subseteq A$. Thus precision and recall are given as $\mathrm{TP} / |{\mathcal{S}}|$ and $\mathrm{TP} / 2^{|A|}$, respectively, where $\mathrm{TP} = |\{J \in {\mathcal{S}}\mid J \subseteq A\}|$.

To examine the efficiency of methods, we also measure the running time needed for enumerating all significant combinations. In the binarization-based method, we exclude the time used for binarization as the process is efficient enough and negligible compared to the pattern mining step.

Results on Synthetic Data
-------------------------

To examine the evaluation criteria, we generate synthetic data, varying the sample size $N$ from $100$ to $100,000$, the number $n$ of features from $10$ to $200$, and setting the class ratio to $r_1 = 0.5$ or $r_1 = 0.1$, i.e., the number $N_1$ of samples in the minor class is $N/2$ or $N/10$. In each dataset, we generate 20% of features that are associated with the class labels. More precisely, first we generate the entire dataset from the uniform distribution from $0$ to $1$ and assign the class label $1$ to the first $N_1$ data point. Then each $j$ of the 20% of associated features, we sort ${\boldsymbol{v}}^j$ in decreasing order and add Gaussian noise with $(\mu, \sigma^2) = (0, 0.2)$. These 20% of correlated features are the ground truth and used to compute the true positives. The other 80% of features are uninformative features.

Results are plotted in Figure \[figure:synth\](**a**) for $r_1 = 0.5$ and Figure \[figure:synth\](**b**) for $r_1 = 0.1$. In the respective figure, we plot results with varying $N$ while fixing $n = 20$ on the left column and those with varying $n$ while fixing $N = 200$ on the right column. In comparison with the median-based binarization method (plotted in green), our method has a clear advantage regarding the quality of detected significant combinations. In all cases, both precision and recall of our method are higher than those of the median-based binarization method. This is why binarization-based method cannot distinguish correlated and uncorrelated combinations as we discussed in the Introduction and illustrated in Figure \[figure:example\].

  ------ ------ ----- ---------------------------------------- ------------------------------------ ----------------------------------------- ------------------------------------- ---------------------------------------- ------------------------------------
    data    $N$   $n$                                                                                                                                                                                                        
                      Ours                                     Binarization                         Ours                                      Binarization                          Ours                                     Binarization
      D1    942    44 [**2.08 $\times$ 10$^{-\text{1}}$**]{}   [5.23 $\times$ 10$^{-\text{2}}$]{}   [**1.85 $\times$ 10$^{-\text{3}}$**]{}    [2.33 $\times$ 10$^{-\text{8}}$]{}    [**9.59 $\times$ 10$^{+\text{1}}$**]{}   [4.54 $\times$ 10$^{+\text{3}}$]{}
      D2    316    54 [**8.89 $\times$ 10$^{-\text{1}}$**]{}   [4.48 $\times$ 10$^{-\text{1}}$]{}   [**1.79 $\times$ 10$^{-\text{7}}$**]{}    [1.27 $\times$ 10$^{-\text{13}}$]{}   [**8.13 $\times$ 10$^{+\text{0}}$**]{}   [1.42 $\times$ 10$^{+\text{1}}$]{}
      D3    252    68 [**9.17 $\times$ 10$^{-\text{1}}$**]{}   [8.93 $\times$ 10$^{-\text{1}}$]{}   [**6.40 $\times$ 10$^{-\text{10}}$**]{}   [1.73 $\times$ 10$^{-\text{17}}$]{}   [**9.14 $\times$ 10$^{+\text{0}}$**]{}   [2.57 $\times$ 10$^{+\text{1}}$]{}
      D4    660    38 [**1.56 $\times$ 10$^{-\text{1}}$**]{}   [8.57 $\times$ 10$^{-\text{2}}$]{}   [**1.06 $\times$ 10$^{-\text{3}}$**]{}    [1.27 $\times$ 10$^{-\text{7}}$]{}    [**2.11 $\times$ 10$^{+\text{1}}$**]{}   [5.49 $\times$ 10$^{+\text{1}}$]{}
      D5   3314    42 [**1.01 $\times$ 10$^{-\text{1}}$**]{}   —                                    [**1.64 $\times$ 10$^{-\text{3}}$**]{}    —                                     [**5.52 $\times$ 10$^{+\text{3}}$**]{}   $\ge$10h
      D6    424    60 [**8.00 $\times$ 10$^{-\text{1}}$**]{}   [1.55 $\times$ 10$^{-\text{1}}$]{}   [**3.27 $\times$ 10$^{-\text{5}}$**]{}    [1.82 $\times$ 10$^{-\text{12}}$]{}   [**6.20 $\times$ 10$^{+\text{1}}$**]{}   [2.06 $\times$ 10$^{+\text{2}}$]{}
  ------ ------ ----- ---------------------------------------- ------------------------------------ ----------------------------------------- ------------------------------------- ---------------------------------------- ------------------------------------

  : Results on real data. Datasets from D1 to D6 means `ctg`, `faults`, `ionosphere`, `segment`, `waveform`, and `wdbc`, respectively. Best scores are denoted in Bold.[]{data-label="table:real_results"}

In both methods, precision drops when the sample size $N$ becomes large: As we gain more and more statistical power for larger $N$, many feature combinations, even those with very small dependence to the class labels, tend to reach statistical significance. In contrast, recall drops when the number $n$ of features becomes large, due to the growing number of testable combinations, resulting in a loss of statistical power in finding significant combinations.

Our approach is also constantly faster than the binarization approach on all datasets, as can be seen in Figure 3 (bottom row). To summarize, we observe that our proposal improves over the competing binarization approach in terms of runtime and precision and recall in detecting truly significant feature combinations.

Results on Real Data
--------------------

We also evaluate our proposal on real-world datasets shown in Table \[table:real\_results\], which are benchmark datasets for binary classification from the UCI repository[^3]. We treat all $n$ features as truly associated features, and add $n$ fake features drawn from a uniform distribution from $0$ to $1$. All datasets are balanced to maximize the statistical power for comparing detected significant combinations, that is, $r_1 = r_0 = 0.5$. If they are not balanced in the original dataset, we randomly subsample data points from the larger class. We summarize results in Table \[table:real\_results\], and we observe the same trend as before. In all dataset, our method shows superior performance: both precision and recall are better and running time is faster than for the binarization-based method. These results confirm the effectiveness of our proposal.

Conclusion {#sec:conclusion}
==========

In this paper, we have addressed the challenge of finding all (multiplicative) combinations of *continuous* features that are significantly associated with an output variable after rigorously controlling for multiple testing. While this problem has been studied for binary features, our approach is the first solution to this problem on continuous features. Key to our approach is the derivation of a lower bound on the $p$-value of association for each feature combination, which can be used to prune unnecessary combinations that can never reach statistical significance. Our experimental results demonstrate that our approach is superior in detecting true feature combinations compared to existing pattern mining approaches that require a prior binarization of the features.

Our work opens the door to many applications of searching significant feature combinations, in which the data is not adequately described by binary features, including large fields such as data analysis for high-throughput technologies in biology and medicine. Our work here addresses the problem of finding continuous *features*. Finding significant combinations associated with *continuous output variables* is an equally challenging and practically relevant problem, that we will tackle in future work.

[19]{} \[1\][\#1]{} \[1\][`#1`]{} urlstyle \[1\][doi: \#1]{}

M. Atzmueller. Subgroup discovery. *Wiley Interdisciplinary Reviews: Data Mining and Knowledge Discovery*, 50 (1):0 35–49, 2015.

C. E. Bonferroni. Teoria statistica delle classi e calcolo delle probabilit[à]{}. *Pubblicazioni del R Istituto Superiore di Scienze Economiche e Commerciali di Firenze*, 8:0 3–62, 1936.

R. A. Fisher. On the mathematical foundations of theoretical statistics. *Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences*, 2220 (594-604):0 309–368, 1922.

H. Grosskreutz and S. R[ü]{}ping. On subgroup discovery in numerical domains. *Data Mining and Knowledge Discovery*, 190 (2):0 210–226, 2009.

I. Guyon and A. Elisseeff. An introduction to variable and feature selection. *Journal of Machine Learning Research*, 3:0 1157–1182, 2003.

F. Herrera, C. J. Carmona, P. Gonz[á]{}lez, and M. J. del Jesus. An overview on subgroup discovery: foundations and applications. *Knowledge and Information Systems*, 290 (3):0 495–525, 2011.

M. Kaytoue, S. O. Kuznetsov, and A. Napoli. Revisiting numerical pattern mining with formal concept analysis. In *Proceedings of the 22nd International Joint Conference on Artificial Intelligence*, pages 1342–1347, 2011.

J. D. Lee, D. L. Sun, Y. Sun, and J. E. Taylor. Exact post-selection inference, with application to the lasso. *Annals of Statistics*, 440 (3):0 907–927, 06 2016.

F. Llinares-L[ó]{}pez, M. Sugiyama, L. Papaxanthos, and K. M. Borgwardt. Fast and memory-efficient significant pattern mining via permutation testing. In *Proceedings of the 21st ACM SIGKDD Conference on Knowledge Discovery and Data Mining*, pages 725–734, 2015.

S. Minato, T. Uno, K. Tsuda, A. Terada, and J. Sese. A fast method of statistical assessment for combinatorial hypotheses based on frequent itemset enumeration. In *Machine Learning and Knowledge Discovery in Databases*, volume 8725 of *LNCS*, pages 422–436. Springer, 2014.

P. K. Novak, N. Lavra[č]{}, and G. I. Webb. Supervised descriptive rule discovery: [A]{} unifying survey of contrast set, emerging pattern and subgroup mining. *The Journal of Machine Learning Research*, 10:0 377–403, 2009.

L. Papaxanthos, F. Llinares-Lopez, D. Bodenham, and K. M. Borgwardt. Finding significant combinations of features in the presence of categorical covariates. In *Advances in Neural Information Processing Systems*, volume 29, pages 2271–2279, 2016.

M. Sugiyama, F. Llinares-L[ó]{}pez, N. Kasenburg, and K. M Borgwardt. Significant subgraph mining with multiple testing correction. In *Proceedings of 2015 SIAM International Conference on Data Mining*, pages 37–45, 2015.

R. E. Tarone. A modified [B]{}onferroni method for discrete data. *Biometrics*, 460 (2):0 515–522, 1990.

N. Tatti. Itemsets for real-valued datasets. In *2013 IEEE 13th International Conference on Data Mining*, pages 717–726, 2013.

J. Taylor and R. J. Tibshirani. Statistical learning and selective inference. *Proceedings of the National Academy of Sciences*, 1120 (25):0 7629–7634, 2015.

A. Terada, M. Okada-Hatakeyama, K. Tsuda, and J. Sese. Statistical significance of combinatorial regulations. *Proc. Natl. Acad. Sci. USA*, 1100 (32):0 12996–13001, 2013.

T. Uno, T. Asai, Y. Uchida, and H. Arimura. An efficient algorithm for enumerating closed patterns in transaction databases. In *Discovery Science*, volume 3245 of *Lecture Notes in Computer Science*, pages 16–31, 2004.

B. Woolf. The log likelihood ratio test (the *G*-test). *Annals of human genetics*, 210 (4):0 397–409, 1957.

[^1]: This setting can be directly generalized to more than two classes when considering one-versus-rest classification.

[^2]: We used the implementation available at  <https://www.bsse.ethz.ch/mlcb/research/machine-learning/wylight.html>

[^3]: <http://archive.ics.uci.edu/ml/>
