Invariant Causal Prediction for Nonlinear Models

TL;DR

This paper develops and evaluates nonlinear, nonparametric methods for Invariant Causal Prediction (ICP), demonstrating a robust residual distribution test approach for identifying causal parents across diverse settings and applying it to fertility rate modeling.

Contribution

It introduces a new nonlinear, nonparametric ICP method called invariant residual distribution test and evaluates its robustness and effectiveness in various simulation and real-world scenarios.

Findings

The residual distribution test approach is robust across many simulations.

Performance depends on the true causal structure, especially with larger parental sets.

Application to fertility modeling confirms causal influence of child mortality rates.

Abstract

An important problem in many domains is to predict how a system will respond to interventions. This task is inherently linked to estimating the system's underlying causal structure. To this end, Invariant Causal Prediction (ICP) (Peters et al., 2016) has been proposed which learns a causal model exploiting the invariance of causal relations using data from different environments. When considering linear models, the implementation of ICP is relatively straightforward. However, the nonlinear case is more challenging due to the difficulty of performing nonparametric tests for conditional independence. In this work, we present and evaluate an array of methods for nonlinear and nonparametric versions of ICP for learning the causal parents of given target variables. We find that an approach which first fits a nonlinear model with data pooled over all environments and then tests for differences between the residual distributions across environments is quite robust across a large variety of simulation settings. We call this procedure "invariant residual distribution test". In general, we observe that the performance of all approaches is critically dependent on the true (unknown) causal structure and it becomes challenging to achieve high power if the parental set includes more than two variables. As a real-world example, we consider fertility rate modelling which is central to world population projections. We explore predicting the effect of hypothetical interventions using the accepted models from nonlinear ICP. The results reaffirm the previously observed central causal role of child mortality rates.

Tables

Table 1. Table 1: Coverage

Coverage guarantee	0.95	0.90	0.8	0.5
Coverage with nonlinear ICP	0.99	0.95	0.88	0.58
Coverage with Random Forest	0.76	0.71	0.61	0.32
Coverage with mean change	0.95	0.88	0.68	0.36

Table 2. Table 2: Overview of implemented test combinations in CondIndTests package

CIT	R function name/Method	Test	$E$	$d$
(A)	KCI() KCI (without GP support)	–	D/C	$\ge 1$
(B)	ResidualPredictionTest() Residual prediction test with – Nyström random features (RBF and polynomial kernel) – Fourier random features – fixed basis expansion	–	D/C	$\ge 1$
(C)	InvariantEnvironmentPrediction() Random forest	${\chi }^2$ test (prop.test()) Wilcoxon test (wilcox.test())	D/C	$1$
(D)	InvariantTargetPrediction() Random forest GAM	F-Test Wilcoxon test (wilcox.test())	D/C	$\ge 1$
(E)	InvariantResidualDistributionTest() Random forest GAM	Kolmogorov-Smirnov (ks.test()) Levene’s test + Wilcoxon test (levene.test(), wilcox.test())	D	$1$
(F)	InvariantConditionalQuantilePrediction() Quantile regression forest	Fisher’s exact test (fisher.test())	D	$1$

Equations

$$Z_1$$

$$Z_2$$

$$Z_q$$

$$Y\leftarrow f(X_{S^{\ast }})+\epsilon ,$$

$$Y\mbox{${}\perp\mkern-11.0mu\perp{}$}E\;|\;X_{S^{*}}.$$

$$H_{0,S}:\qquad Y\mbox{${}\perp\mkern-11.0mu\perp{}$}E\;|\;X_{S}.$$

$$\hat{S}:=\bigcap _{S:H_{0,S}\text{ not rejected}}S.$$

$$P\big{(}H_{0,S^{*}}\mbox{ is rejected at level }\alpha\big{)}\leq\alpha,$$

$$P\big{(}\hat{S}\subseteq S^{*}\big{)}\;\geq\;P\big{(}H_{0,S^{*}}\mbox{ accepted}\big{)}\;\geq\;1-\alpha.$$

$$Y=2X+N,N\mbox{${}\perp\mkern-11.0mu\perp{}$}X,$$

$$Y=2X+M,M\mbox{${}\perp\mkern-11.0mu\perp{}$}X,$$

$$X_1$$

$$X_2$$

$$X_3$$

$$Y$$

$$X_1$$

$$Y$$

$$X_2$$

$$\forall x\in \mathbb{R}:f(x)$$

$$\forall x\in \mathbb{R}:g(x)$$

$$\beta {\gamma }^2-\gamma$$

$$S_1$$

$$S_2$$

$$S_3$$

$$P(S^{\ast }\cap \hat{D}=\varnothing )\le P(H_{0,S^{\ast }}\text{ rejected})\le \alpha .$$

$${\hat{D}}_1$$

$${\hat{D}}_2$$

$${\hat{D}}_3$$

$${\hat{D}}_4$$

$${\hat{D}}_5$$

$${\hat{D}}_6$$

$${\hat{D}}_7$$

$$P(f\in {\hat{\mathcal{F}}}_{S^{\ast }})\ge 1-\alpha .$$

$$\hat{\mathcal{F}}:=\bigcup _{S:H_{0,S}\text{ not rejected}}{\hat{\mathcal{F}}}_S.$$

$$P\big{(}f\in\hat{\mathcal{F}}\big{)}\;\geq\;1-2\alpha,$$

$${\mathit{ACE}}(\tilde{x},x)\;:=\;E\big{(}Y\big{|}\mbox{do}(X=\tilde{x})\big{)}-E\big{(}Y\big{|}\mbox{do}(X={x})\big{)}.$$

$$E\big{(}Y\big{|}\mbox{do}(X=x)\big{)}\;=\;E\big{(}Y\big{|}\mbox{do}(X_{S^{*}}=x_{S^{*}})\big{)},$$

$$E\big{(}Y\big{|}\mbox{do}(X_{S^{*}}=x_{S^{*}})\big{)}\;=\;E\big{(}Y\big{|}X_{S^{*}}=x_{S^{*}}\big{)},$$

$$\widehat{{\mathit{ACE}}}(\tilde{x},x):=\big{[}\inf_{g\in\hat{\mathcal{F}}}(g(\tilde{x})-g(x)),\;\sup_{g\in\hat{\mathcal{F}}}(g(\tilde{x})-g(x))\big{]},$$

Full text

In this work, we present and evaluate an array of methods for nonlinear and nonparametric versions of ICP for learning the causal parents of given target variables. We find that an approach which first fits a nonlinear model with data pooled over all environments and then tests for differences between the residual distributions across environments is quite robust across a large variety of simulation settings. We call this procedure “invariant residual distribution test”. In general, we observe that the performance of all approaches is critically dependent on the true (unknown) causal structure and it becomes challenging to achieve high power if the parental set includes more than two variables.

As a real-world example, we consider fertility rate modeling which is central to world population projections. We explore predicting the effect of hypothetical interventions using the accepted models from nonlinear ICP. The results reaffirm the previously observed central causal role of child mortality rates.

1 Introduction

Invariance based causal discovery (Peters et al., 2016) relies on the observation that the conditional distribution of the target variable $Y$ given its direct causes remains invariant if we intervene on variables other than $Y$. While the proposed methodology in Peters et al. (2016) focuses on linear models, we extend Invariant Causal Prediction to nonlinear settings. We first introduce the considered structural causal models in Section 1.1 and review related approaches to causal discovery in Section 1.2. The invariance approach to causal discovery from Peters et al. (2016) is briefly summarized in Section 1.3 and we outline our contribution in Section 1.4. In Section 1.5 we introduce the problem of fertility rate modeling which we consider as a real-world example throughout this work.

1.1 Structural causal models

Assume an underlying structural causal model (also called structural equation model) (e.g. Pearl, 2009)

[TABLE]

for which the functions $g_{k}$, $k=1,\ldots,q$, as well as the parents $\text{pa}_{k}\subseteq\{1,\ldots,q\}\setminus\{k\}$ of each variable are unknown. Here, we have used the notation $Z_{S}=(Z_{i_{1}},\ldots,Z_{i_{s}})$ for any set $S=\{i_{1},\ldots,i_{s}\}\subseteq\{1,\ldots,q\}$. We assume the corresponding directed graph to be acyclic. We further require the noise variables $\eta_{1},\ldots,\eta_{q}$ to be jointly independent and to have zero mean, i.e. we assume that there are no hidden variables.

Due to its acyclic structure, it is apparent that such a structural causal model induces a joint distribution $P$ over the observed random variables. Interventions on the system are usually modeled by replacing some of the structural assignments (e.g. Pearl, 2009). If one intervenes on variable $Z_{3}$, for example, and sets it to the value $5$, the system again induces a distribution over $Z_{1},\ldots,Z_{q}$, that we denote by $P(\cdot|do(Z_{3}\leftarrow 5))$. It is usually different from the observational distribution $P$. We make no counterfactual assumptions here: we assume a new realization $\eta$ is drawn from the noise distribution as soon as we make an intervention.111The new realization of $\eta$ under an intervention and the realization under observational data can be assumed to be independent. However, such an assumption is untestable since we can never observe realizations under different interventions simultaneously and we do not make statements or assumptions about the joint distribution of observational and interventional settings.

