Analysis of cause-effect inference by comparing regression errors

Problem setting and notation

In this work, we use the framework of structural causal models (Pearl, 2009) with the goal of correctly identifying cause and effect variables of given observations from X and Y. As illustrated in Fig. 2, this can be described by the problem of identifying whether the causal DAG of X → Y or X←Y is true. Throughout this paper, a capital letter denotes a random variable and a lowercase letter denotes values attained by the random variable. Variables X and Y are assumed to be real-valued and to have a joint probability density (with respect to the Lebesgue measure), denoted by p_X,Y. By slightly abusing terminology, we will not further distinguish between a distribution and its density since the Lebesgue measure as a reference is implicitly understood. The notations p_X, p_Y, and p_Y|X are used for the corresponding marginal and conditional densities, respectively. The derivative of a function f is denoted by f′.

General idea

As mentioned before, the general idea of our approach is to simply compare the MSE of regressing Y on X and the MSE of regressing X on Y. If we denote cause and effect by C, E ∈ {X, Y}, respectively, our approach explicitly reads as follows. Let ϕ denote the function that minimizes the expected least squares error when predicting E from C, which implies that ϕ is given by the conditional expectation ϕ(c) = 𝔼[E|c]. Likewise, let ψ be the minimizer of the least squares error for predicting C from E, that is, ψ(e) = 𝔼[C|e]. Then we will postulate assumptions that imply (2)${\displaystyle \mathbb{E}\left[{\left(E-\varphi \left(C\right)\right)}^2\right]\le \mathbb{E}\left[{\left(C-\psi \left(E\right)\right)}^2\right]}$ in the regime of almost deterministic relations. This conclusion certainly relies on some kind of scaling convention. For our theoretical results we will assume that both X and Y attain values between 0 and 1. However, in some applications, we will also scale X and Y to unit variance to deal with unbounded variables. Equation (2) can be rewritten in terms of conditional variance as ${\displaystyle \mathbb{E}\left[Var\left[E|C\right]\right]\le \mathbb{E}\left[Var\left[C|E\right]\right].}$

Assumptions

First, recall that we assume throughout the paper that either X is the cause of Y or vice versa in an unconfounded sense, i.e., there is no common cause. Therefore, the general structural equation is defined as (3)${\displaystyle E=\zeta \left(C,\tilde{N}\right),}$ where $C⫫\tilde{N}$. For our analysis, we first define a function ϕ to be the conditional expectation of the effect given the cause, i.e., ${\displaystyle \varphi \left(c\right):=\mathbb{E}\left[E|c\right]}$ and, accordingly, we define a noise variable N as the residual (4)${\displaystyle N:=E-\varphi \left(C\right).}$ Note that (4) implies that 𝔼[N|c] = 0. The function ϕ is further specified below. Then, to study the limit of an almost deterministic relation in a mathematically precise way, we consider a family of effect variables E_α by (5)${\displaystyle E_{\alpha }:=\varphi \left(C\right)+\alpha N,}$ where α ∈ ℝ⁺ is a parameter controlling the noise level and N is a noise variable that has some (upper bounded) joint density p_N,C with C. Note that N here does not need to be statistically independent of C (in contrast to ANMs), which allows the noise to be non-additive. Therefore, (5) does not, a priori, restrict the set of possible causal relations, because for any pair (C, E) one can always define the noise N as (4) and thus obtain E_α=1 = E for any arbitrary function ϕ³.

For this work, we make use of the following assumptions:

1. Invertible function: ϕ is a strictly monotonically increasing two times differentiable function ϕ:[0, 1] → [0, 1]. For simplicity, we assume that ϕ is monotonically increasing with ϕ(0) = 0 and ϕ(1) = 1 (similar results for monotonically decreasing functions follow by reflection E → 1 − E). We also assume that ϕ^−1′ is bounded.

2. Compact supports: The distribution of C has compact support. Without loss of generality, we assume that 0 and 1 are, respectively, the smallest and the largest values attained by C. We further assume that the distribution of N has compact support and that there exist values n₊ > 0 > n₋ such that for any c, [n₋, n₊] is the smallest interval containing the support of p_N|c. This ensures that we know [αn₋, 1 + αn₊] is the smallest interval containing the support of p_Eα. Then the shifted and rescaled variable (6)${\displaystyle {\tilde{E}}_{\alpha }:=\frac{1}{1+\alpha n_{+}-\alpha n_{-}}\left(E_{\alpha }-\alpha n_{-}\right)}$ attains 0 and 1 as minimum and maximum values and thus is equally scaled as C.

3. Unit noise variance: The expected conditional noise variance is 𝔼[Var[N|C]] = 1 without loss of generality, seeing that we can scale the noise arbitrary by the parameter α and we are only interested in the limit α → 0.

4. Independence postulate: While the above assumptions are just technical, we now state the essential assumption that generates the asymmetry between cause and effect. To this end, we consider the unit interval [0, 1] as probability space with uniform distribution as probability measure. The functions c↦ϕ′(c) and c↦Var[N|c]p_C(c) define random variables on this space, which we postulate to be uncorrelated, formally stated as (7)${\displaystyle Cov\left[{\varphi }^{\prime },Var\left[N|c\right]p_C\right]=0.}$ More explicitly, (7) reads: (8)${\displaystyle {\int }_0^1{\varphi }^{\prime }\left(c\right)Var\left[N|c\right]p_C\left(c\right)dc-{\int }_0^1{\varphi }^{\prime }\left(c\right)dc{\int }_0^{1}Var\left[N|c\right]p_C\left(c\right)dc=0.}$

