Invariant Causal Prediction for Nonlinear Models

1.1 Structural causal models

Assume an underlying structural causal model (also called structural equation model) [e. g. 2]

for which the functions gk, k=1,…,q, as well as the parents pak⊆{1,…,q}∖{k} of each variable are unknown. Here, we have used the notation ZS=(Zi1,…,Zis) for any set S={i1,…,is}⊆{1,…,q}. We assume the corresponding directed graph to be acyclic. We further require the noise variables η1,…,η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. 2]. If one intervenes on variable Z3, for example, and sets it to the value 5, the system again induces a distribution over Z1,…,Zq, that we denote by P(·|do(Z3←5)). It is usually different from the observational distribution P. We make no counterfactual assumptions here: we assume a new realization η is drawn from the noise distribution as soon as we make an intervention.^[1]

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. 2], [e. g. 3], [e. g. 4], [e. g. 5], [e. g. 6], [e. g. 7], [e. g. 8].

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;^[2] (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 learning.^[3]

The PC algorithm [3] 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 [12] 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. [13] and Claassen et al. [14], 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 [7] 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 [15] or equal variances [16] of the errors or non-linearity of the structural equations in additive noise models [17], [6].

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) [1] 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 [18], [19], [20], [2], [21]. 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 II 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.

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 [30] and the more recent Huinink et al. [31]. 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. [32] 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. [1] 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 [1, 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.1 Invariance approach for causal discovery

[1] 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.

where f:R|S∗|→Y. We let F be the function class of f and let FS be the subclass of functions that depend only on the set S⊆{1,…,p} of variables. With this notation we have f∈FS∗.

The assumption of no direct effect of E on Y is analogous to the typical assumptions about instrumental variables [33], [34]. See Section 5 in [1] 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 [1, Proposition 5].

Example (fertility data)

In this work, we analyze a data set provided by the [35]. Here, Y,X and E correspond to the following quantities:

1. Y∈R is the total fertility rate (TFR) in a country in a given year,

2. X∈R9 are potential causal predictor variables for TFR:

   1. IMR – infant mortality rate

   2. Q5 – under-five mortality rate

   3. Education expenditure (% of GNI)

   4. Exports of goods and services (% of GDP)

   5. GDP per capita (constant 2005 US$)

   6. GDP per capita growth (annual %)

   7. Imports of goods and services (% of GDP)

   8. Primary education (% female)

   9. Urban population (% of total)

3. E∈{C1,C2,C3,C4,C5,C6} 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.

Figure 1

Three candidates for a causal DAG with target total fertility rate (TFR) and four potential causal predictor variables. We would like to infer the parents of TFR in the true causal graph. We use the continent as the environment variable E. If the true DAG was one of the two graphs on the left, the environmental variable would have no direct influence on the target variable TFR and ‘Continent’ would be a valid environmental variable, see Definition 1.

The basic yet central insight underlying the invariance approach is the fact that for the true causal parental set S∗:=paY we have the following conditional independence relation under Definition 1 of environmental variables:

(2)Y⫫E|XS∗.

This follows directly from the local Markov condition [e. g. 36]. The goal is to find S∗ by exploiting the above relation (2). Suppose we have a test for the null hypothesis

(3)H0,S:Y⫫E|XS.

It was then proposed in [1] to define an estimate Sˆ for the parental set S∗ by setting

(4)Sˆ:=⋂S:H0,Snot rejectedS.

Here, the intersection runs over all sets S, s.t. E∩S=∅. If the index set is empty, i. e. H0,S is rejected for all sets S, we define Sˆ to be the empty set. If we can test (3) with the correct type I error rate in the sense that

(5)P(H0,S∗is rejected at levelα)≤α,

then we have as immediate consequence the desired statement

P(Sˆ⊆S∗)≥P(H0,S∗accepted)≥1−α.

This follows directly from the fact that S∗ is accepted with probability at least 1−α since H0,S∗ is true; see [1] for details.

In the case of linear models, the method proposed by Peters et al. [1, Eq. (16)] considers a set S as invariant if there exist linear regression coefficients β and error variance σ 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. 1, Eq. (3)]:

Example 1.

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

Y=2X+N,N⫫X,

and inE=2

Y=2X+M,M⫫X,

where M and N have the same mean and variance but differ in higher order moments. In this case, we would haveE⫫̸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 [1]. 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 whichX2andX3are direct causes of Y.

X1←E+ηXX2←3X1+ηX1X3←2X1+ηX2Y←X22−X32+ηY

Due to the nonlinear effect, a linear regression from Y onX2andX3does not yield an invariant model. If we regress Y onX1, 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

X1←E+η1Y←f(X1)+ηYX2←g(X1)+γY+η2

with independent Gaussian error terms. Furthermore, assume that