1.2 Causal discovery

In causal discovery (also called structure learning) one tries to reconstruct the structural causal model or its graphical representation from its joint distribution (e.g. Pearl, 2009; Spirtes et al., 2000; Peters et al., 2017; Chickering, 2002; Peters et al., 2014; Heckerman, 1997; Hauser and Bühlmann, 2015).

Existing methods for causal structure learning can be categorized along a number of dimensions, such as (i) using purely observational data vs. using a combination of interventional and observational data; (ii) score-based vs. constraint-based vs. “other” methods; (iii) allowing vs. precluding the existence of hidden confounders; (iv) requiring vs. not requiring faithfulness;222 A distribution satisfies faithfulness and the global Markov condition with respect to a graph $G$ if the following statement holds for all disjoint sets $A$, $B$, and $C$ of variables: $A$ is independent of $B$, given $C$, if and only if $A$ is $d$-separated (in $G$) from $B$, given $C$. The concept of $d$-separation (Pearl, 1985, 1988) is defined in Peters et al. (2017, Def. 6.1), for example. (v) type of object that the method estimates. Moreover, different methods vary by additional assumptions they require. In the following, we give brief descriptions of the most common methods for causal structure learning333Also see Heinze-Deml et al. (2018) for a review and empirical comparison of recently proposed causal structure learning algorithms..

The PC algorithm (Spirtes et al., 2000) uses observational data only and estimates the Markov equivalence class of the underlying graph structure, based on (conditional) independence tests under a faithfulness assumption. The presence of hidden confounders is not allowed. Based on the PC algorithm, the IDA algorithm (Maathuis et al., 2009) computes bounds on the identifiable causal effects.

The FCI algorithm is a modification of the PC algorithm. It also relies on purely observational data while it allows for hidden confounders. The output of FCI is a partial ancestral graph (PAG), i.e. it estimates the Markov equivalence class of the underlying maximal ancestral graph (MAG). Faster versions, RFCI and FCI+, were proposed by Colombo et al. (2012) and Claassen et al. (2013), respectively.

The PC, FCI, RFCI and FCI+ algorithms are formulated such that they allow for an independence oracle that indicates whether a particular (conditional) independence holds in the distribution. These algorithms are typically applied in the linear Gaussian setting where testing for conditional independence reduces to testing for vanishing partial correlation.

One of the most commonly known score-based methods is greedy equivalence search (GES). Using observational data, it greedily searches over equivalence classes of directed acyclic graphs for the best scoring graph (all graphs within the equivalence class receive the same score) where the score is given by the Bayesian information criterion, for example. Thus, GES is based on an assumed parametric model such as linear Gaussian structural equations or multinomial distributions. The output of GES is the estimated Markov equivalence class of the underlying graph structure. Heckerman (1997) describe a score-based method with a Bayesian score.

Greedy interventional equivalence search (GIES) extends GES to operate on a combination of interventional and observational data. The targets of the interventions need to be known and the output of GIES is the estimated interventional Markov equivalence class. The latter is typically smaller than the Markov equivalence class obtained when using purely observational data.

Another group of methods makes restrictive assumptions which allows for obtaining full identifiability. Such assumptions include non-Gaussianity (Shimizu et al., 2006) or equal variances (Peters and Bühlmann, 2014) of the errors or non-linearity of the structural equations in additive noise models (Hoyer et al., 2009; Peters et al., 2014).

Instead of trying to infer the whole graph, we are here interested in settings, where there is a target variable $Y$ of special interest. The goal is to infer both the parental set $S^{*}$ for the target variable $Y$ and confidence bands for the causal effects.

1.3 Invariance based causal discovery

This work builds on the method of Invariant Causal Prediction (ICP) (Peters et al., 2016) and extends it in several ways. The method’s key observation is that the conditional distribution of the target variable $Y$ given its direct causes remains invariant if we intervene on variables other than $Y$. This follows from an assumption sometimes called autonomy or modularity (Haavelmo, 1944; Aldrich, 1989; Hoover, 1990; Pearl, 2009; Schölkopf et al., 2012). In a linear setting, this implies, for example, that regressing $Y$ on its direct causes yields the same regression coefficients in each environment, provided we have an infinite amount of data. In a nonlinear setting, this can be generalized to a conditional independence between an index variable indicating the interventional setting and $Y$, given $X$; see (3). The method of ICP assumes that we are given data from several environments. It searches for sets of covariates, for which the above property of invariance cannot be rejected. The method then outputs the intersection of all such sets, which can be shown to be a subset of the true set with high probability, see Section 2.1 and Algorithm 1 in Appendix B for more details. Such a coverage guarantee is highly desirable, especially in causal discovery, where information about ground truth is often sparse.

In many real life scenarios, however, relationships are not linear and the above procedure can fail: The true set does not necessarily yield an invariant model and the method may lose its coverage guarantee, see Example 2. Furthermore, environments may not come as a categorical variable but as a continuous variable instead. In this work, we extend the concept of ICP to nonlinear settings and continuous environments. The following paragraph summarizes our contributions.

1.4 Contribution

Our contributions are fivefold.

Conditional independence tests.

We extend the method of ICP to nonlinear settings by considering conditional independence tests. We discuss in Section 3 and in more length in Appendix B several possible nonlinear and nonparametric tests for conditional independence of the type (3) and propose alternatives. There has been some progress towards nonparametric independence tests (Bergsma and Dassios, 2014; Hoeffding, 1948; Blum et al., 1961; Rényi, 1959; Székely et al., 2007; Zhang et al., 2011). However, in the general nonparametric case, no known non-trivial test of conditional independence has (even asymptotically) a type I error rate less than the pre-specified significance level. This stresses the importance of empirical evaluation of conditional independence tests.

Defining sets.

We discuss in Section 2.2 cases of poor identifiability of the causal parents. If there are highly correlated variables in the dataset, we might get an empty estimator if we follow the approach proposed in (Peters et al., 2016). We can, however, extract more information via defining sets. The results are to some extent comparable to similar issues arising in multiple testing (Goeman and Solari, 2011). For example, if we know that the parental set of a variable $Y$ is either $S=\{1,3\}$ or $S=\{2,3\}$, we know that $\{3\}$ has to be a parent of $Y$. Yet we also want to explore the information that one variable out of the set $\{1,2\}$ also has to be causal for the target variable $Y$, even if we do not know which one out of the two.

Confidence bands for causal effects.

Beyond identifying the causal parents, we can provide nonparametric or nonlinear confidence bands for the strength of the causal effects, as shown in Section 2.3.

Prediction under interventions.

Using the accepted models from nonlinear ICP, we are able to forecast the average causal effect of external interventions. We will discuss this at hand of examples in Section 2.4.

Software.

R (R Core Team, 2017) code for nonlinear ICP is provided in the package nonlinearICP. The proposed conditional independence tests are part of the package CondIndTests. Both packages are available from CRAN.

1.5 Fertility rate modeling

At the hand of the example of fertility rate modeling, we shall explore how to exploit the invariance of causal models for causal discovery in the nonlinear case.

Developing countries have a significantly higher fertility rate compared to Western countries. The fertility rate can be predicted well from covariates such as ‘infant mortality rate’ or ‘GDP per capita’. Classical prediction models, however, do not answer whether an active intervention on some of the covariates leads to a change in the fertility rate. This can only be answered by exploiting causal knowledge of the system.

Traditionally, in statistics the methods for establishing causal relations rely on carefully designed randomized studies. Often, however, such experiments cannot be performed. For instance, factors like ‘infant mortality rate’ are highly complex and cannot be changed in isolation. We may still be interested in the effect of a policy that aims at reducing the infant mortality rate but this policy cannot be randomly assigned to different groups of people within a country.

There is a large body of work that is trying to explain changes in fertility; for an interesting overview of different theories see Hirschman (1994) and the more recent Huinink et al. (2015). There is not a single established theory for changes in fertility and we should clarify in the beginning that all models we will be using will have shortcomings, especially the shortcoming that we might not have observed all relevant variables. We would nevertheless like to take the fertility data as an example to establish a methodology that allows data-driven answers; discussing potential shortfalls of the model is encouraged and could be beneficial in further phrasing the right follow-up questions and collecting perhaps more suitable data.

An interesting starting point for us was the work of Raftery et al. (1995) and very helpful discussions with co-author Adrian Raftery. That work tries to distinguish between two different explanatory models for a decline in fertility in Iran. One model argues that economic growth is mainly responsible; another argues that transmission of new ideas is the primary factor (ideation theory). What allows a distinction between these models is that massive economic growth started in 1955 whereas ideational changes occurred mostly 1967 and later. Since the fertility began to drop measurably already in 1959, the demand theory seems more plausible and the authors conclude that reduced child mortality is the key explanatory variable for the reduction in fertility (responsible for at least a quarter of the reduction).

