Fast Conditional Independence Test for Vector Variables with Large Sample Sizes

2018·Arxiv

Abstract

Abstract

We present and evaluate the Fast (conditional) Independence Test (FIT) – a non-parametric conditional independence test. The test is based on the idea that when P(X | Y, Z) = P(X | Y ), Z is not useful as a feature to predict X, as long as Y is also a regressor. On the contrary, if might improve prediction results. FIT applies to thousand-dimensional random variables with a hundred thousand samples in a fraction of the time required by alternative methods. We provide an extensive evaluation that compares FIT to six extant nonparametric independence tests. The evaluation shows that FIT has low probability of making both Type I and Type II errors compared to other tests, especially as the number of available samples grows. Our implementation of FIT is publicly available1.

1 Introduction

Two random variables X and Y are conditionally independent given a third variable Z if and only if P(X, Y | Z) = P(X | Z)P(Y | Z). We denote this relationship by , and its negation, the case of conditional dependence, by . This article develops and evaluates the Fast (conditional) Independence Test (FIT), a nonparametric conditional or unconditional independence test. Given a finite sample from the joint distribution P(X, Y, Z), FIT returns the p-value under the null hypothesis that . FIT applies to datasets with large sample sizes on scalar- or vector-valued variables, and returns in short time. Extensive empirical comparison with existing alternatives shows that FIT is currently the only conditional independence test to achieve this goal.

1.1 Motivation

Independence tests are ubiquitous in scientific inference. They are the main work horse of classical hypothesis testing (see Fisher [1925], Chapter 21), which underlies most experimental and many observational techniques of data analysis. A conditional independence test (CIT) extends this idea by testing for independence between two variables given the value of a third, or more generally, given the values of a set of further variables [Dawid, 1979]. Probabilistic dependence – conditional or not – is often taken as an indicator of informational exchange (see e.g. Cover and Thomas [2012], Chapter 2), causal connection [Dawid, 1979, Spirtes et al., 2000] or simply as tool for prediction and

diagnosis [Koller and Sahami, 1996]. Consequently, knowledge of the (conditional) independence relations among a set of variables provides some of the most basic indicators for scientific relations of interest.

Our motivation comes from the challenge of causal discovery. Many extant causal discovery algorithms use conditional independence tests to infer the underlying causal structure over a set of variables from observational data [Pearl, 2000, Spirtes et al., 2000, Hoyer et al., 2009]. The most general of these methods do not make any specific assumptions about the functional form or parameterization of the underlying causal system. Discovery of independence structure from data, rather than estimation of any specific parameters, allows these methods to establish causal structure. However, such generality can only be maintained if the CITs used to determine the independence structure in the first place are themselves non-parametric. That is, the CITs have to check for general dependence relations among the variables, and not just, say, for a non-zero (partial) correlation.

A further motivation derives from the need for CITs that are both fast to compute and flexible in the variables they can handle. Generally, causal discovery algorithms have to perform a number of tests that is at least polynomial, if not exponential, in the number of variables [Spirtes et al., 2000]. For genetic datasets where thousands of genes are measured, or neuroscience data containing tens of thousands of voxels, this implies millions of independence tests. Moreover, for example, in social network data, the variables may be categorical or continuous, and scalar or vector valued. These facts indicate a clear need for non-parametric tests that are fast and flexible with regard to the nature of the variables being tested.

While our motivation derives from the domain of causal discovery, the test is completely domain general and can be applied to any area in which (conditional) independence tests are used.

1.2 Previous work

For continuous variables, extant non-parametric CITs have largely focused on kernel methods [Schölkopf and Smola, 2002]. The methods differ in how the null-distribution is estimated, which distance measure is used in the reproducing kernel Hilbert space (RKHS) or which test statistic is applied. Here we consider six methods we are aware of for comparison:

The Conditional Hilbert-Schmidt Independence Criterion (CHSIC) [Fukumizu et al., 2008] maps the data into a RKHS. It determines independence using the normalized cross-covariance operator, which can be thought of as an extension of standard correlations to higher order moments.

