Causal Discovery with Cascade Nonlinear Additive Noise Model

1 Introduction

Understanding causal relationships is a fundamental problem in various disciplines of science, and causal direction identification is an essential issue in causality studies. It is well known that using randomized experiments to identify causal influences usually encounters unethical or substantial expense issues. Fortunately, inferring causal relations from pure observations, also known as causal discovery from observational data, has demonstrated its power in empirical studies and has been a focus in causality research.

Various methods have been proposed to infer the causal direction, by exploring properly constrained forms of functional causal models (FCMs). A functional causal model represents the effect Y as a function of its direct causes X and independent noise, i.e., Y � f�X; ϵ�,X ⎡ ϵ. Without constraints on f, then for any two variables one can always express one of them as a function of the other and independent noise [Zhang et al., 2015]. However, it is interesting to note that with properly constrained FCMs, the causal direction between X and Y is identifiable because the independence condition between the noise and cause holds only for the true causal direction and is violated for the wrong direction. Such FCMs include the Linear, Non-Gaussian, Acyclic Model (LiNGAM) [Shimizu et al., 2006], in which Y � a�X � ϵ with linear coefficients a, the nonlinear additive noise model (ANM) [Hoyer et al., 2009], in which Y � f�X��ϵ, and the post-nonlinear (PNL) causal model [Zhang and Hyv¨arinen, 2009], which also considers possible nonlinear sensor or measurement distortion f2 in the causal process: Y � f2�f1�X� � ϵ�. It has been shown that in the generic case, for data generated by the above FCMs, the reverse direction will not admit the same FCM class with independent noise. One can then find causal direction by estimating the FCM followed by testing for independence between the hypothetical cause and estimated noise [Hoyer et al., 2009; Zhang and Hyv¨arinen, 2009].

In reality, we can usually record only a subset of all variable which are causally related. If some variable is the direct cause of only one measured variable and is not measured, it is considered as part of the omitted factors, or noise. If a hidden variable is a direct cause of two measured variables, it is a confounder, and causal discovery in the presence of confounders is challenging, although there exist some methods with asymptotic correctness guarantees, such as the FCI algorithm [Spirtes et al., 2000]. In this paper, we are concerned with unmeasured intermediate causal variables. Suppose X1 � X2 � X3, with X2 unmeasured, and that each direct causal influence can be represented by a FCM in a certain class. If the direct causal relations are linear with additive noise, then the causal influence X1 � X3 still follows a linear model with additive noise. However, if each direct causal influence follows the ANM, the causal influence X1 � X3 does not necessarily follow the same model class. Fig. 1 gives an illustration of this phenomenon of "non-transitivity of nonlinear causal model classes," in which X2 � 2 tanh�5X1��N2, and X3 � �X2♥2�3 �N3, with X1, N2, and N3 mutually independent and following the uniform distribution between �0.5 and 0.5. As seen from the heterogeneity of the noise in X3 relative to X1, given in Fig. 1(c), the causal influence from X1 to X3 clearly does not admit a nonlinear model with additive noise. Hence, even for the correct causal direction, which is from X1 to X3, the independent noise condition is violated, and existing methods for causal direction determination by checking whether regression residual is independent from the hypothetical cause may fail. The PNL is more general than the additive noise model - in this example, if N3 is zero, then X1 �X3 will follow this model. However, the PNL model class is also non-transitive.

This paper deals with such indirect, nonlinear causal relations, which seem to be ubiquitous in practice. Finding causal direction for such causal relations has recently been posed as an open problem [Spirtes and Zhang, 2016]. In particular, we aim to find the causal direction between X and Y that are generated according the process given in Fig. 2, in which there might be a number of unmeasured intermediate causal variables Zi in between and each direct causal influence, e.g., the influence from Z1 and Z2 on Z3, follows the ANM. We name the causal model from X and Y given in Fig. 2 a Cascade Additive Noise Model (CANM). We note that the considered problem is different from causal discovery in the presence of confounders, for which there have been a number of studies, including the FCI [Spirtes et al., 2000], RFCI [Colombo et al., 2012], M3B [Yu et al., 2018] algorithms, and methods relying on stronger assumptions [Janzing et al., 2009; Zhang et al., 2010]. [Kocaoglu et al., 2018] propose an algorithm to search for the latent variable along the path X and Y but they only consider discrete random variables.