Note the way we decide between two potential causal theories for a decline in fertility: if a causal model is valid, it has to be able to explain the decline consistently. In particular, the predictions of the model have to be valid for all time-periods, including the time of 1959 with the onset of the fertility decline. The ideation theory wrongly places the onset of fertility decline later and is thus dismissed as less plausible.

The invariance approach of Peters et al. (2016) we follow here for linear models has a similar basic idea: a causal model has to work consistently. In our case, we choose geographic location instead of time for the example and demand that a causal model has to work consistently across geographic locations or continents. We collect all potential models that show this invariance and know that if the underlying assumption of causal sufficiency is true and we have observed all important causal variables then the causal model will be in the set of retained models. Clearly, there is room for a healthy and interesting debate to what extent the causal sufficiency assumption is violated in the example. It has been argued, however, that missing variables do not allow for any invariant model, which renders the method to remain conservative (Peters et al., 2016, Prop. 5).

We establish a framework for causal discovery in nonlinear models. Incidentally, the approach also identifies reduced child mortality as one of key explanatory variables for a decline in fertility.

2 Nonlinear Invariant Causal Prediction

We first extend the approach of (Peters et al., 2016) to nonlinear models, before discussing defining sets, nonparametric confidence bands and prediction under interventions.

2.1 Invariance approach for causal discovery

Peters et al. (2016) proposed an invariance approach in the context of linear models. We describe the approach here in a notationally slightly different way that will simplify statements and results in the nonlinear case and allow for more general applications. Assume that we are given a structural causal model (SCM) over variables $(Y,X,E)$, where $Y$ is the target variable, $X$ the predictors and $E$ so-called environmental variables.

Definition 1** (Environmental variables)**

We know or assume that the variables $E$ are neither descendants nor parents of $Y$ in the causal DAG of $(Y,X,E)$. If this is the case, we call $E$ environmental variables.

In Peters et al. (2016), the environmental variables were given and non-random. Note that the definition above treats the variables as random but we can in practice condition on the observed values of $E$. The definition above excludes the possibility that there is a direct causal connection between one of the variables in $E$ and $Y$. We will talk in the following about the triple of random variables $(Y,X,E)$, where the variable $X$ of predictor variables is indexed by $X_{1},\ldots,X_{p}$. With a slight abuse of notation, we let $S^{*}\subseteq\{1,\ldots,p\}$ be the indices of $X$ that are causal parents $\text{pa}_{Y}$ of $Y$. Thus, the structural equation for $Y$ can be written as

[TABLE]

where $f:\mathbb{R}^{|S^{*}|}\to\mathcal{Y}$. We let $\mathcal{F}$ be the function class of $f$ and let $\mathcal{F}_{S}$ be the subclass of functions that depend only on the set $S\subseteq\{1,\ldots,p\}$ of variables. With this notation we have $f\in\mathcal{F}_{S^{*}}$.

The assumption of no direct effect of $E$ on $Y$ is analogous to the typical assumptions about instrumental variables (Angrist et al., 1996; Imbens, 2014). See Section 5 in Peters et al. (2016) for a longer discussion on the relation between environmental variables and instrumental variables. The two main distinctions between environmental and instrumental variables are as follows. First, we do not need to test for the “weakness” of instrumental/environmental variables since we do not assume that there is a causal effect from $E$ on the variables in $X$. Second, the approaches are used in different contexts. With instrumental variables, we assume the graph structure to be known typically and want to estimate the strength of the causal connections, whereas the emphasis is here on both causal discovery (what are the parents of a target?) and then also inference for the strength of causal effects. With a single environmental variable, we can identify in some cases multiple causal effects whereas the number of instrumental variables needs to match or exceed the number of variables in instrumental variable regression. The instrumental variable approach, on the other hand, can correct for unobserved confounders between the parents and the target variable if their influence is linear, for example. In these cases, our approach could remain uninformative (Peters et al., 2016, Proposition 5).

Example (Fertility data).

In this work, we analyze a data set provided by the United Nations (2013). Here, $Y,X$ and $E$ correspond to the following quantities:

1. (a)

$Y\in\mathbb{R}$ is the total fertility rate (TFR) in a country in a given year, 2. (b)

$X\in\mathbb{R}^{9}$ are potential causal predictor variables for TFR:

• –

IMR – infant mortality rate

• –

Q5 – under-five mortality rate

• –

Education expenditure (% of GNI)

• –

Exports of goods and services (% of GDP)

• –

GDP per capita (constant 2005 US$)

• –

GDP per capita growth (annual %)

• –

Imports of goods and services (% of GDP)

• –

Primary education (% female)

• –

Urban population (% of total) 3. (c)

$E\in\{C_{1},C_{2},C_{3},C_{4},C_{5},C_{6}\}$ is the continent of the country, divided into the categories Africa, Asia, Europe, North and South America and Oceania. If viewed as a random variable (which one can argue about), the assumption is that the continent is not a descendant of the fertility rate, which seems plausible. For an environmental variable, the additional assumption is that the TFR in a country is only indirectly (that is via one of the other variables) influenced by which continent it is situated on (cf. Figure 1).

Clearly, the choices above are debatable. We might for example also want to include some ideation-based variables in $X$ (which are harder to measure, though) and also take different environmental variables $E$ such as time instead of geographic location. We could even allow for additive effects of the environmental variable on the outcome of interest (such as a constant offset for each continent) but we do not touch this debate much more here as we are primarily interested in the methodological development.

The basic yet central insight underlying the invariance approach is the fact that for the true causal parental set $S^{*}:=\text{pa}_{Y}$ we have the following conditional independence relation under Definition 1 of environmental variables:

[TABLE]

This follows directly from the local Markov condition (e.g. Lauritzen, 1996). The goal is to find $S^{*}$ by exploiting the above relation (2). Suppose we have a test for the null hypothesis

[TABLE]

It was then proposed in Peters et al. (2016) to define an estimate $\hat{S}$ for the parental set $S^{*}$ by setting

[TABLE]

Here, the intersection runs over all sets $S$, s.t. $E\cap S=\emptyset$. If the index set is empty, i.e. $H_{0,S}$ is rejected for all sets $S$, we define $\hat{S}$ to be the empty set. If we can test (3) with the correct type I error rate in the sense that

[TABLE]

then we have as immediate consequence the desired statement

[TABLE]

This follows directly from the fact that $S^{*}$ is accepted with probability at least $1-\alpha$ since $H_{0,S^{*}}$ is true; see Peters et al. (2016) for details.

In the case of linear models, the method proposed by Peters et al. (2016, Eq. (16)) considers a set $S$ as invariant if there exist linear regression coefficients $\beta$ and error variance $\sigma$ which are identical across all environments. We consider the conditional independence relation in (3) as a generalization, even for linear relations. In the following example the regression coefficients are the same in all environments, and the residuals have the same mean and variance, but differ in higher order moments (cf. Peters et al., 2016, Eq. (3)):

Example 1

Consider a discrete environmental variable $E$. If in $E=1$ we have

[TABLE]

and in $E=2$

[TABLE]

where $M$ and $N$ have the same mean and variance but differ in higher order moments. In this case, we would have E\mbox{{}\not!\perp\mkern-11.0mu\perp{}}Y\,|\,X, but the hypothesis “same linear regression coefficients and error variance” cannot be rejected.

The question remains how to test (3). If we assume a linear function $f$ in the structural equation (1), then tests that can guarantee the level as in (5) are available (Peters et al., 2016). The following examples show what could go wrong if the data contain nonlinearities that are not properly taken into account.

Example 2** (Linear model and nonlinear data)**

Consider the following SCM, in which $X_{2}$ and $X_{3}$ are direct causes of $Y$.

[TABLE]

Due to the nonlinear effect, a linear regression from $Y$ on $X_{2}$ and $X_{3}$ does not yield an invariant model. If we regress $Y$ on $X_{1}$, however, we obtain invariant prediction and independent residuals. In this sense, the linear version of ICP fails but it still chooses a set of ancestors of $Y$ (it can be argued that this failure is not too severe).

Example 3** (Linear model and nonlinear data)**

In this example, the model misspecification leads to a wrong set that includes a descendant of $Y$. Consider the following SCM

[TABLE]

with independent Gaussian error terms. Furthermore, assume that

[TABLE]

for some $\alpha,\beta$ and $h:\mathbb{R}\rightarrow\mathbb{R}$. Then, in the limit of an infinite sample size, the set $\{X_{1},X_{2}\}$ is the only set that, after a linear regression, yields residuals that are independent of $E$. (To see this write $Y=f(X_{1})+\eta_{Y}$ as a linear function in $X_{1}$, $X_{2}$ and show that the covariance between the residuals and $X_{2}$ is zero.) Here, the functions have to be “fine-tuned” in order to make the conditional $Y|X_{1},X_{2}$ linear in $X_{1}$ and $X_{2}$.444This example is motivated by theory that combines linear and nonlinear models with additive noise (Rothenhäusler et al., 2018). As an example, one may choose $Y\leftarrow X_{1}+0.5X_{1}^{2}+\eta_{Y}$ and $X_{2}\leftarrow 0.5X_{1}^{2}-X_{1}+Y+\eta_{2}$ and $\eta_{1},\eta_{Y},\eta_{2}$ i.i.d. with distribution $\mathcal{N}(0,\sigma^{2}=0.5)$.

