Causal Inference on Discrete Data using Additive Noise Models
Jonas Peters, Dominik Janzing, Bernhard Schölkopf
Introduction
Inferring causal relations by analyzing statistical dependences among observed random variables is a challenging task if no controlled randomized experiments are available. So-called constraint-based approaches to causal discovery (Pearl, 2000; Spirtes et al., 1993) select among all directed acyclic graphs (DAGs) those that satisfy the Markov condition and the faithfulness assumption, i.e., those for which the observed independences are imposed by the structure rather than being a result of specific choices of parameters of the Bayesian network. These approaches are unable to distinguish among causal DAGs that impose the same independences. In particular, it is impossible to distinguish between and . More recently, several methods have been suggested that use not only conditional independences, but also more sophisticated properties of the joint distribution. For simplicity, we explain the ideas for the two variable setting since this case is particularly challenging. Kano & Shimizu (2003) use models
where is a linear function and is additive noise that is independent of the hypothetical cause . This is an example for an additive noise model from to . Apart from trivial cases, can only admit such a model from to and from to in the bivariate Gaussian case. Hoyer et al. (2009) generalize the method to non-linear functions and showed that generic models of this form generate joint distributions that do not admit such an additive noise model from to . Zhang & Hyvarinen (2009) augment the model by applying a non-linear function to the rhs of eq. (1) and still obtain identifiability for generic cases. Peters et al. (2009) use independent linear additive noise models in order to detect whether a sample of a time series has been reversed. Their positive results further support this way of causal reasoning. All these proposals, however, were only designed for real-valued variables and .
For discrete variables, Sun et al. (2008) propose a method to measure the complexity of causal models via a Hilbert space norm of the logarithm of conditional densities and prefer models that induce smaller norms. Sun et al. (2006) fit joint distributions of cause and effect with conditional densities whose logarithm is a second order polynomial (up to the log-partition function) and show that this often makes causal directions identifiable when some or all variables are discrete. For discrete variables, several Bayesian approaches (Heckerman et al., 1999) are also applicable, but the construction of good priors are challenging and often the latter are designed such that Markov equivalent DAGs still remain indistinguishable.
The main idea of the causal inference method we propose goes as follows: If such an additive noise model exists in one direction but not in the other, we prefer the former based on Occam’s Razor and infer it to be the causal direction.
Such a procedure is only sensible if there are only few instances, in which there is an additive noise models in both directions. If, for example, all additive noise models from to also allow an additive noise model from to , we could not draw any causal conclusions at all. We will show that reversible cases are very rare and thereby answer this theoretical question.
For a practical causal inference method we have to test whether the data admits an additive noise model and thus have to perform a discrete regression. But since in the discrete case regularization of the regression function is not necessary, in principle we would have to check all possible functions and test whether they result in independent residuals. This is highly intractable, of course, and we therefore propose an efficient heuristic procedure that proved to work very well in practice.
In section 2 we extend the concept of additive noise models to discrete random variables and show the corresponding identifiability results for generic cases in section 3. In section 4 we introduce an efficient algorithm for causal inference on finite data, for which we show experimental results in section 5. We conclude in section 6.
Additive Noise Models for Discrete Variables
As it has been proposed for the continuous case by Shimizu et al. (2006); Hoyer et al. (2009); Zhang & Hyvarinen (2009) we assume the following causal principle to hold throughout the remainder of this article:
Causal Inference Principle (for discrete random variables) Whenever satisfies an additive noise model with respect to and not vice versa then is a cause for , and we write .
Janzing & Steudel (2009) give further theoretical support for this principle using the concept of Kolmogorov complexity and Peters et al. (2009) use this way of reasoning for detecting the arrow of time.
Note that whenever there is no additive noise model in any direction (which may well happen) we do not draw any causal conclusions and other causal inference methods should be tried.
Furthermore we require for all . This does not restrict the model class, but is due to a freedom we have in choosing and : If Y=f(X)+N,\,N\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}X, then we can always construct a new function , such that Y=f_{j}(X)+N_{j},\,N_{j}\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}X by choosing and .
Such an additive noise model is called reversible if there is also an additive noise model from to , i.e. if it satisfies an additive noise model in both directions.
2 Cyclic Constraint
We can extend additive noise models to random variables which inherit a cyclic structure and therefore take values in a periodic domain. Random variables are usually defined as measurable maps from a probability space into the real numbers. We thus make the following definition
Let be a probability space. A function is called an -cyclic random variable if . All other concepts of probability theory (like distributions and expectations) can be constructed analogous to the well-known case, in which takes values .
3 Relations
The following two remarks are essential in order to understand the relationship between integer and cyclic constraints:
(1) The difference between these two models manifests in the target domain. If we consider an ANM from to it is important whether we put integer or cyclic constraints on (and thus on ). It does not make a difference, however, whether we consider the regressor to be cyclic (with a cycle larger than the support of ) or not. The independence constraint remains the same.
(2) In the finite case additive noise models with cyclic constraints are more general than the ones with integer constraints: Assume there is an ANM , where all variables are taken to be non-cyclic and takes values between and , say. Then we still have an ANM if we regard to be -cyclic because remains independent of . It is possible, however, that N\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}\diagup\hskip 3.98337ptX, but N\mod(l-k+1)\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}X (see Example 2 in Section 3.2.1).
Identifiability
Whether an additive noise model is allowed depends on the form of the joint distribution . Let be the subset of the set of all possible joint distributions that allow an additive noise model from to in the “forward direction”, whereas allows an additive noise model in the backward direction from to . Some trivial examples like and immediately show that there are joint distributions allowing additive noise models in both directions, meaning (see Figure 1). But how large is this intersection? Our method would not be useful if we find out that and are almost the same sets. Then most additive noise models can be fit in either both directions or in none. For additive noise models with integer constraints and with cyclic constraints we identify the intersection and show that it is indeed a very small set. If we are unlucky and the data generating process we consider happens to be in , our method does not give wrong results, but answers “I do not know the answer”. In all other situations the method identifies the correct direction given that we observe enough data. The proofs are provided in the appendix.
Figure 3 shows a (rather non-generic) example that allows an ANM in both directions if we choose for and for .
The s are shifted versions of each other
The probability distributions on the s are shifted and scaled versions of each other with the same shift constant as above: For
As for the other theorems of this section the proof is provided in the appendix. Its main point is based on the asymmetric effects of the “corners” of the joint distribution. In order to allow for an infinite support of (or ) the proof generalizes the concept of “corners”.
1.2 X𝑋X and Y𝑌Y have infinite support
Consider an additive noise model where both and have infinite support. We distinguish between two cases
Note that the first case is again a complete characterization of all instances of a joint distribution, an ANM in both directions is conform with. The second case does not yield a complete characterization, but shows how restricted we are for a given function and noise in order to choose a distribution that yields a reversible ANM.
2 Cyclic Constraint
Note that the model is reversible if and only if there is a function , such that
First, we give three examples of additive noise models that are not identifiable. This restricts the class of situations in which identifiability can be expected. Figure 4 shows instances of all three examples.
If for a bijective and affine and uniformly distributed , then the model is reversible.
2.2 Identifiability Results
Motivated by the counter examples we now make the assumptions that is not constant (Example 1(i)), is not uniformly distributed (Example 1(ii)) and that not both and are uniformly distributed (Example 2). Without proof we state the conjecture that this is already enough to ensure identifiability meaning an ANM can only hold in one direction.
Assume and are random variables that are not uniformly distributed and non-degenerate (that is they do not only take one value). Assume further that we have an additive noise model Y=f(X)+N,\,N\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}X with non-constant . Then there is no additive noise model from to .
3 Mixed Constraints
With the results developed in the last two section we can cover even models with mixed constraints (for the precise conditions of “usually” see Section 3.2):
Practical Method for Causal Inference
Based on our theoretical findings in Section 3 we propose the following method for causal inference (see Hoyer et al. (2009) for the continuous case):
Given: iid data of the joint distribution .
Regularization of the regression function is at least in principle not necessary in the discrete case. Since we may observe many different values of for one specific value there is no risk in overfitting. This introduces further difficulties compared to continuous regression since in principle we now have to try all possible functions from to and compare the corresponding values of the loss function.
(2) Even if all values of the true function are one of the considered values, the problem of checking all possible functions is not tractable: If and there are possible functions (the amount of particles in the universe is estimated to be ). We thus propose the following heuristic but very efficient procedure (the experimental results will show that it works very reliably in practice):
2 Independence Test and Dependence Measure
Assume we are given joint iid samples of the discrete variables and and we want to test whether and are independent. In our implementation we only use Person’s test (e.g. E. L. Lehmann (2005)), which is most commonly used. It computes the difference between observed frequencies and expected frequencies in the contingency table. The test statistic is known to converge towards a distribution, which is taken as an approximation even in the finite sample case. For very few samples, however, this approximation and therefore the test will usually fail. It has been suggested (e.g. RProject (2009)) that instead of a test, Fisher’s exact test (E. L. Lehmann, 2005) could be used if not more than of the expected counts are larger than 5 (“Cochran’s condition”).
For a dependence measure we simply use the -value of the independence test or the test statistic if the -value is too small (in a computer system the -value is sometimes regarded to be zero).
Experiments
We first investigate the performance of our method on synthetic data sets. Therefore we simulate data from additive noise models and check whether the method is able to rediscover the true model. We showed in Section 3 that only very few examples allow a reversible ANM.
Data sets 1a and 1b support these theoretical results. We simulate a large amount of data (1000–2000 data points) from many randomly chosen models. All models that allow an ANM in both directions are instances of our examples from above (without exception).
Data sets 2a and 2b show how well our method performs for small data size and models that are close to non-identifiable.
Data set 3a investigates empirically the run-time performance of our regression method and compares it with a brute-force search.
Data set 3b shows that the method does not favor one direction if the supports of and are of different size.
The results given in Table 1 show that the methods works well on almost all simulated data sets. The algorithm outputs “bad fit in both directions” in roughly of all cases, which corresponds to the chosen test level. The model is non-identifiable only in very few cases, which are shown in Table 1. All of these cases are instances of the counter examples from above. This experiment further supports our proposition that the model is identifiable in the generic case.
Data set 2a (close to non-identifiable). For this data set we sample from the model with
Depending on the parameter we sample from
For each value of the parameter ranging between we use 100 different samples, each of which has the size 400.
In Theorem 2 we proved that the ANM is reversible if and only if . Figure 5 shows that the algorithm identifies the correct direction for . Again, the test level of introduces indecisiveness of roughly the same size, which can be seen for . The number of such cases can be reduced by decreasing (but would lead to some wrongly accepted backward models, too).
Each time we simulated a uniformly distributed with values between and for . For each noise/regressor distribution we simulated 100 data sets.
For and , for example, there are possible functions in total and functions with positive empirical support. Our method only checked functions before termination. The full results are shown in Figure 6.
1.2 Cyclic Constraints
The results given in Table 2 show that the method works well on almost all simulated data sets. The algorithm outputs “bad fit in both directions” in roughly of all cases, which corresponds to the chosen test level. The model is non-identifiable only in very few cases, which are shown in Table 3. All of these cases are instances of the counter examples from above. This experiment further supports our theoretical result that the model is identifiable in the generic case.
Example 1 and the fact that (see proof of Proposition 7) show that the model is not identifiable if and only if the noise distribution is uniform, i.e. if and only if . The further is away from , the more the noise differs from a uniform distribution and the easier it should be for our method to detect the true direction. For each value of the parameter we use 100 different samples, each of which has size 200. Figure 8 shows the results.
For and (indicated by the arrows in Figure 8) we further investigate the dependence on the data size. Clearly, results in a model that is still very close to non-identifiability and thus we need more data to perform well (see Figure 8). Note that non-identifiable models lead to very few, but not to wrong decisions.
2 Real Data.
Data set 4 (abalone). We also applied our method to the abalone data set (Nash et al., 1994) from the UCI Machine Learning Repository (Asuncion & Newman, 2007). We tested the sex of the abalone (male (1), female (2) or infant (0)) against length , diameter and height , which are all measured in mm, and have and different values, respectively. Compared to the number of samples (up to 4177) we treat this data as being discrete. Because we do not have information about the underlying continuous length we have to assume that the data structure has not been destroyed by the user-specific discretization. We regard , and as being the ground truth, since the sex is probably causing the size of the abalone, but not vice versa.
Clearly, the variables do not have a cyclic structure. For the sex variable, however, the most natural model would be a structureless set which is contained in the cyclic constraints; for comparison we try both models for . Our method with integer constraints is able to identify all 3 directions correctly. Since may be cyclic we also try to fit an ANM from to with the cyclic constraints. Again, these models are rejected (see Table 5 and Figure 9). We used and the first 1000 samples of the data set.
As mentioned earlier it is not surprising that we would accept an ANM from to even if we put cyclic constraints on (which are certainly violated for this data set). We would obtain the following -values: and .
For this data set the method proposed by (Sun et al., 2006) returns a slightly higher likelihood for the true causal directions than for the false directions, but this difference is so small, that the algorithm does not consider it to be significant.
The abalone data set also shows that working with -values requires some carefulness. For synthetic data sets that we simulate from one fixed model the -values do not depend on the data size. In real world data, however, this often is the case. If the data generating process does not exactly follow the model we assume, but is reasonable close to it, we get good fits for moderate data sizes. Only including more and more data reveals the small difference between process and model, which therefore leads to small -values. Figure 10 shows how the -values vary if we include the first data points of the abalone data set (in total: 4177). One can see that although the -values for the correct direction decrease they are clearly preferable to the -values of the wrong direction. This is a well-known problem in applied statistics that also has to be considered using our method.
For data points both directions are rejected (, ). Figure 11 shows, however, that again the are decreasing much slower than thus using other criteria than simple -values we still may prefer this direction and propose it as the true one.
Conclusions and Future Work
We proposed a method that is able to infer the cause-effect relationship between two discrete random variables. We proved that for generic choices the direction of a discrete additive noise model is identifiable in the population case and we developed an efficient algorithm that is able to infer the causal relationship between two variables for a finite amount of data. Since it is known that fails for small data sizes, changing the independence test for those cases may lead to an even higher performance of the algorithm.
Our method can be generalized in two directions: (1) handling more than two variables is straightforward from a practical point of view (although one may have to introduce regularization to make the regression computationally feasible) and (2) it should be investigated how our procedure can be applied to the case, where one variable is discrete and the other continuous. Corresponding identifiability results remain to be shown.
In future work additive noise models should be tested on more real world data sets in order to support (or disprove) additive noise models as a principle in causal inference. Furthermore we hope that more fundamental and general principles for identifying causal relationships will be developed that cover additive noise models as a special case. Nevertheless we regard our work as an important step towards understanding the difference between cause and effect.
Now we consider the other case, namely that has finite support. Then we define to be disjoint sets, such that is constant on each of them: . This time, it does not matter which of these sets is called . Since
The rest of the proof is valid for both cases (either or has finite support): Consider for any . According to the assumption that an additive noise model holds we have
with . Thus (including ).
For (which implies ) we have
In order to show that we have an reversible ANM we define the function as follows
B Proof of Theorem 3
and , else.
Start with any arbitrary and define
This implies and .
We have that either or : Otherwise and satisfy
Because of this contradicts the existence of an backward additive noise model. Wlog we therefore assume . Then we even have ,
(Otherwise we can use the same argument as above with and .) Define further
Since , but , such a value must exist. Again, we can define in the same way as above.
over the sequence . But since we computed the boxes in a deterministic way, the same satisfies
Note that a corresponding equation with the same constant holds for the opposite direction. This leads to a contradiction, since there is no probability distribution for with infinite support that can fulfill this condition (no matter if is greater, equal or smaller than ).
This direction is proved in exactly the same way as in Theorem 2.
C Proof of Theorem 5
regarding : Each distribution has to have the same support (up to an additive shift) and thus the same number of elements with probability greater than :
For we now consider 3 different cases and show necessary conditions for reversibility each.
Assume Y=f(X)+N,\,N\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}X for bijective and . If the model is reversible with a bijective , then and are uniformly distributed.
Proof. Since is bijective we have that . It follows from (2)
Since is bijective it follows that . This holds for all and thus and are uniformly distributed.
is not injective. Assume . From (2) it follows that
which imply equality constraints on . To determine the number of constraints we define a function that maps the arguments of the numerator to those of the denominator
Assume Y=f(X)+N,\,N\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}X and . Assume further that the model is reversible with a non-injective .
introduces a functional relationship between and .
is not injective. Assume . In a slight abuse of notation we write
Assume Y=f(X)+N,\,N\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 2.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 2.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 2.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 2.0mu{\scriptscriptstyle\perp}}}X, is not injective and . Assume further that the model is reversible for a function .
The rest follows analogously to the proof of Proposition 7.
Note that these three cases are sufficient since and injective implies and and bijective.