Finding Statistically Significant Interactions between Continuous Features

TL;DR

This paper introduces a novel algorithm for efficiently identifying statistically significant interactions among continuous features, overcoming computational challenges and improving detection power in high-dimensional data.

Contribution

It presents the first method to detect significant interactions among continuous features by deriving p-value bounds for effective pruning, enhancing computational efficiency and statistical power.

Findings

Successfully detects all significant interactions in synthetic datasets

Efficiently handles high-dimensional continuous feature interactions

Outperforms existing methods in computational speed and accuracy

Abstract

The search for higher-order feature interactions that are statistically significantly associated with a class variable is of high relevance in fields such as Genetics or Healthcare, but the combinatorial explosion of the candidate space makes this problem extremely challenging in terms of computational efficiency and proper correction for multiple testing. While recent progress has been made regarding this challenge for binary features, we here present the first solution for continuous features. We propose an algorithm which overcomes the combinatorial explosion of the search space of higher-order interactions by deriving a lower bound on the p-value for each interaction, which enables us to massively prune interactions that can never reach significance and to thereby gain more statistical power. In our experiments, our approach efficiently detects all significant interactions in a variety of synthetic and real-world datasets.

Tables

Table 1. Table 1: Contingency tables.

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

Table 2. (a) Expected distribution.

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

Table 3. (b) Observed distribution.

	$X_J=1$	$X_J=0$	Total
$Y=1$	${\eta }_1\left(J\right)$	$r_1-{\eta }_1\left(J\right)$	$r_1$
$Y=0$	${\eta }_0\left(J\right)$	$r_0-{\eta }_0\left(J\right)$	$r_0$
Total	$\eta \left(J\right)$	$1-\eta \left(J\right)$	$1$

Table 4. Table 2: Statistics of real data.

Data	$N$	$d$	# candidate combinations
			(search space)
ctg	942	22	4,194,304
faults	316	27	134,217,728
ijcnn	9,706	22	4,194,304
magic	13,376	10	1,024
segment	660	19	524,288
transfusion	356	4	16
waveform	3,314	21	2,097,152
wdbc	424	30	1,073,741,824
wine	3,198	11	2,048
yacht	308	6	64

Equations

$$x_{Ji}={\prod }_{j\in J}\pi (v_i^j)$$

$$\eta (J)=\frac{1}{N}\sum _{i=1}^N{x}_{Ji}=\frac{1}{N}\sum _{i=1}^N\prod _{j\in J}\pi (v_i^j).$$

$$\displaystyle\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{)}.$$

$$\displaystyle\boldsymbol{p}_{\mathrm{E}}=\bigl{(}\,\eta(J)r_{1},\,\eta(J)r_{0},\,r_{1}-\eta(J)r_{1},\,r_{0}-\eta(J)r_{0}\,\bigr{)},$$

$$\displaystyle\boldsymbol{p}_{\mathrm{O}}=\bigl{(}\,\eta_{1}(J),\,\eta_{0}(J),\,r_{1}-\eta_{1}(J),\,r_{0}-\eta_{0}(J)\,\bigr{)},$$

$$D_{\mathrm{KL}}(\mathbf{p},{\mathbf{p}}_{\mathrm{E}})<a\log \frac{1}{b}$$

$$\displaystyle f(x)=\$$

$$\frac{{\partial }^{2}f}{\partial x^2}=\frac{1}{1-a-b+x}+\frac{1}{b-x}+\frac{1}{a-x}+\frac{1}{x},$$

$$\underset{x\to 0}{\lim }f(x)$$

$$\underset{x\to a}{\lim }f(x)$$

$$\delta \ge -b\log b+b^{\prime }\log b^{\prime }+b^{\prime }\log b^{\prime }-b^{\prime }\log b^{\prime }=0.$$

$$\psi (J_1)\le \psi (J_2)\le \psi (J_3)\le \cdots \le \psi (J_{2^d})$$

$$m\cdot \psi (J_m)<\alpha \text{and}(m+1)\cdot \psi (J_{m+1})\ge \alpha .$$