The examples show that ICP loses its coverage guarantee if we assume linear relationships for testing (3) while the true data generating process is nonlinear.

In the general nonlinear and nonparametric case, however, it becomes more difficult to guarantee the type I error rate when testing the conditional independence (3) (Shah and Peters, 2018). This in contrast to nonparametric tests for (unconditional) independence (Bergsma and Dassios, 2014; Székely et al., 2007). In a nonlinear conditional independence test setting, where we know an appropriate parametric basis expansion for the causal effect of the variables we condition on, we can of course revert back to unconditional independence testing. Apart from such special circumstances, we have to find tests that guarantee the type I error rate in (5) as closely as possible under a wide range of scenarios. We describe some methods that test (3) in Section 3 but for now let us assume that we are given such a test. We can then apply the method of nonlinear ICP (4) to the example of fertility data.

Example (Fertility data).

The following sets were accepted at the level $\alpha=0.1$ when using nonlinear ICP with invariant conditional quantile prediction (see Appendix B for details) as a conditional independence test:

[TABLE]

As the intersection of $S_{1},\ldots,S_{3}$ is empty, we have $\hat{S}=\emptyset$. This motivates the concept of defining sets.

2.2 Defining sets

It is often impossible to distinguish between highly correlated variables. For example, infant mortality IMR and under-five mortality Q5 are highly correlated in the data and can often be substituted for each other. We accept sets that contain either of these variables. When taking the intersection as in (4), this leads to exclusion of both variables in $\hat{S}$ and potentially to an altogether empty set $\hat{S}$. We can instead ask for the defining sets (Goeman and Solari, 2011), where a defining set $\hat{D}\subseteq\{1,\ldots,p\}$ has the properties

1. (i)

$S\cap\hat{D}\neq\emptyset$ for all $S$ such that $H_{0,S}$ is accepted. 2. (ii)

there exists no strictly smaller set $D^{\prime}$ with $D^{\prime}\subset\hat{D}$ for which property (i) is true.

In words, we are looking for subsets $\hat{D}$, such that each accepted set $S$ has at least one element that also appears in $\hat{D}$. If the intersection $\hat{S}$ (4) is non-empty, any subset of $\hat{S}$ that contains only one variable is a defining set. Defining set are especially useful, however, in cases where the intersection $\hat{S}$ is empty. We still know that, with high probability, at least one of the variables in the defining set $\hat{D}$ has to be a parent. Defining sets are not necessarily unique. Given a defining set $\hat{D}$, we thus know that

[TABLE]

That is, a) at least one of the variables in the defining set $\hat{D}$ is a parent of the target, and b) the data do not allow to resolve it on a finer scale.

Example (Fertility data).

We obtain seven defining sets:

[TABLE]

Thus the highly-correlated variables infant mortality IMR and under-five mortality Q5 indeed form one of the defining sets in this example in the sense that we know at least one of the two is a causal parent for fertility but we cannot resolve which one it is or whether both of them are parents.

2.3 Confidence bands

For a given set $S$, we can in general construct a $(1-\alpha)$-confidence band $\hat{\mathcal{F}}_{S}$ for the regression function when predicting $Y$ with the variables $X_{S}$. Note that if $f$ is the regression function when regressing $Y$ on the true set of causal variables $X_{S^{*}}$ and hence, then, with probability $1-\alpha$, we have

[TABLE]

Furthermore, from Section 2.1 we know that $H_{0,S^{*}}$ is accepted with probability $1-\alpha$. We can hence construct a confidence band for the causal effects as

[TABLE]

Using a Bonferroni correction, we have the guarantee that

[TABLE]

where the coverage guarantee is point-wise or uniform, depending on the coverage guarantee of the underlying estimators $\hat{\mathcal{F}}_{S}$ for all given $S\subseteq\{1,\ldots,p\}$.

2.4 Average causal effects

The confidence bands $\hat{\mathcal{F}}$ themselves can be difficult to interpret. Interpretability can be guided by looking at the average causal effect in the sense that we compare the expected response at $\tilde{x}$ and $x$:

[TABLE]

For the fertility data, this would involve a hypothetical scenario where we fix the variables to be equal to ${x}$ for a country in the second term and, for the first term, we set the variables to $\tilde{x}$, which might differ from $x$ just in one or a few coordinates. Eq. (7) then compares the average expected fertility between these two scenarios. Note that the expected response under a do-operation is just a function of the causal variables $S^{*}\subseteq\{1,\ldots,p\}$. That is—in the absence of hidden variables—we have

[TABLE]

and the latter is then equal to

[TABLE]

that is it does not matter whether we set the causal variables to a specific value $x_{S^{*}}$ or whether they were observed in this state.

Once we have a confidence band as defined in (6), we can bound the average causal effect (7) by the interval

[TABLE]

with the immediate guarantee that

[TABLE]

where the factor $2\alpha$ is guarding, by a Bonferroni correction, against both a probability $\alpha$ that $S^{*}$ will not be accepted—and hence $\hat{S}\subseteq S^{*}$ is not necessarily true—and another probability $\alpha$ that the confidence bands will not provide coverage for the parental set $S^{*}$.

Example (Fertility data).

The confidence bands $\hat{\mathcal{F}}$, required for the computation of $\widehat{{\mathit{ACE}}}(\tilde{x},x)$, are obtained by a time series bootstrap (Künsch, 1989) as the fertility data contain temporal dependencies. The time series bootstrap procedure is described in Appendix A. We use a level of $\alpha=0.1$ which implies a coverage guarantee of 80% as per (8). In the examples below, we set $x$ to an observed data point and vary only $\tilde{x}$.

In the first example, we consider the observed covariates for Nigeria in 1993 as ${x}$. The point of comparison $\tilde{x}$ is set equal to ${x}$, except for the variables in the defining set $\hat{D}_{1}=\{\text{IMR, Q5}\}$. In Figures 2LABEL:sub@fig:def_sets_ind1 and LABEL:sub@fig:def_sets_ind2, these are varied individually over their respective quantiles. The overall confidence interval $\hat{\mathcal{F}}$ consists of the union of the shown confidence intervals $\hat{\mathcal{F}}_{S}$. If $x=\tilde{x}$ (shown by the vertical lines), the average causal effect is zero, of course. In neither of the two scenarios shown in Figures 2LABEL:sub@fig:def_sets_ind1 and LABEL:sub@fig:def_sets_ind2, we observe consistent effects different from zero as some of the accepted models do not contain IMR and some do not contain Q5. However, when varying the variables $\hat{D}_{1}=\{\text{IMR, Q5}\}$ jointly (see Figure 2LABEL:sub@fig:def_sets_jointly), we see that all accepted models predict an increase in expected $\log(\text{TFR})$ as IMR and Q5 increase.

In the second example, we compare the expected fertility rate between countries where all covariates are set to the value ${x}$, which is here chosen to be equal to the observed values of all African countries in 2013. The expected value of log-fertility under this value ${x}$ of covariates is compared to the scenario where we take as $\tilde{x}$ the same value but set the values of the child-mortality variables IMR and Q5 to their respective European averages. The union of intervals in Figure 3LABEL:sub@fig:ace (depicted by the horizontal line segments) correspond to $\widehat{{\mathit{ACE}}}(\tilde{x},x)$ for each country under nonlinear ICP with invariant conditional quantile prediction. The accepted models make largely coherent predictions for the effect associated with this comparison. For most countries, the difference is negative, meaning that the average expected fertility declines if the child mortality rate in a country decreases to European levels. The countries where $\widehat{{\mathit{ACE}}}_{\text{IMR}+\text{Q5}}(\tilde{x},x)$ contains 0 typically have a child mortality rate that is close to European levels, meaning that there is no substantial difference between the two points $\tilde{x},x$ of comparison.

For comparison, in Figure 3LABEL:sub@fig:rf, we show the equivalent computation as in Figure 3LABEL:sub@fig:ace when all covariates are assumed to have a direct causal effect on the target and a Random Forest is used for estimation (Breiman, 2001). We observe that while the resulting regression bootstrap confidence intervals often overlap with $\widehat{{\mathit{ACE}}}_{\text{IMR}+\text{Q5}}(\tilde{x},x)$, they are typically much smaller. This implies that if the regression model containing all covariates was—wrongly—used as a surrogate for the causal model, the uncertainty of the prediction would be underestimated. Furthermore, such an approach ignoring the causal structure can lead to a significant bias in the prediction of causal effects when we consider interventions on descendants of the target variable, for example.

Lastly, we consider a cross validation scheme over time to assess the coverage properties of nonlinear ICP. We leave out the data corresponding to one entire continent and run nonlinear ICP with invariant conditional quantile prediction using the data from the remaining five continents.