In contrast, the Kernel-based Conditional Independence test (KCIT) [Zhang et al., 2012] derives a test statistic based on the traces of the kernel matrices and approximates the null-distribution using a Gamma distribution. This avoids the permutations required for CHSIC, but comes at the expense of significant matrix operations.

The Kernel Conditional Independence Permutation test (KCIPT) [Doran et al., 2014], in many ways a variant of CHSIC, replaces the random sampling of permutations for the null-distribution with a specifically learned permutation that satisfies, among other criteria, that it is representative of conditional independence. KCIPT also replaces the normalized cross-covariance operator with the maximum mean discrepancy (MMD) test statistic developed by Gretton et al. [2012].

Kernel-based tests do not scale well to large sample sizes due to the matrix operations required to kernelize the data and embed it in the RKHS. To address this limitation, the Randomized Conditional Independence test (RCIT) [Strobl et al., 2017], uses a fast Fourier transform to speed up the matrix operations. For similar reasons, the Conditional Correlation Independence test (CCI) [Ramsey, 2014] determines independence of X and Y by only checking whether cov(f(X), g(Y )) = 0 for different choices of f, g taken from an (in practice, truncated) set of basis functions. For conditional independence of , the same test is performed on the non-parametric residuals of X and Y each regressed on Z using a uniform kernel.

In the case of categorical data, two variables are independent conditional on a third categorical variable Z just in case they are independent given any specific categorical value of Z = z. Consequently, the simplest way to handle categorical data is to repeat whatever unconditional test is available for each subset of the data corresponding to a fixed value of Z = z. Unless some further structure among the categories (such as a particular order) holds important information about the probabilistic relations between variables, there is nothing further for a conditional independence test to exploit. Consequently, tests conditioning on categorical variables are in the general case no different than unconditional tests, and inevitably sample intensive. However, a particular – and common – challenge arises when the conditioning variable Z is continuous and at least one of X or Y is categorical. In this case we have a "hybrid" situation. The methods mentioned above do not apply due to the categorical nature of X or Y , and for a continuous conditioning variable Z, by definition one cannot explicitly check conditional independence by checking independence for each value of Z.

The development of the above tests has enabled the construction of genuinely non-parametric methods of discovery, but has left several questions of practical importance open. First, it remains unclear how to choose between the tests given that there is no available systematic comparison. Second, the evaluation of tests has focused on a small set of specialized datasets. Third, there has been no evaluation of how these tests perform for high dimensional variables and large sample sizes.

In our evaluation in Sec. 4 we have attempted to faithfully implement all the above tests2 and systematically evaluated them on all the previous data sets, in each case reporting the actual behavior of the resulting p-values. In addition, we vastly extended the range of parameters that were explored to generate some of the existing datasets, and added datasets that explore the hybrid case with categorical and continuous variables described above.

1.3 Contributions

We address the challenge of testing for (conditional) independence when the variables are high-dimensional and sample sizes are high, but when the relation among the variables cannot be well described by a parametric model. Moreover, in order to enable applications for which it is necessary to run tens of thousands of conditional independence tests – for example, establishing the causal graph over hundreds of variables – the test has to run in a short amount of time.

To our knowledge, none of the existing CITs applies with reasonable runtime in such situations. The main contributions of this article are as follows:

1. The Fast Independence Test (FIT), which can process hundred-dimensional conditional independence queries with samples in less than 60s.

2. An extensive evaluation and comparison of FIT and alternatives (CHSIC, KCIT, KCIPT, RCIT, CCI) on a wide range of datasets.

3. A parallelized Python implementation of FIT, easily installable through Python Package Index3.

2 Fast (conditional) Independence Test (FIT)

The intuition behind FIT is that if , then prediction of Y using both X and Z as covariates should be more accurate than prediction of Y when only Z is used as the covariate. On the other hand, if , then the accuracy of the prediction of Y should not change whether just Z, or X and Z are used as covariates. While this is not always true (see Sec 5.3), it is a relatively weak assumption to make compared to, for example, assuming a specific parametric form of the involved distributions or linearity of the functional relationships. Alg. 1 contains pseudocode for the implementation of FIT. Our implementation uses decision tree regression (DTR, Breiman et al. [1984])4 to predict Y using both X, Z, and also using Z only. he rationale behind this choice and the alternatives are discussed in Sec. 5.2.