The justification of (7) is not obvious at all. For the special case where the conditional variance Var[N|c] is a constant in c (e.g., for ANMs), (7) reduces to (9)${\displaystyle Cov\left[{\varphi }^{\prime },p_C\right]=0,}$ which is an independence condition for deterministic relations stated in Schölkopf et al. (2013). Conditions of similar type as (9) have been discussed and justified in Janzing et al. (2012). They are based on the idea that ϕ contains no information about p_C. This, in turn, relies on the idea that the conditional p_E|C contains no information about p_C.

To discuss the justification of (8), observe first that it cannot be justified as stating some kind of ‘independence’ between p_C and p_E|C. To see this, note that (8) states an uncorrelatedness of the two functions c↦ϕ′(c) and c↦Var[N|c]p_C(c). While ϕ′ depends only on the conditional p_E|C and not on p_C, the second function depends on both p_C|E and p_E, since Var[N|c] is a property of p_E|C. Nevertheless, to justify (8) we assume that the function ϕ represents a law of nature that persists when p_C and N change due to changing background conditions. From this perspective, it becomes unlikely that they are related to the background condition at hand. This idea follows the general spirit of ‘modularity and autonomy’ in structural equation modeling, that some structural equations may remain unchanged when other parts of a system change (see Chapter 2 in Peters, Janzing & Schölkopf (2017) for a literature review)⁴. To further justify (7), one could think of a scenario where someone changes ϕ independently of p_N,C, which then results in vanishing correlations. Typically, this assumption would be violated if ϕ is adjusted to p_N,C or vice versa. This could happen due to an intelligent design by, for instance, first observing p_N,C and then defining ϕ or due to a long adaption process in nature (see Janzing et al. (2012) for further discussions of possible violations in a deterministic setting).

A simple implication of (8) reads (10)${\displaystyle {\int }_0^1{\varphi }^{\prime }\left(c\right)Var\left[N|c\right]p_C\left(c\right)dc=1,}$ due to ${\int }_0^1{\varphi }^{\prime }\left(c\right)dc=1$ and ${\int }_0^{1}Var\left[N|c\right]p_C\left(c\right)dc=\mathbb{E}\left[Var\left[N|C\right]\right]=1$.

In the following, the term independence postulate is used to refer to the aforementioned postulate and the term independence to a statistical independence, which should generally become clear from the context.

Error asymmetry theorem

For our main theorem, we first need an important lemma:

Lemma 1 (Limit of variance ratio) Let the assumptions 1–3 in ‘Assumptions’ hold. Then the following limit holds: (11)${\displaystyle \underset{\alpha \to 0}{lim}\frac{\mathbb{E}\left[Var\left[C|{\tilde{E}}_{\alpha }\right]\right]}{\mathbb{E}\left[Var\left[{\tilde{E}}_{\alpha }|C\right]\right]}={\int }_0^1\frac{1}{{\varphi }^{\prime }{\left(c\right)}^2}Var\left[N|c\right]p_C\left(c\right)dc}$

Proof: We first give some reminders of the definition of the conditional variance and some properties. For two random variables Z and Q the conditional variance of Z, given q is defined by ${\displaystyle Var\left[Z|q\right]:=\mathbb{E}\left[{\left(Z-\mathbb{E}\left[Z|q\right]\right)}^2|q\right],}$ while Var[Z|Q] is the random variable attaining the value Var[Z|q] when Q attains the value q. Its expectation reads ${\displaystyle \mathbb{E}\left[Var\left[Z|Q\right]\right]:=\int Var\left[Z|q\right]p_Q\left(q\right)dq.}$ For any a ∈ ℝ, we have ${\displaystyle Var\left[\frac{Z}{a}\left|\right.q\right]=\frac{Var\left[Z|q\right]}{a^2}.}$ For any function h, we have ${\displaystyle Var\left[h\left(Q\right)+Z|q\right]=Var\left[Z|q\right],}$ which implies Var[h(Q)|q] = 0. Moreover, we have ${\displaystyle Var\left[Z|h\left(q\right)\right]=Var\left[Z|q\right],}$ if h is invertible.

To begin the main part of the proof, we first observe (12)${\displaystyle \mathbb{E}\left[Var\left[E_{\alpha }|C\right]\right]=\mathbb{E}\left[Var\left[\varphi \left(C\right)+\alpha N|C\right]\right]={\alpha }^2\underset{=1\text{(Assumpt. 3))}}{\underbrace{\mathbb{E}\left[Var\left[N|C\right]\right]}}={\alpha }^2.}$ Moreover, one easily verifies that (13)${\displaystyle \underset{\alpha \to 0}{lim}\frac{\mathbb{E}\left[Var\left[C|{\tilde{E}}_{\alpha }\right]\right]}{\mathbb{E}\left[Var\left[{\tilde{E}}_{\alpha }|C\right]\right]}=\underset{\alpha \to 0}{lim}\frac{\mathbb{E}\left[Var\left[C|E_{\alpha }\right]\right]}{\mathbb{E}\left[Var\left[E_{\alpha }|C\right]\right]},}$ due to (6) provided that these limits exist. Combining Eqs. (12) and (13) yields (14)${\displaystyle \underset{\alpha \to 0}{lim}\frac{\mathbb{E}\left[Var\left[C|{\tilde{E}}_{\alpha }\right]\right]}{\mathbb{E}\left[Var\left[{\tilde{E}}_{\alpha }|C\right]\right]}=\underset{\alpha \to 0}{lim}\frac{\mathbb{E}\left[Var\left[C|E_{\alpha }\right]\right]}{{\alpha }^2}=\underset{\alpha \to 0}{lim}\mathbb{E}\left[Var\left[\frac{C}{\alpha }\left|\right.E_{\alpha }\right]\right].}$