We perform this leave-one-continent-out scheme for different values of $\alpha$. For each value of $\alpha$, we then compute the predicted change in the response $\log{(\text{TFR})}$ from 1973 – 2008 for each country belonging to the continent that was left out during the estimation procedure. The predictions are obtained by using the respective accepted models.555Their number differs according to $\alpha$: for a smaller $\alpha$, additional models can be accepted compared to using a larger value of $\alpha$. In other words, the accepted models for $\alpha_{2}$ where $\alpha_{1}<\alpha_{2}$ are a subset of the accepted models for $\alpha_{1}$. We then compare the union of the associated confidence intervals with the real, observed change in $\log{(\text{TFR})}$. This allows us to compute the coverage statistics shown in Table 1. We observe that nonlinear ICP typically achieves more accurate coverage compared to (i) a Random Forest regression model including all variables and (ii) a baseline where the predicted change in $\log{(\text{TFR})}$ for a country is the observed mean change in $\log{(\text{TFR})}$ across all continents other than the country’s own continent. Figures 4 and 5 show the confidence intervals and the observed values for all African countries (Figure 4) and all countries (Figure 5) with observed $\log{(\text{TFR})}$ in 1973 and 2008.

Recall that one advantage of a causal model is that, in the absence of hidden variables, it does not matter whether certain variables have been intervened on or whether they were observed in this state – the resulting prediction remains correct in any of these cases. On the contrary, the predictions of a non-causal model can become drastically incorrect under interventions. This may be one reason for the less accurate coverage statistics of the Random Forest regression model—in this example, it seems plausible that some of the predictors were subject to different external ‘interventions’ across continents and countries.

3 Conditional Independence Tests

We present and evaluate an array of methods for testing conditional independence in a nonlinear setting, many of which exploit the invariance of causal models across different environments. Here, we briefly sketch the main ideas of the considered tests, their respective assumptions and further details are provided in Appendix B. All methods (A) – (F) are available in the package CondIndTests for the R language. Table 2 in Appendix B.7 shows the supported methods and options. An experimental comparison of the corresponding power and type I error rates of these tests can be found in Section 4.

1. (A)

Kernel conditional independence test. Use a kernel conditional independence test for Y\mbox{{}\perp\mkern-11.0mu\perp{}}E\;|\;X_{S} (Fukumizu et al., 2008; Zhang et al., 2011). See Appendix B.1 for further details. 2. (B)

Residual prediction test. Perform a nonlinear regression from $Y$ on $X_{S}$, using an appropriate basis expansion, and apply a variant of a Residual Prediction (RP) test (Shah and Bühlmann, 2018). The main idea is to scale the residuals of the regression such that the resulting test statistic is not a function of the unknown noise variance. This allows for a straight-forward test for dependence between the residuals and ($E$, $X_{S}$). In cases where a suitable basis expansion is unknown, random features (Williams and Seeger, 2001; Rahimi and Recht, 2008) can be used as an approximation. See Appendix B.2 for further details. 3. (C)

Invariant environment prediction. Predict the environment $E$, once with a model that uses $X_{S}$ as predictors only and once with a model that uses $(X_{S},Y)$ as predictors. If the null is true and we find the optimal model in both cases, then the out-of-sample performance of both models is statistically indistinguishable. See Appendix B.3 for further details. 4. (D)

Invariant target prediction. Predict the target $Y$, once with a model that uses $X_{S}$ as predictors only and once with a model that uses $(X_{S},E)$ as predictors. If the null is true and we find the optimal model in both cases, then the out-of-sample performance of both models is statistically indistinguishable. See Appendix B.4 for further details. 5. (E)

Invariant residual distribution test. Pool the data across all environments and predict the response $Y$ with variables $X_{S}$. Then test whether the distribution of the residuals is identical in all environments $E$. See Appendix B.5 for further details. 6. (F)

Invariant conditional quantile prediction. Predict a $1-\beta$ quantile of the conditional distribution of $Y$, given $X_{S}$, by pooling the data over all environments. Then test whether the exceedance of the conditional quantiles is independent of the environment variable. Repeat for a number of quantiles and aggregate the resulting individual $p$-values by Bonferroni correction. See Appendix B.6 for further details.

Another interesting possibility for future work would be to devise a conditional independence test based on model-based recursive partitioning (Zeileis et al., 2008; Hothorn and Zeileis, 2015).

Non-trivial, assumption-free conditional independence tests with a valid level do not exist (Shah and Peters, 2018). It is therefore not surprising that all of the above tests assume the dependence on the conditioning variable to be “simple” in one form or the other. Some of the above tests require the noise variable in (1) to be additive in the sense that we do not expect the respective test to have the correct level when the noise is not additive. As additive noise is also used in Sections 2.3 and 2.4, we have written the structural equations above in an additive form.

One of the inherent difficulties with these tests is that the estimation bias when conditioning on potential parents in (3) can potentially lead to a more frequent rejection of a true null hypothesis than the nominal level suggests. In approaches (C) and (D), we also need to test whether the predictive accuracy is identical under both models and in approaches (E) and (F) we need to test whether univariate distributions remain invariant across environments. While these additional tests are relatively straightforward, a choice has to be made.

Discussion of power.

Conditional independence testing is a statistically challenging problem. For the setting where we condition on a continuous random variable, we are not aware of any conditional independence test that holds the correct level and still has (asymptotic) power against a wide range of alternatives. Here, we want to briefly mention some power properties of the tests we have discussed above.

Invariant target prediction (D), for example, has no power to detect if the noise variance is a function of $E$, as shown by the following example

Example 4

Assume that the distribution is entailed by the following model

[TABLE]

where $\eta_{E},\eta_{X},\eta_{Y}\overset{\text{i.i.d.}}{\sim}\mathcal{N}(0,1)$. Then, any regression from $Y$ on $X$ and $E$ yields the same results as regressing $Y$ on $X$ only. That is,

[TABLE]

although

[TABLE]

The invariant residual distribution test (E), in contrast, assumes homoscedasticity and might have wrong coverage if this assumption is violated. Furthermore, two different linear models do not necessarily yield different distributions of the residuals when performing a regression on the pooled data set.

Example 5

Consider the following data generating process

[TABLE]

where the input variables $X^{e=1}$ and $X^{e=2}$ and the noise variables $N^{e=1}$ and $N^{e=2}$ have the same distribution in each environment, respectively. Then, approach (E) will accept the null hypothesis of invariant prediction.

It is possible to reject the null hypothesis of invariant prediction in Example 5 by testing whether in each environment the residuals are uncorrelated from the input.

Invariant conditional quantile prediction (F) assumes neither homoscedasticity nor does it suffer from the same issue of (D), i.e. no power against an alternative where the noise variance $\sigma$ is a function of $E$. However, it is possible to construct examples where (F) will have no power if the noise variance is a function of both $E$ and the causal variables $X_{S^{*}}$. Even then, though, the noise level would have to be carefully balanced to reduce the power to 0 with approach (F) as the exceedance probabilities of various quantiles (a function of $X_{S^{*}}$) would have to remain constant if we condition on various values of $E$.

4 Simulation Study

For the simulations, we generate data from different nonlinear additive noise causal models and compare the performance of the proposed conditional independence tests. The structural equations are of the form $Z_{k}\;\leftarrow\;g_{k}(Z_{\text{pa}_{k}})+\eta_{k},$ where the structure of the DAG is shown in Figure 6 and kept constant throughout the simulations for ease of comparison. We vary the nonlinearities used, the target, the type and strength of interventions, the noise tail behavior and whether parental contributions are multiplicative or additive. The simulation settings are described in Appendix C in detail.

We apply all the conditional independence tests (CITs) that we have introduced in Section 3, implemented with the following methods and tests as subroutines:

[TABLE]

As a disclaimer we have to note that KCI is implemented without Gaussian process estimation. The KCI results shown below might improve if the latter is added to the algorithm.

Baselines.

We compare against a number of baselines. Importantly, most of these methods contain various model misspecifications when applied in the considered problem setting. Therefore, they would not be suitable in practice. However, it is instructive to study the effect of the model misspecifications on performance.

The method of Causal Additive Models (CAM) (Bühlmann et al., 2014) identifies graph structure based on nonlinear additive noise models (Peters et al., 2014). Here, we apply the method in the following way. We run CAM separately in each environment and output the intersection of the causal parents that were retrieved in each environment. Note that the method’s underlying assumption of Gaussian noise is violated. 2. 2.

We run the PC algorithm (Spirtes et al., 2000) in two different variants. We consider a variable to be the parent of the target if a directed edge between them is retrieved; we discard undirected edges. In the first variant of PC we consider, the environment variable is part of the input; conditional independence testing within the PC algorithm is performed with KCI, for unconditional independence testing we use HSIC (Gretton et al., 2005, 2007), using the implementation from Pfister et al. (2017) (denoted with ‘PC(i)’ in the figures). In the second variant, we run the PC algorithm on the pooled data (ignoring the environment information), testing for zero partial correlations (denoted with ‘PC(ii)’ in the figures). Here, the model misspecification is the assumed linearity of the structural equations. 3. 3.

We compare against linear ICP (Peters et al., 2014) where the model misspecification is the assumed linearity of the structural equations. 4. 4.