We measure the accuracy of the prediction in terms of the mean squared error (MSE). Central to our approach is the simplifying assumption that if and only if the mean squared error (MSE) of the algorithm trained using both X and Z is not smaller than the MSE of the algorithm trained using Z only.

2.1 Obtaining the p-value

Our implementation runs DTR n_perm times to learn both the function as well as . The train/test (or in-sample/out-of-sample) split is resampled each of the n_perm training repetitions. We store the n_perm MSEs resulting from regressing Y on X, Z in a list mses_x, and the n_perm MSEs resulting from regressing Y on Z in a list mses_nox.

The remaining challenge is to determine the p-value of the null hypothesis that the first list contains on average values equal to or larger than the second list. To do this, we run a one-tailed t-test Student [1908], with the null hypothesis that the corresponding values in the two arrays are on average equal. A small p-value rejects the null hypothesis, meaning that mses_x are significantly smaller than . As an alternative to using the t-test (which assumes the distribution of MSEs under the null hypothesis is Gaussian), we experimented with using bootstrap to estimate the null distribution [Efron, 1992]. Since the results were nearly identical across our evaluation tasks, for simplicity of presentation we use the t-test in the final method.

2.2 Speed of the Implementation

Two details of the implementation are worth mentioning:

1. Both the cross-validation procedure (Lines 13, 17) and the training loop (Lines 23- 37) are parallelized to automatically take up all the available CPUs. This provides a multiplicative speedup for larger datasets.

2. Using our algorithm as an unconditional independence test (testing whether ) requires only a simple modification, which our implementation uses automatically when Z is not given. Instead of training using only Z as input (Lines 17 and 36, 37) we train using X permuted randomly across the sampling dimension.

Asymptotically the algorithm’s speed scales as . Sec. 4 shows actual runtimes on an 8 Intel(R) Core(TM) i7-6700K CPU 4.00GHz machine. The runtimes are bounded below by about 1s due to parallelization overhead, and grow to about 100s for a 1000-dimensional dataset with 100K samples.

Algorithm 1 Fast (conditional) Independence Test

3 Evaluating Statistical Test Performance

An ideal nonparametric CIT should be fast, accurate, and generally applicable. To provide a thorough analysis of these factors we found it necessary to deviate from previous work in the evaluation method, as described in this section. We discuss alternative criteria and the problems we found with them in Sec. 5.1. Our evaluation criteria are as follows:

Power: The test should have high power. Whenever , the p-value should be small and should tend to 0 as sample sizes grows.

Size: The size of the test should be equal to or lower than its level . That is, whenever the p-value of the test should not be small for any sample size. Note that in some applications it might be important for the size of the test to be equal (neither lower nor higher) to its level. We discuss this further in Sec. 5.1.

Speed: The test should be fast. This is motivated by the use cases discussed in Sec. 1.1.

Table 1: Summary of evaluation settings. Dimensionality is summed over x, y and z. Type differentiates between continuous (cont) and categorical variables. Complexity is a non-rigorous quantification of the complexity of the conditional distributions that generate the data. See Sec. 3.2 for details.

3.1 Evaluation method

In order to understand the speed, accuracy and breadth of FIT and its competitors, we first chose four settings, described in Sec. 3.2. These four settings cover a range of possible data: low- and high-dimensional data, simple and complex functional relationships between the variables, continuous and mixed discrete-continuous systems5 – Table 1 summarizes the characteristics of each setting. We designed each setting in a way that enables us to vary its dimensionality or complexity. For convenience of plotting the results, we chose nine instantiations of the dimensionality/complexity parameters in each setting, and for each such instantiation, we generated a pair of datasets, one version where and one where

We applied each candidate test described in Sec. 1.2 to each dataset, no matter the dimensionality or complexity. To further understand the change in speed and performance of the tests as the number of samples grows, we ran each test on each dataset for n_samples (number of samples) varying between and logarithmically. For each n_samples we plot the p-value corresponding to the dependent and independent versions of the data, as well as the runtime of each test on each of the datasets and n_samples. This results in plots such as those shown in Fig. 6- 11. These figures provide a very detailed picture of the tests’ performance, which we discuss in Sec. 4. In addition, Fig. 5 reduces the results to answer the question “If I give each test at most 1 minute time, what is the best Type I and Type II error it can achieve on each dataset?”. This question follows naturally from our need to run large numbers of CITs for causal discovery (Sec. 1.1).