Now, we can rewrite (14) as (15)${\displaystyle \underset{\alpha \to 0}{lim}\mathbb{E}\left[Var\left[\frac{C}{\alpha }\left|\right.E_{\alpha }\right]\right]=\underset{\alpha \to 0}{lim}\mathbb{E}\left[Var\left[\frac{{\varphi }^{-1}\left(E_{\alpha }-\alpha N\right)}{\alpha }\left|\right.E_{\alpha }\right]\right]=\underset{\alpha \to 0}{lim}{\int }_{\varphi \left(0\right)+\alpha n_{-}}^{\varphi \left(1\right)+\alpha n_{+}}Var\left[\frac{{\varphi }^{-1}\left(e-\alpha N\right)}{\alpha }\left|\right.e\right]p_{E_{\alpha }}\left(e\right)de=\underset{\alpha \to 0}{lim}{\int }_{\varphi \left(0\right)}^{\varphi \left(1\right)}Var\left[\frac{{\varphi }^{-1}\left(e-\alpha N\right)}{\alpha }\left|\right.e\right]p_{E_{\alpha }}\left(e\right)de.}$ In the latter step, αn₊ and −αn₋ vanishes in the limit seeing that the function ${\displaystyle e\mapsto Var\left[{\varphi }^{-1}\left(e-\alpha N\right)∕\alpha \left|\right.e\right]p_{E_{\alpha }}\left(e\right)}$ is uniformly bounded in α. This is firstly, because ϕ⁻¹ attains only values in [0, 1], and hence the variance is bounded by 1. Secondly, p_Eα(e) is uniformly bounded due to ${\displaystyle p_{E_{\alpha }}\left(e\right)={\int }_{n_{-}}^{n_{+}}{p}_{\varphi \left(C\right),N}\left(e-\alpha n,n\right)dn={\int }_{n_{-}}^{n_{+}}{p}_{C,N}\left({\varphi }^{-1}\left(e-\alpha n\right),n\right){\varphi }^{{-1}^{\prime }}\left(e-\alpha n\right)dn\le \parallel {\varphi }^{{-1}^{\prime }}{\parallel }_{\infty }\parallel p_{C,N}{\parallel }_{\infty }\left(n_{+}-n_{-}\right).}$ Accordingly, the bounded convergence theorem states ${\displaystyle \underset{\alpha \to 0}{lim}{\int }_{\varphi \left(0\right)}^{\varphi \left(1\right)}Var\left[\frac{{\varphi }^{-1}\left(e-\alpha N\right)}{\alpha }\left|\right.e\right]p_{E_{\alpha }}\left(e\right)de={\int }_{\varphi \left(0\right)}^{\varphi \left(1\right)}\underset{\alpha \to 0}{lim}\left(Var\left[\frac{{\varphi }^{-1}\left(e-\alpha N\right)}{\alpha }\left|\right.e\right]p_{E_{\alpha }}\left(e\right)\right)de.}$

To compute the limit of ${\displaystyle Var\left[\frac{{\varphi }^{-1}\left(e-\alpha N\right)}{\alpha }\left|\right.e\right],}$ we use Taylor’s theorem to obtain (16)${\displaystyle {\varphi }^{-1}\left(e-\alpha n\right)={\varphi }^{-1}\left(e\right)-\alpha n{{\varphi }^{-1}}^{\prime }\left(e\right)-\frac{{\alpha }^2{n}^2{{\varphi }^{-1}}^{\prime \prime }\left(E_2\left(n,e\right)\right)}{2},}$ where E₂(n, e) is a real number in the interval (e − αn, e). Since (16) holds for every $n\in \left[-\frac{e}{\alpha },\frac{1-e}{\alpha }\right]$ (note that ϕ and ϕ⁻¹ are bijections of [0, 1], thus e − αn lies in [0, 1]) it also holds for the random variable N if E₂(n, e) is replaced with the random variable E₂(N, e) (here, we have implicitly assumed that the map n↦e₂(n, e) is measurable). Therefore, we see that (17)${\displaystyle \underset{\alpha \to 0}{lim}Var\left[\frac{{\varphi }^{-1}\left(e-\alpha N\right)}{\alpha }\left|\right.e\right]=\underset{\alpha \to 0}{lim}Var\left[-N{{\varphi }^{-1}}^{\prime }\left(e\right)-\frac{\alpha N^2{{\varphi }^{-1}}^{\prime \prime }\left(E_2\left(N,e\right)\right)}{2}\left|\right.e\right]={\varphi }^{{-1}^{\prime }}{\left(e\right)}^{2}Var\left[N|e\right].}$ Moreover, we have (18)${\displaystyle \underset{\alpha \to 0}{lim}{p}_{E_{\alpha }}\left(e\right)=p_{E_0}\left(e\right).}$ Inserting Eqs. (18) and (17) into (15) yields ${\displaystyle \underset{\alpha \to 0}{lim}\mathbb{E}\left[Var\left[\frac{C}{\alpha }\left|\right.E_0\right]\right]={\int }_{\varphi \left(0\right)}^{\varphi \left(1\right)}{\varphi }^{{-1}^{\prime }}{\left(e\right)}^{2}Var\left[N|e\right]p_{E_0}\left(e\right)de={\int }_0^1{\varphi }^{{-1}^{\prime }}{\left(\varphi \left(c\right)\right)}^{2}Var\left[N|\varphi \left(c\right)\right]p_C\left(c\right)dc={\int }_0^1\frac{1}{{\varphi }^{\prime }{\left(c\right)}^2}Var\left[N|c\right]p_C\left(c\right)dc,}$ where the second equality is a variable substitution using the deterministic relation E₀ = ϕ(C) (which implies p_E0(ϕ(c)) = p_C(c)∕ϕ′(c) or, equivalently, the simple symbolic equation p_E0(e)de = p_C(c)dc). This completes the proof due to (14). □