We compare against LiNGAM (Shimizu et al., 2006), run on the pooled data without taking the environment information into account. Here, the model misspecifications are the assumed linearity of the structural equations and the i.i.d. assumption which does not hold. 5. 5.

We also show the outcome of a random selection of the parents that adheres to the FWER-limit by selecting the empty set ($\hat{S}=\emptyset$) with probability $1-\alpha$ and setting $\hat{S}=\{k\}$ for $k$ randomly and uniformly picked from $\{1,\ldots,p\}\setminus k^{\prime}$ with probability $\alpha$, where $k^{\prime}$ is the index of the current target variable. The random selection is guaranteed to maintain FWER at or below $1-\alpha$.

Thus, all considered baseline models in 1. – 4. —except for ‘PC(i)‘—contain at least slight model misspecifications.

Metrics.

Error rates and power are measured in the following by

1. (i)

Type I errors are measured by the family-wise error rate (FWER), the probability of making one or more erroneous selections

[TABLE] 2. (ii)

Power is measured by the Jaccard similarity, the ratio between the size of the intersection and the size of the union of the estimated set $\hat{S}$ and the true set $S^{*}$. It is defined as 1 if both $S^{*}=\hat{S}=\emptyset$ and otherwise as

[TABLE]

The Jaccard similarity is thus between 0 and 1 and the optimal value 1 is attained if and only if $\hat{S}=S^{*}$.

Type-I-error rate of conditional independence tests.

Figure 7 shows the average FWER on the $x$-axis (and the average Jaccard similarity on the $y$-axis) for all methods. The FWER is close but below the nominal FWER rate of $\alpha=0.05$ for all conditional independence tests, that is $P(\hat{S}\subseteq S^{*})\geq 1-\alpha$. The same holds for the baselines linear ICP and random selection. Notably, the average Jaccard similarity of the random selection baseline is on average not much lower than for the other methods. The reason is mostly a large variation in average Jaccard similarity across the different target variables, as discussed further below and as will be evident from Figure 8 (top right plot). In fact, as can be seen from Figure 9, random guessing is much worse than the optimal methods on each target variable. The FWER of the remaining baselines CAM, LiNGAM, PC(i) and PC(ii) lies well above $\alpha$.

A caveat of the FWER control seen in Figure 7 is that while the FWER is maintained at the desired level, the test $H_{0,S^{*}}$ might be rejected more often than with probability $\alpha$. The error control rests on the fact that $H_{0,S^{*}}$ is accepted with probability higher than $1-\alpha$ (since the null is true for $S^{*}$). However, if a mistake is made and $H_{0,S^{*}}$ is falsely rejected, then we might still have $\hat{S}\subseteq S^{*}$ because either all other sets are rejected, too, in which case $\hat{S}=\emptyset$, or because other sets (such as the empty set) are accepted and the intersection of all accepted sets is—by accident—again a subset of $S^{*}$. In other words: some mistakes might cancel each other out but overall the FWER is very close to the nominal level, even if we stratify according to sample size, target, type of nonlinearity and other parameters, as can be seen from Figure 8.

Power.

Figures 7 shows on the $y$-axis the average Jaccard similarity for all methods. The optimal value is 1 and is attained if and only if $\hat{S}=S^{*}$. A value 0 corresponds to disjoint sets $\hat{S}$ and $S^{*}$. The average Jaccard similarity is around 0.4 for most methods and not clearly dependent on the type I errors of the methods. Figure 8 shows the average FWER and Jaccard similarities stratified according to various parameters.

One of the most important determinants of success (or the most important) is the target, that is the variable for which we would like to infer the causal parents; see top right panel in Figure 8. Variables 1 and 5 as targets have a relatively high average Jaccard similarity when trying to recover the parental set. These two variables have an empty parental set ($S^{*}=\emptyset$) and the average Jaccard similarity thus always exceeds $1-\alpha$ if the level of the procedure is maintained as then $\hat{S}=\emptyset=S^{*}$ with probability at least $1-\alpha$ and the Jaccard similarity is 1 if both $\hat{S}$ and $S^{*}$ are empty. As testing for the true parental set corresponds to an unconditional independence test in this case, maintaining the level of the test procedure is much easier than for the other variables.

Figure 9 shows the same plot as Figure 7 for each of the more difficult target variables 2, 3, 4, and 6 separately. As can be seen from the graph in Figure 6 and the detailed description of the simulations in Appendix C, the parents of target variable 3 are difficult to estimate as the paths $1\rightarrow 2\rightarrow 3$ and $1\rightarrow 3$ cancel each other exactly in the linear setting (and approximately for nonlinear data), thus creating a non-faithful distribution. The cancellation of effects holds true if interventions occur on variable 1 and not on variable 2. A local violation of faithfulness leaves type I error rate control intact but can hurt power as many other sets besides the true $S^{*}$ can get accepted, especially the empty set, thus yielding $\hat{S}=\emptyset$ when taking the intersection across all accepted sets to compute the estimate $\hat{S}$ in (4). Variable 4, on the other hand, has only a single parent, namely $S^{*}=\{3\}$, and the recovery of the single parent is much easier, with average Jaccard similarity up to a third. Variable 6 finally again has average Jaccard similarity of up to around a tenth only. It does not suffer from a local violation of faithfulness as variable 3 but the size of the parental set is now three, which again hurts the power of the procedure, as often already a subset of the true parents will be accepted and hence $\hat{S}$ in (4) will not be equal to $S^{*}$ any longer but just a subset. For instance, when variable 5 is not intervened on in any environment it cannot be identified as a causal parent of variable 6, as it is then indistinguishable from the noise term. Similarly, in the linear setting, merely variable 3 can be identified as a parent of variable 6 if the interventions act on variables 1 and/or 2 only.

The baselines LiNGAM and PC show a larger Jaccard similarity for target variables 3, 4 (only LiNGAM), and 6 at the price of large FWER values.

In Appendix D, Figures 10 – 13 show the equivalent to Figure 8, separately for target variables 2, 3, 4 and 6. For the sample size $n$, we observe that increasing it from 2000 to 5000 decreases power in case of target variables 4. This behavior can be explained by the fact that when testing $S^{*}$ in Eq. (3), the null is rejected too often as the bias in the estimation performed as part of the conditional independence test yields deviations from the null that become significant with increasing sample size. For the nonlinearity, we find that the function $f_{4}(x)=\sin(2\pi x)$ is the most challenging one among the nonlinearities considered. It is associated with very low Jaccard similarity values for the target variables that do have parents. For the intervention type, it may seem surprising that ‘all’ does not yield the largest power. A possible explanation is that intervening on all variables except for the target yields more similar intervention settings—the intervention targets do not differ between environments 2 and 3, even though the strength of the interventions is different. So more heterogeneity between the intervention environments, i.e. also having different intervention targets, seems to improve performance in terms of Jaccard similarity. Lastly, we see that power is often higher for additive parental contributions than for multiplicative ones.

In summary, all tests (A) – (F) seem to maintain the desired type I error, chosen here as the family-wise error rate, while the power varies considerably. An invariant residual distribution test using GAM with Levene’s test and Wilcoxon test produces results here that are constantly as good or nearly as good as the optimal methods for a range of different settings. However, it is only applicable for categorical environmental variables. For continuous environmental variables, the results suggest that the residual prediction test with random features might be a good choice.

5 Discussion and future work

Causal structure learning with the invariance principle was proposed Peters et al. (2016). However, the assumption of linear models in Peters et al. (2016) is unrealistic in many applications. In this work, we have shown how the framework can be extended to nonlinear and nonparametric models by using suitable nonlinear and nonparametric conditional independence tests. The properties of these conditional independence tests are critically important for the power of the resulting causal discovery procedure. We evaluated many different test empirically in the given context and highlighted approaches that seem to work robustly in different settings. In particular we find that fitting a nonlinear model with pooled data and then testing for differences between the residual distributions across environments results in desired coverage and high power if compared against a wide range of alternatives.

Our approach allowed us to model how several interventions may affect the total fertility rate of a country, using historical data about decline and rise of fertilities across different continents. In particular, we provided bounds on the average causal effect under certain (hypothetical) interventions such as a reduction in child mortality rates. We showed that the causal prediction intervals for hold-out data have better coverage than various baseline methods. The importance of infant mortality rate and under-five mortality rate on fertility rates is highlighted, reconfirming previous studies that have shown or hypothesized these factors to be important (Hirschman, 1994; Raftery et al., 1995). We stress that the results rely on causal sufficiency of the used variables, an assumption that can and should be debated for this particular example.

We also introduced the notion of ‘defining sets’ in the causal discovery context that helps in situations where the signal is weak or variables are highly correlated by returning sets of variables of which we know that at least one variable (but not necessarily all) in this set are causal for the target variable in question.

Finally, we provide software in the R (R Core Team, 2017) package nonlinearICP. A collection of the discussed conditional independence tests are part of the package CondIndTests and are hopefully of independent interest.