$$\displaystyle\Bigl{(}\,v_{i}^{\leq v_{1}^{j}},v_{i}^{\leq v_{2}^{j}},\dots,v_{i}^{\leq v_{N}^{j}},v_{i}^{\geq v_{1}^{j}},v_{i}^{\geq v_{2}^{j}},\dots,v_{i}^{\geq v_{N}^{j}}\,\Bigr{)},$$

$$\displaystyle\Bigl{(}\,v_{i}^{(v_{1}^{j},v_{2}^{j}]},v_{i}^{(v_{1}^{j},v_{3}^{j}]},\dots,v_{i}^{(v_{1}^{j},v_{N}^{j}]},v_{i}^{(v_{2}^{j},v_{3}^{j}]},\dots,v_{i}^{(v_{N-1}^{j},v_{N}^{j}]}\,\Bigr{)}.$$

Full text

1 Introduction

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 testing 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 Guyon and Elisseeff (2003) has been extensively studied in statistics and machine learning, including the recent advances in selective inference Taylor and Tibshirani (2015), a technique that can assess the statistical significance of features selected by linear models such as the Lasso Lee et al. (2016).

However, current approaches have a crucial limitation: They can find only single features or linear interactions of features, but it is still an open problem to find patterns, that is, multiplicative (potentially higher-order) interactions between features. A relevant line of research towards this goal is significant (discriminative) pattern mining Terada et al. (2013a); Llinares-López et al. (2015); Papaxanthos et al. (2016); Pellegrina and Vandin (2018), which tries to find statistically associated feature interactions while controlling the family-wise error rate (FWER), 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 existing significant pattern mining approaches, a binarization-based method may not be able to distinguish (un)correlated features (see Figure 1).

To date, there is no method that can find all higher-order interactions 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:

How to assess the significance for a multiplicative interaction of continuous features? 2. 2.

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 positives? 3. 3.

How to manage the combinatorial explosion of the candidate space, where the number of possible interactions is $2^{d}$ for $d$ features?

The second problem and the third problem are related with 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 an extensive body of work in statistics on multiple testing correction including the FDR Hochberg (1988); Benjamini and Hochberg (1995) and also several studies on approaches to interaction detection Bogdan et al. (2015); Su and Candès (2016), none of these approaches addresses the problem of dealing with the combinatorial explosion of the $2^{d}$-dimensional search space when trying to finding higher-order significant multiplicative interactions.

Our goal in this paper is to present the first method, called C-Tarone, that can find all higher-order interactions between 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 joint occurrence of each feature combination from continuous data, which is known as copula support Tatti (2013), and apply a likelihood ratio test to assess the significance of association between feature interactions and the class label, which solves the problem 1. We present the tight lower bound on $p$-values of association, which enables us to prune unnecessary interactions that can never be significant through the notion of testability proposed by Tarone (1990), and can solve both problems 2 and 3.

This paper is organized as follows: We introduce our method C-Tarone in Section 2. We introduce a likelihood ratio test as a statistical association test for interactions in Section 2.1, analyze multiple testing correction using the testability in Section 2.2, and present an algorithm in Section 2.3. We experimentally validate our method in Section 3 and summarize our findings in Section 4.

2 The Proposed Method: C-Tarone

Given 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 $d$-dimensional vector $\boldsymbol{v}_{i}=(v_{i}^{1},v_{i}^{2},\dots,v_{i}^{d})\in\mathbb{R}^{d}$ and a binary class label $y_{i}\in\{0,1\}$. We denote the set of features by $V=\{1,2,\dots,d\}$ and the power set of features by $2^{V}$. For each feature $j\in\{1,2,\dots,d\}$, 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$.

Our goal is to find every multiplicative feature interaction that is associated with class labels. To tackle the problem, first we measure the joint occurrence probability $\eta(J)\in\mathbb{R}$ of a feature combination $J\in 2^{V}$ in a dataset. For each combination $J\in 2^{V}$, the size $|J|$ corresponds to the order of an interaction, and we find arbitrary-order interactions in the form of combinations. If data is not real-valued but binary, that is, $v_{i}^{j}\in\{0,1\}$, this problem is easily solved by measuring the support used in frequent itemset mining Agrawal and Srikant (1994); Aggarwal and Han (2014) . Each feature combination $J\in 2^{V}$ is called an itemset, and the support of $J$ is defined as $\eta(J)=(1/N)\sum_{i=1}^{N}\prod_{j\in J}v_{i}^{j}$, which corresponds to joint probability of $J$. The support is always from [math] to $1$, and we can introduce a binary random variable $X_{J}$ corresponding to the joint occurrence of $J$, where $X_{J}=1$ if features in $J$ jointly occurs and $X_{J}=0$ otherwise, $\eta(J)$ corresponds to the empirical estimate of the probability $\Pr(X_{J}=1)$.