While the formal proof is a bit technical, the intuition behind this idea is quite simple: just think of the scatter plot of an almost deterministic relation as a thick line. Then Var[E_α|c] and Var[C|E_α = ϕ(c)] are roughly the squared widths of the line at some point (c, ϕ(c)) measured in vertical and horizontal direction, respectively. The quotient of the widths in vertical and horizontal direction is then given by the slope. This intuition yields the following approximate identity for small α: (19)${\displaystyle Var\left[C|{\tilde{E}}_{\alpha }=\varphi \left(c\right)\right]\approx \frac{1}{{\left({\varphi }^{\prime }\left(c\right)\right)}^2}Var\left[{\tilde{E}}_{\alpha }|C=c\right]={\alpha }^2\frac{1}{{\left({\varphi }^{\prime }\left(c\right)\right)}^2}Var\left[N|c\right].}$ Taking the expectation of (19) over C and recalling that Assumption 3 implies $\mathbb{E}\left[Var\left[{\tilde{E}}_{\alpha }|C\right]\right]={\alpha }^2\mathbb{E}\left[Var\left[N|C\right]\right]={\alpha }^2$ already yields (11).

With the help of Lemma 1, we can now formulate the core theorem of this paper:

Theorem 1 (Error Asymmetry) Let the assumptions 1–4 in ‘Assumptions’ hold. Then the following limit always holds ${\displaystyle \underset{\alpha \to 0}{lim}\frac{\mathbb{E}\left[Var\left[C|{\tilde{E}}_{\alpha }\right]\right]}{\mathbb{E}\left[Var\left[{\tilde{E}}_{\alpha }|C\right]\right]}\ge 1,}$ with equality only if the function stated in Assumption 1 is linear.

Proof: We first recall that Lemma 1 states ${\displaystyle \underset{\alpha \to 0}{lim}\frac{\mathbb{E}\left[Var\left[C|{\tilde{E}}_{\alpha }\right]\right]}{\mathbb{E}\left[Var\left[{\tilde{E}}_{\alpha }|C\right]\right]}={\int }_0^1\frac{1}{{\varphi }^{\prime }{\left(c\right)}^2}Var\left[N|c\right]p_C\left(c\right)dc.}$ We then have ${\displaystyle {\int }_0^1\frac{1}{{\varphi }^{\prime }{\left(c\right)}^2}Var\left[N|c\right]p_C\left(c\right)dc}$ ${\displaystyle ={\int }_0^1\frac{1}{{\varphi }^{\prime }{\left(c\right)}^2}Var\left[N|c\right]p_C\left(c\right)dc\cdot \underset{=1\text{(Assumpt. 3)}}{\underbrace{{\int }_0^{1}Var\left[N|c\right]p_C\left(c\right)dc}}}$ ${\displaystyle ={\int }_0^1{\sqrt{{\left(\frac{1}{{\varphi }^{\prime }\left(c\right)}\right)}^{2}Var\left[N|c\right]}}^2{p}_C\left(c\right)dc\cdot {\int }_0^1{\sqrt{Var\left[N|c\right]}}^2{p}_C\left(c\right)dc}$ ${\displaystyle \ge {\left({\int }_0^1\sqrt{{\left(\frac{1}{{\varphi }^{\prime }\left(c\right)}\right)}^{2}Var\left[N|c\right]}\sqrt{Var\left[N|c\right]}{p}_C\left(c\right)dc\right)}^2}$ (20)${\displaystyle ={\left({\int }_0^1\frac{1}{{\varphi }^{\prime }\left(c\right)}Var\left[N|c\right]p_C\left(c\right)dc\right)}^2,}$ where the inequality is just the Cauchy Schwarz inequality applied to the bilinear form f, g↦∫f(c)g(c)p_C(c)dc for the space of functions f for which ∫f²(c)p_c(c)dc exists. Note that if ϕ is linear, (20) becomes 1, since ϕ′ = 1 according to Assumpt. 1. We can make a statement about (20) in a similar way by using (10) implied by the independence postulate and using Cauchy Schwarz: ${\displaystyle {\int }_0^1\frac{1}{{\varphi }^{\prime }\left(c\right)}Var\left[N|c\right]p_C\left(c\right)dc}$ ${\displaystyle ={\int }_0^1\frac{1}{{\varphi }^{\prime }\left(c\right)}Var\left[N|c\right]p_C\left(c\right)dc\cdot \underset{=1\left(10\right)}{\underbrace{{\int }_0^1{\varphi }^{\prime }\left(c\right)Var\left[N|c\right]p_C\left(c\right)dc}}}$ ${\displaystyle ={\int }_0^1{\sqrt{\frac{1}{{\varphi }^{\prime }\left(c\right)}Var\left[N|c\right]}}^2{p}_C\left(c\right)dc\cdot {\int }_0^1{\sqrt{{\varphi }^{\prime }\left(c\right)Var\left[N|c\right]}}^2{p}_C\left(c\right)dc}$ ${\displaystyle \ge {\left({\int }_0^1\sqrt{\frac{1}{{\varphi }^{\prime }\left(c\right)}Var\left[N|c\right]}\sqrt{{\varphi }^{\prime }\left(c\right)Var\left[N|c\right]}{p}_C\left(c\right)dc\right)}^2}$ (21)${\displaystyle ={\left(\underset{=1\text{(Assumpt. 3)}}{\underbrace{{\int }_0^{1}Var\left[N|c\right]p_C\left(c\right)dc}}\right)}^2=1.}$ Combining Eqs. (20) and (21) with Lemma 1 completes the proof. □