3.2 Datasets

To evaluate FIT and compare it to alternatives, we used four settings to generate synthetic datasets, LINGAUSS, CHAOS, HYBRID and PNL, that cover a broad variety of use cases – see Table 1. We here provide descriptions of these settings in some detail to give a better sense of what the data can look like. Figures 1-4 show low dimensional examples from these settings in order to provide at least some illustration of their difficulty and diversity.

LINGAUSS: Any reasonable statistical independence test should work on linear Gaussian models. Surprisingly, previous work was not evaluated on this basic case. For the version of the data, given dim (the desired data dimensionality), we sample n_samples samples from the linear-Gaussian graphical model

identity matrix. 2. Sample independent standard normal variables and arrange them into a

-"edge coefficient" matrix). 3. Sample independent standard normal variables and arrange them into a

Figure 1: An example LINGAUSS dataset. For both the dependent and independent (1-dimensional) instantiation of the dataset we created random samples, and plotted x vs y that correspond to the cluster of z’s between the 50th and 51st percentile of z-values.

For the case, given dim (the desired data dimensionality), we sample n_samples samples from the linear-Gaussian graphical model

1. Let Z again be a dim-dimensional independent standard Gaussian. 2. Sample independent standard normal variables and arrange them into a

-"edge coefficient" matrix). 3. Sample independent standard normal variables and arrange them into a

Note that we deliberately did not create the (conditionally) dependent dataset by adding an extra "edge" between X and Y to the conditionally independent "common cause model" in order to ensure that the dimensionalities of the relation between X and Y conditional on Z remain the same between the independent and dependent case.

We evaluated the CITs using LINGAUSS with dim = 1, 2, 4,... , 256.

CHAOS: Doran et al. [2014] used this dataset to evaluate their CIT. The data is sampled from a highly nonlinear, chaotic dynamical system (see Fig. 2). The dimensionality of the data is fixed: X and Y both have dim = 4 and Z has dim = 2. A scalar parameter controls the complexity

To create the data, take

To create the data, take

Figure 2: Example CHAOS dataset when . a) The generative model for the conditionally independent version of the dataset. Only three timesteps are shown and the third and fourth dimensions of X and Y are omitted, since they are independent noise variables. b) The generative model for the conditionally dependent version of the dataset. c) For the independent version of the dataset we created random samples. We then clustered the z’s into clusters of size roughly 1000 using K-Means, and plotted vs. corresponding to a randomly chosen cluster of z’s (since Z is multivariate continuous, the distribution cluster of z values) is an approximation of P(X, Y | z) for one value of z.) d) Same as c) but for the dependent data.

Finally, to make the task more difficult, append two independent Gaussian variables with N(0, 0.5) as the third and fourth dimension of X and Y . To make the samples independent, we generated samples and chose a random subset of the desired size to form our datasets. Fig. 2a,b illustrate the choices of X, Y, Z in the generative model.

We evaluated the CITs using CHAOS with = .01, .04, .16, .32, .5, .68, .84, .96, .99.

HYBRID: We created this hybrid continuous-categorical dataset as we have not seen any previous methods evaluated on hybrid data (see Sec. 1.1). We deliberately constructed categorical variables, as opposed to (e.g. numerical) discrete variables, whose values still have some order.

Given a count parameter , each sample for a HYBRID()-dataset, say , was constructed using the following pseudocode:

1. Let be a a dim-dimensional continuous vector sampled from the uniform Dirichlet distribution with

2. Draw two independent samples from the multinomial distribution with probability parameter and count parameter

3. Since we want and to represent multi-dimensional categorical variables, convert each dimension of and to the one-hot encoding. For example, take dim=3 and . Suppose for a particular sample. Then might be (0, 2, 0), corresponding to zero draws in the first and third bucket and two draws in the second bucket of the sample from the -multinomial. Then the one-hot encoding of makes . The one-hot encoding enforces the discrete metric on the X-space.