This approach can be generalized to continuous (real-valued) data using the copula support, which is the prominent result given by Tatti (2013) in the context of frequent pattern mining from continuous data. The copula support allows us to define the binary random variable $X_{J}$ of joint occurrence of $J$ and estimate $\eta(J)=\Pr(X_{J}=1)$ from continuous data. The key idea is to use the normalized ranking for each feature $j\in V$, which is 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}$. Then we convert a dataset $D$ with respect to features in a combination $J\in 2^{V}$ into a single $N$-dimensional vector $\boldsymbol{x}_{J}=(x_{J1},x_{J2},\dots,x_{JN})\in[0,1]^{N}$ by

[TABLE]

for each data point $i\in\{1,2,\dots,N\}$. Tatti (2013) showed that the empirical estimate of the probability $\Pr(X_{J}=1)$ of the joint occurrence of a combination $J$ is obtained by the copula support defined as

[TABLE]

Intuitively, joint occurrence of $J$ means that rankings among features in $J$ match with each other. The definition in Equation (2) is analogue to the support for binary data, where the only difference here is $\pi(v_{i}^{j})$ instead of binary value $v_{i}^{j}$.

The copula support always satisfies $\eta(J)\leq 0.5$ by definition and has the monotonicity with respect to the inclusion relationship of combinations: $\eta(J)\geq\eta(J\cup\{j\})$ for all $J\in 2^{V}$ and $j\in V\setminus J$. Note that, although Tatti (2013) 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 and the FWER was not controlled.

2.1 Statistical Testing

Now we can formulate our problem of finding feature interactions, or itemset mining on continuous features, that are associated with class labels as follows: Let $Y$ be an output binary variable of which class labels are realizations. The task is to find all feature combinations $J\in 2^{V}$ such that the null hypothesis $X_{J}\perp\!\!\!\perp Y$, that is, $X_{J}$ and $Y$ are statistically independent, 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$.

As a statistical test, we propose to use a likelihood ratio test Fisher (1922), which is a generalized $\chi^{2}$-test and often called G-test Woolf (1957). Although Fisher’s exact test has been used as the standard statistical test in recent studies Llinares-López et al. (2015); Sugiyama et al. (2015); Terada et al. (2013a), it can be applied to only discrete test statistics and cannot be used in our setting.

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^{\prime})$ with $l,l^{\prime}\in\{0,1\}$ and can be described as a four-dimensional probability vector $\boldsymbol{p}$:

[TABLE]