Remark

Theorem 1 states that the inequality holds for all values of α smaller than a certain finite threshold. Whether this threshold is small or whether the asymmetry with respect to regression errors already occurs for large noise cannot be concluded from the theoretical insights. Presumably, this depends on the features of ϕ, p_C, p_N|C in a complicated way. However, the experiments in ‘Experiments’ suggest that the asymmetry often appears even for realistic noise levels.

If the function ϕ is non-invertible, there is an information loss in anticausal direction, since multiple possible values can be assigned to the same input. Therefore, we can expect that the error difference becomes even higher in these cases, which is supported by the experiments in ‘Simulated cause–effect pairs with strong dependent noise’.

Causal inference methods for comparison

In the following, we briefly discuss and compare the causal inference methods which we used for the evaluations.

Artificial data

For experiments with artificial data, we performed evaluations with simulated cause–effect pairs generated for a benchmark comparison in Mooij et al. (2016). Further, we generated additional pairs with linear, nonlinear invertible and nonlinear non-invertible functions where input and noise are strongly dependent.

Simulated benchmark cause–effect pairs

The work of Mooij et al. (2016) provides simulated cause–effect pairs with randomly generated distributions and functional relationships under different conditions. As pointed out by Mooij et al. (2016), the scatter plots of these simulated data look similar to those of real-world data. We took the same data sets as used in Mooij et al. (2016) and extend the reported results with an evaluation with SLOPE, CURE, LiNGAM, RECI and further provide results in the preprocessed data.

The data sets are categorized into four different categories:

• SIM: Pairs without confounders. The results are shown in Figs. 3A–3B

• SIM-c: A similar scenario as SIM, but with one additional confounder. The results are shown in Figs. 3C–3D

• SIM-ln: Pairs with low noise level without confounder. The results are shown in Figs. 4A-4B

• SIM-G: Pairs where the distributions of C and N are almost Gaussian without confounder. The results are shown in Figs. 4C–4D.

The general form of the data generation process without confounder but with measurement noise⁵ is ${\displaystyle C^{\prime }\sim p_C,N\sim p_N}$ ${\displaystyle N_C\sim \mathcal{N}\left(0,{\sigma }_C\right),N_E\sim \mathcal{N}\left(0,{\sigma }_E\right)}$ ${\displaystyle C=C^{\prime }+N_C}$ ${\displaystyle E=f_E\left(C^{\prime },N\right)+N_E}$ and with confounder ${\displaystyle C^{\prime }\sim p_C,N\sim p_N,Z\sim p_Z}$ ${\displaystyle C^{\prime \prime }=f_C\left(C^{\prime },Z\right)}$ ${\displaystyle N_C\sim \mathcal{N}\left(0,{\sigma }_C\right),N_E\sim \mathcal{N}\left(0,{\sigma }_E\right)}$ ${\displaystyle C=C^{\prime \prime }+N_C}$ ${\displaystyle E=f_E\left(C^{\prime \prime },Z,N\right)+N_E,}$ where N_C, N_E represent independent observational Gaussian noise and the variances σ_C and σ_E are chosen randomly with respect to the setting. Note that only N_E is Gaussian, while the regression residual is non-Gaussian due to the nonlinearity of f_E and non-Gaussianity of N, Z. Thus, the noise in SIM, SIM-c and SIM-G is non-Gaussian. More details can be found in Appendix C of Mooij et al. (2016).

Generally, ANM performs the best in all data sets. However, the difference between ANM and RECI, depending on the regression model, becomes smaller in the preprocessed data where isolated points were removed. According to the observed performances, removing these points seems to often improve the accuracy of RECI, but decrease the accuracy of ANM. In case of the preprocessed SIM-G data sets, the accuracy of RECI is decreased seeing that in these nearly Gaussian data the removal of low-density points leads to an underestimation of the noise distribution.