Figure 3: Example HYBRID dataset with and dim= 32. For both the dependent and independent version of the dataset we created random samples, and picked the x, y values that correspond to . The plot shows the heatmap (more saturated color = more samples) of the number of samples falling into given x, y bin.

4. (If creating the dependent version of the dataset), toss an unbiased coin. If heads, then overwrite the value of with the value of

Since Z is a sufficient statistic of the distributions of X and Y , we have in Steps 1-3, and Step 4 is used to create an obvious dependency (which nevertheless turns out to be non-trivial to detect by many of the algorithms, see Sec. 4).

We evaluated the CITs using HYBRID with = (2, 2), (2, 8), (2, 32), (8, 2), (8, 8), (8, 32), (32, 2), (32, 8), (32, 32).

PNL: Various versions of this dataset have been used before in [Zhang et al., 2012, Doran et al., 2014, Strobl et al., 2017]. We used the version from the most recent article [Strobl et al., 2017]. Since this version offers only scalar X, Y, Z, we extended it to the multi-dimensional case in a straightforward way. To create the generating mode for the dim-dimensional PNL dataset:

1. Let -dimensional matrix of independent standard normal scalars. 2. Let Z be a dim-dimensional Gaussian with a covariance and mean 0. 3. Let X be a random function of the first coordinate of Z (see below). 4. Let Y be a (re-sampled) random function of the first coordinate of Z (see below). 5. (To create the dependent version of the dataset,) add a Gaussian variable with mean zero and standard deviation .5 to both X and Y .

This procedure uses “random functions”. To obtain such functions, following Strobl et al. [2017], we sample with equal probability from the set of functions

Note that in this dataset, X and Y are always one-dimensional. We evaluated the CITs using PNL with dim= 1, 2, 4, ..., 256.

Figure 4: Example PNL dataset. For both the dependent and independent (one-dimensional) version of the dataset we created random samples, and plotted x vs. y that correspond to the cluster of z’s between the 50th and 51st percentile of z-values.

4 Empirical Results and Comparison with Previous Work

For each out of the 72 dataset versions (4 settings times 9 parameter instantiations for both the dependent and independent version), we obtain p-values of FIT and each of the five conditional independence tests described in Sec. 1.2. In each case, we varied the number of samples between and logarithmically. In addition, we stopped any test that exceeded the time threshold of about 100s. In this section, for each CIT we picked the results plot that best summarizes the broad behavior of the test across the datasets. In addition, Figure 5 provides a summary view of the full results. Due to their sheer volume we leave the detailed full results to Appendix A. Nevertheless, the full results plots form an important part of the evidence that motivates the discussion in this section.

Figure 5: Type I and II errors at significance level . For each dataset setting (x-axis) and each dataset difficulty (narrow columns within each x-axis block), we computed the errors each method (y-axis) returns on the largest number of samples it can process within 60s (colors indicate error magnitude). Black regions correspond to the cases where an algorithm was not able to process even 100 samples within 100s. Black vertical lines separate the four settings, and for each setting the columns are arranged in increasing difficulty, e.g. for LINGAUSS the first column corresponds to dim=1, last column to dim=256. See Figures 12-35 for full results plots, covering all sample sizes and all p-values.

For all methods, as expected, lower dataset difficulty results in smaller Type I and higher Type II errors. The average of the two error types, shown in the right column of Figure 5, shows that FIT achieves the best performance across the board. The remaining methods either have high errors (CHSIC, RCIT, CCI) or are too slow to process a large proportion of the datasets (KCIT, KCIPT). More detailed remarks on each method’s performance follow.

4.1 FIT

Figure 6 illustrates typical behavior of the test, in this case applied to the PNL dataset. The only data that causes problems for the FIT algorithm is the high-dimensional, sparse versions of the HYBRID dataset. None of the other methods succeed here either (see full result plots in Appendix A). The results suggest the following considerations for the FIT test:

1. The test can fail for small sample sizes. It becomes reliable once the sample size reaches

2. The test is fast even for a large number of dimensions and sample sizes, most often returning in less than 10s.

3. Type I and Type II errors decrease consistently as the number of samples increases.