Let $\boldsymbol{p}_{\mathrm{E}}$ be the probability vector under the null hypothesis $X_{J}\perp\!\!\!\perp Y$ and $\boldsymbol{p}_{\mathrm{O}}$ be the empirical 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$ is translated into the condition $D_{\mathrm{KL}}(\boldsymbol{p}_{\mathrm{O}},\boldsymbol{p}_{\mathrm{E}})=0$. In the G-test, which is a special case of likelihood ratio test, the test statistic is given as $\lambda=2ND_{\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 probability vector $\boldsymbol{p}_{\mathrm{E}}$ under the null is given as

[TABLE]

where $\eta(J)$ is the copula support of a feature combination $J$ and $r_{l}$ is the ratio of the label $l\in\{0,1\}$ in the dataset. In contrast, the observed probability vector $\boldsymbol{p}_{\mathrm{O}}$ is given as

[TABLE]

where $\eta_{1}(J)=(1/N)\boldsymbol{x}_{J}\cdot\boldsymbol{y}$ and $\eta_{0}(J)=(1/N)\boldsymbol{x}_{J}\cdot(1-\boldsymbol{y})$ with $\boldsymbol{x}_{J}$ defined in Equation (1). These two distributions are shown in Table 1. While it is 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, such situation does not usually occur since $\eta(J)$ is not too large as $\eta(J)\leq 0.5$ by definition and not too small as such combinations are not testable, which will be shown in Section 2.3.

In the following, we write the set of all possible probability vectors by $\boldsymbol{\mathcal{P}}=\{\,\boldsymbol{p}\mid p_{i}\geq 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\,\}$.

2.2 Multiple Testing Correction

Since we have $2^{d}$ hypotheses as each feature combination translated into a hypothesis, we need to perform multiple testing correction to control the FWER (family-wise error rate), otherwise we will find massive false positive combinations. The most popular multiple testing correction method is Bonferroni correction Bonferroni (1936). 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^{d}$ in our case. Each hypothesis $J$ is declared as significant only if $\text{\emph{p}-value}(J)<\delta_{\mathrm{Bon}}=\alpha/2^{d}$. Then the resulting FWER $<\alpha$ is guaranteed. However, it is well known that Bonferroni correction is too conservative. In particular in our case, the correction factor $2^{d}$ 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 Tarone (1990) and widely used in the literature of significant pattern mining Llinares-López et al. (2015); Papaxanthos et al. (2016); Sugiyama et al. (2015); Terada et al. (2013a), which allows us to prune unnecessary feature combinations that can never be significant while controlling the FWER at the same level. Hence Tarone’s testability always offers better results than Bonferroni correction if it is applicable. The only requirement of Tarone’s testability is the existence of the lower bound of the $p$-value given the marginals of the contingency table, which are $\eta(J)$ and $r_{0}$ (or $r_{1}$) in our setting. In the following, we prove that we can analytically obtain the tight upper bound of the KL divergence, which immediately leads to the tight lower bound of the $p$-value.

Theorem 1**.**

(Tight upper bound of the KL divergence)* For $a,b\in[0,1]$ with $a\leq b\leq 1/2$ and a probability vector $\boldsymbol{p}_{\mathrm{E}}=(ab,a(1-b),(1-a)b,(1-a)(1-b))$,*

[TABLE]

for all $\boldsymbol{p}\in\boldsymbol{\mathcal{P}}(a,b)$ and this is tight (see Appendix for its proof).

Proof.   Let $f(x)=D_{\mathrm{KL}}(\boldsymbol{p},\boldsymbol{p}_{\mathrm{E}})$ with $\boldsymbol{p}=(x,a-x,b-x,(1-b)-(a-x))$ and $a^{\prime}=1-a$, $b^{\prime}=1-b$. We have

[TABLE]

The second derivative of $f(x)$ is given as

[TABLE]

which is always positive from the constraint $0<x<\min\{a,b\}=a$. Thus $f(x)$ is maximized when $x$ goes to [math] or $a$. The limit of $f(x)$ is obtained as

[TABLE]

To check which is larger, let $\delta$ be the difference $\lim_{x\to a}f(x)-\lim_{x\to 0}f(x)$. Then it is obtained as $\delta=-b\log b+b^{\prime}\log b^{\prime}+(b-a)\log(b-a)-(b^{\prime}-a)\log(b^{\prime}-a)$. The partial derivative of $\delta$ with respect to $a$ is $\partial\delta/\partial a=\log(1-a-b)-\log(b-a)$, which is always positive as $(1-a-b)-(b-a)=1-2b\geq 0$ with the condition $b\leq 1/2$. Hence the difference $\delta$ takes the minimum value at $a=0$ and we obtain

[TABLE]

Thus $f(x)$ is the tight upper bound when $x\to a$.$\hfill\square$

Here we formally introduce how to prune unnecessary hypotheses by Tarone’s testability. Let $\psi(J)$ be the lower bound of the $p$-value of $J$ obtained from the upper bound of the KL-divergence proved in Theorem 1. Suppose that $J_{1},J_{2},\dots,J_{2^{d}}$ be the sorted sequence of all combinations such that

[TABLE]

is satisfied. Let $m$ be the threshold such that

[TABLE]

Tarone (1990) 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 $\text{\emph{p}-value}(J)<\delta_{\mathrm{Tar}}=\alpha/m$. On the other hand, combinations $J_{m+1},\dots,J_{2^{d}}$ are called untestable as they can never be significant. Since $m\ll 2^{d}$ usually holds, we can expect to obtain higher statistical power in Tarone’s method compared to Bonferroni method. Moreover, we present C-Tarone in the next subsection, which can enumerate such testable combinations without seeing untestable ones, hence we overcome combinatorial explosion of the search space of combinations.

Let us denote by $B(a,b)$ the upper bound provided in Equation (3). We analyze the behavior of the bound $B(a,b)$ as a function of $a$ with fixed $b$. This is a typical situation in our setting as $a$ corresponds to the copula support $\eta(J)$, which varied across combinations, while $b$ corresponds to the class ratio $r_{1}$, which is fixed in each analysis. Assume that $b\leq 1/2$. When $a<b$, we have $\partial B(a,b)/\partial a=\log(1-a)/(b-a)>0$, hence it is monotonically increases as $a$ increases. When $b<a<1/2$, we have $\partial B(a,b)/\partial a=\log(a-b)/a<0$, 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 4 in Appendix.

2.3 Algorithm

We present the algorithm of C-Tarone that efficiently finds testable combinations $J_{1},J_{2},\dots$, $J_{m}$ such that $\psi(J_{1})\leq\psi(J_{2})\leq\dots\leq\psi(J_{m})$, where $m$ satisfies the condition (4). We summarize C-Tarone in Algorithm 1, which performs depth-first search to find $J_{1},J_{2},\dots$, $J_{m}$ such that $\psi(J_{1})\leq\psi(J_{2})\leq\dots\leq\psi(J_{m})$, where $m$ satisfies the condition (4). Suppose that $r_{1}\leq 1/2\leq r_{0}$. Since the lower bound of the $p$-value $\psi(J)$ takes the minimum value when $\eta(J)=r_{1}$ (see Figure 4 in Appendix by letting $a=\eta(J)$ and $b=r_{1}$) and is monotonically decreasing as $\eta(J)$ decreases, for any $\eta(J)\leq r_{1}$, $\mathcal{C}\supseteq\mathcal{T}$ with $\mathcal{C}=\{J\in 2^{V}\mid\eta(J)\geq\sigma\}$ and $\mathcal{T}=\{J\in 2^{V}\mid\psi(J)\leq B(\sigma,r_{1})\}$ is always guaranteed, where $\sigma$ is a threshold for copula supports. Thus if the condition (4) is satisfied for some $m\leq|\mathcal{T}|$, the $m$ smallest combinations in $\mathcal{T}$ are the testable combinations.

Moreover, since $\eta(J)$ has the monotonicity with respect to the inclusion relationship, that is, $\eta(J)\geq\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)\geq\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 the threshold $\sigma$ to enumerate testable combinations while pruning massive unnecessary combinations. 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 (4). If $|\mathcal{T}|B(\sigma,r_{1})>\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}|$.