In applications where it is unclear whether the underlying models are linear or not, we suggest the following. While our proposed methods also hold the significance level if the underlying models are linear, we expect the linear version of ICP to have more power. Therefore, it is advisable to use the linear version of ICP if one has strong reasons to believe that the underlying model is indeed linear. In practice, one might first apply ICP with linear models and apply a nonlinear version if, for example, all linear models are rejected. One would then need to correct for multiple testing by a factor of 2.

Acknowledgements

We thank Jan Ernest and Adrian Raftery for helpful discussions and an AE and two referees for very helpful comments on an earlier version of the manuscript.

A Time series bootstrap procedure

In the time series bootstrap procedure used to obtain the confidence bands $\hat{\mathcal{F}}$ in Section 2.4, $B$ bootstrap samples of the response $Y$ are generated by first fitting the model on all data points. We then use the fitted values and residuals from this model. Each bootstrap sample is generated by resampling the residuals of this fit block- and country-wise. In more detail, we define the block-length $l_{b}$ of residuals that should be sampled consecutively (we use $l_{b}=3$) and we sample a number of time points $t_{s_{1}},\ldots,t_{s_{k}}$ from which the residuals are resampled. For a country $a$ and the first time point $t_{1}$, consider the fitted values at point $t_{1}$ and the fitted values for the $l_{b}-1$ consecutive observations. We then sample a country $b$ and add country $b$’s residuals from time points $t_{s_{1}},t_{s_{1}+1},\ldots t_{s_{1}+l_{b}-1}$ to the fitted values of country $a$ for the considered period $t_{1},\ldots,t_{l_{b}}$. We then proceed with the next $l_{b}$ consecutive fitted values for country $a$ and add country $b$’s residuals from observations $t_{s_{2}},t_{s_{2}+1},\ldots t_{s_{2}+l_{b}-1}$, until all fitted values of country $a$ are covered. This procedure is applied to each country. Finally, to obtain the confidence intervals, we fit the model on each of the $B$ bootstrap samples $(Y^{b},X)$, consisting of the response $Y^{b}$ generated from the fitted values and the resampled residuals, and the observations $X$ which have not been modified.

B Conditional independence tests

For completeness, we first restate the generic method for Invariant Causal Prediction from Peters et al. (2016):

The conditional independence tests discussed in this work can be used to perform the test in Step 2 of Algorithm 1. Therefore, the inputs to these tests consist of an i.i.d. sample of $(Y,X_{S},E)$ and $\alpha$ where $X_{S}$ contains the variables corresponding to $S\subseteq\{1,\ldots,p\}$, i.e. the subset to be tested. Additionally, some test specific parameters might need to be specified. The return value of the tests is the respective test’s decision about $H_{0,S}$.

For most tests, $E\in\mathbb{R}^{d}$ can be either discrete or continuous. As all empirical results in this work use an environment variable that is discrete and one-dimensional, the descriptions below focus on this setting. We then denote the index set of different environments with ${\mathcal{E}}$. We will comment on the required changes for the continuous and higher-dimensional case in the respective sections. Whenever applying the test for environmental variables $E\in\mathbb{R}^{d}$ with $d>1$ is infeasible with the method, each test can be applied separately for each variable in $E$. The overall $p$-value is obtained by multiplying the minimum of the individual $p$-values by $d$, i.e. by applying a Bonferroni correction for the number of environmental variables. When applying the function CondIndTest() from the R package CondIndTests with a conditional independence test that does not support a multidimensional environment variable, the described Bonferroni correction is applied.

B.1 Kernel conditional independence test

Setting and assumptions.

We use the kernel conditional independence test proposed in Zhang et al. (2011). When $E$ is discrete, we use a delta kernel for $E$, and otherwise an RBF kernel. The test is also applicable when $E$ contains more than one environmental variable as the inputs can be sets of random variables.

B.2 Residual Prediction test

Setting and assumptions.

We do not expect this test to have the correct level when the noise in Eq. (1) is not additive. The described procedure does not need to be modified for higher-dimensional and/or continuous environmental variables $E$.

We consider a version of a Residual Prediction test as proposed in Shah and Bühlmann (2018) to determine whether $H_{0,S}$ holds at level $\alpha$ for a particular set of variables $S$. The main idea is to find a suitable basis expansion of $f$ that allows us to regress $Y$ on $X_{S}$ by reverting back to the linear case. Given an appropriate basis expansion, the scaled ordinary least squares residuals can then be tested for possible remaining nonlinear dependencies between the scaled residuals and $(E,X_{S})$. The scaling ensures that the resulting test statistic is not a function of the noise variance. Under the null, the scaled residuals are expected to behave roughly like the noise term. In other words, there should be no dependence between the scaled residuals and the environmental variables and $X_{S}$, so there should be no signal left in the residuals that could be fitted by a nonparametric method like a random forest using $E$ and $X_{S}$ as predictors. This necessitates to make an assumption on the noise distribution $F_{\varepsilon}$, e.g. $\varepsilon\sim\mathcal{N}(0,1)$.

In order to generalize the method to settings where an appropriate basis expansion of $f$ is unknown, we look at ways to find such a suitable basis expansion automatically by using random features (Williams and Seeger, 2001; Rahimi and Recht, 2008).

Step 1.

The choice of $h_{m}(X_{S}),m=1,\ldots,M$ can be based on domain knowledge, e.g. when the nonlinearity in Eq. (1) is known to be a polynomial of a given order. If such domain knowledge is not available, the linear basis expansion can be approximated by random features, e.g. using the Nyström method or by random Fourier features. For these methods, the kernel function needs to be chosen as well as the kernel parameters and the number of random features to be generated.

Step 3.

For instance, a random forest can be used for the estimation. If the residuals only differ in the second moments, predicting the expectation of the residuals is not sufficient as the predictors $E$ have no discriminative power for this task. In such a setting, the absolute value of the residuals can be predicted to exploit the heterogeneity in the second moments across environments.

Step 4.

For instance, the mean squared error can be used here.

Step 5.

If the error term is non-Gaussian, the appropriate distribution can be used at this stage to accommodate non-Gaussianity of the noise.

Parameter settings used in simulations.

In the simulations, we use $B=250$ and $\varepsilon\sim\mathcal{N}(0,1)$. In step 1, we consider the following options: (a) Fourier random features (approach (B)(i) in Section 4), (b) Nyström random features and RBF kernel ((B)(ii)), (c) Nyström random features and polynomial kernel of random degree ((B)(iii)), (d) polynomial basis of random degree ((B)(iv)). The number of random features in (c) and (d) is chosen to be $\lceil n/4\rceil$. In step 7, we predict the mean as well as the absolute value of the residuals and aggregate the results using a Bonferroni correction.

B.3 Invariant environment prediction

Setting and assumptions.

The described procedure does not need to be modified for continuous environmental variables $E$. For higher-dimensional $E$ the test would need to be applied for each variable separately and the resulting $p$-values would need to be aggregated with a Bonferroni correction.

Step 3.

When a random forest is used to predict the environment variable, one can also use $X_{S}$ and a permutation of $Y$ as predictors to ensure the random forest fits are based on the same number of predictor variables. As the number of variables considered for each split in the random forest estimation procedure is a function of the total number of predictor variables, this helps to mitigate differences between the prediction accuracies that are only due to artefacts of the estimation procedure. This is especially relevant for small sets $S$.

Step 5.

For instance, a $\chi^{2}$ test can be used here. If the null is true and we find the optimal model in both cases, then the out-of-sample performance of both models is statistically indistinguishable as $Y$ is independent of $E$ given $X_{S}$. If the null is not true, we expect the model containing the response to perform better as $Y$ contains additional information in this case (since $Y$ is not independent of $E$ given $X_{S}$).

Parameter settings used in simulations.

In step 1, we use $2/3$ of the data points for training and $1/3$ for testing. In step 3, we use a random forest to predict the environment variable and use $X_{S}$ and a permutation of $Y$ as predictors. In step 4, we use the $\chi^{2}$ test implemented in prop.test() (Wilson, 1927) from the stats package in R.

B.4 Invariant target prediction

Setting and assumptions.

The described procedure does not need to be modified for continuous and/or higher-dimensional environmental variables $E$.

Step 3.

When a random forest is used, one can also use $X_{S}$ and a permutation of $E$ as predictors to ensure the random forest fit is based on the same number of predictor variables. As the number of variables considered for each split in the random forest estimation procedure is a function of the total number of predictor variables, this helps to mitigate differences between the prediction accuracies that are only due to artefacts of the estimation procedure. This is especially relevant for small sets $S$. As an alternative to using a random forest, one could use GAMs as the estimation procedure, implying the implicit assumption that the components in $f(X)$ in Eq. (1) are additive.

Step 5.

For instance, an F-test can be used here. Another option is a Wilcoxon test using the difference between the absolute residuals. If the null is true and we find the optimal model in both cases, then the out-of-sample performance of both models is statistically indistinguishable as $Y$ is independent of $E$ given $X_{S}$. If the null is not true, we expect the model additionally containing $E$ to perform better as $E$ contains additional information in this case (since $Y$ is not independent of $E$ given $X_{S}$).

Parameter settings used in simulations.