Figure 6: Typical FIT results: all the p-values plotted against sample size on log-log scales. Here we show as example the results on PNL data, one plot for each dimensionality setting. Black dots correspond to dataset versions where . To achieve low Type I error, they should be close to 1. White dots correspond to dataset versions where . To achieve low Type II error, they should be close to 0. (Circled black dots result when the independent and dependent version of the dataset returned the same p-value.) The values are clipped at . The horizontal lines mark the p-values of 0.05 and 0.01, respectively. Lack of datapoints indicates that the method ran out of time (e.g. for high sample sizes in the highest dimensionality datasets). The last plot shows the runtime for different sample sizes on each of the datasets. We stopped the methods after 100s.

4.2 CHSIC

CHSIC performs well on low-dimensional, low-sample-size and low-difficulty data, but has trouble detecting conditional dependence in the more difficult cases. Fig. 7 illustrates this using the LINGAUSS datasets. Overall, our results suggests that CHSIC:

1. Achieves very good Type I and Type II error levels for most datasets settings, but only for the lowest difficulty settings.

2. Fails to return meaningful results for sample sizes over 1000 and/or large difficulty data.

Figure 7: Typical CHSIC results (shown here on the LINGAUSS datasets). Black dots correspond to dataset versions where . To achieve low Type I error, they should be close to 1. White dots correspond to dataset versions where . To achieve low Type II error, they should be close to 0. The values are clipped at . Lack of datapoints indicates that the method ran out of time. Note the comparison to the runtime plot in Fig. 6.

4.3 KCIT

KCIT does well in most cases when it returns, but is often too slow to approach the harder datasets – as illustrated in Fig. 8, again using the LINGAUSS datasets.

1. KCIT achieves great Type I and Type II errors for low dimensionality / low difficulty datasets.

2. Type I error of the method increases as dataset difficulty increases, while Type II error stays low, across all the settings.

3. The method’s speed only allows it to tackle low-dimensionality data with sample size smaller than 1000.

Figure 8: Typical KCIT results (shown here on the LINGAUSS datasets). Black dots correspond to dataset versions where . To achieve low Type I error, they should be close to 1. White dots correspond to dataset versions where . To achieve low Type II error, they should be close to 0. The values are clipped at . Lack of datapoints indicates that the method ran out of time – almost always when number of samples was over 1000.

4.4 RCIT

In almost all cases (with the only exception of 1-dimensional LINGAUSS and 1-dimensional PNL data), RCIT returns very low p-values as the number of samples increases, thus yielding high Type I errors. Fig. 9 illustrates this behavior.

Figure 9: Typical RCIT results (shown here on the CHAOS datasets). Black dots correspond to dataset versions where . To achieve low Type I error, they should be close to 1. White dots correspond to dataset versions where . To achieve low Type II error, they should be close to 0. The values are clipped at . Lack of datapoints indicates that the method ran out of time.

4.5 KCIPT

KCIPT fails on a majority of the datasets, and is too slow to process any but the least-dimensional and lowest-sample cases – as shown in Fig. 10.

Figure 10: Typical KCIPT results (shown here on LINGAUSS data). Black dots correspond to dataset versions where . To achieve low Type I error, they should be close to 1. White dots correspond to dataset versions where . To achieve low Type II error, they should be close to 0. The values are clipped at . Lack of datapoints indicates that the method ran out of time.

4.6 CCI

CCI fails on most datasets. In addition, instead of outputting a p-value it outputs a binary decision: either 1 (indicating the data is conditionally independent) or 0 (for conditionally dependent data). Unfortunately, in our evaluation it has a tendency to output 0 for most of the continuous datasets, and 1 for the HYBRID dataset (as shown in Fig. 11).

Figure 11: Typical CCI results (shown here on HYBRID data). Black dots correspond to dataset versions where . To achieve low Type I error, they should be close to 1. White dots correspond to dataset versions where . To achieve low Type II error, they should be close to 0. The values are clipped at . Lack of datapoints indicates that the method ran out of time. CCI returns only a binary result of 0 (for dependence) or 1 (for independence).

5 Discussion