Since Algorithm 1 always outputs the set of testable combinations $\mathcal{T}$, which directly follows from $\mathcal{C}\supseteq\mathcal{T}$ and the monotonicity of $\eta(J)$, the SignificanceTesting function finds all significant combinations with the FWER $\leq\alpha$. Thus Algorithm 1 always finds all significant feature combinations with controlling the $FWER$ under $\alpha$. Moreover, C-Tarone is independent of the feature ordering and the above completeness with respect to the set of significant combinations is always satisfied. The time complexity of C-Tarone is $O(|\mathcal{T}|)$.

3 Experiments

We examine the effectiveness and the efficiency of C-Tarone using synthetic and real-world datasets. We used Amazon Linux AMI release 2017.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.5. The FWER level $\alpha=0.05$ throughout experiments.

There exists no existing method that can enumerate significant feature combinations from continuous data with multiple testing correction. Thus we compare C-Tarone 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}^{\leq\mathop{\mathrm{med}}(j)},v_{i}^{>\mathop{\mathrm{med}}(j)})$, where $v_{i}^{\leq\mathop{\mathrm{med}}(j)}=\mathtt{1}$ if $v_{i}^{j}\leq\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 $2d$ features. We use the state-of-the-art significant itemset mining algorithm LAMP ver.2 Minato et al. (2014) and employ implementation provided by Llinares-López et al. (2018), which incorporates the fastest frequent pattern mining algorithm LCM Uno et al. (2004) to enumerate testable feature combinations from binary data. Note that, in terms of runtime, comparison with the brute-force approach of testing all the $2^{d}$ combinations with the Bonferroni correction is not valid as the resulting FWER is different between the brute-force and the proposed method.