In step 1, we use $2/3$ of the data points for training and $1/3$ for testing. In step 3, to predict $Y$ we use a GAM or a random forest. In step 5, we use an F-test or a Wilcoxon test (wilcox.test() from the stats package in R). These combinations yield approaches (D)(i) – (iv) in Section 4. When using a random forest in step 3, we use $X_{S}$ and a permutation of $E$ as predictors.

B.5 Invariant residual distribution test

Setting and assumptions.

We do not expect this test to have the correct level when the noise in Eq. (1) is not additive. It is only applicable to discrete environmental variables. For higher-dimensional $E$ the test would need to be applied for each variable separately and the resulting $p$-values would need to be aggregated with a Bonferroni correction.

Step 1.

For instance, one could use a random forest or a GAM as the estimation procedure. The latter implicitly assumes that the components in $f$ in Eq. (1) are additive.

Step 4.

For instance, a nonparametric test such as Kolmogorov-Smirnov can be used here. Alternatively, we can limit the test to assess equality of first and second moments by first using a Wilcoxon test for the expectation with an one-vs-all scheme as described in the algorithm. Subsequently, Levene’s test for homogeneity of variance across groups can be used to test for equality of the second moments of the residual distributions. In this case, the final $p$-value would be twice the minimum of (a) the Bonferroni-corrected $p$-value from the one-vs-all Wilcoxon test and (b) the $p$-value from Levene’s test.

Parameter settings used in simulations.

In step 1, we use a GAM or a random forest. In step 4, we use both approaches described above, using (a) ks.test() from the stats package in R (Conover, 1971) and (b) wilcox.test() and levene.test() (the latter being contained in the lawstat package in R (Levene, 1960; Gastwirth et al., 2015)). These combinations yield approaches (E)(i) – (iv) in Section 4.

B.6 Invariant conditional quantile prediction

Setting and assumptions.

For continuous and/or higher-dimensional environmental variables $E$ the test described in Steps 4 – 11 which assesses whether \text{Exceedance}\mbox{{}\perp\mkern-11.0mu\perp{}}E would need to be modified according to the structure of $E$.

Step 3.

For instance, a Quantile Regression Forest (Meinshausen, 2006) can be used here.

Step 7.

For instance, Fisher’s exact test can be used here by computing the $2\times 2$ contingency table of the exceedance of the residuals for the quantile $1-\beta$ for $I=0$ and for $I=1$.

Parameter settings used in simulations.

In step 3, we use a quantile regression forest for $\mathcal{B}=\{0.1,0.5,0.9\}$. In step 7, we use fisher.test() from the stats package in R.

B.7 Overview of conditional independence tests in CondIndTests package

The described conditional independence tests are available in the R package CondIndTests. A wrapper function CondIndTest() is provided which takes the respective test as the argument method. The package supports the estimation procedures, subroutines and statistical tests shown in Table 2. The column $E$ indicates whether the environmental variables can be discrete (’D’), continuous (’C’), or both; the column $d$ shows the supported dimensionality of $E$.

As described at the beginning of Appendix B, a Bonferroni correction is applied when calling the function CondIndTest() with a conditional independence test that does not support a multidimensional environment variable. Similarly, a Bonferroni correction is applied when the first input argument Y to the respective test is multidimensional and if the specified test does not support this internally.

C Experimental settings for numerical studies

For each simulation, we compare the performance of all methods and conditional independence tests while choosing the following parameters randomly (but keeping them constant for one simulation): In total, there are 27478 simulations from 1240 distinct settings that are evaluated for each of the 22 considered methods.

1. (i)

Sample size. Sample size ‘n’ is chosen randomly in the set $\{100,200,500,2000,5000\}$. 2. (ii)

Target variable. We sample one of the variables in the graph in Figure 6 at random as a target variable (variable ‘target’ is chosen uniformly from $\{1,2,\ldots,6\}$ in other words). 3. (iii)

Tail behavior of the noise. The noise $\eta_{k}$ for $k=1,\ldots,6$ is sampled from a t-distribution and the degrees of freedom are chosen at random from $\mathrm{df}\in\{2,3,5,10,20,50,100\}$, where the latter is very close to a Gaussian distribution. 4. (iv)

Multiplicative or additive effects. For each simulation setting, we determine whether $g_{k}(\cdot)$ has additive or multiplicative components. We sample additive components of the form $g_{k}(Z_{\text{pa}_{k}})=\sum_{j\in pa_{k}}f(\epsilon_{j,k}\cdot Z_{j})$ ($\mathit{multiplic}=\mathrm{FALSE}$) and multiplicative components of the form $g_{k}(Z_{\text{pa}_{k}})=\prod_{j\in pa_{k}}f(\epsilon_{j,k}\cdot Z_{j})$ ($\mathit{multiplic}=\mathrm{TRUE}$) with equal probability, where the signs $\epsilon_{j,k}\in\{-1,1\}$ are as shown in Figure 6 along the relevant arrows. 5. (v)

Shift- or do-Interventions. The variable ‘shift’ is set with equal probability to either $\mathrm{TRUE}$ (shift-interventions) or $\mathrm{FALSE}$ (do-interventions). For do-interventions we replace the structural equation of the intervened variable $k\in\{1,\ldots,q\}$ by

[TABLE]

where $e_{k}$ is the randomly chosen intervention value which is sampled i.i.d for each observation as described under (vi). For shift-interventions, the value $e_{k}$ is added as

[TABLE]

See for example Section 5 of Peters et al. (2016) for a more detailed discussion of shift interventions. 6. (vi)

Strength of interventions The intervention values $e_{k}$ are chosen independently for all variables from a t-distribution with ‘df’ degrees of freedom, shifted by a constant ‘meanshift’ (chosen uniformly at random in $\{0,0.1,0.2,0.5,1,2,5,10\}$, and scaled by a factor ‘strength’, chosen uniformly at random in $\{0,0.1,0.2,0.5,1,2,5,10\}$). 7. (vii)

Non-linearities. For the functions $f=f_{\mathrm{id}}$ we consider the following four nonlinear functions, where the index ‘id’ is sampled uniformly from $\{1,2,3,4\}$ and the same nonlinearity is used throughout the graph:

[TABLE] 8. (viii)

Location of interventions. Each sample is independently assigned into one environment in ${\mathcal{E}}=\{1,2,3\}$, where $\{1\}$ corresponds to observational data, that is all samples in environment $\{1\}$ are sampled as observational data, where samples in environments $\{2,3\}$ are intervention data. The intervention targets and strengths for samples in environment $\{2\}$ are drawn as per the description below and kept constant for all samples in environment $\{2\}$ and then analogously for environment $\{3\}$, where intervention targets are drawn independently and identically to environment $\{2\}$. The variable ‘interv’ is set uniformly at random to one of the values $\{$‘all’,‘rand’,‘close’$\}$. If it is equal to ‘all’, then interventions in environments $\{2,3\}$ occur at all variables except for the target variable. If it is equal to ‘rand’, then interventions occur at one ancestor, chosen uniformly at random, and one descendant of the target variable, chosen again uniformly at random (the empty set is chosen in case there are no ancestors or descendants). Finally, if it is equal to ‘close’, interventions occur at one parent, chosen uniformly at random, and one child of the target variable, chosen again uniformly at random (and again no interventions occur if these sets are empty).

D Additional experimental results

Figures 10 – 13 show the equivalent to Figure 8, separately for target variables 2, 3, 4 and 6. For the sample size $n$, we observe that increasing it from 2000 to 5000 decreases power in case of target variable 4. This behavior can be explained by the fact that when testing $S^{*}$ in Eq. (3), the null is rejected too often as the bias in the estimation performed as part of the conditional independence test yields deviations from the null that become significant with increasing sample size. For the nonlinearity, we find that the function $f_{4}(x)=\sin(2\pi x)$ is the most challenging one among the nonlinearities considered. It is associated with very low Jaccard similarity values for the target variables that do have parents. For the intervention type, it may seem surprising that ‘all’ does not yield the largest power. A possible explanation is that intervening on all variables except for the target yields more similar intervention settings—the intervention targets do not differ between environments 2 and 3, even though the strength of the interventions is different. So more heterogeneity between the intervention environments, i.e. also having different intervention targets, seems to improve performance in terms of Jaccard similarity. Lastly, we see that power is often higher for additive parental contributions than for multiplicative ones.

E Example

Here we illustrate the methods presented in this manuscript by considering a causal DAG $X_{1}\rightarrow X_{2}\rightarrow X_{3}$. Figure 14 visualizes the generated data. There are six environments with shift interventions. The latter act on $X_{1}$ in two environments (green, yellow) and on $X_{3}$ in four environments (green, cyan, blue, magenta). The red environment consists of observational data. We run the proposed approaches (A) – (F) to retrieve the parents of $X_{2}$, i.e. $S^{*}=\{X_{1}\}$. Below we give an overview of which sets were accepted by the respective methods with $\alpha=0.05$. We see that approaches (A), (B)(i)+(ii), (E)(i)-(iii) and (F) retrieve $S^{*}$ correctly, while the other approaches return the empty set.

[TABLE]