In all data sets, except for SIM-G, RECI always outperforms SLOPE, IGCI, CURE and LiNGAM if a simple logistic or polynomial function is utilized for the regression. However, in the SIM-G data set, our approach performs comparably poor, which could be explained by the violation of the assumption of a compact support. In this nearly Gaussian setting, IGCI performs the best with a Gaussian reference measure. However, we also evaluated RECI with standardized data in the SIM-G data sets, which is equivalent to a Gaussian reference measure for IGCI. A summary of the results can be found in Figures 1(c)–1(d) in the supplements and more detailed results in Table 4 and Table 5 in the supplements. These results are significantly better than normalizing the data in this case. However, although our theorem only justifies a normalization, a different scaling, such as standardization, might be a reasonable alternative.

Even though Theorem 1 does not exclude cases of a high noise level, it makes a clear statement about low noise level. Therefore, as expected, RECI performs the best in SIM-ln, where the noise level is low. In all cases, LiNGAM performs very poorly due to the violations of its core assumptions. Surprisingly, although PNL is a generalization of ANM, we found that PNL performs generally worse than ANM, but better than SLOPE, CURE and LiNGAM.

ANM and RECI require a least-squares regression, but ANM additionally depends on an independence test, which can have a high computational cost and a big influence on the performance. Therefore, even though RECI does not outperform ANM, it represents a competitive alternative with a lower computational cost, depending on the regression model and MSE estimation. Also, seeing that RECI explicitly allows both cases, a dependency and an independency between C and N and ANM only the latter, it can be expected that RECI performs significantly better than ANM in cases where the dependency between C and N is strong. This is evaluated in ‘Simulated cause–effect pairs with strong dependent noise’. In comparison with PNL, SLOPE, IGCI, LiNGAM and CURE, RECI outperforms in almost all data sets. Note that Mooij et al. (2016) performed more extensive experiments and showed more comparisons with ANM and IGCI in these data sets, where additional parameter configurations were tested. However, they reported no results for the preprocessed data.

Simulated cause–effect pairs with strong dependent noise

Since the data sets of the evaluations in ‘Simulated benchmark cause–effect pairs’ are generated by structural equations with independent noise variables, we additionally performed evaluations with artificial data sets where the input distribution and the noise distribution are strongly dependent. For this, we considered a similar data generation process as described in the work by Daniušis et al. (2010). We generated data with various cause and noise distributions, different functions and varying values for α ∈ [0, 1]. In order to ensure a dependency between C and N, we additionally introduced two unobserved source variables S₁ and S₂ that are randomly sampled from different distributions. Variables C and N then consist of a randomly weighted linear combination of S₁ and S₂. The general causal structure of these data sets is illustrated in Fig. 5. Note that S₁ and S₂ can be seen as hidden confounders affecting both C and E.

Apart from rather simple functions for ϕ, Daniušis et al. (2010) proposed to generate more general functions in the form of convex combinations of mixtures of cumulative Gaussian distribution functions ψ(C|μ_i, σ_i): ${\displaystyle s_n\left(C\right)=\sum _{i=1}^n{\beta }_i\psi \left(C|{\mu }_i,{\sigma }_i\right),}$ where β_i, μ_i ∈ [0, 1] and σ_i ∈ [0, 0.1]. For the experiments, we set n = 5 and chose the parameters of s₅(C) randomly according to the uniform distribution. Note that ψ(C|μ_i, σ_i) is always monotonically increasing and thus s₅(C) can have an arbitrary random shape while being monotonically increasing.

Cause-effect pairs were then generated in the following way ${\displaystyle w_1,w_2\sim U\left(0,1\right)}$ ${\displaystyle S_1,S_2\sim p_S}$ ${\displaystyle S_1=S_1-\mathbb{E}\left[S_1\right]}$ ${\displaystyle S_2=S_2-\mathbb{E}\left[S_2\right]}$ ${\displaystyle C^{\prime }=w_1\cdot f_1\left(S_1\right)+\left(1-w_1\right)\cdot f_2\left(S_2\right)}$ ${\displaystyle N^{\prime }=w_2\cdot f_3\left(S_1\right)+\left(1-w_2\right)\cdot f_4\left(S_2\right)}$ ${\displaystyle C=normalize\left(C^{\prime }\right)}$ ${\displaystyle N=\alpha \cdot standardize\left(N^{\prime }\right)}$ ${\displaystyle E=\varphi \left(C\right)+N,}$ where the distributions of S₁ and S₂ and the functions f₁, f₂, f₃, f₄ were chosen randomly from p_S and f in Table 1, respectively. Note that S₁ and S₂ can follow different distributions. The choice of ϕ depends on the data set, where we differentiated between three data sets:

Table 1:

All distributions p_S, functions ϕ and functions f that were used for the generation of the Linear, Non-invertible and Invertible data sets.

In case of the functions for Non-invertible, rescale (X,  − n, n) denotes a rescaling of the input data X on [ − n, n]. GM_μ,σ denotes a Gaussian mixture distribution with density $p_{G{M}_{\mu ,\sigma }}\left(c\right)=\frac{1}{2}\left(\phi \left(c|{\mu }_1,{\sigma }_1\right)+\phi \left(c|{\mu }_2,{\sigma }_2\right)\right)$ and Gaussian pdf φ(c|μ, σ).