To the best of our knowledge, this is a first study as to finding causal direction between indirectly and nonlinearly related variables. The considered causal model can be seen as a cascade of processes, each of which follows the ANM, and the intermediate variables are unmeasured. Intuitively, the independence between the noise and cause is still helpful in finding causal direction-the wrong direction will not follow the independence noise condition in the generic case, allowing us to correctly identify causal direction. This will be supported by our theoretical studies and empirical results in subsequent sections.

2 Cascade Additive Noise Model

Without loss of generality, let X be the cause of effect Y (X � Y ), with unmeasured intermediate variables Zi between them, as shown in Fig. 2. We further assume there is no confounder in the mechanism and the data generation follows the nonlinear additive noise assumption. Then, such an indirect causal mechanism can be formalized by the CANM in the following definition.

Definition 1. A CANM for cause X and effect Y is that there exists a sequence of unmeasured intermediate variables between X and Y such that no variable in the latter is the cause of the former one:

where X, Ni, and ϵ are mutually independent, T denotes depth of the chain, and Zpa�t�,Zpa�y� denote parents of the Zt and y, respectively. To ensure the cascade structure, the causal relations among Zi are recursive. Let f � ρf1, f2, ..., fTξ and N � ρN1,N2, ...,NTξ denote a set of nonlinear functions and the corresponding additive noises at each depth in the chain, respectively. Naturally, here the direct cause and the noises are independent from each other.

We are given a set of data D � ρx�i� , y�i� ξmi�1. Let θ be the parameters of the causal mechanism. Combing all the independence relations of CANM, we can derive its marginal log-likelihood as follows:

(2)

Eq. (2) first decomposes the joint likelihood based on the Markov condition [Spirtes et al., 2000], then applies the independence property between the cause and the noise in the second equality, i.e., p�Zt∂Zpa�t�� � p�Nt � Zt � ft�Zpa�t��∂Zpa�t�� Zpa�t�⎡Nt ���������� p�Nt � Zt � ft�Zpa�t���. At the same time, we replace dz with dn and rewrite function fT�1�Zpa�y�� as f�X,N�, because the last unobserved direct cause ZT Λ Zpa�t� contains all the information of the noise N and cause X relative to Y .

In the above derivation, we used the transformation from X and noises to Y . The property of the transformation helps study identifiability and find a practical solution. In light of the independence property of the noises, below we propose a variational approach to approximating the marginal loglikelihood as well as identifying the causal direction.

2.1 Variational Solution of CANM

The variational solution to estimation of CANM consists of two steps. First, we take advantage of the independence property in CANM to replace the latent variable Z with N. Second, we find an amortized inference distribution qφ�N∂X, Y � with respect to the parameter φ to approximate the intractable posterior pθ�N∂X, Y � and jointly optimize a variational lower bound of the marginal log-likelihood. Note that, different from the vanilla VAE, Y can be seen as a function of X and N and, as a result, N is a function of X and Y and we need to recover N from both X and Y . According to Eq. (2), the (log) marginal likelihood, as the sum over of the marginal likelihoods over individual data points:

where L�θ, φ; x�i� , y�i� � be the lower bound at data point �x�i� , y�i� �, resulting from approximating an intractable posterior pθ�n∂x�i� , y�i� � by qφ�n∂x�i� , y�i� �. Under the framework of CANM, the lower bound of the total marginal likelihood can be further estimated as follows:

The details of derivation can be found in Supplementary A in [Cai et al., 2019]. As shown in Eq. (3), the lower bound L is tight at KL�qφ�n∂x�i� , y�i� �⏐pθ�n∂x�i� , y�i� �� � 0. That is, when qφ�n∂x�i� , y�i� � � pθ�n∂x�i� , y�i� �, the marginal loglikelihood is equal to the lower bound. Below we will maximize the variational lower bound.

Here, we assume the distributions of noise N can be factorized as pθ�N� � 4Tt�1 pθ�Nt�. Note that if N is an empty set, the above lower bound is equivalent to the log-likelihood of the standard additive noise model.

2.2 Variational Auto-encoder

The design of the variational auto-encoder (VAE) generally follows the typical configuration in [Kingma and Welling, 2014]. We denote qφ�n∂x�i� , y�i� � as encoder and pθ�y�i� ∂n, x�i� � as decoder, using a multilayer perceptron (MLP) as an universal approximator for this two functions.