We also tried other binarization approaches, interordinal scaling used in numerical pattern mining Kaytoue et al. (2011) and interval binarization in subgroup discovery Grosskreutz and Rüping (2009) as a preprocessing of significant itemset mining. However, both preprocessing generate too many binarized dense features, resulting in the lack of scalability in the itemset mining step in the enumeration of testable combinations. Hence we do not employ them as comparison partners. Details of these binarization techniques are summarized in Appendix.

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

To examine the quality of detected combinations, we compute precision, recall, and the F-measure by comparing such combinations with those obtained by the standard decision tree method CART Breiman et al. (1984), which obtains multiplicative combinations of features in the form of binary trees. We used the rpart function in R with its default parameter setting, where the Gini index is used for splitting and the minimum number of data points that must exist in a node is $20$. We apply the decision tree to each dataset and retrieve all the paths from the root to leaf nodes of the learned tree. In each path, we use the collection of features used in the path as a positive feature combination of the ground truth, that is, a feature combination found by C-Tarone (or binarization method) is deemed to be true positive if and only if it constitutes one of full paths from the root to a leaf of the learned decision tree. Note that the FWER is always controlled under $\alpha$ in both of C-Tarone and binarization method, and our aim is to empirically examine the quality of the feature combinations compared to those selected by a standard decision tree. We used not forests but a single tree as the ground truth depends on the number of trees if we use forests, resulting in arbitrary results.

Please note that, due to the fact that there are up to $2^{d}$ feature combinations and we only count an exact match between a retrieved pattern and a true pattern as a hit, precision and F-measure will be close to 0. Still, we compute them to allow for a comparison of the relative performance of the different approaches. Evaluation criteria that take overlaps between retrieved patterns and true patterns (partial matches) into account would lead to higher levels of precision, but they are a topic of ongoing research and not a focus of this work.

Results on Synthetic Data. First we evaluate C-Tarone on synthetic data with varying the sample size $N$ from $1,000$ to $200,000$, the number $d$ of features from $20$ to $100$, and setting the class ratio to $r_{1}=0.5$ or $r_{1}=0.2$, i.e., the number $N_{1}$ of samples in the minor class is $N/2$ or $N/5$. 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 [math] to $1$ and assign the class label $1$ to the first $N_{1}$ data point. Then, for the $N_{1}$ data points in the class $1$, we pick up one of the 20% of associated features and copy it to every associated feature with adding Gaussian noise with $(\mu,\sigma^{2})=(0,0.1)$. Hence there is no correlation in the class [math] across all features and there are positive correlations among such 20% of features in the class $1$. The other 80% are uninformative features.

Results are plotted in Figure 2 for $r_{1}=0.5$ (classes are balanced). See Figure 5 in Appendix for $r_{1}=0.2$ (classes are imbalanced). In the figure, we plot results with varying $N$ while fixing $d=20$ on the left column and those with varying $d$ while fixing $N=1,000$ on the right column.

In comparison with the median-based binarization method (plotted in blue), C-Tarone (plotted in red) has a clear advantage regarding the precision, which is several orders of magnitude higher than the binarization method in every case. This is why binarization method cannot distinguish correlated and uncorrelated combinations as we discussed in Introduction and illustrated in Figure 1, resulting in including uncorrelated features into significant combinations. Although recall is competitive across various $N$ and $d$, in all cases, the F-measure of C-Tarone are higher than those of the median-based binarization method. 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 and not used in the decision tree, tend to reach statistical significance.

Although the algorithm LCM (itemset mining) used in the binarization approach is highly optimized with respect to the efficiency, C-Tarone is competitive with it on all datasets, as can be seen in Figure 2 (bottom row). To summarize, we observe that C-Tarone improves over the competing binarization approach in terms of the F-measure in detecting higher quality feature combinations in classification.