Data set	ϕ(C)
Linear	C	pS	f(X)
Invertible	s5(C)	U(0, 1)	X
Non-invertible	rescale(C,  − 2, 2)2	$\mathcal{N}\left(0,{\sigma }^2\right)$	exp(X)
	rescale(C,  − 2, 2)4	$\mathcal{N}\left(0.5,{\sigma }^2\right)$	s5(X)
	sin(rescale(C,  − 2⋅π, 2⋅π))	$\mathcal{N}\left(1,{\sigma }^2\right)$
		GM[0.3,0.7]T,[0.1,0.1]T

DOI: 10.7717/peerjcs.169/table-1

• Linear: Only the identity function ϕ(C) = C

• Invertible: Arbitrary invertible functions ϕ(C) = s₅(C)

• Non-invertible: Functions that are not invertible on the respective domain

In total, we generated 100 data sets for each value of parameter α, which controls the amount of noise in the data. In each generated data set, we randomly chose different distributions and functions. For Linear and Non-invertible the step size of α is 0.1 and for Invertible 0.025. Here, we only performed one repetition on each data set for all algorithms. Figs. 6A–6C summarize all results and Table 6 in the supplements shows the best performing functions and parameters of the different causal inference methods. Note that we omitted experiments with CURE in these data sets due to the high computational cost.

Linear: As expected, ANM, PNL, SLOPE, IGCI and RECI perform very poorly, since they require nonlinear data. In case of RECI, Theorem 1 states an equality of the MSE if the functional relation is linear and, thus, the causal direction can not be inferred. While LiNGAM performs well for α = 0.1, and probably for smaller values of α too, the performance drops if α increases. The poor performances of LiNGAM and ANM can also be explained by the violation of its core assumption of an independence between cause and input.

Invertible: In this data set, IGCI performs quite well for small α, since all assumptions approximately hold, but the performance decreases when the noise becomes stronger and violates the assumption of a deterministic relation. In case of RECI, we made a similar observation, but it performs much better than IGCI if α < 0.5. Aside from the assumption of linear data for LiNGAM, the expected poor performance of LiNGAM and ANM can be explained by the violation of the independence assumption between C and N. In contrast to the previous results, PNL performs significantly better in this setting than ANM, although, likewise LiNGAM and ANM, the independence assumption is violated.

Non-invertible: These results seem very interesting, since it supports the argument that the error asymmetry becomes even clearer if the function is not invertible due to an information loss of regressing in anticausal direction. Here, IGCI and SLOPE perform reasonably well, while ANM and LiNGAM perform even worse than a baseline of just guessing. Comparing ANM and PNL, PNL has a clear advantage, although the overall performance is only slightly around 60% in average. The constant results of each method can be explained by the rather simple and similar choice of data generating functions.

While the cause and noise also have a dependency in the SIM-c data sets, the performance gap between ANM and RECI is vastly greater in Invertible and Non-invertible than in SIM-c due to a strong violation of the independent noise assumption. Therefore, RECI might perform better than ANM in cases with a strong dependency between cause and noise.

Real-world data

In real-world data, the true causal relationship generally requires expert knowledge and can still remain unclear in cases where randomized controlled experiments are not possible. For our evaluations, we considered the commonly used cause–effect pairs (CEP) benchmark data sets for inferring the causal direction in a bivariate setting. These benchmark data sets provided, at the time of these evaluations, 106 data sets with given cause and effect variables and can be found on https://webdav.tuebingen.mpg.de/cause-effect/. However, since we only consider a two variable problem, we omit six multivariate data sets, which leaves 100 data sets for the evaluations. These data sets consist of a wide range of different scenarios, such as data from time dependent and independent physical processes, sociological and demographic studies or biological and medical observations. An extensive analysis and discussion about the causal relationship of the first 100 data sets can be found in the work by Mooij et al. (2016). Each data set comes with a corresponding weight determined by expert knowledge. This is because several data sets are too similar to consider them as independent examples, hence they get lower weights. Therefore, the weight w_m in Eq. (22) depends on the corresponding data set. The evaluation setup is the same as for the artificial data sets, but we doubled the number of internal repetition of CURE to eight times in order to provide the same conditions as in Sgouritsa et al. (2015).

Figures 7A and 7B shows the results of the evaluations in the original and preprocessed data, respectively. In all cases, SLOPE and RECI perform significantly better than ANM, PNL, IGCI, CURE and LiNGAM. While SLOPE performs slightly better in the original data sets than RECI, RECI performs better overall in the preprocessed data⁶. Further, the performance gap even increases in the preprocessed data with removed low-density points. Surprisingly, as Table 1 in the supplements indicates, the simple shifted monomial function ax² + c performs the best, even though it is very unlikely that this function is able to capture the true function ϕ. We obtained similar observations in the artificial data sets, where the simple logistic function oftentimes performs the best.

In order to show that RECI still performs reasonably well under a different scaling, we also evaluated RECI in the real-world data set with standardized data. These results can be found summarized in Figures 1(a)–1(b) in the supplements and more detailed in Table 3 and Table 4 in the supplements. While standardizing the data improves the performance in the SIM-G data, it slightly decreases the performance in the real-world data as compared to a normalization of the data, but still performs reasonably well. This shows some robustness with respect to a different scaling.