∀x∈R:f(x)=αx+βh(x)∀x∈R:g(x)=h(x)−γf(x)βγ2−γ=−βvar(η2)/var(ηY)

for someα,βandh:R→R. Then, in the limit of an infinite sample size, the set{X1,X2}is the only set that, after a linear regression, yields residuals that are independent of E. (To see this writeY=f(X1)+ηYas a linear function inX1,X2and show that the covariance between the residuals andX2is zero.) Here, the functions have to be “fine-tuned” in order to make the conditionalY|X1,X2linear inX1andX2.^[4]As an example, one may chooseY←X1+0.5X12+ηYandX2←0.5X12−X1+Y+η2andη1,ηY,η2i.i.d. with distributionN(0,σ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) [38]. This in contrast to nonparametric tests for (unconditional) independence [22], [26]. 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.

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 Sˆ and potentially to an altogether empty set Sˆ. We can instead ask for the defining sets [28], where a defining set Dˆ⊆{1,…,p} has the properties

1. S∩Dˆ≠∅ for all S such that H0,S is accepted.

2. there exists no strictly smaller set D′ with D′⊂Dˆ for which property (i) is true.

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

2.3 Confidence bands

For a given set S, we can in general construct a (1−α)-confidence band FˆS for the regression function when predicting Y with the variables XS. Note that if f is the regression function when regressing Y on the true set of causal variables XS∗ and hence, then, with probability 1−α, we have

Furthermore, from Section 2.1 we know that H0,S∗ is accepted with probability 1−α. We can hence construct a confidence band for the causal effects as

Using a Bonferroni correction, we have the guarantee that

where the coverage guarantee is point-wise or uniform, depending on the coverage guarantee of the underlying estimators FˆS for all given S⊆{1,…,p}.

2.4 Average causal effects

The confidence bands 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 x˜ and x:

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 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∗⊆{1,…,p}. That is—in the absence of hidden variables—we have

and the latter is then equal to

that is it does not matter whether we set the causal variables to a specific value xS∗ 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

with the immediate guarantee that

where the factor 2α is guarding, by a Bonferroni correction, against both a probability α that S∗ will not be accepted—and hence Sˆ⊆S∗ is not necessarily true—and another probability α that the confidence bands will not provide coverage for the parental set S∗.

Figure 2

Data for Nigeria in 1993: The union of the confidence bands FˆS, denoted by Fˆ, bounds the average causal effect of varying the variables in the defining set Dˆ1={IMR, Q5} on the target log(TFR). IMR and Q5 have been varied individually, see panels (a) and (b), as well as jointly, see panel (c), over their respective quantiles. In panels (a) and (b), we do not 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 Dˆ1={IMR, Q5} jointly (see panel (c)), we see that all accepted models predict an increase in expected log(TFR) as IMR and Q5 increase.

Example (fertility data)

The confidence bands Fˆ, required for the computation of ACEˆ(x˜,x), are obtained by a time series bootstrap [39] as the fertility data contain temporal dependencies. The time series bootstrap procedure is described in Appendix I. We use a level of α=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 x˜.

In the first example, we consider the observed covariates for Nigeria in 1993 as x. The point of comparison x˜ is set equal to x, except for the variables in the defining set Dˆ1={IMR, Q5}. In Figures 2(a) and (b), these are varied individually over their respective quantiles. The overall confidence interval Fˆ consists of the union of the shown confidence intervals FˆS. If x=x˜ (shown by the vertical lines), the average causal effect is zero, of course. In neither of the two scenarios shown in Figures 2(a) and (b), 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 Dˆ1={IMR, Q5} jointly (see Figure 2(c)), we see that all accepted models predict an increase in expected log(TFR) as IMR and Q5 increase.

Figure 3