Results on Real Data. We also evaluate C-Tarone on real-world datasets shown in Table 2 in Appendix, which are benchmark datasets for binary classification from the UCI repository Lichman (2013). To clarify the exponentially large search space, we also show the number $2^{d}$ of candidate combinations for $d$ features in the table. All datasets are balanced to maximize the statistical power for comparing detected significant combinations, i.e. $r_{1}=0.5$. If they are not balanced in the original dataset, we randomly subsample data points from the larger class.

We summarize results in Figure 3. Again, C-Tarone shows higher precision in all datasets than the binarization method and better or competitive recall, resulting in higher F-measure scores in all datasets. In addition, running time is competitive with the binarization method, which means that C-Tarone can successfully prune the massive candidate space for significant feature combinations. These results demonstrate the effectiveness of C-Tarone.

4 Conclusion

In this paper, we have proposed a solution to the open problem of finding all multiplicative feature interactions between continuous features that are significantly associated with an output variable after rigorously controlling for multiple testing. While interaction detection with multiple testing has been studied before, our approach, called C-Tarone, is the first to overcome the problem of detecting all higher-order interactions from the enormous search space $2^{d}$ for $d$ 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.

Appendix A Related Work

Since the only line of work that tries to find multiplicative feature combinations while controlling the FWER is significant pattern mining, we provide an overview of this field in the following.

Significant pattern mining introduces statistical significance into the task of contrast (or discriminative) pattern mining Dong and Bailey (2013), where the objective is to find discriminative patterns with respect to class partitioning of a dataset. After early work on multiple testing correction in association rule mining Hämäläinen (2012); Webb (2007), Terada et al. (2013a) were the first to achieve control of the FWER in itemset mining, by successfully combining a pattern mining algorithm and Tarone’s trick Tarone (1990). The enumeration algorithm has been improved in LAMP ver. 2 Minato et al. (2014) and significant subgraph mining Sugiyama et al. (2015). To date, significant pattern mining has been extended to various types of tests and data, including a Westfall-Young permutation test to treat independencies among patterns Llinares-López et al. (2015); Terada et al. (2013b), logistic regression Terada et al. (2016) or a Cochran–Mantel–Haenszel (CMH) test Llinares-López et al. (2017); Papaxanthos et al. (2016) for categorical covariates, hypothesis streams Webb and Petitjean (2016), and top-$K$ significant patterns Pellegrina and Vandin (2018). However, none of the above studies succeeded to directly perform significant pattern mining on continuous variables without prior binarization.

Although the field of subgroup discovery Atzmueller (2015); Novak et al. (2009); Herrera et al. (2011) also considers measures of statistical dependence for finding multiplicative feature combinations, e.g. Grosskreutz and Rüping (2009); Mampaey et al. (2015); van Leeuwen and Ukkonen (2016), none of these methods accounts for multiple testing by controlling the FWER.

Appendix B Additional Binalization Methods

In interordinal scaling, each binarized feature is in the form of “$\leq a$” or “$\geq a$”, where endpoints $a$ are from a dataset, that is, $a\in\{v_{1}^{j},v_{2}^{j},\dots,v_{N}^{j}\}$ for a feature $j\in\{1,2,\dots,d\}$. Thus, for an $d$-dimensional real-valued vector $\boldsymbol{v}_{i}\in D$, each element $v_{i}^{j}$ is expanded as the $2N$-dimensional binary vector such that

[TABLE]

where each value for the binarized feature $v_{i}^{\leq v_{k}^{j}}=\mathtt{1}$ if $v_{i}^{j}\leq v_{k}^{j}$ and 0 otherwise. As a result, the dataset $D$ is converted into the binary dataset with $2dN$ features. In interval binarization, each binarized feature is in the form of “$(a,b]$”, where endpoints $a$, $b$ are from data, and each element $v_{i}^{j}$ of an $d$-dimensional vector $\boldsymbol{v}_{i}$ is expanded as the $(N(N-1)/2)$-dimensional binary vector such that

[TABLE]

Thus a dataset $D$ is converted into the binary dataset with $dN(N-1)/2$ features. Both interordinal scaling and interval binarization could finish their computation for a tiny dataset with $(N,d)=(50,5)$ in approximately 24 hours, but did not finish after 48 hours for $(N,d)=(100,10)$, and they exceeded the memory limit (2.0 TB) for larger datasets.