Error ratio as rejection criterion

 
__________________________________________________________________________________________ 
Algorithm 2 Causal inference algorithm that uses Eq. (23) as rejection criterion. 
__________________________________________________________________________________________ 
  function RECI(X, Y, t)  ⊳ X and Y are the observed data and t ∈ [0,1] is the confidence 
  threshold for rejecting a decision. 
      (X,Y) ← RescaleData(X,Y) 
      f ← FitModel(X,Y)                           ⊳ Fit regression model f:X → Y 
      g ← FitModel(Y,X)                           ⊳ Fit regression model g:Y → X 
      MSEY|X ← MeanSquaredError(f,X,Y) 
      MSEX|Y  ← MeanSquaredError(g,Y,X) 
      ξ ← 1 − min(MSEX|Y,MSEY|X) 
      max(MSEX|Y,MSEY|X) 
      if ξ ≥ t then 
           if MSEY|X < MSEX|Y  then 
               return X causes Y 
           else 
               return Y causes X 
           end if 
      else 
           return No decision 
      end if 
  end function 
___________________________________________________________________________________________    

It is not clear how a confidence measure for the decision of RECI can be defined. However, since Theorem 1 states that the correct causal direction has a smaller error, we evaluated the idea of utilizing the error ratio for a confidence measure in terms of: (23)${\displaystyle \text{confidence}=1-\frac{min\left(\mathbb{E}\left[Var\left[X|Y\right]\right],\mathbb{E}\left[Var\left[Y|X\right]\right]\right)}{max\left(\mathbb{E}\left[Var\left[X|Y\right]\right],\mathbb{E}\left[Var\left[Y|X\right]\right]\right)},}$ The idea is that, the smaller the error ratio, the higher the confidence of the decision, due to the large error difference. Note that we formulated the ratio inverse to Theorem 1 in order to get a value on [0, 1]. Algorithm 1 can be modified in a straight forward manner to utilize this confidence measure. The modification is summarized in Algorithm 2⁷.

We re-evaluated the obtained results by considering only data sets where Algorithm 2 returns a decision with respect to a certain confidence threshold. Figs. 8A–8D show some examples of the performance of RECI if we use Eq. (23) to rank the confidence of the decisions. A decision rate of 20%, for instance, indicates the performance when we only force a decision on 20% of the data sets with the highest confidence. In this sense, we can get an idea of how useful the error ratio is as rejection criterion. While Figs. 8A, 8C and 8D support the intuition that the smaller the error ratio, the higher the chance of a correct decision, Fig. 8B has a contradictive behavior. In Fig. 8B, it seems that a small error ratio (big error difference) is rather an indicator for an uncertain decision (probably caused by over- or underfitting issues). Therefore, we can not generally conclude that Eq. (23) is a reliable confidence measure in all cases, but it seems to be a good heuristic approach in the majority of the cases. More plots can be found in the supplements, where Fig. 2 shows the plots in the original data, Fig. 3 in the preprocessed data and Fig. 4 in the standardized data. Note that a deeper analysis of how to define an appropriate confidence measure for all settings is beyond the scope of this paper and we rather aim to provide some insights of utilizing the error ratio for this purpose.

Run-time comparison

In order to have a brief overview and comparison of the run-times, we measured the execution durations of each method in the evaluation of the original CEP data sets. All measures were performed with a Intel Xeon E5-2695 v2 processor. Table 2 summarizes the measured run-times, where we stopped the time measurement of CURE after 9999 s. As the table indicates, IGCI is the fastest method, followed by LiNGAM and RECI. The ranking of ANM, PNL, SLOPE and RECI is not surprising; ANM and PNL need to evaluate the independence between input and residual on top of fitting a model. In case of SLOPE, multiple regression models need to be fitted depending on a certain criterion that requires to be evaluated. Therefore, by construction, RECI can be expected to be faster than ANM, PNL and SLOPE.

Discussion

Due to the greatly varying behavior and the choice of various optimization parameters, a clear rule of which regression function is the best choice for RECI remains an unclear and difficult problem. Overall, it seems that simple functions are better in capturing the error asymmetries than complex models. However, a clear explanation for this is still lacking. A possible reason for this might be that simple functions in causal direction already achieve a small error, while in anticausal direction, more complex models are required to achieve a small error. To justify speculative remarks of this kind raises deep questions about the foundations of causal inference. According to Eq. (1), it is possible to conclude that the joint distribution has a simpler description in causal direction than in anticausal direction. Seeing this, a model selection based on the regression performance and model complexity considered in a dependent manner might further improve RECI’s practical applicability. Regarding the removal of low-density points, the performance of methods that are based on the Gaussiantiy assumption, such as FN and IGCI with Gaussian reference measure, seems not to be influenced by the removal. On the other hand, the performance of HSIC, ENT, and IGCI with uniform measure is negatively affected, while the performance of LiNGAM and RECI increases. In case of RECI, this can be explained by a better estimation of the true MSE with respect to the regression function class.

Regarding the computational cost, we want to emphasize again that RECI, depending on the implementation details, can have a significantly lower computational cost than ANM, SLOPE and CURE, while providing comparable or even better results. Further, it can be easily implemented and applied.