Finding Statistically Significant Interactions between Continuous Features

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 combinations 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;

∗Contact Author

Llinares-L´opez 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:

1. How to assess the significance for a multiplicative interaction of continuous features?
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. How to manage the combinatorial explosion of the candidate space, where the number of possible interactions is 2d 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`es, 2016], none of these approaches addresses the problem of dealing with the combinatorial explosion of the 2d -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 ProposedMethod: C-Tarone

Given a supervised dataset D = { (v1, y1), (v2, y2), . . . , (vN, yN) }, where each data point is a pair of an ddimensional vector vi = (v1 i , v2 i , . . . , vd i ) ∈ Rd and a binary class label yi ∈ {0, 1}. We denote the set of features by V = {1, 2, . . . , d} and the power set of features by 2V . For each feature j ∈ {1, 2, . . . , d}, we write vj = (vj1 , vj2 , . . . , vjN ), which is the N-dimensional vector composed of the jth 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 η(J) ∈ R of a feature combination J ∈ 2V in a dataset. For each combination J ∈ 2V , 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, vj i ∈ {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 ∈ 2V is called an itemset, and the support of J is defined as η(J) = (1/N) ∑Ni=1 ∏j∈J vj i , which corresponds to joint probability of J. The support is always from 0 to 1, and we can introduce a binary random variable XJ corresponding to the joint occurrence of J, where XJ = 1 if features in J jointly occurs and XJ = 0 otherwise, η(J) corresponds to the empirical estimate of the probability Pr(XJ = 1).

This approach can be generalized to continuous (realvalued) 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 XJ of joint occurrence of J and estimate η(J) = Pr(XJ = 1) from continuous data. The key idea is to use the normalized ranking for each feature j ∈ V , which is defined as π(vj i ) = (k - 1)/(N - 1) if vj i is the kth smallest value among N values vj1 , vj2 , . . . , vjN . Then we convert a dataset D with respect to features in a combination J ∈ 2V into a single N-dimensional vector xJ = (xJ1, xJ2, . . . , xJN) ∈ [0, 1]N by

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

Intuitively, joint occurrence of J means that rankings among features in J match with each other. The definition in Equation (2) is analogous to the support for binary data, where the only difference here is π(vj i ) instead of binary value vj i .

The copula support always satisfies η(J) ≤ 0.5 by definition and has the monotonicity with respect to the inclusion relationship of combinations: η(J) ≥ η(J∪{j}) for all J ∈ 2V and j ∈ V \ 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 ∈ 2V such that the null hypothesis XJ ⊥⊥ Y , that is, XJ 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 α.

As a statistical test, we propose to use a likelihood ratio test [Fisher, 1922], which is a generalized χ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´opez et al., 2015; Sugiyama et al., 2015;

Table 1: Contingency tables.

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) = rl for each class l ∈ {0, 1}. From two binary variables XJ and Y , we obtain a 2 × 2 contingency table, where each cell denotes the joint probability Pr(XJ = l, Y = l′ ) with l, l′ ∈ {0, 1} and can be described as a four-dimensional probability vector p:

Let pE be the probability vector under the null hypothesis XJ ⊥⊥ Y and pO be the empirical vector obtained from N observations. The difference between two distributions pO and pE can be measured by the Kullback-Leibler (KL) divergence DKL(pO, pE) = ∑i pOi log(pOi/pEi), and the independence XJ ⊥⊥ Y is translated into the condition DKL(pO, pE) = 0. In the G-test, which is a special case of likelihood ratio test, the test statistic is given as λ = 2NDKL(pO, pE), which follows the χ2 -distribution with the degree of freedom 1.

In our case, for each combination J ∈ 2V , the probability vector pE under the null is given as

where η(J) is the copula support of a feature combination J and rl is the ratio of the label l ∈ {0, 1} in the dataset. In contrast, the observed probability vector pO is given as

where η1(J) = (1/N)xJ · y and η0(J) = (1/N)xJ · (1-y) with xJ defined in Equation (1). These two distributions are shown in Table 1. While it is known that the statistic λ does not exactly follow the χ2 -distribution if one of components of pE is too small, such situation does not usually occur since η(J) is not too large as η(J) ≤ 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 P = {p | pi ≥ 0, ∑ pi = 1}, and its subset given marginals a and b as P(a, b) = { p ∈ P | p1 + p2 = a, p1 + p3 = b }.

2.2 Multiple Testing Correction

Since we have 2d 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 α (e.g. α = 0.05) is corrected as δBon = α/m, where m is called a correction factor and is the number of hypotheses m = 2d in our case. Each hypothesis J is declared as significant only if p-value(J) < δBon = α/2d . Then the resulting FWER < α is guaranteed. However, it is well known that Bonferroni correction is too conservative. In particular in our case, the correction factor 2d 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´opez 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 η(J) and r0 (or r1) 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 ∈ [0, 1] with a ≤ b ≤ 1/2 and a probability vector pE = (ab, a(1 - b), (1 - a)b, (1 - a)(1 - b)),

for all p ∈ P(a, b) and this is tight.

Here we formally introduce how to prune unnecessary hypotheses by Tarone's testability. Let ψ(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 J1, J2, . . . , J2d be the sorted sequence of all combinations such that

is satisfied. Let m be the threshold such that

Tarone [1990] showed that the FWER is controlled under α with the correction factor m. The set T = {J1, J2, . . . , Jm} is called testable combinations, and each J ∈ T is significant if p-value(J) < δTar = α/m. On the other hand, combinations Jm+1, . . . , J2d are called untestable as they can never be significant. Since m ≪ 2d 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 η(J), which varied across combinations, while b corresponds to the class ratio r1, which is fixed in each analysis. Assume that b ≤ 1/2. When a < b, we have ∂B(a, b)/∂a = log(1 - a)/(b - a) > 0, hence it is monotonically increases as a increases. When b < a < 1/2, we have ∂B(a, b)/∂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 J1, J2, . . . , Jm such that ψ(J1) ≤ ψ(J2) ≤ · · · ≤ ψ(Jm), where m satisfies the condition (4). We summarize C-Tarone in Algorithm 1, which performs depth-first search to find J1, J2, . . . , Jm such that ψ(J1) ≤ ψ(J2) ≤ · · · ≤ ψ(Jm), where m satisfies the condition (4). Suppose that r1 ≤ 1/2 ≤ r0. Since the lower bound of the p-value ψ(J) takes the minimum value when η(J) = r1 (see Figure 5 in Appendix by letting a = η(J) and b = r1) and is monotonically decreasing as η(J) decreases, for any η(J) ≤ r1, C ⊇ T with C = {J ∈ 2V | η(J) ≥ σ} and T = {J ∈ 2V | ψ(J) ≤ B(σ, r1)} is always guaranteed, where σ is a threshold for copula supports. Thus if the condition (4) is satisfied for some m ≤ |T |, the m smallest combinations in T are the testable combinations.

Moreover, since η(J) has the monotonicity with respect to the inclusion relationship, that is, η(J) ≥ η(J ∪ {j}) for all j ∈ V \ J, finding the set C is simply achieved by DFS (depth-first search): Starting from the smallest combination ∅ with assuming η(∅) = 1, if η(J) ≥ σ for a combination J, we update J by adding j ∈ V \ J and recursively check J. Otherwise if η(J) < σ, we prune all combinations K ⊇ J as η(K) < σ holds.

Our algorithm dynamically updates the threshold σ to enumerate testable combinations while pruning massive unnecessary combinations. First we set σ = 0, which means that all combinations are testable. Whenever we find a new combination J, we update C and T and check the condition (4). If |T |B(σ, r1) > α holds, σ is too low and T will become too large with the current σ, hence we update σ to minJ∈C η(J) and remove the combination argminJ∈C η(J). Finally, when the algorithm stops, it is clear that the set T coincides with the set of testable combinations, hence each combination J ∈ T is significant if the p-value of J is smaller than δTar = α/|T |.

Since Algorithm 1 always outputs the set of testable combinations T , which directly follows from C ⊇ T and the monotonicity of η(J), the SIGNIFICANCETESTING function finds all significant combinations with the FWER ≤ α. Thus Algorithm 1 always finds all significant feature combinations with

Algorithm 1: C-Tarone.

1 C-TARONE(D, α) 2 σ ← 0; C ← ∅; // σ is a global variable 3 DFS(∅, 0, C, D, α); 4 T ← {K ∈ C | ψ(K) ≤ B(σ, r1)}; // The set of testable combinations 5 SIGNIFICANCETESTING(T , α); 6 DFS(J, jprev, C, D, α) 7 foreach j ∈ {jprev + 1, . . . , d} do 8 J ← J ∪ {j}; 9 Compute η(J) by Equation (2); 10 if η(J) > σ then 11 Compute ψ(J); C ← C ∪ {J}; 12 T ← {K ∈ C | ψ(K) ≤ B(σ, r1)}; 13 while |T |B(σ, r1) ≥ α do 14 Jmin ← argminK∈C η(K); 15 σ ← η(Jmin); C ← C \ {Jmin}; 16 T ← T \ {Jmin}; 17 DFS(J, j, C, D, α); 18 J ← J \ {j}; 19 SIGNIFICANCETESTING(T , α) 20 foreach J ∈ T do 21 if p-value(J) < α / |T | then output J;

controlling the FWER under α. 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(|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 α = 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 ∈ V , we pick up the median of vj = (vj1 , vj2 , . . . , vjN ), denoted by med(j), and binarize each value vj i as a pair (v≤med(j) i , v>med(j) i ), where v≤med(j) i = 1 if vj i ≤ med(j) and 0 otherwise, and v>med(j) i = 1 if vj i > med(j) and 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´opez 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 2d combinations with the Bonferroni correction is not valid as the resulting FWER is different between the bruteforce 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¨uping, 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 Section C 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 α 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 2d 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 r1 = 0.5 or r1 = 0.2, i.e., the number N1 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 0 to 1 and assign the class label 1 to the first N1 data point. Then, for the N1 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 (µ, σ2 ) = (0, 0.1). Hence there is no correlation in the class 0 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 r1 = 0.5 (classes are balanced). See Figure 5 in Appendix for r1 = 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 medianbased 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 2d 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. r1 = 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 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 2d 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.

Acknowledgments

This work was supported by JSPS KAKENHI Grant Numbers JP16K16115, JP16H02870, and JST, PRESTO Grant Number JPMJPR1855, Japan (M.S.). This work was funded in part by the SNSF Starting Grant "Significant Pattern Mining" (K.B.).

A Proof of Theorem 1

Let f (x) = DKL(p, pE) with p = (x, a - x, b - x, (1 - b) - (a - x)) and a′ = 1 - a, b′ = 1 - b. We have

The second derivative of f(x) is given as

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

To check which is larger, let δ be the difference limx→a f(x) - limx→0 f(x). Then it is obtained as δ = -b log b+b′ log b′ +(b-a) log(b-a)-(b′ -a) log(b′ -a). The partial derivative of δ with respect to a is ∂δ/∂a = log(1 - a - b) - log(b - a), which is always positive as (1 - a - b) - (b - a) = 1 - 2b ≥ 0 with the condition b ≤ 1/2. Hence the difference δ takes the minimum value at a = 0 and we obtain

Thus f(x) is the tight upper bound when x → a. □

B RelatedWork

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¨am¨al¨ainen, 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´opez et al., 2015; Terada et al., 2013b], logistic regression [Terada et al., 2016] or a Cochran-Mantel-Haenszel (CMH) test [Llinares-L´opez et al., 2017; Papaxanthos et al., 2016] for categorical covariates, hypothesis streams [Webb and Petitjean, 2016], and topK 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¨uping, 2009; Mampaey et al., 2015; van Leeuwen and Ukkonen, 2016], none of these methods accounts for multiple testing by controlling the FWER.

C Additional BinarizationMethods

In interordinal scaling, each binarized feature is in the form of "≤ a" or "≥ a", where endpoints a are from a dataset, that is, a ∈ {vj1 , vj2 , . . . , vjN } for a feature j ∈ {1, 2, . . . , d}. Thus, for an d-dimensional real-valued vector vi ∈ D, each element vj i is expanded as the 2N-dimensional binary vector such that

where each value for the binarized feature v≤vjk i = 1 if vj i ≤ vjk 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 vj i of an d-dimensional vector vi is expanded as the (N(N - 1)/2)- dimensional binary vector such that

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.