We showed that FIT is a fast and accurate conditional independence test. Our evaluation covered a broad range of dataset sizes, dimensionalities and complexities. In this evaluation, FIT performs as well as or better than the alternatives. In addition, it is the only algorithm that applies successfully to large sample sizes. In this section, we discuss alternatives to our evaluation strategy, and clarify some aspects of FIT.

5.1 Alternative evaluation metrics

Previous work used a variety of methods to evaluate the performance of CITs. Fukumizu et al. [2008] and Zhang et al. [2012] report Type I and Type II errors on a number of datasets. This is similar to our condensed presentation in Fig. 5. However, their evaluation considered a fixed, small sample size of 200 and 400 datapoints (as the evaluated methods are too slow to attack a larger dataset). In contrast, we allow each algorithm to process as many datapoints as possible with a fixed time limit.

As Doran et al. [2014] point out however, Type I and Type II errors can bias the results. Some algorithms might perform particularly well at one significance level, and some applications might require non-standard significance levels. Thus, Doran et al. [2014] proposed to condense the performance of a CIT across all significance levels to two statistics. First, their “area under the power curve” (AUPC) is the area under the cumulative distribution function of the p-values the algorithm returns on a conditionally dependent dataset. Second, they perform a test on the distribution of p-values: their “KS-test p-value” is the p-value of the Kolmogorov-Smirnoff test [Stephens, 1974] for the null hypothesis that the p-values of a CIT are uniformly distributed on a conditionally independent dataset. A good test has low AUPC and a non-low KS-test p-value.

The problem with both of these approaches is that they do not consider how the behavior of a CIT changes as the number of samples grows. To us, an important requirement is that for any type of data, the test is “correct” in the limit of an infinite number of samples. Looking at a limited range of sample sizes can be deceptive. Consider Fig. 31 with dim=128. RCIT seems to be doing very well for sample sizes between 100 and 1000. However, it fails for any larger number of samples. In reality, both for the dependent and independent version of the dataset RCIT’s p-values follow a decaying trend as the number of samples grows, which is not a correct behavior.

Evaluating a CIT on a fixed sample size is deceptive in a similar way to evaluating it at a fixed significance level. The best (although energy-consuming) way we found to understand the behavior of a CIT and compare it with other algorithms is to carefully examine how its p-value changes with a changing sample size on a variety of datasets. Figures 12- 35 show these p-value plots. We encourage any reader interested in a thorough comparison of the CIT methods to take the time to examine these plots.

5.2 Why use a Decision Tree Regression?

Algorithm 1 uses decision tree regression (DTR) as its machine learning component. Other regression algorithms could serve as a substitute for DTR. There are two crucial requirements for such algorithms (apart from good predictive power): 1) computational efficiency and 2) ability to handle structured multi-dimensional inputs. Initially, we implemented FIT using a neural network as the back-end, efficiently implemented in Tensorflow [Abadi et al., 2016], and trained on a GPU (Titan Xp). However, given that we want FIT to apply to both small and large, high- and low-dimensional data, the hyperparameter choice for the neural network regressor becomes a complex problem – for example Jaderberg et al. [2017] propose an simple and effective hyperparameter search method that nevertheless assumes having access to a cluster of a hundred GPU machines as a reasonable setting. DTRs are fast and require us to choose the value of only one hyperparameter.

5.3 Failure cases

Our Fast (Conditional) Independence Test assumes that if , then the MSE obtained after regressing Y on X, Z is smaller than that obtained after regressing Y on Z only. This assumption doesn’t always hold. Take the system where X, Y, Z are one-dimensional variables and N(Z, X). In other words, Y is a Gaussian with mean Z and variance X. Then, whether X is known or not, the true regression line is simply Y = Z, and the MSE does not change. However, Y clearly depends on X, even when Z is given. This represents just one possible failure mode; we do not attempt to characterize all the cases in which nonlinear regression does not benefit from a conditionally dependent regressor and simply assume them away.

5.4 Which test should I use?

Our experimental results suggest that in general, FIT is the best choice. The only situation where an alternative should be considered is low-dimensional data where it i impossible to obtain more than 1000 samples. In such cases, CHSIC or KCIT are likely to outperform FIT.