In the encoder phase, the noises of CANM are inferred by an encoder network with a reparameterization trick. That is, reparameterize the random variable n � qφ�n∂x, y� with a differentiable transformation hφ�x, y, u� such that n � hφ�x, y, u� with u � p�u�. Then we can estimate the expectation of lower bound En�qφ�n∂x,y� �p�ϵ�i� � y�i� � f �x�i� , n; θ��� using Monte Carlo with the reparameterization trick.

In the decoder phase, we estimate the ϵ�i� by subtracting sample y�i� from the reconstruction value of decoder f �x�i� , hφ�x�i� , y�i� , u�l� �; θ�, where u�l� � p�u�. Then, alternatively, processes the encoder and decoder phases, we can optimize the lower bound until it converges.

Fig. 3 shows a toy example of the structure of the CANM variational auto-encoder with qφ�n∂x�i� , y�i� � � N�n; µφ �x�i� , y�i� � , σφ �x�i� , y�i� � I�, where µφ and σφ are deterministic function with parameter φ. In the encoder phase, we encode the samples into the noises using a reparameterization trick n�l� � µφ �x�i� , y�i� � � σφ �x�i� , y�i� �u�l� where u�l� � N�0, 1�. In the decoder phase, the sample y�i� is reconstructed by the decoder y←�i� � f �x�i� , n�l� ; θ�.

2.3 Practical Algorithm

Finally, we propose a general principle that makes use of the VAE to estimate the marginal log-likelihood as well as identify the causal direction.

Algorithm 1 Inferring causal direction with CANM

Input: Data samples ρ�x�i� , y�i� �ξmi�1.

Output: The causal direction.

• 1: Split the data into training and test sets;
• 2: Choose the best number of latent variables by optimizing the variational lower bound (Eq. (3)) on the training set and evaluating the performance on the test set;
• 3: Optimize the lower bound in both directions with the best number of latent variables on the full dataset, obtaining LX�Y and LY�X (see Eq. (4)), respectively.
• 4: if LX�Y % LY�X � δ, where δ is a pre-asigned small positive number, then,
• 5: Infer X � Y
• 6: else if LX�Y ∃ LY�X � δ, then
• 7: Infer Y �X
• 8: else
• 9: Non-identifiable
• 10: end if

Algorithm 1 consists of two phases; the first is model selection, selecting the best number of latent noises, and the second is to identify the causal direction. In phase 1, by splitting the data into training and testing sets, the best number of noises is selected based on the performance on the test set (Line 1-2). In phase 2, we use the number of the latent noises determined in phase 1 to optimize the variational lower bound on the full dataset and then identify causal direction according to the likelihood for both directions (Line 3-10).

3 Identifiability

In this section, we investigate whether there exist any CANMs whose generated data also admit a CANM in the reverse (anti-causal) direction. In the following theorem, we propose a way to derive the noise distribution for the reverse direction p�ˆ ϵ� by making use of the theory of Fourier transform [Bracewell and Bracewell, 1986]. The causal direction is unidentifiable according to the CANM if ˆ ϵ is independent from Y and ˆN (i.e., the marginal likelihoods for both directions are equal).

Theorem 1. Let X � Y follow the cascade additive noise model, while there exists a backward model following the same form, i.e.

then the noise distribution of the reverse direction pˆ ϵ must be

(6)

where f, g denote the function implied by the cascade process.

Roughly speaking, regardless of the linear case, Theorem 1 implies that the noise distribution in the reverse direction is generally coherent with y. To ensure such noise is independent from Y , one strict condition must holds, i.e., ˆ ϵ should be independent from Y in the sense that ⎯y1, y2, ∆ e2πiˆ ϵ�ν ∆∆ p�x�p�n�pϵ�y1�f�x,n��e�2πix�ν dndx p�y1� ∆ p�ˆn�e�2πig�y1,ˆn��ν dˆn dν � ∆ e2πiˆ ϵ�ν ∆∆ p�x�p�n�pϵ�y2�f�x,n��e�2πix�ν dndx p�y2� ∆ p�ˆn�e�2πig�y2,ˆn��ν dˆn dν. However, in general, it seems that such a condition holds only in restrictive cases. Therefore, in most cases, after the latent noise is recovered, we can identify the causal direction by using the independence property for �X,N, ϵ�.