(a) Bounds for the average causal effect of setting the variables IMR and Q5 in the African countries in 2013 to European levels, that is x˜ differs from the country-specific observed values x in that the child mortality rates IMR and Q5 have been set to their respective European average. The implied coverage guarantee is 80 % as we chose α=0.1. (b) Random Forest regression model using all covariates as input. The (non-causal) regression effect coverage is again set to 80 %. We will argue below that the confidence intervals obtained by random forest are too small, see Table 1 and Figure 5.

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 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 3(a) (depicted by the horizontal line segments) correspond to ACEˆ(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 ACEˆIMR+Q5(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 x˜,x of comparison.

For comparison, in Figure 3(b), we show the equivalent computation as in Figure 3(a) when all covariates are assumed to have a direct causal effect on the target and a Random Forest is used for estimation [40]. We observe that while the resulting regression bootstrap confidence intervals often overlap with ACEˆIMR+Q5(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.

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

Figure 4

The confidence intervals show the predicted change in log(TFR) from 1973–2008 for all African countries when not using their data in the nonlinear ICP estimation (using invariant conditional quantile prediction with α=0.1). In other words, only the data from the remaining five continents was used during training. The horizontal line segments mark the union over the accepted models’ intervals for the predicted change; the blue squares show the true change in log(TFR) from 1973–2008. Only those countries are displayed for which the response was not missing in the data, i. e. where log(TFR) in 1973 and in 2008 were recorded. The coverage is 25/26≈0.96.

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 α. For each value of α, we then compute the predicted change in the response log(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.^[5] We then compare the union of the associated confidence intervals with the real, observed change in log(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(TFR) for a country is the observed mean change in log(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(TFR) in 1973 and 2008.

Figure 5

The confidence intervals show the predicted change in log(TFR) from 1973–2008 for all countries when not using the data of the country’s continent in the estimation (with implied coverage guarantee 80 %). Only those countries are displayed for which log(TFR) in 1973 and 2008 were not missing in the data. For nonlinear ICP the shown intervals are the union over the accepted models’ intervals.

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.

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

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.

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 σ 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 XS∗. 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 XS∗) would have to remain constant if we condition on various values of E.

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.

1. The method of Causal Additive Models (CAM) [47] identifies graph structure based on nonlinear additive noise models [6]. 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. We run the PC algorithm [3] 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 [48], [49], using the implementation from [50] (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. We compare against linear ICP [6] where the model misspecification is the assumed linearity of the structural equations.

4. We compare against LiNGAM [15], 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. We also show the outcome of a random selection of the parents that adheres to the FWER-limit by selecting the empty set (Sˆ=∅) with probability 1−α and setting Sˆ={k} for k randomly and uniformly picked from {1,…,p}∖k′ with probability α, where k′ is the index of the current target variable. The random selection is guaranteed to maintain FWER at or below 1−α.

Metrics

Error rates and power are measured in the following by

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

   P(Sˆ⊈S∗).

2. 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 Sˆ and the true set S∗. It is defined as 1 if both S∗=Sˆ=∅ and otherwise as

   |Sˆ∩S∗||Sˆ∪S∗|.

   The Jaccard similarity is thus between 0 and 1 and the optimal value 1 is attained if and only if 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 α=0.05 for all conditional independence tests, that is P(Sˆ⊆S∗)≥1−α. 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 α.

Figure 7

Average Jaccard similarity (y-axis) against average FWER (x-axis), stratified according to which conditional independence test (A)–(F) or baseline method has been used. The nominal level is α=0.05, illustrated by the vertical dotted line. The shown results are averaged over all target variables. Since the empty set is the correct solution for target variable 1 and 5, methods that mostly return the empty set (such as random or linear ICP) perform still quite well in terms of average Jaccard similarity. Since all variables are highly predictive for the target variable Y, see Figure 6, classical variable selection techniques as LASSO have a FWER that lies far beyond α. Importantly, the considered baselines are not suitable for the considered problem setting due to various model misspecifications. We show their performance for comparison to illustrate the influence of these misspecifications.

Figure 8

Average Jaccard similarity over the conditional independence tests (A)–(F) (y-axis) against average FWER (x-axis), stratified according to various parameters (from top left to bottom right): sample size ‘n’, type of nonlinearity ‘id’, ‘target variable’, intervention location ‘interv’, multiplicative effects indicator ‘multiplic’, ‘strength’ of interventions, mean value of interventions ‘meanshift’, shift intervention indicator ‘shift’ and degrees of freedom for t-distributed noise ‘df’. For details, see the description in Appendix III. The FWER is within the nominal level in general for all conditional independence tests. The average Jaccard similarity is mostly determined by the target variable under consideration, see top right panel.

Figure 9

The identical plot to Figure 7 separately for target variables 2, 3, 4 and 6. For all target variables, method (E)(ii)—an invariant residual distribution test using GAM with Levene’s test + Wilcoxon test—performs constantly as good or nearly as good as the optimal method among the considered tests.

A caveat of the FWER control seen in Figure 7 is that while the FWER is maintained at the desired level, the test H0,S∗ might be rejected more often than with probability α. The error control rests on the fact that H0,S∗ is accepted with probability higher than 1−α (since the null is true for S∗). However, if a mistake is made and H0,S∗ is falsely rejected, then we might still have Sˆ⊆S∗ because either all other sets are rejected, too, in which case Sˆ=∅, 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 Sˆ=S∗. A value 0 corresponds to disjoint sets 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∗=∅) and the average Jaccard similarity thus always exceeds 1−α if the level of the procedure is maintained as then Sˆ=∅=S∗ with probability at least 1−α and the Jaccard similarity is 1 if both 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 III, the parents of target variable 3 are difficult to estimate as the paths 1→2→3 and 1→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 Sˆ=∅ when taking the intersection across all accepted sets to compute the estimate 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 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 IV, 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 f4(x)=sin(2π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.