RCIT and CCI do not work well in our evaluation. KCIPT is too slow to justify using it over FIT, KCIT or CHSIC.

References

Martín Abadi, Ashish Agarwal, Paul Barham, Eugene Brevdo, Zhifeng Chen, Craig Citro, Greg S. Corrado, Andy Davis, Jeffrey Dean, and Matthieu Devin. Tensorflow: Large-scale machine learning on heterogeneous distributed systems. arXiv preprint arXiv:1603.04467, 2016.

Leo Breiman, Jerome Friedman, Charles J. Stone, and Richard A. Olshen. Classification and regression trees. CRC press, 1984.

Thomas M Cover and Joy A Thomas. Elements of information theory. John Wiley & Sons, 2012.

Alexander Philip Dawid. Conditional independence in statistical theory. Journal of the Royal Statistical Society. Series B (Methodological), pages 1–31, 1979.

Gary Doran, Krikamol Muandet, Kun Zhang, and Bernhard Schölkopf. A Permutation-Based Kernel Conditional Independence Test. In UAI, pages 132–141, 2014.

Bradley Efron. Bootstrap methods: another look at the jackknife. Springer, 1992.

Ronald Aylmer Fisher. Statistical methods for research workers. In Breakthroughs in Statistics, pages 66–70. Springer, 1925.

Kenji Fukumizu, Arthur Gretton, Xiaohai Sun, and Bernhard Schölkopf. Kernel Measures of Conditional Dependence. In Twenty-First Annual Conference on Neural Information Processing Systems (NIPS 2007), pages 489–496, 2008.

Arthur Gretton, Karsten M. Borgwardt, Malte J. Rasch, Bernhard Schölkopf, and Alexander Smola. A kernel two-sample test. Journal of Machine Learning Research, 13(Mar):723–773, 2012.

Patrik Hoyer, Dominik Janzing, Joris Mooij, Jonas Peters, and Bernhard Schölkopf. Nonlinear causal discovery with additive noise models. In Advances in Neural Information Processing Systems, pages 689–696, 2009.

Max Jaderberg, Valentin Dalibard, Simon Osindero, Wojciech M. Czarnecki, Jeff Donahue, Ali Razavi, Oriol Vinyals, Tim Green, Iain Dunning, Karen Simonyan, Chrisantha Fernando, and Koray Kavukcuoglu. Population Based Training of Neural Networks. arXiv:1711.09846 [cs], November 2017. URL http://arxiv.org/abs/1711.09846. arXiv: 1711.09846.

Daphne Koller and Mehran Sahami. Toward optimal feature selection. In International Conference on Machine Learning (ICML), 1996.

Judea Pearl. Causality. Cambridge university press, 2000.

Joseph D. Ramsey. A scalable conditional independence test for nonlinear, non-gaussian data. arXiv preprint arXiv:1401.5031, 2014.

Bernhard Schölkopf and Alexander Smola. Learning with kernels: support vector machines, regularization, optimization, and beyond. MIT press, 2002.

Peter Spirtes, Clark N. Glymour, and Richard Scheines. Causation, prediction, and search. MIT press, 2000.

Michael A. Stephens. EDF statistics for goodness of fit and some comparisons. Journal of the American statistical Association, 69(347):730–737, 1974.

Eric V. Strobl, Kun Zhang, and Shyam Visweswaran. Approximate Kernel-based Conditional Independence Tests for Fast Non-Parametric Causal Discovery. arXiv preprint arXiv:1702.03877, 2017.

Student. The probable error of a mean. Biometrika, pages 1–25, 1908.

Kun Zhang, Jonas Peters, Dominik Janzing, and Bernhard Schölkopf. Kernel-based conditional independence test and application in causal discovery. arXiv preprint arXiv:1202.3775, 2012.

A Full Results

Below, we present the full results plots for all the datasets and methods. Fig. 5 compresses (not losslessly!) information presented here. In addition, the full results support the conclusions of Sec. 5. Figures 12- 35 show plots of p-values obtained by each algorithm as the dataset and number of samples varies.

A.1 LINGAUSS datasets

A.2 CHAOS datasets

A.3 HYBRID datasets

A.4 PNL datasets

designed for accessibility and to further open science