To further illustrate the implication of Theorem 1, we provide two special cases in the following corollaries. In Corollary 1, we show that CANM is unidentifiable if the generation process is linear Gaussian. In Corollary 2, we show the connection with ANM when there is no unmeasured intermediate variables, and shows a generic choices of f , pX�x�, and pϵ�ϵ� for the identification of the model. Those two special cases are consistent with the previous results.

Corollary 1. Assume that CANM is linear Gaussian, i.e.,

where X, N, ϵ � N�0, 1�, then their exist a backward CANM

where ˆN, ˆ ϵ � N�0, 1� and ˆ ϵ is independent of Y and ˆN.

Corollary 2. Suppose that there is no unmeasured intermediate noises in CANM, if the solution of Eq. (6) exists, then the triple �f, pX, pϵ� must satisfy the differential equation from ANM [Hoyer et al., 2009, Theorem 1] for all x, y with ν ←← �y � f�x��f ← �x� ϕ 0:

where ν � log pϵ, ξ � log pX

Proof. See Supplementary D in [Cai et al., 2019].

4 Experiments

4.1 Synthetic Data

In this section, we design three experiments with known ground truth, with the depth � ρ0, 1, 2, 3, 4, 5ξ, sample size � ρ250, 500, 1000, 2000, 3000, 4000, 5000, 6000ξ, and with different sample sizes for some fix structure. The default setting is marked in bold. All experimental results are averaged over 1000 randomly generated causal pairs by the cascade additive noise model. Code for CANM is available online 1 .

To make the synthetic data general enough, in each depth, we randomly generate an additive noise model and then stack it together to obtain the cascade additive noise model. In detail, the cause �X� is sampled from a random Gaussian Mixture model of 3 components p�xi∂θ� �

1 https://github.com/DMIRLAB-Group/CANM

<3k�1 πkN�xi∂µk, σk� where µk � N�0, 1�, σk � Super � Gaussian. For each layer xt � ft�xt�1� � nt where nt � N�0, 1� and ft is generated from a cubic spline interpolation using a 6-dimensional grid from min�xt�1� to max�xt�1� as input with respect to 6 random generated points as knots for the interpolation; the generated points are sampled from N�0, 1� and the number of knots is used to control nonlinearity of the function. Such generative process follows the instrument given in [Prestwich et al., 2016].

The following four algorithms are taken as baseline methods: ANM [Hoyer et al., 2009], CAM [B¨uhlmann et al., 2014], IGCI [Janzing et al., 2012], and LiNGAM [Shimizu et al., 2006]. We also improve the implementation for ANM by using the XGBoost [Chen and Guestrin, 2016] for regression and the Hilbert-Schmidt independence criterion (HSIC) [Gretton et al., 2008] as the independence test. Therefore, ANM can be evaluated in two ways. First, we compare the HSIC statistic to determine the direction and second, we select the best significance level (p � 0.01) range from 0.01 to 1 to determine the causal direction. At the same time, the best parameter setting of IGCI is chosen. For the other baseline methods, we use the parameter settings in their original papers. The implementation and the parameter settings of LiNGAM and CAM are based on the CompareCausalNetworks packages in R [Heinze-Deml et al., 2018].

Sensitivity to Depth

Fig. 4 shows the accuracy with different depths in 3000 samples. First, when the depth equals 0 (the original ANM), all CANM, ANM, and CAM achieve a high accuracy. Note that CANM still has a similar performance comparing with ANM even though CANM assume that there might exist unmeasured intermediate variables, which demonstrates the robustness of our method. Second, as the depth increases, the accuracy of CANM is stable and around 90% accuracy with a slight decrease, while the performance of the rest methods decreases rapidly as the depth grows. In particular, the ANM with the significance level of 0.01 gives almost random decisions when the cascade structure exists.

Sensitivity to Sample Size

Fig. 5 shows the accuracy with different sample sizes while the depth is fixed at 3. The result shows that even in the small sample size, CANM still outperforms the other methods. As the sample size increases, the accuracy of CANM grows faster than the other methods. Thus, large samples are beneficial to CANM, because of the variational auto-encoder framework employed in CANM. A similar result also can be observed in ANM and CAM while the other methods are less sensitive to the sample size due to the model restriction.

Sensitivity to Sample Size in a Fixed Structure

Fig. 6 shows the accuracy with different numbers of samples while we use a fixed causal mechanism, which was randomly generated with depth=3. When the sample size is small, the variance of the likelihood is large; however, the asymmetry in the causal direction is still clear. As the sample size increases, the variance of the likelihood decreases and the accuracy increases, which implies the effectiveness and robustness of CANM as the sample size grows.

4.2 Real World Data

Electricity Consumption Dataset

The electricity consumption dataset [Prestwich et al., 2016] has 9504-hour measurements from the energy industry, containing the hour of data, outside temperature and the electricity load on the power station. The causal mechanism among the three variables are hour of day � temperature and temperature � electricity load. The first pair is generally caused by the heating of sunlight and the second pair is base on the fact that the usage of heating or air condition depends upon temperature. We are interested to know whether we can identify the hour of day �X� is the cause of the electricity load �Y � and what intermediate variable will be inferred via CANM.

In general, we successfully identify the correct causal direction with average score LX�Y � �2.62 % LY�X � �2.67 while ANM fails on this pair (the p-value � 0 on both directions). The prediction of electricity is given in Fig. 7. It is interesting to note that there might exist more than one unmeasured variable, e.g., season, causing a different electricity load at the same hour of day. Such unmeasured variables are successfully captured by CANM as the prediction separating into both upper and lower parts. Furthermore, the intermediate variable inferred by our method has rather high correlation (ρ � �0.35) with the temperature as shown in Fig. 8, which means that CANM not only recovers the information of the season but also the information of the temperature.

Stock Market Dataset

The stock market dataset is collected by T¨ubingen causal effect benchmark (https://webdav.tuebingen.mpg.de/ cause-effect/) as pairs 66-67. It contains the stock return of Hutchison, Cheung Kong and Sun Hung Kai with the causal relationship: Hutchison � Cheung Kong and Cheung Kong � Sun Hung Kai. The reason for the first pair is that Cheung Kong owns about 50% of Hutchison. For the second pair, Sun Hung Kai Prop., a typical stock in the Hang Seng Property subindex, is believed to depend on the major stock Cheung Kong. Similarly to the previous experiment, we are interested to know whether we can identify the Hutchison �X� is the cause of the Sun Hung Kai �Y �.

Since these three stocks form a causal chain that Hutchison � Cheung Kong � Sun Hung Kai, using CANM, we successfully identify the indirect causal direction with average score LX�Y � �2.49 % LY �X � �2.51 while ANM fails on this pair (the p-value � 0.006 ∃ 0.05 on the causal direction and p-value � 0.29 % 0.05 on the reverse direction). Fig. 9 shows the prediction of the stock return of the Sun Hung Kai. We also find that the fitted intermediate variable has a high correction (ρ � �0.54) with the stock return of Cheung Kong as shown in Fig. 10.

5 Conclusion

In this paper, we proposed the cascade nonlinear additive noise model, as an extension of the nonlinear additive noise model, to represent indirect causal influences, which result from unmeasured intermediate causal variables. We have demonstrated that, the independence between the noise and cause is still generally helpful to determine causal direction between two variables, as long as the cascade additive noise process holds. We propose to estimate the model as well as the intermediate causal variables with the variational autoencoder framework, and the produced likelihood indicates the asymmetry between cause and effect. As supported by our theoretical and empirical results, the proposed approach provides an effective method for causal direction determination from data generated by nonlinear, indirect causal relations.

Acknowledgments

This research was supported in part by NSFC-Guangdong Joint Found (U1501254), Natural Science Foundation of China (61876043), Natural Science Foundation of Guangdong (2014A030306004, 2014A030308008), Guangdong High-level Personnel of Special Support Program (2015TQ01X140) and Pearl River S&T Nova Program of Guangzhou (201610010101). KZ would like to acknowledge the support by National Institutes of Health (NIH) under Contract No. NIH-1R01EB022858-01, FAINR01EB022858, NIH-1R01LM012087, NIH-5U54HG008540-02, and FAINU54HG008540, by the United States Air Force under Contract No. FA8650-17-C-7715, and by National Science Foundation (NSF) EAGER Grant No. IIS-1829681. The NIH, the U.S. Air Force, and the NSF are not responsible for the views reported here. We appreciate the comments from anonymous reviewers, which greatly helped to improve the paper.