Journal of Machine Learning Research 21 (2020) 1-108

Submitted 3/17; Revised 1/20; Published 3/20

Joint Causal Inference from Multiple Contexts
Joris M. Mooij∗

j.m.mooij@uva.nl

Korteweg-De Vries Institute, University of Amsterdam
Postbox 94248, 1090 GE Amsterdam, The Netherlands

Sara Magliacanesara.magliacane@ibm.com

MIT-IBM Watson AI Lab, IBM Research
75 Binney St, Cambridge, MA 02142, USA

Tom Claassen

tomc@cs.ru.nl

Institute for Computing and Information Sciences, Radboud University Nijmegen
Postbox 9010, 6500 GL Nijmegen, The Netherlands

Editor: Peter Spirtes

Abstract
The gold standard for discovering causal relations is by means of experimentation. Overthe last decades, alternative methods have been proposed that can infer causal relations
between variables from certain statistical patterns in purely observational data. We introduce Joint Causal Inference (JCI), a novel approach to causal discovery from multiple
data sets from different contexts that elegantly unifies both approaches. JCI is a causal
modeling framework rather than a specific algorithm, and it can be implemented using any
causal discovery algorithm that can take into account certain background knowledge. JCI
can deal with different types of interventions (e.g., perfect, imperfect, stochastic, etc.) ina unified fashion, and does not require knowledge of intervention targets or types in caseof interventional data. We explain how several well-known causal discovery algorithms
can be seen as addressing special cases of the JCI framework, and we also propose novel
implementations that extend existing causal discovery methods for purely observational
data to the JCI setting. We evaluate different JCI implementations on synthetic data andon flow cytometry protein expression data and conclude that JCI implementations canconsiderably outperform state-of-the-art causal discovery algorithms.
Keywords: causal discovery, causal modeling, causal inference, observational and experimental data, interventions, randomized controlled trials

1. Introduction
The aim of causal discovery is to learn the causal relations between variables of a system ofinterest from data. As a simple example, suppose a researcher wants to find out whether
playing violent computer games causes aggressive behavior. She gathers observational databy taking a sample from pupils at several high schools in different countries and observes asignificant correlation between the daily amount of hours spent on playing violent computer
games, and aggressive behavior at school (see also Figure 1). This in itself does not yet imply
a causal relation between the two in either direction. Indeed, an alternative explanation of∗. Part of this work was done while the authors were with the Informatics Institute of the University of
Amsterdam.
c 2020 Joris M. Mooij, Sara Magliacane and Tom Claassen.
License: CC-BY 4.0, see https://creativecommons.org/licenses/by/4.0/. Attribution requirements are provided
at http://jmlr.org/papers/v21/17-123.html.

Mooij, Magliacane and Claassen

the observed correlation could be the presence of a confounder (a latent common cause),for example, a genetic predisposition towards violence that makes the carrier particularly
enjoy such games and also make him behave more aggressively. The most reliable way toestablish whether playing violent computer games causes aggressive behavior, is by meansof experimentation, for example by a randomized controlled trial (Fisher, 1935). This would
imply assigning each pupil to one out of two groups randomly, where the pupils in one groupare forced to play violent computer games for several hours a day, while the pupils in theother group are forced to abstain from playing those games. After several months, theaggressive behavior in both groups is measured. If a significant correlation between groupand outcome is observed (or equivalently, the outcome is significantly different between thetwo groups), it can then be concluded that playing violent computer games indeed causes
aggressive behavior.
Given the ethical and practical problems that such an experiment would involve, one
might wonder whether there are alternative ways to answer this question. One such alternative is to combine data from different contexts. For example, in some countries thegovernment may have decided to forbid certain ultra-violent games from being sold. In
addition, some schools may have introduced certain measures to discourage aggressive behavior. By combining the data from these different contexts in an appropriate way, one may
be able to identify the presence or absence of a causal effect of playing violent computer
games on aggressive behavior. For example, in the setting of Figure 1(c), the causal relationship between the two variables of interest turns out to be identifiable from conditional
independence relationships in pooled data from all the contexts. In particular, in that casethe observed correlation between playing violent computer games and aggressive behavior
could be unambiguously attributed to a causal effect of one on the other, just from combining multiple readily available data sets, without the need for an impractical experiment.1
In this paper, we propose a simple and general way to combine and analyze data sets fromdifferent contexts that enables one to draw such strong causal conclusions.
While experimentation is still the gold standard to establish causal relationships, researchers realized in the early nineties that there are other methods that require only purely
observational data (Spirtes et al., 2000; Pearl, 2009). Many methods for causal discovery
from purely observational data have been proposed over the last decades, relying on different
assumptions. These can be roughly divided into constraint-based causal discovery methods,
such as the PC (Spirtes et al., 2000), IC (Pearl, 2009) and FCI algorithms (Spirtes et al.,
1999; Zhang, 2008a), score-based causal discovery methods (e.g., Cooper and Herskovits,
1992; Heckerman et al., 1995; Chickering, 2002; Koivisto and Sood, 2004), and methods
exploiting other statistical patterns in the joint distribution (e.g., Mooij et al., 2016; Peters
et al., 2017). Originally, these methods were designed to estimate the causal graph of thesystem from a single data set corresponding to a single (purely observational) context.
More recently, various causal discovery methods have been proposed that extend thesetechniques to deal with multiple data sets from different contexts. As an example, the data
sets may correspond with a baseline of purely observational data consisting of measurements
concerning the “natural” state of the system, and data consisting of measurements under
1. One can show that the conditional dependence Cα ⊥6 ⊥ X2 | Cβ and conditional independence
Cα ⊥⊥ X2 | {X1 , Cβ } in the pooled data that are entailed by the causal graph, together with the assumption that neither Cα nor Cβ is caused by X1 or X2 , suffice to arrive at this conclusion.

2

Joint Causal Inference from Multiple Contexts

(a)

(b)
X1

X2

(c)
X1

Cα

Cβ

X1

X2

X2

Figure 1: Different causal graphs relating X1 , the daily amount of hours spent on playing
violent computer games, and X2 , a measure of aggressive behavior. (a) Playing
violent computer games causes aggressive behavior; (b) The observed correlation
between X1 and X2 is explained by a latent confounder, e.g., a genetic predisposition towards violence. (c) Hypothetical causal graph also involving context
variables Cα , which indicates whether ultra-violent games have been banned by
the government, and Cβ , which represents school interventions to stimulate social
behavior. Without considering contexts, it is not possible to distinguish between(a) and (b) based on conditional independences in the data. In scenario (c),JCI allows one to infer from conditional independences in the pooled data that
X1 causes X2 and that X1 and X2 are not confounded (assuming that context
variables Cα and Cβ are not caused by system variables X1 and X2 ).

different perturbations of the system due to external interventions on the system.2 Moregenerally, they can correspond to measurements of the system in different environments.
These methods can be divided into two main approaches:(a) methods that obtain statistics or constraints from each context separately and thenconstruct a single context-independent causal graph by combining these statistics, butnever directly compare data from different contexts (Claassen and Heskes, 2010; Tillman
and Spirtes, 2011; Hyttinen et al., 2012, 2014; Triantafillou and Tsamardinos, 2015;
Rothenhäusler et al., 2015; Forré and Mooij, 2018);(b) methods that pool all data and construct a single context-independent causal graph directly from the pooled data (Cooper, 1997; Cooper and Yoo, 1999; Tian and Pearl, 2001;
Sachs et al., 2005; Eaton and Murphy, 2007; Chen et al., 2007; Hauser and Bühlmann,
2012; Mooij and Heskes, 2013; Peters et al., 2016; Oates et al., 2016a; Zhang et al.,2017).
In this paper, we propose Joint Causal Inference (JCI), a framework for causal modelingof a system in different contexts and for causal discovery from multiple data sets consistingof measurements obtained in different contexts, which takes the latter approach. As will
be discussed in more detail in Section 4.3, JCI is the most generally applicable of those
approaches—for example, it allows for the presence of latent confounders and cyclic causal
2. In certain parts of the causal discovery literature, the word “intervention” has become synonymous to“perfect intervention” (i.e., an intervention that precisely sets a variable or set of variables to a certain
value without directly affecting any other variables in the system), but in this work we use it in the moregeneral meaning of any external perturbation of the system.

3

Mooij, Magliacane and Claassen

relationships—and also offers most flexibility in terms of its implementation. While theingredients of the JCI framework are not novel, the added value of the framework is that on
the one hand it arrives at a unifying description of a diverse spectrum of existing approaches,
while on the other hand it serves to inspire new implementations, such as the adaptations of
FCI that we propose in this work. Technically, this is achieved by formulating the problem
in terms of a (standard) Structural Causal Model that considers system and environment
as subsystems of one joint system, rather than other types of representations in which thesystem is modeled conditionally on its environment (Dawid, 2002; Bareinboim and Pearl,2013; Oates et al., 2016a; Yang et al., 2018; Forré and Mooij, 2019). This allows us toapply the standard notion of statistical independence in the same ways as is commonly
done in the purely observational setting. As we observed in our experiments (that arereported in Section 5), the novel algorithms proposed in this work compare favorably withthe state-of-the-art in causal discovery on synthetic data in many settings.
The key idea of JCI is to (i) consider auxiliary context variables that describe the context
of each data set, (ii) pool all the data from different contexts, including the values of thecontext variables, into a single data set, and finally (iii) apply standard causal discovery
methods to the pooled data, incorporating appropriate background knowledge on the causal
relationships involving the context variables. The framework is simple and very generally
applicable as it allows one to deal with latent confounding and cycles (if the causal discovery
method supports this) and various types of interventions in a unified way. It does not require
background knowledge on the intervention types and targets, making it very suitable to theapplication on complex systems in which the effects of certain interventions are not knowna priori, a situation that often occurs in practice. On the other hand, if such backgroundknowledge is available, it can be exploited.
JCI can be implemented using any causal discovery method that can incorporate theappropriate background knowledge on the relationships between context and system variables. This allows one to benefit from the availability of sophisticated and powerful causal
discovery methods that have been primarily designed for a single data set from a single
context by extending their application domain to the setting of multiple data sets frommultiple contexts. For example, we will show in this work how FCI (Spirtes et al., 1999;
Zhang, 2008a) can easily be adapted to the JCI setting. At the same time, JCI accommodates various well-known causal discovery methods as special cases, such as the standard
randomized controlled trial setting (Fisher, 1935), Local Causal Discovery (LCD) (Cooper,
1997) and Invariant Causal Prediction (ICP) (Peters et al., 2016). By explicitly introducing the context variables and treating them analogously to the system variables (but withadditional background knowledge about their causal relations with the system variables),
JCI makes it possible to elegantly combine the principles of causal discovery from experimentation with those of causal discovery from purely observational data to achieve a causal
discovery framework that is more powerful than either of the two separately.
This paper is structured as follows. In Section 2 we describe the relevant causal modelingand discovery concepts and define terminology and notation. In Section 3 we introduce
the JCI framework and modeling assumptions. In Section 4, we show how JCI can be
implemented using various causal discovery methods, and compare it with related work. In
Section 5 we report experimental results on synthetic and flow cytometry data. We conclude
in Section 6 with some promising directions for future developments.
4

Joint Causal Inference from Multiple Contexts

2. BackgroundIn this section, we present the background material on which we will base our exposition.
We start in Section 2.1 with a brief subsection stating the basic definitions and resultsin the field of graphical causal modeling that we will use in this paper. In addition tocovering material that is standard in the field, we review more recent extensions to thecyclic setting (Bongers et al., 2020). Because the cyclic setting is quite similar to theacyclic one that is mostly considered in the literature, we decided to present both cases inparallel rather than first explaining the acyclic setting and then explaining how everything
generalizes to the cyclic setting.3 In Section 2.2, we discuss the key idea of causal discovery
from experimentation (in the setting of a randomized controlled trial, or A/B-testing) inthese terms. We finish with Section 2.3 that briefly illustrates the basic idea underlying
constraint-based causal discovery from purely observational data in a simple setting.
2.1. Graphical Causal Modeling
We briefly summarize some basic definitions and results in the field of graphical causal
modeling. For more details, we refer the reader to Pearl (2009) and Bongers et al. (2020).
2.1.1. Directed Mixed Graphs
A Directed Mixed Graph (DMG) is a graph G = hV, E, Fi with nodes V and two types ofedges: directed edges E ⊆ V 2 , and bidirected edges F ⊆ {{i, j} : i, j ∈ V, i 6= j}. We will
denote a directed edge (i, j) ∈ E as i → j or j ← i, and call i a parent of j and j a child of i.
We denote all parents of j in the graph G as paG (j) := {i ∈ V : i → j ∈ E}, and all childrenof i in G as chG (i) := {j ∈ V : i → j ∈ E}. We allow for self-cycles i → i, so a variable can beits own parent and child. We will denote a bidirected edge {i, j} ∈ F as i ↔ j or j ↔ i, andcall i and j spouses. Two nodes i, j ∈ V are called adjacent in G if they are connected by anedge (or multiple edges), i.e., if i → j ∈ E or i ← j ∈ E or i ↔ j ∈ F. For a subset of nodes
W ⊆ V, we define the induced subgraph GW := (W, E ∩ W 2 , F ∩ {{i, j} : i, j ∈ W, i 6= j}),i.e., with nodes W and exactly those edges of G that connect nodes in W.
A walk between i, j ∈ V is a tuple hi0 , e1 , i1 , e2 , i3 , . . . , en , in i of alternating nodes andedges in G (n ≥ 0), such that all i0 , . . . , in ∈ V, all e1 , . . . , en ∈ E ∪ F, starting with
node i0 = i and ending with node in = j, and such that for all k = 1, . . . , n, the edge
ek connects the two nodes ik−1 and ik in G. If the walk contains each node at mostonce, it is called a path. A trivial walk (path) consists just of a single node and zeroedges. A directed walk (path) from i ∈ V to j ∈ V is a walk (path) between i and jsuch that every edge ek on the walk (path) is of the form ik−1 → ik , i.e., every edge isdirected and points away from i. By repeatedly taking parents, we obtain the ancestors of
j: anG (j) := {i ∈ V : i = i0 → i1 → i2 → · · · → in = j in G}. Similarly, we define thedescendants of i: deG (i) := {j ∈ V : i = i0 → i1 → i2 → · · · → in = j in G}. In particular,
each node is ancestor and descendant of itself. A directed cycle is a directed path from ito j such that in addition, j → i ∈ E. An almost directed cycle is a directed path from
3. The disadvantage is that our notation and definitions deviate somewhat from those commonly used in theacyclic causal discovery literature. Therefore, we recommend reading this section also to those readersthat are already familiar with the theory of acyclic structural causal models.

5

Mooij, Magliacane and Claassen

i to j such that in addition, j ↔ i ∈ F. All nodes on directed cycles passing throughi ∈ V together form the strongly-connected component scG (i) := anG (i) ∩ deG (i) of i. We
extend the definitions to sets I ⊆ V by setting anG (I) := ∪i∈I anG (i), and similarly fordeG (I) and scG (I). A directed mixed graph G is acyclic if it does not contain any directed
cycle, in which case it is known as an Acyclic Directed Mixed Graph (ADMG). A directed
mixed graph that does not contain bidirected edges is known as a Directed Graph (DG). If a
directed mixed graph does not contain bidirected edges and is acyclic, it is called a Directed
Acyclic Graph (DAG).
A node ik on a walk (path) π = hi0 , e1 , i1 , e2 , i3 , . . . , en , in i in G is said to form a collider
on π if it is a non-endpoint node (1 ≤ k < n) and the two edges ek , ek+1 meet head-to-headon their shared node ik (i.e., if the two subsequent edges are of the form ik−1 → ik ← ik+1 ,ik−1 ↔ ik ← ik+1 , ik−1 → ik ↔ ik+1 , or ik−1 ↔ ik ↔ ik+1 ). Otherwise (that is, if
it is an endpoint node, i.e., k = 0 or k = n, or if the two subsequent edges are of theform ik−1 → ik → ik+1 , ik−1 ← ik ← ik+1 , ik−1 ← ik → ik+1 , ik−1 ↔ ik → ik+1 , orik−1 ← ik ↔ ik+1 ), ik is called a non-collider on π. We will denote the colliders on a walk
π as col(π) and the non-colliders on π (including the endpoints of π) as ncol(π). A tripleof nodes hi, j, ki in G is called an unshielded triple if i is adjacent to j, j is adjacent to k
and i is not adjacent to k in G.
2.1.2. Structural Causal Models
Directed Mixed Graphs form a convenient graphical representation for variables (labelledby the nodes) and their functional relations (expressed by the edges) in a Structural Causal
Model (SCM) (Pearl, 2009), also known as a (non-parametric) Structural Equation Model(SEM) (Wright, 1921). Several slightly different definitions of SCMs have been proposed in
the literature, which all have their (dis)advantages. Here we use a variant of the definition
in Bongers et al. (2020) that is most convenient for our purposes. The reason we use SCMs
to formulate JCI (rather than for example the more well-known causal Bayesian networks)
is that SCMs are expressive enough to model both latent common causes and cyclic causal
relationships.
Definition 1 A Structural Causal Model (SCM) is a tuple M = hI, J , H, X , E, f , PE i of:(i) a finite index set I for the endogenous variables in the model;(ii) a finite index set J for the latent exogenous variables in the model (disjoint from I);(iii) a directed graph H with nodes I Q∪ J , and directed edges pointing from I ∪ J to I;
4(iv) a product of Borel spaces X = i∈I Xi , which define the domains of the endogenous
variables;Q(v) a product of Borel spaces E = j∈J Ej , which define the domains of the exogenous
variables;Q(vi) a product probability measure PE = j∈J PEj on E specifying the exogenous distribution;(vii) a measurable function f : X × E → X , the causal mechanism, such that each of its
components fi only depends on a particular subset of the variables, as specified by the
4. A Borel space is both a measurable and a topological space, such that the sigma-algebra is generated by
the open sets. Most spaces that one encounters in applications as the domain of a random variable are(isomorphic to) Borel spaces.

6

Joint Causal Inference from Multiple Contexts

Intervened SCM

Interventional Distribution

Simple SCM
Marginal SCM
Augmented Graph

Observational Distribution
Faithfulness

Latent Confounders

Graph

Direct Causes

d/σ-separations

Causal Relations

Markov
Property(Conditional) Independences
Figure 2: Relationships between various representations of simple SCMs. Directed edges
represent mappings. Intervened and marginal SCMs are always defined and arealso simple.

directed graph H:
fi : X paH (i)∩I × E paH (i)∩J → Xi ,i ∈ I.

In discussing the concepts and properties of SCMs, the graphical representation of various
objects and their relations in Figure 2 may be helpful. This shows how the SCM is thebasic object containing all information, and how other representations can be derived fromthe SCM. In the rest of this section, we will discuss this in more detail.
We refer to the graph H in Definition 1(iii) as the augmented graph of M. In contrast, thegraph of M, denoted G(M), is the directed mixed graph with nodes I, directed edges i1 → i2
iff i1 → i2 ∈ H, and bidirected edges i1 ↔ i2 iff there exists j ∈ paH (i1 )∩paH (i2 )∩J .5 While
the augmented graph H shows in detail the functional dependence of endogenous variableson the (independent) exogenous variables, the graph G(M) provides an abstraction by notincluding the exogenous variables explicitly, but using bidirected edges to represent any
shared dependence of pairs of endogenous variables on a common exogenous parent. If
G(M) is acyclic, we call the SCM M acyclic, otherwise we call the SCM cyclic. If G(M)contains no bidirected edges, we call the endogenous variables in the SCM M causally
sufficient.
A pair of random variables (X, E) is called a solution of the SCM M if X = (Xi )i∈I
with Xi ∈ Xi for all i ∈ I, E = (Ej )j∈J with Ej ∈ Ej for all j ∈ J , the distribution P(E)
5. This definition of graph makes a slight simplification: a more precise definition would leave out edges
that are redundant. For example, if the structural equation for X2 reads X2 = 0 · X1 + X3 it could bethat 1 → 2 ∈ H, but this edge would not appear in G(M). For the rigorous version of this definition,
see Bongers et al. (2020).

7

Mooij, Magliacane and Claassen

is equal to the exogenous distribution PE , and the structural equations:
Xi = fi (XpaH (i)∩I , EpaH (i)∩J )a.s.

hold for all i ∈ I. An SCM is often specified informally by specifying only the structural
equations and the density6 of the exogenous distribution with respect to some productmeasure, for example:(
Xi = fi (XpaH (i)∩I , EpaH (i)∩J ),i ∈ I,
M:
Q
p(E) = j∈J p(Ej ).
For acyclic SCMs, solutions exist and have a unique distribution that is determined by
the SCM. This is not generally the case in cyclic SCMs, as these could have no solution atall, or could have multiple solutions with different distributions (Bongers et al., 2020).
Definition 2 An SCM M is said to be uniquely solvable w.r.t. O ⊆ I if there exists a
measurable mapping gO : X (paH (O)\O)∩I × E paH (O)∩J → X O such that for PE -almost every
e for all x ∈ X :xO = gO (x(paH (O)\O)∩I , epaH (O)∩J )⇐⇒xO = fO (x, e) .

(Loosely speaking: the structural equations for O have a unique solution for XO in termsof the other variables appearing in those equations.)
If M is uniquely solvable with respect to I (in particular, this holds if M is acyclic), thenit induces a unique observational distribution PM (X).
Given an SCM that models a certain system, we can model the system after an idealizedintervention in which an external influence enforces a subset of endogenous variables to takeon certain values, while leaving the rest of the system untouched.
Definition 3 Let M be an SCM. The perfect intervention with target I ⊆ I and valueξI ∈ X I induces the intervened SCM Mdo(I,ξI ) obtained by copying M, but letting H̃ be H
without the edges {j → i ∈ H : j ∈ I ∪ J , i ∈ I}, and modifying the causal mechanism intof˜ such that(ξii∈If˜i (x, e) =fi (x, e) i ∈/ I.
The interpretation is that the causal mechanisms that normally determine the values of thecomponents i ∈ I are replaced by mechanisms that assign the values ξi . Other types ofinterventions are possible as well (see also Section 3.3). If the intervened SCM M
 do(I,ξI )induces a unique observational distribution, this is denoted as PM X | do(I, ξI ) and referred to as the interventional distribution of M under the perfect intervention do(I, ξI ).
Pearl (2009) derived the do-calculus for acyclic SCMs, consisting of three rules that express
relationships between interventional distributions of an SCM.
6. We denote a probability measure (or distribution) of a random variable X by P(X), and a density of X
with respect to some fixed product measure by p(X).

8

Joint Causal Inference from Multiple Contexts

2.1.3. Simple Structural Causal Models
The theory of general cyclic Structural Causal Models is rather involved (Bongers et al.,2020). In this work, for simplicity of exposition, we will focus on a certain subclass of SCMs
that has many convenient properties and for which the theory simplifies considerably:
Definition 4 An SCM M is called simple if it is uniquely solvable with respect to any
subset O ⊆ I.
All acyclic SCMs are simple. Simple SCMs provide a special case of the more general classof modular SCMs (Forré and Mooij, 2017). The class of simple SCMs can be thought of asa generalization of acyclic SCMs that allows for (weak) cyclic causal relations, but preserves
many of the convenient properties that acyclic SCMs have.
Indeed, a simple SCM induces a unique observational distribution. Its marginalizations
are always defined (Bongers et al., 2020), and are also simple; in other words, the classof simple SCMs is closed under marginalizations. The class of simple SCMs is also closed
under perfect interventions, and hence, all perfect interventional distributions of a simple
SCM are uniquely defined. Without loss of generality, one can assume that simple SCMs
have no self-cycles. The causal interpretation of the graph of an SCM with cycles and/or
bidirected edges can be rather subtle in general. However, for graphs of simple SCMs there
is a straightforward causal interpretation:
Definition 5 Let M be a simple SCM. If i → j ∈ G(M) we call i a direct cause of jaccording to M. If there exists a directed path i → · · · → j ∈ G(M), i.e., if i ∈ anG(M) (j),then we call i a cause of j according to M. If there exists a bidirected edge i ↔ j ∈ G(M),then we call i and j confounded according to M.
We conclude that the graph G(M) of a simple SCM can be interpreted as its causal
graph. In the next subsection, we will discuss how the same graph G(M) of a simple
SCM M also represents the conditional independences that must hold in the observational
distribution of M.
2.1.4. Structural Causal Models: Markov Properties
Under certain conditions, the graph G(M) of an SCM M can be interpreted as a statisticalgraphical model, i.e., it allows one to read off conditional independences that must holdin the observational distribution PM (X). One of the most common formulations of such
Markov properties involves the following notion of d-separation, first proposed by Pearl(1986) in the context of DAGs, and later shown to be more generally applicable:7
Definition 6 (d-separation) We say that a walk hi0 . . . in i in DMG G = hV, E, Fi is dblocked by C ⊆ V if:(i) its first node i0 ∈ C or its last node in ∈ C, or(ii) it contains a collider ik ∈/ anG (C), or(iii) it contains a non-collider ik ∈ C.
If all paths in G between any node in set A ⊆ V and any node in set B ⊆ V are d-blocked
by a set C ⊆ V, we say that A is d-separated from B by C, and we write A ⊥dG B | C.
7. It is also sometimes called “m-separation” in the ADMG literature.

9

Mooij, Magliacane and Claassen

In the general cyclic case, however, the notion of d-separation is too strong, as was alreadypointed out by Spirtes (1994). A solution is to replace it with a non-trivial generalizationof d-separation, known as σ-separation (Forré and Mooij, 2017):
Definition 7 (σ-separation) We say that a walk hi0 . . . in i in DMG G = hV, E, Fi is σblocked by C ⊆ V if:(i) its first node i0 ∈ C or its last node in ∈ C, or(ii) it contains a collider ik ∈/ anG (C), or(iii) it contains a non-collider ik ∈ C that points to a neighboring node on the walk inanother strongly-connected component (i.e., ik−1 → ik → ik+1 or ik−1 ↔ ik → ik+1
with ik+1 ∈/ scG (ik ), ik−1 ← ik ← ik+1 or ik−1 ← ik ↔ ik+1 with ik−1 ∈/ scG (ik ), or
ik−1 ← ik → ik+1 with ik−1 ∈/ scG (ik ) or ik+1 ∈/ scG (ik )).
If all paths in G between any node in set A ⊆ V and any node in set B ⊆ V are σ-blocked
by a set C ⊆ V, we say that A is σ-separated from B by C, and we write A ⊥σG B | C.
Forré and Mooij (2017) proved the following fundamental result for modular SCMs, which
we formulate here only for the special case of simple SCMs:
Theorem 8 (Generalized Directed Global Markov Property) Any solution (X, E)
of a simple SCM M obeys the Generalized Directed Global Markov Property with respect
to the graph G(M):σA ⊥ B | C =⇒ XA
G(M)

⊥
⊥
PM (X)

XB | XC∀A, B, C ⊆ I.

The following stronger Markov properties, in which σ-separation is replaced by the morefamiliar notion of d-separation, have been derived for special cases by Forré and Mooij(2017) (where again we consider only the special case of simple SCMs):
Theorem 9 (Directed Global Markov Property) Let M = hI, J , H, X , E, f , PE i bea simple SCM. If M satisfies at least one of the following three conditions:(i) M is acyclic;(ii) all endogenous spaces Xi are discrete;(iii) M is linear (i.e., Xi = R for each i ∈ I, Ej = R for each j ∈ J , and each causal
mechanism fi : X paH (i)∩I × EpaH (i)∩J → Xi is linear), each causal mechanism fi
depends non-trivially on some exogenous variable(s), and its exogenous distribution
has a density p(E) with respect to Lebesgue measure;then any solution (X, E) of M obeys the Directed Global Markov Property with respect tothe graph G(M):dA ⊥ B | C =⇒ XA
G(M)

⊥
⊥
PM (X)

XB | XC∀A, B, C ⊆ I.

Of these cases, the acyclic and linear cases are well-known.8
8. The acyclic case was first shown in the context of linear-Gaussian structural equation models (Spirtes
et al., 1998; Koster, 1999). The discrete case fixes the erroneous theorem by Pearl and Dechter (1996),
for which a counterexample was found by Neal (2000), by adding the unique solvability condition, andextends it to allow for latent common causes. The linear case extends existing results for the linearGaussian setting without latent common causes (Spirtes, 1994, 1995; Koster, 1996) to a linear (possibly
non-Gaussian) setting with latent common causes.

10

Joint Causal Inference from Multiple Contexts

We conclude that simple SCMs also have convenient Markov properties. A simple
SCM induces a unique observational distribution that satisfies the Generalized Directed
Global Markov Property; under additional conditions, it satisfies even the Directed Global
Markov Property. Similarly, for any perfect intervention, a simple SCM induces a uniqueinterventional distribution that satisfies the (Generalized) Directed Global Markov Property
with respect to the intervened graph. We conclude that the graph G(M) of a simple SCM
has two interpretations: it expresses both the causal structure between the variables as wellas the conditional independence structure of the solutions. These two interpretations of thegraph G(M) of a simple SCM can be combined into a causal do-calculus (Forré and Mooij,
2019) that extends the acyclic do-calculus of Pearl (2009) to the class of simple (or moregenerally, modular) SCMs.
The starting point for constraint-based approaches to causal discovery from observational data is to assume that the data is modelled by an (unknown) SCM M, such thatits observational distribution PM (X) exists and satisfies a Markov property with respect
to its graph G(M). In addition, one usually assumes the faithfulness assumption to hold
(Spirtes et al., 2000; Pearl, 2009), i.e., that the graph explains all conditional independences
present in the observational distribution. For the cases in which the d-separation criterion
Theorem 9 applies, this amounts to assuming the following implication:dA ⊥ B | C ⇐= XA
G(M)

⊥
⊥
PM (X)

XB | XC∀A, B, C ⊆ V.

Meek (1995) has shown completeness properties of d-separation. More specifically, Meek
(1995) showed that faithfulness holds generically for DAGs if (i) all variable domains are
finite, or (ii) if all variables are real-valued, linearly related and have a multivariate Gaussian
distribution. This in particular provides some justification for assuming faithfulness. On
the other hand, no completeness results are known yet for the general cyclic case in whichthe σ-separation criterion Theorem 8 applies. Nevertheless, we believe that such results canbe shown, and we will assume for simple SCMs a similar faithfulness assumption as for thed-separation case:σA ⊥ B | C ⇐= XA
G(M)

⊥
⊥
PM (X)

XB | XC∀A, B, C ⊆ V.

2.2. Causal Discovery by Experimentation
The gold standard for causal discovery is by means of experimentation. For example,
randomized controlled trials (Fisher, 1935) form the foundation of modern evidence-based
medicine. In engineering, A/B-testing is a common protocol to optimize certain causal
effects of an engineered system. Toddlers learn causal representations of the world throughplayful experimentation.
We will discuss here the simplest randomized controlled trial setting by formulating it interms of the graphical causal terminology introduced in the last section. The experimental
procedure is as follows. Consider two variables, “treatment” C and “outcome” X. In thesimplest setting, one considers a binary treatment variable, where C = 1 corresponds to“treat with drug” and C = 0 corresponds to “treat with placebo”. For example, the drug
could be aspirin, and outcome could be the severity of headache perceived two hours later.
11

Mooij, Magliacane and Claassen(a) Two separate data sets:
Placebo(C = 0):
X
-0.2
0.6
-1.7
...

(b) Pooled data:
C
0
0
0
0
1
1
1
1

Drug
(C = 1):
X
-0.3
1.8
-0.1
...

X
-0.2
0.6
-1.7
...
-0.3
1.8
-0.1
...

Figure 3: Illustration of the data from an example randomized controlled trial. The data
can either be interpreted as (a) two separate data sets, one for the treatment
and one for the control group, or (b) as a single data set including a context
variable indicating treatment/control. Note that in this particular example, C isdependent on X in the pooled data (or equivalently, the distribution of X differs
between contexts C = 0 and C = 1), which implies that C is a cause of X.

Patients are split into two groups, the treatment and the control group, by means of acoin flip that assigns a value of C to every patient.9 Patients are treated depending on
the assigned value of C, i.e., patients in the treatment group are treated with the drug
and patients in the control group are treated with a placebo. Some time after treatment,
the outcome X is measured for each patient. This yields a data set (Cn , Xn )N
n=1 with twothmeasurements (Cn , Xn ) for the n patient. If the distribution of outcome X significantly
differs between the two groups, one concludes that treatment is a cause of outcome.
The important underlying causal assumptions that ensure the validity of the conclusionare:(i) outcome X is not a cause of treatment C (which is commonly deemed justified if theoutcome is an event that occurs later in time than the treatment event);(ii) there is no latent confounder of treatment and outcome (this is where the randomization comes in: if treatment is decided solely by a proper coin flip, then it seems
reasonable to assume that there cannot be any latent common cause of the coin flip
C and the outcome X that is not just a combination of two statistically independent
separate causes of C and X),(iii) no selection bias is present in the data (in other words, no data is missing; for example,
if only those patients that did not suffer from certain treatment side effects are includedin the data set, then this assumption will be violated).
Under these assumptions, one can show that if the distribution of the outcome X differs
between the two groups of patients (“treatment group” with C = 1 vs. “control group”
with C = 0), then treatment must be a cause of outcome, at least in this population ofpatients (see Proposition 10). There are two conceptually slightly different ways of testing
9. Usually this is done in a double-blind way, so that neither the patient nor the doctor knows which groupa patient has been assigned to.

12

Joint Causal Inference from Multiple Contexts

this in the data, depending on whether we treat the data as a single pooled data set, or
rather as two separate data sets (each one corresponding to a particular patient group), seealso Figure 3. If we consider the data about outcome X in the two groups as two separate
data sets (corresponding to the same variable X, but measured in different contexts C),
then the question is whether the distribution of X is statistically different in the two data
sets. This can be tested with a two-sample test, for example, a t-test or a Wilcoxon
test. The other alternative is to consider the data as a single pooled data set (by poolingthe data for the two groups), and let the value of C indicate the context of each sample(treatment or control). The question now becomes whether the conditional distribution of
X given C = 0 differs from the conditional distribution of X given C = 1, i.e., whether
P(X | C = 0) 6= P(X | C = 1). In other words, we have to test whether there is a statisticallysignificant dependence C ⊥6 ⊥ X in the pooled data between treatment C and outcome X;
if there is, it must be due to the treatment C causing the outcome X, as the following
proposition shows:
Proposition 10 Suppose that the data-generating process on context variable C and outcome variable X can be modeled by a simple SCM M and no selection bias is present.10
Under the randomized controlled trial assumptions:(i) C ← X ∈/ G(M) (“outcome X is not a cause of treatment C”)(ii) C ↔ X ∈/ G(M) (“there is no latent confounder of treatment C and outcome X”),
a dependence C ⊥6 ⊥ X in the joint distribution P(C, X) implies that C causes X. Furthermore, the causal effect of C on X is given by:

PM X | do(C = c) = PM (X | C = c).
(1)
Proof Out of the eight possible graphs G(M), only two satisfy the assumptions:

C

X

C

X

C

X

C

X

X

C

X

C

X

C

X

C ⊥
⊥X

C

By the Markov property (Theorem 8), if the edge C → X were absent in G(M), then C
would be independent of X. Therefore, if C ⊥6 ⊥ X, the edge C → X must be in G(M). In
both cases, the causal do-calculus applied to G(M) yields the identity (1).
Of course, in this straightforward example the equivalence between the two approaches(differences between two separate data sets vs. properties of a single pooled data set) is
trivial, and the reader may wonder why we emphasize it. The reason is that the key
10. The context variable C is here considered as an endogenous variable in the SCM, as explained in Section 3.1.

13

Mooij, Magliacane and Claassen

idea of our approach is precisely this: reducing an apparently complicated causal discovery
problem with multiple data sets to a more standard causal discovery problem involving asingle pooled data set. The Joint Causal Inference framework that we propose in this paper
can be considered as an extension of this randomized controlled trial setting to multiple
treatment and outcome variables.
It is important to realize that the simple causal reasoning for the RCT cannot bemade when looking at the two data sets in isolation (i.e., by considering only properties ofP(X | C = 0) and P(X | C = 1) separately, and not using in addition any other properties of
the joint distribution P(X, C)). The latter approach is commonly used by constraint-based
methods for causal discovery from multiple data sets (e.g., Tillman, 2009; Claassen and Heskes, 2010; Tillman and Spirtes, 2011; Hyttinen et al., 2014; Triantafillou and Tsamardinos,
2015; Rothenhäusler et al., 2015; Forré and Mooij, 2018). Under the assumptions made,the crucial (and possibly very strong) signal in the data that allows one to draw the conclusion that C causes X is the dependence C ⊥6 ⊥ X that can only be seen in the pooled data.
Methods that only test for conditional independences within each context and subsequently
combine these into a single context-independent causal model will not yield any conclusionin this setting. The approach taken by JCI, on the other hand, is to analyze the pooleddata jointly, so that informative signals like these can be taken into account.
2.3. Causal Discovery from Purely Observational Data
In the previous section, we discussed the current gold standard for discovering causal relations. Over the last two decades, alternative methods have been proposed to performcausal discovery from purely observational data. This is intriguing and of high relevance,since experiments may be impossible, infeasible, impractical, unethical or too expensive
to perform. These causal discovery methods can be divided into constraint-based causal
discovery methods, such as the PC (Spirtes et al., 2000), IC (Pearl, 2009) and FCI algorithms (Spirtes et al., 1999; Zhang, 2008a), and score-based causal discovery methods (e.g.,
Heckerman et al., 1995; Chickering, 2002; Koivisto and Sood, 2004). The PC and IC algorithms and most score-based methods assume causal sufficiency (i.e., the absence of latent
confounders), while the FCI algorithm and other modern constraint-based algorithms allow
for latent confounders and selection bias. Originally, these methods have been designed toestimate the causal graph of the system from a single data set corresponding to a single(purely observational) context.
All these methods try to infer causal relationships on the basis of subtle statistical
patterns in the data. The most important of these patterns are conditional independences
between variables. These are exploited by most constraint-based methods, and implicitly,
by score-based methods. Other patterns, such as “Verma constraints” (Shpitser et al.,2014), algebraic constraints in the linear-Gaussian case (van Ommen and Mooij, 2017),
non-Gaussianity in linear models (Kano and Shimizu, 2003), and non-additivity of noise
in nonlinear models (Peters et al., 2014) can also be exploited. Another class of methods
that has become popular more recently are methods that try to infer the causal direction(A → B vs. B → A) from purely observational data of variable pairs (see e.g., Mooij et al.,2016).
14

Joint Causal Inference from Multiple Contexts

Since our main goal is to enable constraint-based causal discovery from multiple contexts,
we will focus on this approach here, while noting that the JCI framework that we propose in
the next section is compatible with all approaches to causal discovery from purely observational data that allow for multiple variables and can handle certain background knowledge(to be made precise in Section 3.4).
As discussed in detail by Spirtes et al. (2000), causal discovery from conditional independence patterns in purely observational data becomes possible under strong assumptions.
The simplest example of how certain patterns of conditional independences in the observational distribution can lead to conclusions about the causal relations of the variables isgiven by the “Y-structure” pattern (Mani, 2006), which is illustrated in Figure 4. We show
here that the Y-structure pattern also generalizes to the cyclic case.
Proposition 11 Suppose that the data-generating process on four variables X1 , X2 , X3 , X4
can be modeled by a simple SCM M. Assume that the sampling procedure is not subject
to selection bias, and that faithfulness holds. If the following conditional (in)dependencies
hold in the observational distribution PM (X):
X1 ⊥6 ⊥ X4 ,X2 ⊥6 ⊥ X4 ,X1 ⊥⊥ X2 ,X1 ⊥⊥ X4 | X3 , X2 ⊥⊥ X4 | X3 , X1 ⊥6 ⊥ X2 | X3 ,then X3 is a direct cause of X4 according to M. Furthermore, X3 and X4 are unconfounded
according to M and the causal effect of X3 on X4 is given by:
PM X4 | do(X3 = x3 ) = PM (X4 | X3 = x3 ).
(2)
Proof By the assumed Markov and faithfulness properties, one can check that the only
(cyclic or acyclic) graphs that are compatible with the observed conditional independences
are the ones in Figure 4 (left), where X1 must be adjacent to X3 via at least one of the twodashed edges, and similarly, X2 must be adjacent to X3 via at least one of the two dashededges. Hence, X3 is a direct cause of X4 according to M, but X4 is not a direct cause of
X3 according to M. Also, X3 and X4 cannot be confounded according to M. By applying
the causal do-calculus, we arrive at (2).
This example illustrates how conditional independence patterns in the observational distribution allow one to infer certain features of the underlying causal model. This principle
is exploited more generally by constraint-based methods, and implicitly, by score-based
methods that optimize a penalized likelihood over (equivalence classes of) causal graphs.
Typically, the graph cannot be completely identified from purely observational data. For
example, in the Y-structure case, the conditional independences in the observational data
do not allow to conclude whether the dependence between X1 and X3 is explained by X1
being a cause of X3 , or by X1 and X3 having a latent confounder, or both. However, under
the assumption of faithfulness, one can deduce the Markov equivalence class of the graphfrom the conditional independences in the observational data, i.e., the class of all DMGs
that induce the same separations. Another disadvantage of causal discovery methods frompurely observational data is that they typically need very large sample sizes and strong
assumptions in order to work reliably. These are some of the motivations to combine theseideas with those of causal discovery by experimentation, as we will do in the next section.
15

Mooij, Magliacane and Claassen

X1

X4

X3

X4

X2
X2

X1

X1

X2

X4
Figure 4: Left: Causal graphs satisfying the “Y-structure” pattern on four variables(X1 ,X2 ,X3 ,X4 ). Right: Scatter plots illustrating the Y-structure pattern inpurely observational data, where X3 is discrete-valued and its value is indicatedby color (red/blue).

3. Joint Causal Inference
In this section we present Joint Causal Inference (JCI), a novel framework for causal discovery from multiple data sets corresponding to measurements that have been performed
in different contexts. JCI combines the existing approaches towards causal discovery that
we discussed in Sections 2.2 and 2.3.
3.1. The Distinction between System and Context
Henceforth, we will distinguish system variables (Xi )i∈I describing the system of interest,and context variables (Ck )k∈K describing the context in which the system has been observed.
An observation that will turn out to be crucial in what follows is that the decision of what
to consider part of the “system” and what to consider part of its “context” does not reflect
an objective property of nature, but is a choice of the modeler.
While the system variables are treated as endogenous variables of the system of interest, we usually (but not necessarily) think of the context variables as observed exogenous
variables for the system of interest. In particular, context variables could describe which interventions have been performed on the system (or more specifically, how these interventions
have been performed), in which case we will also refer to them as intervention variables.
The possible interventions are not limited to the perfect interventions modeled by the dooperator of Pearl (2009), but can also be more general types of interventions that appear
in practice, like mechanism changes (Tian and Pearl, 2001), soft interventions (Markowetz
et al., 2005), fat-hand interventions (Eaton and Murphy, 2007), activity interventions (Mooij
and Heskes, 2013), and stochastic versions of all these. This will be discussed in more detailin Section 3.3. Even more generally, a context variable could describe any property of theenvironment of the system, including those properties that one would not normally think
about as an intervention. Examples are the lab in which measurements have been done, thetime of the day, the patient population, variables like “gender” or “age”, etc. Like systemvariables, context variables can be discrete or continuous (or more generally, take values insome Borel space).
16

Joint Causal Inference from Multiple Contexts

context
meta-systemsystemFigure 5: JCI reduces modeling a system in its environment to modeling the meta-system
consisting of the system and its environment.

The idea of explicitly considering context variables is not novel: they have been discussed
in the literature under various names, such as “policy variables” (Spirtes et al., 2000), “force
variables” (Pearl, 1993), “decision variables” in influence diagrams (Dawid, 2002), “regime
indicators” (Didelez et al., 2006), “selection variables” in selection diagrams (Bareinboim
and Pearl, 2013), and “environment variable” (Peters et al., 2016). Their use for causal
discovery was already suggested by Cooper and Yoo (1999). Formal aspects in how thesevariables are treated vary across accounts, however. For example, Dawid (2002) treatssystem variables as random variables and chooses to not treat context (“decision”) variables
as random variables. In this work we simply consider context variables as random variables
with added background knowledge on their causal relations, which expresses their assumed
exogeneity with respect to the system.
Conceptually, context variables provide a more general notion than intervention variables, since every intervention can be seen as a change of context, but not every change ofcontext is naturally thought of as an intervention. For example, the causal effect of somedrug on a certain health outcome may differ for males and females. Taking “gender” as acontext variable that just encodes the specific subpopulation of patients we are considering
is more natural than considering it to be an intervention variable that encodes the result ofa gender-changing operation on the patient. Furthermore, interventions usually come withan “observational baseline” of “doing nothing”, but this is not always naturally available
for more general context variables (e.g., “male” and “female” could both qualify as a baseline, while neither of the two would provide a more natural “observational” baseline thanthe other). When considering context variables, we do not have to specify such a baseline,
whereas if we consider them as intervention variables, one can always ask “which value ofthe variable corresponds with no intervention?”. Ultimately, though, both interpretations
can be treated equally from a mathematical modeling perspective. Henceforth, we will use
the term “context variable” in general, but “intervention variable” specifically for context
variables that model an external intervention on the system.
That being said, the approach we take in JCI is simple (see also Figure 5): rather
than considering a causal model of the system alone (i.e., modeling only the endogenous
system variables), we broaden its scope to include relevant parts of the environment of thesystem (i.e., we include the context variables as additional endogenous variables). Thereby,
we “internalize” parts of the environment of the system, which makes the meta-system(consisting of both system and its environment) amenable to formal causal modeling. The
17

Mooij, Magliacane and Claassen

meta-system can now formally be considered as occurring in just a single (meta)-context,
and thereby we have reduced the problem of how to deal with multiple contexts to one ofdealing with a single context only. We will formalise this idea in the next subsection.
3.2. Joint Causal Modeling of Multiple Contexts
Different approaches to modeling multiple contexts can be taken, e.g., using influence diagrams (Dawid, 2002), using selection diagrams (Bareinboim and Pearl, 2013), considering
only conditional models (i.e., for the conditional probability of the system given the context)(Eaton and Murphy, 2007; Mooij and Heskes, 2013), or using ioSCMs (Forré and Mooij,
2019). Here, we will take what is perhaps the simplest approach: we treat both context
and system variables as endogenous variables in an SCM.
We will use a simple SCM to model the meta-system (i.e., the system and its contexts)
causally. The endogenous
Q variables of the SCM consist of the system variables X = (Xi )i∈I
withvaluesx∈
X=
i∈I Xi and the context variables C = (Ck )k∈K with values c ∈ C =
Q
C
.
The
latentexogenous
variables of the SCM are denoted E = (Ej )j∈J with valuesk∈K k Qe ∈ E = j∈J Ej . The SCM modeling the meta-system is then assumed to be of thefollowing form:Ck = fk (XpaH (k)∩I , CpaH (k)∩K , EpaH (k)∩J ),M : Xi = fi (XpaH (i)∩I , CpaH (i)∩K , EpaH (i)∩J ),
Q
P(E) = j∈J P(Ej ).

k ∈ K,i ∈ I,(3)

The system variables X and context variables C are all treated as endogenous variables ofthe meta-system, and the exogenous variables E are independent latent variables that are
assumed not to be caused by the system variables X or the context variables C.11 The
augmented graph H has nodes I ∪ J ∪ K and directed edges corresponding to the functional
dependencies of the causal mechanisms on the variables. The graph G(M) has only nodes
I ∪ K, and may contain both directed and bidirected edges between the nodes, expressing
direct causal relations and latent confounders.
Note that the most general way to use SCMs to model multiple contexts would be touse separate SCMs, one for each context. In that approach, we could have a different graphfor each context. Representing the contexts jointly, as in (3), we simply obtain the union of
those graphs. In particular, even if within each context, the system is acyclic, it could bethat the mixture of systems in different contexts has a cyclic graph. As a simple example,
consider a system with two system variables X1 and X2 , and consider two different contexts,
where in the first context X1 causes X2 (but not vice versa), and in the second context, X2
causes X1 (but not vice versa); see also Figure 6. As a more concrete example, the engine
drives the wheels of a car when going uphill, but when going downhill, the rotation of thewheels drives the engine. Modeling this in a joint SCM as in (3) requires a cyclic graph.
The model (3) imposes a probability distribution P(C) on the context variables, thecontext distribution. The context distribution will reflect the empirical distribution of thecontext variables in the pooled data P̂(C), by using as the probability of a context the
11. At this stage, we have not yet incorporated the assumption that context variables are exogenous to thesystem, and they are still treated equally to system variables in (3).

18

Joint Causal Inference from Multiple Contexts

(a)

(b)

(c)

C = 0:X1

C = 1:

X2

X1

C

X2

X1

X2

Figure 6: The graph of a mixture of two acyclic SCMs can be cyclic. (a) X1 causes X2 incontext C = 0; (b) X2 causes X1 in context C = 1; (c) X1 and X2 cause eachother in the joint model.

fraction of the total number of samples that have been measured in that context. In casethe context variables are used to model interventions, for example, the context distribution
is determined by the experimental design. One might object that this makes the model
very specific to the particular setting, since it also specifies the relative numbers of samples
in each data set, but as it turns out, the conclusions of the causal discovery procedure do
not depend on these details under reasonable assumptions, and therefore generalize to othercontext distributions. In other words, the behavior of the system is invariant of the context
distribution.
Because the context variables are treated as endogenous variables (similarly to the system variables), we have “internalized” them. The main advantage of our modeling approachover alternative approaches is that in (3), context variables are formally treated in exactly
the same way as the system variables. This implies in particular that all standard definitions and terminology of Section 2.1, and all causal discovery methods that are applicablein that setting, can be directly applied.
3.3. Modeling Interventions as Context Changes
The causal model in (3) allows one to model a perfect intervention in the usual way (Pearl,2009). Specifically, the perfect intervention that forces
Q XI to take on the value ξI (“do(XI =ξI )”) for some subset I ⊆ I and some value ξI ∈ i∈I Xi can be modeled by replacing thestructural equations for the system variables in (3) by:
(
ξi
Xi =
fi (Xpa(i)∩I , Cpa(i)∩K , Epa(i)∩J )i∈Ii ∈ I \ I,while leaving the rest of the model invariant.12
Alternatively, the context variables can be used to model interventions. For example,
the same perfect intervention could be modeled by introducingQa context variable Ck thathas ch(k) = I, no parents or spouses, and domain Ck = {∅} ∪ i∈I Xi , by taking fI to be
12. For brevity, we dropped the subscript H of paH (·).

19

Mooij, Magliacane and Claassenof the following form:(f˜i (Xpa(i)∩I , Cpa(i)∩K\{k} , Epa(i)∩J )fi (Xpa(i)∩I , Cpa(i)∩K , Epa(i)∩J ) =(Ck )iCk = ∅
Q
Ck ∈ i∈I Xi(4)for i ∈ I. Here, Ck = ∅ corresponds to no intervention (i.e., the observational baseline).
Modeling a perfect intervention in this way is similar to the concept of “force variables”
introduced by Pearl (1993). The observational distribution of the system variables is thengiven by the conditional distribution P(X | Ck = ∅), the interventional distribution corresponding to the perfect intervention do(XI = ξI ) is given by the conditional distribution
P(X | Ck = ξI ), and the marginal distribution P(X) represents a mixture of those. This is
illustrated in Figure 7.
More general types of interventions such as mechanism changes (Tian and Pearl, 2001)
can be modeled in a similar way, simply by not enforcing the dependence on Ck to beof the form (4), but allowing more general forms of functional dependence. For example,
switching the causal mechanism of system variable Xi from mechanism A to mechanism
B can be modeled as follows by introducing a context variable Ck with ch(k) = {i} anddomain Ck = {A, B}:(f˜iA (Xpa(i)∩I , Cpa(i)∩K\{k} , Epa(i)∩J ) Ck = A
fi (Xpa(i)∩I , Cpa(i)∩K , Epa(i)∩J ) =f˜iB (Xpa(i)∩I , Cpa(i)∩K\{k} , Epa(i)∩J ) Ck = B.
As another example, a stochastic perfect intervention on Xi that is only successful with acertain probability can be modeled by having one of the latent exogenous variables Ej withj ∈ pa(i) determine whether the intervention was successful:fi (Xpa(i)∩I , Cpa(i)∩K , Epa(i)∩J )(f˜i (Xpa(i)∩I , Cpa(i)∩K\{k} , Epa(i)∩J \{j} ) Ck = ∅ or Ej = 0=
Ck
Ck ∈ Xi and Ej = 1.
This approach of modeling interventions by means of context variables is very general,
as it allows to treat various types of interventions in a unified way. For example, it can deal
with perfect interventions (Pearl, 2009), mechanism changes (Tian and Pearl, 2001), softinterventions (Markowetz et al., 2005), fat-hand interventions (Eaton and Murphy, 2007),
activity interventions (Mooij and Heskes, 2013), and stochastic versions of all these. Incase the context variables are used to model interventions in this way, we also refer to thecontext distribution P(C) (the probability for each context to occur) as the experimental
design.
3.4. JCI Assumptions
In this subsection, we discuss additional background knowledge on the causal relationshipsof context variables that one may often have in practice, and that can be very helpful forcausal discovery.
20

Joint Causal Inference from Multiple Contexts

(a)

(b)

Observational:(Cα = 0)

X1

X2

X3

Interventional:(Cα = 1)

X1

X2

X3

Observational:(Cα = 0)X1
X2 X3
-0.2 -0.4 0.6
0.6
0.8 1.3-1.7 0.1 0.3

Jointly:

Interventional:(Cα = 1)
X1
X2
X3
-0.3 1.8 -0.1
1.8 -2.2 -0.2

Cα

X1

X2

X3

Cα
0
0
0
1
1

Pooled:X1
X2-0.2 -0.4
0.6
0.8-1.7 0.1-0.3 1.8
1.8 -2.2

X3
0.6
1.3
0.3
-0.1
-0.2

Figure 7: Two ways of representing interventions, either through modeling contexts separately (left), or by modeling system and context jointly (right). In this example,
we consider a perfect intervention on X2 , though the same idea applies to othertypes of interventions. (a) shows the corresponding causal graphs, as separate
graphs for each context (left), or as a single joint graph that includes a context
variable (right); (b) shows different ways of grouping the data: as separate data
sets for each context (left), or as a single joint data set after pooling (right).

3.4.1. JCI Assumption 0
First, we restate formally our basic modeling assumption:
Assumption 0 (“Joint SCM”) The data-generating mechanism is described by a simple
SCM M of the form:Ck = fk (XpaH (k)∩I , CpaH (k)∩K , EpaH (k)∩J ),M : Xi = fi (XpaH (i)∩I , CpaH (i)∩K , EpaH (i)∩J ),
Q
P(E) = j∈J P(Ej ),k ∈ K,i ∈ I,(5)that jointly models the system and the context. Its graph G(M) has nodes I ∪ K (corresponding to system variables {Xi }i∈I and context variables {Ck }k∈K ).
Whereas we will always make this assumption in order to facilitate the formulation of JCI,
the following three assumptions that we discuss are optional, and their applicability has tobe decided based on a case-by-case basis.
3.4.2. JCI Assumption 1
Typically, when a modeler decides to distinguish a system from its context, the modelerpossesses background knowledge that expresses that the context is exogenous to the system:
21

Mooij, Magliacane and Claassen

Assumption 1 (“Exogeneity”, optional) No system variable causes any context variable,i.e.,∀k ∈ K, ∀i ∈ I : i → k ∈/ G(M).
This exogeneity assumption is often easy to justify, for example if context is gender or age.
Another common case is that the context encodes interventions on the system that have been
decided and performed on the system before measurements on the system are performed:
this already rules out any causal influence of system variables on the intervention (context)
variables if time travel is not deemed possible. Of course, one can imagine settings in whicha system variable was measured before an intervention was performed on the system. For
example, a doctor typically first diagnoses a patient before deciding on treatment. For
system variables containing the results of the medical examination used for the diagnosis
and intervention variables describing the treatment that was decided after —and based
upon—the medical examination, JCI Assumption 1 would not apply.
3.4.3. JCI Assumption 2
The second JCI assumption generalizes the randomization assumption for randomized controlled trials:
Assumption 2 (“Complete randomized context”, optional) No context variable is confounded with a system variable, i.e.,∀k ∈ K, ∀i ∈ I :i↔k∈/ G(M).

This assumption is often harder to justify in practice. It is justifiable in experimental
protocols in which the decision of which intervention to perform on the system does notdepend on anything else that might also affect the system of interest, and in which theobserved context variables provide a complete description of the context. This is ensured
for example in case of proper randomization in a double-blind randomized trial setting, i.e.,
in which neither the patient nor the physician knows whether the patient was assigned adrug or a placebo.
Many experimental protocols that do not involve explicit coin flips or random number
generators are implicitly performing randomization. For example, in the experimental procedure described by Sachs et al. (2005) (see also Section 5.8), one starts with a collection ofhuman immune system cells. These are divided into batches randomly, without taking into
account any property of the cells. When done carefully, the experimenter tries to ensurethat for example the size of a cell cannot influence the batch it ends up in, by stirring
the liquid that contains the cells before pipetting. Then, after randomly assigning cells tobatches, interventions are performed on each batch separately, by adding some chemicalcompound to the batch of cells. Finally, properties of each individual cell within each batchare measured. If the system variables reflect the measured properties of the individual cells,
and the context variables encode the batch ID, this experimental procedure justifies JCI
Assumption 2.
However, one should be careful not to jump to the conclusion that the chemical compound administered to the batch is what actually causes the observed system behavior, as
there may be other factors that vary across batches due to unintentional side effects of the22

Joint Causal Inference from Multiple Contexts

experimental procedure. For example, the lab assistant that carries out the experiment fora particular batch of cells might influence the outcome, because slightly different experimental procedures are used by different lab assistants. Another example is that the time of the
day may affect the measurements, and also correlate with batch ID. In situations like those,
identifying the batch ID with the chemical compound administered to that batch could be
misleading, and could lead one to incorrectly attribute the inferred causal relation between
batch ID and a certain system variable to the causal effect of the intended intervention
corresponding to that batch on the system variable. This is a subtle type of error thatthe causal modeler should beware of. Even though we have good reasons to assume thatproper randomization was performed for batch ID in the Sachs et al. (2005) experiment,
it is questionable whether the interpretation of the context variables as concerning solely
the addition of certain chemical compounds (and not any other factors that actually varied
across batches) is appropriate.
The issue can also be understood by noting that JCI Assumption 2 may not be preserved
when marginalizing out context variables, as illustrated in Figure 8. The following example
describes a situation in which this phenomenon may occur.
Example 1 Consider a randomized trial setup for establishing whether sugar causes plants
to grow. Context variable Cα denotes the coin flip result, Cβ indicates whether sugar isadministered to the plant, and Cγ indicates whether water is administered to the plant. The
experimenter decided to use an experimental design with two groups, and assigning plants
to groups with a coin flip. One group of plants was administered a solution consisting ofsugar dissolved in water on a daily basis, the other (control) group was not treated in any
way. The growth rate X1 of the plants was measured for both groups. Suppose the following
experimental design was used:
P(C = c)12
12Cα (coin flip)01

Cβ (sugar)
0
1

Cγ (water)0
1

If one would only take context variable Cβ (did the plant get sugar?) into account and would
treat Cα and Cγ as latent, as in Figure 8(b), and would make JCI Assumptions 1 and 2,
one would arrive at the (wrong) conclusion that sugar causes plants to grow. However, if
one would take all three context variables into account, and make JCI Assumptions 1 and
2, one would obtain the right conclusion that at least one of the three context variables mustcause plants to grow.
A simple remedy to avoid the wrong conclusion if only Cβ is observed would be to drop JCI
Assumption 2: then it is no longer identifiable whether Cβ causes X1 , or whether Cβ andX1 are just confounded.
3.4.4. JCI Assumption 3
We have seen that JCI Assumption 1 is often easily justifiable, but the applicability of JCI
Assumption 2 may be less obvious in practice. We will now state JCI Assumption 3, which
can be useful whenever both JCI Assumptions 1 and 2 have been made as well.
23

Mooij, Magliacane and Claassen

(a)

(b)

Cα
Cβ

Cβ

Cγ

X1

X1

Figure 8: Confounding between system and context variables due to unobserved context
variables. (a) If all three context variables Cα , Cβ , Cγ are observed, JCI Assumption 2 would be valid. (b) After marginalizing out Cα and Cγ , leaving only context
variable Cβ as observed, JCI Assumption 2 is no longer valid.

Assumption 3 (“Generic context model”, optional) The context graph13 G(M)K is of thefollowing special form:∀k 6= k 0 ∈ K :k ↔ k 0 ∈ G(M)∧k → k0 ∈/ G(M).

In Figure 9(b), this assumption is satisfied, while in Figure 9(a), it is not. We will show thatJCI Assumption 3 seems stronger than it is, since it can be made without loss of generality
in many cases occurring in practice.
In order to precisely formulate and prove that claim, the following definition is needed.
Definition 12 Given an SCM M satisfying JCI Assumption 0, define the conditional
system graph G(M)do(K) as the DMG with context nodes K and system nodes I, and as
directed and bidirected edges those edges in G(M) that contain at least one system node in
I (i.e., excluding edges between context nodes). We will graphically represent the systemnodes I of G(M)do(K) by ellipses and the context nodes K of G(M)do(K) by squares.
Figure 9(c) shows the common conditional system graph for SCMs with graph as given
in Figure 9(a) and for SCMs with graph as given in Figure 9(b). The conditional systemgraph provides a particular graphical representation for an SCM with context and systemvariables that is less expressive than its graph. This representation is useful when we are
not interested in describing relationships between context variables, but only in describing
the relationships between system variables and how the context affects the system.
The following key result essentially states that when one is only interested in modelingthe causal relations involving the system variables (under JCI Assumptions 1 and 2), one
does not need to care about the causal relations between the context variables, as long as
one correctly models the context distribution.

13. Remember that G(M)K denotes the subgraph on the context variables K induced by the causal graph
G(M).

24

Joint Causal Inference from Multiple Contexts

Theorem 13 Assume that JCI Assumptions 0, 1 and 2 hold for SCM M:k ∈ K,
Ck = fk (CpaH (k)∩K , EpaH (k)∩J ),M : Xi = fi (XpaH (i)∩I , CpaH (i)∩K , EpaH (i)∩J ),i ∈ I,
QP(E) = j∈J P(Ej ),For any other SCM M̃ satisfying JCI Assumptions 0, 1 and 2 that is the same as M exceptthat it models the context differently, i.e., of the form˜Ck = fk (CpaH̃ (k)∩K , EpaH̃ (k)∩J˜ ),M̃ : Xi = fi (XpaH (i)∩I , CpaH (i)∩K , EpaH (i)∩J ),
Q
P(E) = j∈J˜ P(Ej ),k ∈ K,i ∈ I,

with J ⊆ J˜ and paH (i) = paH̃ (i) for all i ∈ I, we have that(i) the conditional system graphs coincide: G(M)do(K) = G(M̃)do(K) ;(ii) if M̃ and M induce the same context distribution, i.e., PM (C) = PM̃ (C), then forany perfect intervention on the system variables do(I, ξI ) with I ⊆ I (including thenon-intervention I = ∅), M̃do(I,ξI ) is observationally equivalent to Mdo(I,ξI ) .
(iii) if the context graphs G(M̃)K and G(M)K induce the same separations, then also G(M̃)
and G(M) induce the same separations (where “separations” can refer to either dseparations or σ-separations).
Proof See Appendix A.
The following corollary of Theorem 13 states that JCI Assumption 3 can be madewithout loss of generality for the purposes of constraint-based causal discovery if the context
distribution contains no conditional independences:
Corollary 14 Assume that JCI Assumptions 0, 1 and 2 hold for SCM M. Then there
exists an SCM M̃ that satisfies JCI Assumptions 0, 1 and 2 and 3, such that(i) the conditional system graphs coincide: G(M)do(K) = G(M̃)do(K) ;(ii) for any perfect intervention on the system variables do(I, ξI ) with I ⊆ I (including
the non-intervention I = ∅), M̃do(I,ξI ) is observationally equivalent to Mdo(I,ξI ) ;(iii) if the context distribution PM (C) contains no conditional or marginal independences,then the same σ-separations hold in G(M̃) as in G(M); if in addition, the Directed
Global Markov Property holds for M, then also the same d-separations hold in G(M̃)
as in G(M).
Proof This follows from Theorem 13 by showing that there exists an M̃ that satisfies all
requirements in Theorem 13 and JCI Assumption 3 by construction, and that induces thesame context distribution as M does. For a detailed proof, see Appendix A.
25

Mooij, Magliacane and Claassen

(a)

(b)
Cα

Cβ

X1

X2

Cγ

X3

(c)
Cα

Cβ

X1

X2

Cγ

X3

Cα

Cβ

X1

X2

Cγ

X3

Figure 9: Example graphs of (a) a true SCM M and (b) the modified SCM M̃ constructedin the proof of Corollary 14 that satisfies JCI Assumption 3, and (c) their corresponding conditional system graph G(M)do(K) . Corollary 14 gives sufficientconditions for M̃ and M to be equivalent for our purposes.

(a)

(b)
Cα

Cβ

Cα

Cβ

X1

X2

X1

X2

Figure 10: Example that illustrates that the genericity assumption in statement (iii) of
Corollary 14 is necessary. Graphs of (a) the true SCM M and (b) the modified
SCM M̃ constructed in the proof of Corollary 14 that are not Markov equivalent. The graph in (a) is identifiable under JCI Assumptions 0–2. The joint
distribution PM (X, C) is not faithful with respect to the graph in (b), whichis the minimal one that also satisfies JCI Assumption 3 such that PM (X, C) is
Markov with respect to it.

An example illustrating this corollary is provided in Figure 9.
JCI Assumption 3 is typically made for convenience. When our aim is not to model
the causal relations between the context variables, but just to use the context variables as
an aid to model the causal relations between system variables and between context andsystem variables, Corollary 14 shows that we may assume JCI Assumption 3 without lossof generality if JCI Assumptions 1 and 2 are made and the context distribution contains no(conditional) independences. The causal discovery algorithm then does not need to waste
time on learning the causal relations between context variables but can focus directly on
learning the causal relations involving the system variables.
Note that the genericity assumption in statement (iii) of Corollary 14 (i.e., PM (C)containing no conditional independences) is necessary, as the simple counterexample in
Figure 10 shows. Depending on how well the causal discovery algorithm can handle faithfulness violations, model misspecification due to incorrectly assuming JCI Assumption 3even though P(C) contains conditional independences might prevent successful identifica26

Joint Causal Inference from Multiple Contexts

Cα
0
1
0
0
0
0

Cβ
0
0
1
0
0
0

Cγ
0
0
0
1
0
0

Cδ
0
0
0
0
1
0

C
0
0
0
0
0
1

possible interpretation
observational
intervention αintervention βintervention γ
intervention δintervention 

Table 1: Example of a diagonal design with 5 context variables. If the context variables
are indicators of interventions, the context with Ck = 0 for all k ∈ K correspondswith the purely observational setting, and the other contexts in which one Ck = 1
and the other Cl = 0 for l 6= k correspond with a particular intervention each.

tion of the causal relationships between system variables. Therefore, it is prudent to checkthat the empirical context distribution P̂(C) indeed contains no conditional independences
before making JCI Assumption 3.
An example of a common situation in which the context distribution contains no conditional independences is what we refer to as a diagonal design (see also Table 1). This is a
simple experimental design that is often used to discover the effects of single interventions
when one is not interested in understanding the interactions that multiple interventions
might have. Note that two non-constant binary variables X, Y can only be independent if
P(X = 1, Y = 1) > 0. Even more, they can only be conditionally independent given a third
discrete variable Z if P(X = 1, Y = 1 | Z = z) > 0 for all z with P(Z = z) > 0. Therefore,
each pair of context variables is dependent in a diagonal design (as there is no context inwhich a pair of context variables simultaneously obtains the value 1), even conditionally on
any subset of the other context variables. In other words, the context distribution P(C)corresponding to any such diagonal design (with non-zero probability for each context)
contains no conditional independences.
JCI Assumption 3 can easily be modified for situations in which the context distribution does contain conditional independences. For example, in the extreme case in which
all context variables are jointly independent, one would simply assume that G(M) contains
no directed and no bidirected edges between context variables. Such situations may occur
for symmetric experimental designs in which all context variables are jointly independent
by design (for example, factorial designs with equal sample sizes in each experimental context). However, we believe that this occurs less often in practice than the generic case inwhich all context variables are (conditionally) dependent, because resource constraints often
lead experimenters to deviate from completely symmetric experimental designs. Therefore,
rather than assuming the context variables to be jointly independent as a default, we have
opted here for the more generic default of assuming that no conditional independences holdbetween context variables in the context distribution.
More generally, one could replace JCI Assumption 3 by assuming that G(M)K equalsa certain graph that expresses the known conditional independences in the experimental
design. Theorem 13 can be applied to these more general situations as well and shows thatfor the purpose of constraint-based causal discovery, any context graph that implies the27

Mooij, Magliacane and Claassen

(a)

(b)
Cα

Cβ

Cα

Cβ

X1

X2

X1

X2

Figure 11: If JCI Assumption 2 does not apply, the causal relations between context variables have testable consequences for the conditional independences in the joint
distribution. (a) X1 ⊥6 ⊥ X2 | {Cα , Cβ }; (b) X1 ⊥⊥ X2 | {Cα , Cβ }.

observed conditional independences (i.e., any graph that is Markov equivalent to the truecontext graph) works.
Another alternative is to omit JCI Assumption 3 and instead try to infer the context
subgraph G(M)K from the data. This would typically be computationally more expensive,
but in our experience does not seem to make much of a difference in terms of accuracy in
our experiments (as we report in Section 5).
JCI Assumption 3 only makes sense when both JCI Assumptions 1 and 2 are made.
If we would not make JCI Assumption 1 or 2, the causal relations between the observed
context variables will have testable consequences in the joint distribution in general. For
an example of this,14 see Figure 11. Here, Cα could be “lab”, and Cβ could be the “temperature” at which an experiment is performed. In this case, we get different conditional
independences in the joint distribution P(X1 , X2 , Cα , Cβ ) if lab causes temperature thanwhen they are confounded (for example, by geographical location). Something similar canhappen if context variables are caused by system variables.
3.4.5. Summary of JCI Assumptions and Other Background Knowledge
Summarizing, the JCI framework rests on different assumptions, one of which is required,
whereas the others are all optional. The basic assumption that is required is JCI Assumption 0, which states that the meta-system consisting of context and system can be describedby a simple SCM. This is just the standard assumption made throughout the causal discovery literature, but now applied to the meta-system rather than to the system only. In
addition, assumptions about the causal relationships of the context variables can be made,which are all optional and can be decided on a case-by-case basis. In most cases, we would
expect JCI Assumption 1 (no system variable causes any context variable) to apply. In
some cases, also JCI Assumption 2 (no system variable is confounded with any context
variables) applies. If both apply, one can assume JCI Assumption 3 for convenience if thecontext distribution contains no (conditional) independences. More generally, we only needto model the observed conditional independences in the context distribution, not necessarily
their causal relations, when our interest is in modeling the causal relations involving systemvariables only.
14. We are grateful to Thijs van Ommen for pointing this out.

28

Joint Causal Inference from Multiple Contexts

The reader may wonder when one can ever be sure in practice that JCI Assumption 2
applies. There is one very common scenario in which JCI Assumption 2 holds. This is in ascientific experiment in which, in chronological order:(i) an ensemble of systems is prepared in an initial state;(ii) the systems are randomly permuted and randomly divided into batches;(iii) for each batch, all systems in the batch are intervened upon simultaneously in thesame way (following an experimental protocol determined in advance);(iv) measurements of the system variables are performed.
The experimental protocol specifies the intended interventions for each batch, which should
be completely encoded as context variables. Since the intended interventions have been
decided before system variables are measured, the intended interventions cannot be causedby the system variables. Because the systems were randomly permuted and divided into
batches, the assigned batch cannot be caused by prior values of the system variables, or byanything else that may also have an effect on the system variables. Because the intended
interventions for each system in each batch are determined completely by the batch, this
implies that intended interventions and system variables cannot be confounded. As long as
the context variables provide a complete encoding of the intended interventions (i.e., theintended interventions are in one-to-one correspondence to values of the context variables),JCI Assumptions 1 and 2 then apply to the context variables.15 If additionally, no (conditional) independences hold in the empirical context distribution, we can also make use of
JCI Assumption 3 to simplify and speed up the causal discovery procedure.
In more general scenarios, such as the example in the introduction (concerning thequestion whether playing violent computer games causes aggressive behavior), the validityof JCI Assumption 2 (no confounding between context and system) should not be taken for
granted. For example, it could be that precisely the schools with a more violent populationof pupils see themselves forced to actively take measures to promote social behavior. In
that case, the level of violence in the past would confound Cβ (does a school take measures
to stimulate social behavior) and X2 (how violently do the pupils of the school behave).
Thus, in scenarios like these, it seems safer not to rely on JCI Assumption 2 as incorrectly
assuming it might lead to wrong conclusions (although we do not currently understand theprecise impact of such model misspecification).
For causal discovery in the JCI framework, knowledge of the intervention targets (ormore generally, which system variables are affected directly by which context variables) is
not necessary, but it is certainly helpful and can be exploited similarly to other available
background knowledge, depending on the algorithm used to implement JCI. When applying
JCI on a combination of different interventional data sets, intervention targets can be learnt
15. If the experimenter sticks to the experimental protocol that was fixed before the experiment was performed, any possible influence of the system variables on the performed interventions is excluded, andtherefore the performed interventions will equal the intended interventions. This means that the JCI
modeling framework (with JCI Assumptions 1 and 2) applies also when interpreting the context variables
as the performed interventions. This may explain why it is considered good scientific practice to perform
an experiment according to an experimental protocol that was fixed beforehand, and not deviate fromit in case of unexpected measurement outcomes, for example.

29

Mooij, Magliacane and Claassenfrom data when they are not known (as the direct effects of intervention variables), similarly
to how the effects of system variables can be learnt. One main advantage of the JCI
framework is that it offers a unified way to deal with different types of interventions, as
discussed in Section 3.3. Therefore, knowledge of intervention types (e.g., is it a perfectintervention, or a mechanism change?) is also not necessary, but can still be helpful as itprovides additional background knowledge that may be exploited for causal discovery.
In concluding this subsection, we observe that the JCI framework generalizes and combines the ideas of causal discovery from purely observational data and of causal discovery
by means of randomized controlled trials. Indeed, note that if JCI is applied to a single
context (i.e., 0 context variables), it reduces to the standard setting of causal discovery
from purely observational data described in Section 2.3. If JCI is applied to a setting witha single context variable and a single system variable, JCI (with Assumptions 1 and 2)
reduces to the randomized controlled trial setting described in Section 2.2. Therefore, theJoint Causal Inference framework truly generalizes both these special cases.

4. Causal Discovery from Multiple Contexts with JCI
In this section, we discuss how causal discovery from multiple contexts can be performedin the Joint Causal Inference framework. Our starting point is the assumption that somemodel of the form (5) is an appropriate causal model for the system and its context, andwe have obtained samples of all system variables in multiple contexts.16 Suppose that theexact model M and in particular, its causal graph G(M), are unknown to us. The goalof causal discovery is to infer as much as possible about the causal graph G(M) from theavailable data and from available background knowledge about context and system.
Nc(c)
Let us denote the data set for context c ∈ C as D(c) = (xin )i∈I n=1, and for simplicity,
assume that no values are missing. The number of samples in each context, given by Nc ,is allowed to depend on the context. As a first step,Pwe pool the data, thereby representing
it as a single data set D = (xn , cn )N
n=1 where N =c∈C Nc . We then assume that D is an
i.i.d. sample of PM (X, C), where (X, C, E) is a solution of the SCM M of the form (5).17
In setting up the problem, we have made the simplifying assumptions that the measurement procedure is not subject to selection bias, nor to (independent) measurement error(Blom et al., 2018). We will assume that the data has been generated by an SCM in ac16. An interesting problem setting considered by several researchers (Claassen and Heskes, 2010; Tillman
and Spirtes, 2011; Triantafillou and Tsamardinos, 2015; Hyttinen et al., 2014; Forré and Mooij, 2018)that we do not consider here would be to allow for each context a (possibly context-dependent) subsetof system variables to remain unobserved.
17. Although this may sound as an innocuous assumption, it is not necessarily satisfied by the data generatingprocess. For example, suppose that in a randomized controlled trial, it is decided a priori that a certain
number N0 of patients will be assigned to the control group, and a number N1 of patients to the treatmentgroup, but which patients end up in which group is completely randomized. The resulting pooled data
is not i.i.d.; indeed, if we repeat this procedure, we will always end up with the same number of patientsin each group, whereas for an i.i.d. sample, the numbers would fluctuate around their expected values.
Nevertheless, this assumption can be made here without losing much generality. In particular, for the caseof binary treatment and binary outcome, Wasserman (2004, Section 15.5) shows that for independence
tests based on the (log) odds ratio the i.i.d. assumption can be weakened accordingly. Alternatively, ina bootstrapping procedure (which we will use in practice for most implementations of JCI in Section 5),
the resampled pooled data is i.i.d. by construction.

30

Joint Causal Inference from Multiple Contexts

cordance with JCI Assumption 0, and optionally, a subset of JCI Assumptions 1, 2 and3. To enable constraint-based causal discovery, we will assume that the joint distribution
PM (X, C) is faithful with respect to the graph G(M), using the appropriate separation
criterion (σ-separation in general, or d-separation for specific cases, as discussed in Section 4.1). We will discuss the ramifications of the faithfulness assumption in more detail in
Section 4.1.
Definition 15 We say that a particular feature of G(M) is identifiable from PM (X, C)
and background knowledge if the feature is present in the graph G(M̃) of any SCM M̃ withPM̃ (X, C) = PM (X, C) that incorporates the background knowledge.
“Feature” could refer to the presence or absence of a direct edge, a directed path, a bidirectededge, arbitrary subgraphs, or even the complete graph. The task of causal discovery is thento identify as many features of G(M) as possible from the data, the i.i.d. sample D of
PM (X, C), and the available background knowledge.
The key insight of the Joint Causal Inference framework that allows one to deal withdata from multiple contexts (D(c) )c∈C is that by incorporating the context variables explicitly, and pooling the data, we have now reduced the causal discovery problem to one that ismathematically equivalent to causal discovery from purely observational data D and applicable background knowledge on the causal relations between context and system variables(a subset of JCI Assumptions 1 and 2). If applicable, JCI Assumption 3 can be made to
reduce the computational effort.
This trick also allows us to easily learn intervention targets from data in a similar way
as we usually learn causal effects between variables from data: the intervention targets are
simply encoded as the direct effects of the intervention variables.
After discussing the faithfulness assumption in more detail, we will give a few suggestionsof how JCI can be implemented in Section 4.2.
4.1. Faithfulness Assumption
In this subsection we will discuss the subtleties of the faithfulness assumption in the JCI
setting and compare it with alternative faithfulness assumptions that have been made in
the literature.
Given a simple SCM M of the form (5) with graph G := G(M), the joint distribution
PM (X, C) induced by the SCM satisfies the Generalized Directed Global Markov Property(Theorem 8) with respect to the graph G of the SCM, i.e., any σ-separation U ⊥σG V | W
between sets of nodes U, V, W ⊆ I ∪ K in the graph G implies a conditional independenceX̃U ⊥⊥ PM (X,C) X̃V | X̃W , where we write X̃ := (X, C). Under the additional assumptionsof Theorem 9, the stronger Directed Global Markov Property holds (i.e., the d-separation
criterion).
For constraint-based causal discovery, some type of faithfulness assumption is usually
made. For simplicity, the faithfulness assumption that we make in this work is the standard
one, but we apply it to the combination of system and its environment: we assume thatthe joint distribution PM (X, C) is faithful with respect to the graph G(M) of M. In otherwords, any conditional independence X̃U ⊥⊥ PM (X,C) X̃V | X̃W for sets of nodes U, V, W ⊆I ∪K is due to the σ-separation U ⊥σG V | W in G (or d-separation, if applicable), and no other31

Mooij, Magliacane and Claassen

conditional independences in PM (X, C) exist. In particular, this assumption rules out any
conditional independence between context variables in case JCI Assumption 3 is made. If
JCI Assumption 3 is not made, conditional independences in the context distribution thatare faithfully described by any DMG are allowed.
This faithfulness assumption allows us to deal with different types of interventions,
including perfect interventions. For example, for the perfect intervention on X2 illustrated(c)in Figure 7 the causal graphs GI (restricted to the system variables {Xi }i∈I ) depend on
the context c ∈ C: in the observational context (Cα = 0), X1 → X2 , whereas in theinterventional context (Cα = 1), this direct causal relation is no longer present (as it has
been overruled by the perfect intervention). This does not invalidate the faithfulness ofthe joint distribution P(Cα , X1 , X2 , X3 ) with respect to the joint causal graph. Indeed,even though X1 ⊥⊥ X2 | Cα = 0, we still have X1 ⊥6 ⊥ X2 | Cα because X1 ⊥6 ⊥ X2 | Cα = 1.18
In other words, the fact that P(X | Cα = 1) is not faithful to the system subgraph GI(i.e., the induced subgraph of the causal graph G on the system nodes I) does not lead toany problem as long as we are not going to test for independences in the subset of data
corresponding to context Cα = 1 separately, but restrict ourselves to testing independences
only in the pooled data set that combines all contexts.
Causal discovery methods that analyze data from each context separately (e.g., Hauser
and Bühlmann, 2012; Triantafillou and Tsamardinos, 2015; Hyttinen et al., 2014) typically make another faithfulness assumption. In our notation, such approaches assume that(c)P(X | C = c) is faithful w.r.t. a causal subgraph GI that may be context-dependent, andmust then reason about how these context-dependent subgraphs are related, explicitly relying on knowledge about the type of interventions (typically assuming that the interventions
are perfect interventions with known targets). This faithfulness assumption is to a certain
extent stronger than ours because it requires faithfulness of the system within each context.
On the other hand, it is to a certain extent weaker than ours because ours implies restrictions on the context distribution (PM (C) must be faithful to some DMG, at the very least)that the alternative does not have. We consider the extension of the applicability of JCI(because no knowledge of intervention types or targets is needed) due to the faithfulness
assumption we chose here to outweigh the limitations in applicability (because not every
context distribution can be handled). In the rest of this section, we will discuss some simple
workarounds that can be applied in practice when dealing with faithfulness violations in
the context distribution PM (C).
Under JCI Assumption 3, the faithfulness assumption implies that the context distribution P(C) does not contain any conditional independences. In case the empirical context
distribution P̂(C) does contain conditional independences, one has several options. The
first (assuming also JCI Assumptions 1 and 2 are made) is to modify the assumed graph ofthe context variables in JCI Assumption 3 such that the context distribution is faithful toit. For example, if all context variables are jointly independent, one could simply assume
that the context graph G(M)K has no (directed or bidirected) edges at all. The second
option is to omit JCI Assumption 3 or some analogue of it completely. In that case, thefaithfulness assumption still imposes the restriction that the context distribution can be
18. Note that for a discrete context domain C, we have that A ⊥⊥ B | C if and only if A ⊥⊥ B | C = c for allc with p(c) > 0. More generally, A ⊥⊥ B | C if and only if A ⊥⊥ B | C = c for almost all c.

32

Joint Causal Inference from Multiple Contexts

faithfully modeled by some DMG. The third option is to use only data corresponding toa certain subset of the contexts and to ignore data from other contexts. In addition, one
can sometimes work around conditional independences in the context distribution by partitioning the set of context variables into groups of context variables, and using combined
context variables instead of the original context variables. This will be illustrated in thenext paragraph. Finally, one could ignore the faithfulness violations in the context distribution and hope that the causal discovery algorithm will handle them well. This last approach
usually means that it will be harder to guarantee consistency of the approach. Note thatthe faithfulness assumption for the context variables is actually testable, since the empiricalcontext distribution is available, and can be directly tested for conditional independences.
The faithfulness assumption also rules out deterministic relations between the variables
that lead to faithfulness violations. In particular, there could be deterministic relations
between context variables. For example, in the experimental design of the experiments
in Sachs et al. (2005) described in Tables 2 and 3, Cα is a (deterministic) function of Cθ
and Cι : Cα = ¬(Cθ ∨ Cι ). One might naı̈vely believe that this could be dealt with bysimply removing context variable Cα from consideration, leaving only context variables
Cβ , . . . , Cι as observed context variables, none of which is a (deterministic) function of theothers. However, marginalizing out context variables may give rise to violations of JCI
Assumption 2, as we have seen in Section 3.4.3. An operation that is generally allowed isgrouping context variables together. In the case of the Sachs et al. (2005) experimental
design, we can combine Cα , Cθ and Cι together into a single context variable given by the
triple (Cα , Cθ , Cι ). Accidentally, in this case this would be mathematically equivalent to
the pair (Cθ , Cι ), but the interpretation of (Cα , Cθ , Cι ) is different from that of (Cθ , Cι ).
Another option would be to ignore a subset of the contexts. In this case, one could exclude
the two contexts with Cθ = 1 or Cι = 1. Then, Cα becomes a constant, and constants canbe safely ignored (or, trivially combined with any other context variable).19
To wrap up: under JCI Assumption 3, it is advisable to check that there are indeed no
conditional independences between context variables in the empirical context distribution.
More generally, one should check whether the conditional independences in the context
distribution can be faithfully described by a DMG. If not, one can try to work aroundby applying the tricks discussed in the last paragraph (grouping context variables, andomitting certain contexts). When grouping context variables, one should note that theinferred causal relations from context variables to system variables may no longer be easilyinterpretable. A simple example that illustrates this is to consider two interventions thatare always performed together: when drug A is prescribed, also drug B is prescribed, andvice versa. In that case we cannot be sure whether the effect on outcome is due to drug
A or to drug B (or to both combined). Nonetheless, we can still use the inferred causal
relations between system variables.

19. In an earlier draft of this work (Magliacane et al., 2016a), we proposed to handle deterministic relations
between context variables by using the notion of D-separation, first presented in Geiger et al. (1990) andlater extended in Spirtes et al. (2000). However, this notion does not provide a complete characterization of conditional independences due to a combination of graph structure and deterministic relations.
Therefore, in this work we use the simpler techniques of grouping context variables and ignoring certain
contexts to deal with faithfulness violations due to deterministic relations between context variables.

33

Mooij, Magliacane and Claassen

Cα
1
1
1
1
1
1
1
0
0
1
1
1
1
1

Cβ
0
1
0
0
0
0
0
0
0
1
1
1
1
1

Cγ
0
0
1
0
0
0
0
0
0
1
0
0
0
0

Cδ
0
0
0
1
0
0
0
0
0
0
1
0
0
0

C
0
0
0
0
1
0
0
0
0
0
0
1
0
0

Cζ
0
0
0
0
0
1
0
0
0
0
0
0
1
0

Cη
0
0
0
0
0
0
1
0
0
0
0
0
0
1

Cθ
0
0
0
0
0
0
0
1
0
0
0
0
0
0

Cι
0
0
0
0
0
0
0
0
1
0
0
0
0
0

Nc
853
902
911
723
810
799
848
913
707
899
753
868
759
927

Cβ
0
1
0
0
0
0
0
0
0
1
1
1
1
1

Cγ
0
0
1
0
0
0
0
0
0
1
0
0
0
0

Cδ
0
0
0
1
0
0
0
0
0
0
1
0
0
0

C
0
0
0
0
1
0
0
0
0
0
0
1
0
0

Cζ
0
0
0
0
0
1
0
0
0
0
0
0
1
0

Cη
0
0
0
0
0
0
1
0
0
0
0
0
0
1

(Cα , Cθ , Cι )(1,0,0)(1,0,0)(1,0,0)(1,0,0)(1,0,0)(1,0,0)(1,0,0)(0,1,0)(0,0,1)(1,0,0)(1,0,0)(1,0,0)(1,0,0)(1,0,0)

NC
853
902
911
723
810
799
848
913
707
899
753
868
759
927

Table 2: Left: Experimental design used by Sachs et al. (2005). Nc is the number of data
samples in context C. Interpretation of context variables is provided in Table 3.
Right: Different choice of context variables: Cα , Cθ and Cι have been grouped
together into a single combined context variable (Cα , Cθ , Cι ) in order to deal withthe deterministic relation Cα = ¬(Cθ ∨ Cι ) (see main text for details).

Cα
Cβ
Cγ
Cδ
C
Cζ
Cη
Cθ
Cι

Reagent
α-CD3, α-CD28
ICAM-2
AKT inhibitor
G0076
Psitectorigenin
U0126
LY294002
PMA
β2CAMPIntervention
global activator
global activator
activity of AKT
activity of PKC
abundance of PIP2
MEK activity
PIP2, PIP3 mechanism change
PKC activity
PKA activity

Table 3: For each context variable in Table 2: reagents used in this experimental setting,and expected intervention type and targets as based on (our interpretation of)
biological background knowledge described in Sachs et al. (2005).

4.2. Implementing JCI
Any causal discovery method that is applicable under the assumptions described in Section 4 can be used for Joint Causal Inference. In this section we will describe some concrete
examples of JCI implementations. Identifiability may greatly benefit from taking into account the available background knowledge on the causal graph stemming from the applicable
JCI assumptions as discussed in Section 3.4. In addition, taking into account backgroundknowledge on targets of intervention variables may help considerably. Some logic-based
causal discovery methods (e.g., Hyttinen et al., 2014; Triantafillou and Tsamardinos, 2015;
Forré and Mooij, 2018), are ideally suited to exploit such background knowledge. For othermethods, e.g., FCI (Spirtes et al., 1999; Zhang, 2008a), or methods that focus on ancestral
relations, e.g., ACI (Magliacane et al., 2016b), incorporating all background knowledge isless straightforward and as far as we know cannot be done with off-the-shelf implementations. Often, simple adaptations of off-the-shelf implementations can be made that do allow
34

Joint Causal Inference from Multiple Contexts

one to benefit from all JCI background knowledge. For example, in Section 4.2.4 we willpropose a simple adaptation of FCI that does so.
Given a causal discovery algorithm for purely observational data that can exploit the
JCI background knowledge, we can implement JCI in a straightforward fashion:(i) introduce context variables, if not already provided;(ii) pool all data sets, including the values of the context variables;(iii) handle faithfulness violations between context variables by grouping context variables
and/or leaving out certain contexts, if necessary;(iv) apply the causal discovery algorithm on the pooled data, taking into account theappropriate JCI background knowledge.
Any soundness, completeness and consistency results for the causal discovery algorithm that
hold for the algorithm in the purely observational setting, but including the backgroundknowledge, directly apply to the JCI setting, as long as there is no model misspecification(i.e., if the assumed JCI assumptions do hold for the true model). If we use JCI Assumption
3 (or something similar), we can use Theorem 13 to show that what the algorithm concludes
about the causal relations concerning system variables is still correct.
In the remainder of this subsection, we discuss four JCI implementations. We will
first describe two existing algorithms (LCD and ICP) that can be used off-the-shelf forcausal discovery within the JCI framework. Then we propose two adaptations of existingalgorithms (ASD and FCI) so that they can be used for causal discovery within the JCI
framework. In Section 5, we will provide empirical results on the properties of those fouralgorithms.
4.2.1. Local Causal Discovery (LCD)
Perhaps the first implementation of JCI (apart from randomized controlled trials) is provided by the LCD algorithm by Cooper (1997). LCD is a very simple constraint-based
causal discovery algorithm that can be used for the purely observational causal discovery
setting where certain background knowledge is available, and in particular, in the JCI setting. The basic idea behind the LCD algorithm is the following result (which we generalized
to allow for cycles):Proposition 16 Suppose that the data-generating process on three variables X1 , X2 , X3
can be represented by a faithful, simple SCM M with I = {1, 2, 3} and that the sampling
procedure is not subject to selection bias. If X2 is not a cause of X1 according to M, thefollowing conditional (in)dependencies in the observational distribution PM (X1 , X2 , X3 )X1 ⊥6 ⊥ X2 ,X2 ⊥6 ⊥ X3 ,X1 ⊥⊥ X3 | X2

imply that the graph G(M) must be one of the three DMGs in Figure 12. In particular,(i) X3 is not a cause of X2 according to M;(ii) X2 is a direct cause of X3 according to M;(iii) X2 and X3 are not confounded according to M;
35

Mooij, Magliacane and Claassen

X1

X2

X3

X1

X2

X3

X1

X2

X3

Figure 12: All possible DMGs detected by LCD.
(iv) the causal effect of X2 on X3 is given by:
PM X3 | do(X2 = x2 ) = PM (X3 | X2 = x2 ).

(6)

Proof The proof proceeds by enumerating all (possibly cyclic) DMGs on three variables
and ruling out the ones that do not satisfy the assumptions. The assumption that X2 is nota cause of X1 implies that there is no directed edge X2 → X1 in the graph G(M). If there
were an edge between X1 and X3 , X1 ⊥⊥ X3 | X2 would not hold (faithfulness). Also, sinceX1 ⊥6 ⊥ X2 , X1 and X2 must be adjacent (Markov property). Similarly, X2 and X3 must be
adjacent. X2 cannot be a collider on any path between X1 and X3 (faithfulness). Since theonly possible edges between X1 and X2 are X1 → X2 and X1 ↔ X2 (both of which havean arrowhead at X2 ), this means that there must be a directed edge X2 → X3 , but there
cannot be a bidirected edge X2 ↔ X3 or directed edge X2 ← X3 . In other words, the only
three possible graphs are the ones in Figure 12. The causal do-calculus applied to G(M)yields (6).
In a JCI setting where JCI Assumption 1 is made, we can directly apply LCD for causal
discovery on tuples hCk , Xi , Xi0 i with k ∈ K, i 6= i0 ∈ I on the pooled data.
A conservative version of LCD has been applied by Triantafillou et al. (2017) to thetask of inferring signaling networks from mass-cytometry data. A high-dimensional version
of LCD has been shown to be successful in predicting the effects of gene knockouts on
gene expression levels (Versteeg and Mooij, 2019) from large-scale interventional yeast gene
expression data (Kemmeren et al., 2014). An algorithm closely related to LCD, named“Trigger”, has been applied on genomics data (Chen et al., 2007). Chen et al. (2007)
motivate the JCI assumptions in the setting of learning the causal relations between gene
expression levels using single nucleotide polymorphisms (SNPs) as context variables. Since
the DNA content cannot be caused by gene expression levels, JCI Assumption 1 is satisfied.
Chen et al. (2007) then argue that Mendelian randomization justifies JCI Assumption 2.
Finally, a single conditional independence in the pooled data (as in LCD) provides thedesired evidence for an unconfounded causal relation between two gene expression levels.
4.2.2. Invariant Causal Prediction (ICP)
ICP exploits invariance of the conditional distribution of a target variable given its directcauses across multiple contexts, assuming that none of the contexts corresponds with anintervention that targets the target variable (Peters et al., 2016). The implementation
described in Peters et al. (2016) handles linear relationships, arbitrary interventions (aslong as they do not change the conditional distribution of the effect variable given its
direct causes), assumes the absence of latent confounders between target variable and its
direct causes, and the absence of cycles involving the target variable. One of the main
advantages of this method over others is that it provides (conservative) confidence intervals
36

Joint Causal Inference from Multiple Contexts

on direct causal relationships that do not require the faithfulness assumption to be made;however, that only works under the assumption of causal sufficiency and acyclicity. The
authors discuss several possible extensions to broaden the scope of the method, but do
not address this in all generality. A nonlinear extension of the method has been proposed
recently (Heinze-Deml et al., 2018). ICP has been successfully applied to predict the effects
of gene knockouts on gene expression levels (Meinshausen et al., 2016) from large-scaleinterventional yeast gene expression data (Kemmeren et al., 2014).
ICP can be interpreted as a particular implementation of the JCI framework, even in
the general setting with nonlinear relations between variables and with latent confounders
and cycles present (although faithfulness is then required). The following result broadens
the conditions under which ICP identifies (possibly indirect) causal relations, strengthening
results of Peters et al. (2016):
Corollary 17 Consider the JCI setting with a single context variable C and multiple systemvariables {Xi }i∈I . Under JCI Assumptions 0 and 1 and faithfulness, the ICP estimator fortarget i ∈ I:\
Ji∗ := {I ⊆ I \ {i} : C⊥⊥ Xi | XI }PM (C,X)

satisfies Ji∗ ⊆ anG(M) (i), i.e., the set Ji∗ consists only of (possibly indirect) causes of i.
Proof Follows immediately from Proposition 28 in Appendix A.2.
Note that this means that asymptotically, ICP outputs a subset of ancestors of the target
variable, even in the presence of confounders and linear or nonlinear cycles. Hence, ICP canbe interpreted as a particular causal discovery algorithm implementing the JCI framework.
4.2.3. ASD-JCI
Here we introduce a novel JCI implementation that builds on the algorithm by Hyttinen
et al. (2014) and its generalization to σ-separation (Forré and Mooij, 2018) and some of theextensions proposed by Magliacane et al. (2016b). Since adapting this algorithm to the JCI
setting is straightforward, and the algorithm itself has been described in detail in the cited
papers, we here only provide a brief description of how it works. More details can be foundin Appendix C.
Hyttinen et al. (2014) proposed formulating causal discovery as an optimization problem over possible causal graphs, where the loss function sums the weights of all the conditional (in)dependencies present in the data that would be violated for a certain underlying
causal graph, assuming Markov and faithfulness properties. The input consists of a list of
weighted conditional independence statements. The weights λ encode the confidence in theconditional (in)dependence, where a weight of λ = ∞ corresponds to a “hard constraint”(absolute certainty) and a weight of λ = 0 corresponds to “no evidence at all”. Hyttinen
et al. (2014) provide an encoding of the notion of d-separation in Answer Set Programming (ASP), a declarative programming language that can be used amongst others forsolving discrete optimization problems. Forré and Mooij (2018) generalize the encoding to
σ-separation. The optimization problem is solved by making use of an off-the-shelf ASP
solver.
There may be multiple optimal solutions to the optimization problem, because theunderlying causal graph may not be identifiable from the inputs. Nonetheless, some of37

Mooij, Magliacane and Claassen

the features of the causal graph (e.g., the presence or absence of a certain directed edge)
may still be identifiable. We employ the method proposed by Magliacane et al. (2016b) for
scoring the confidence that a certain feature is present or absent by calculating the difference
between the optimal losses under the additional hard constraints that the feature is present
vs. that the feature is absent. Magliacane et al. (2016b) showed that this algorithm forscoring features is sound for oracle inputs and asymptotically consistent under reasonableassumptions.
We will make use of the weights proposed in Magliacane et al. (2016b): λj = log pj −log α, where pj is the p-value of a statistical test for the j th conditional independence
statement, with independence as null hypothesis, and α is a significance level (e.g., 1%)that should decrease with sample size at a suitable rate. These weights have the desirable
property that independences get a lower weight than strong dependencies.
As we will need an acronym for this algorithm later, we will henceforth refer to it as
ASD (Accounting for Strong Dependencies), as it essentially tries to explain the observed
dependencies in the data, taking into account the statistical strength of these dependencies.
This is fundamentally different from other constraint-based algorithms such as PC or FCI,
which give priority to observed independences and do not take into account the strength ofdependencies.
Taking into account the JCI background knowledge (a subset of JCI Assumptions 1, 2
and 3), and possible background knowledge on intervention targets, is trivial thanks to theexpressive power of ASP, and can be done with a few lines of ASP code. The resulting
algorithm is very accurate but scales only up to a few variables due to the combinatoric
explosion. Incorporating JCI Assumption 3 considerably reduces computation time, as it
removes the need to learn the causal relations of the context variables.
Since the JCI background knowledge can be completely and exactly encoded as constraints on the possible causal graphs, we can directly extend the known soundness, completeness and consistency results for ASD:
Theorem 18 ASD-JCI is sound and complete for oracle inputs. It is asymptotically consistent if the weights are asymptotically consistent.
Proof For the precise meaning of these statements, we refer the reader to Appendix C,
and in particular, to Theorem 43 and 44.

4.2.4. FCI-JCI
Here we introduce an adaptation of the constraint-based causal discovery algorithm FCI(Spirtes et al., 1999; Zhang, 2008a) that can be used in a JCI setting. The FCI algorithm
was designed to work under the assumption that the data was generated by an acyclic SCM.
While FCI can deal with selection bias, we will assume here for simplicity that no selection
bias is present.
The FCI algorithm consists of two main phases: an adjacency search phase leading to
the skeleton, followed by an edge orientation phase. In the adjacency search the algorithm
searches for conditional independences to eliminate edges from the graph. The subsequentorientation stage consists of a set of graphical rules that allow invariant edge marks, signify38

Joint Causal Inference from Multiple Contexts

ing either causal (tail marks) or non-causal (arrowhead marks) relations, to be added to theskeleton. For a single observational data set the final result is a so-called Partial Ancestral
Graph (PAG) that is a concise representation of ancestral relations and conditional independences. The PAG represents a set of Maximal Ancestral Graphs (MAGs) (Richardson andSpirtes, 2002), and each MAG represents a set of ADMGs (Triantafillou and Tsamardinos,2015), and each ADMG represents an infinite set of DAGs (with arbitrary number of latentvariables). In the case of purely observational data, the PAG output by FCI provides acomplete description of the Markov equivalence class (Zhang, 2008a; Ali et al., 2009). Fora more detailed discussion of MAGs and PAGs, we refer the reader to Appendix B.1.
Extending FCI such that it can take into account the additional JCI background knowledge on the adjacency and causal relations between the combined set of context and systemvariables (see also Section 3.4) is straightforward:• If JCI Assumption 3 is made, all context variables are connected by bidirected edges,
and the adjacency phase of FCI is adapted accordingly by not removing any edges
between context variables; afterwards, all edges between context variables are oriented
as k ↔ k 0 for k 6= k 0 ∈ K. In the subsequent phase of orienting unshielded triples,
only system variables can take on the role of the collider.
• If JCI Assumption 1 is made, then (since we are assuming no selection bias) any
adjacent pair of a context variable k ∈ K and a system variable i ∈ I must beconnected in the MAG by an edge with an arrowhead at i. Therefore, after theadjacency phase, all edges between a context and a system variable are oriented ask ∗→ i, with an arrowhead at the system variable i ∈ I.
• If both JCI Assumptions 1 and 2 are made, any adjacent pair of a context variablek ∈ K and a system variable i ∈ I must be connected by a directed edge k → i (since
we are assuming no selection bias) in the MAG. Hence, after the adjacency phase, alledges between a context and a system variable are oriented as k → i, pointing fromthe context variable k ∈ K to the system variable i ∈ I.
The subsequent orientation phase of the FCI algorithm does not need to be adapted.
We will refer to this adaptation of the FCI algorithm as FCI-JCI. In particular, wedistinguish three variants: FCI-JCI0 (which only makes JCI Assumption 0), FCI-JCI1(also JCI Assumption 1), and FCI-JCI123 (also JCI Assumptions 1, 2, 3). In Appendix B.4
we discuss how one can read off the identified causal and non-causal relations from the PAG
output by FCI or FCI-JCI. Furthermore, in Appendix B.5 we discuss how one can read off
the direct targets and direct non-targets of interventions represented by context nodes fromthe PAG output by FCI-JCI123.
4.2.5. Speeding up FCI-JCI123
Adding the context nodes may make FCI-JCI considerably slower than FCI on a single
context. In this subsection, we propose a further adaptation of FCI-JCI123. By exploiting
the following observation, we can achieve a considerable speedup:

39

Mooij, Magliacane and Claassen

Lemma 19 Let M be an SCM that satisfies JCI Assumptions 0, 1, 2. Then for X ⊆ I
and Y, Z ⊆ I ∪ K:
X ⊥ Y | Z =⇒ X ⊥ Y | Z ∪ (K \ Y )
G(M)

G(M)

(for both d-separation as well as for σ-separation).
Proof By contradiction: Suppose there exists a path hv0 , e1 , v1 , e2 , v3 , . . . , en−1 , vn i between
X and Y that is open given Z ∪ (K \ Y ) and contains no non-endpoint nodes in X ∪ Y ,but is closed given Z. That can only happen if the path contains a collider vi (0 < i < n)that is G(M)-ancestor of K \ (Y ∪ Z). Then vi ∈ K (because no system node is ancestorof K by JCI Assumption 1). Let 1 ≤ j < i be the lowest index such that vl ∈ K for alll ∈ {j, j + 1, . . . , i}. Then vj−1 ← vj on the path (by JCI Assumptions 1, 2), where vj−1
is in another strongly-connected component than vj , and therefore vj ∈ K \ Y blocks thepath, which is a contradiction.
This implies that in the skeleton search of FCI-JCI12 and FCI-JCI123, the search spaces
for finding separating sets between pairs of nodes can be reduced. Indeed, instead of testingfor each subset B of A ⊆ I ∪ K whether B d-separates node v from node w, one can testfor each subset B of A \ K whether B ∪ K \ {v, w} d-separates v from w.20 For the nextstages of the FCI algorithm, it does not matter which separating set is found (as long as
any separating set is found if there is one), as follows from Zhang (2006, Lemma 3.2.1 and3.2.2). This modification to the skeleton phase reduces the worst-case number of conditional
independence tests by a factor that is exponential in the number of context variables. We
will refer to the adapted version of FCI-JCI123 that implements this modified version of
the skeleton search as FCI-JCI123r.
4.2.6. Soundness, Completeness and Consistency Results for FCI-JCIThe FCI algorithm was shown to be sound and complete (Zhang, 2008a) for oracle inputs.
In Appendix B, we prove:
Theorem 20 Let M be an acyclic SCM that satisfies JCI Assumption 0. Assume that its
distribution PM (X, C) is faithful w.r.t. the graph G(M). Then, with input PM (X, C):• FCI-JCI0 is sound and complete;• FCI-JCI1 is sound if M also satisfies JCI Assumption 1;• FCI-JCI12 is sound if M also satisfies JCI Assumptions 1 and 2;• FCI-JCI123 is sound and complete if M also satisfies JCI Assumptions 1, 2, 3.
Here, sound means that the output of the algorithm is a DPAG that contains the true
DMAG(M), and complete means that all edge marks of the true DMAG(M) that canbe identified from the (conditional) independences in PM (X, C) and the JCI backgroundknowledge have been oriented in the DPAG output by the algorithm.
20. Of course, the power of the conditional independence test may be reduced when conditioning on many
variables. However, if the number of contexts in the pooled data is small, one can design conditional
independence tests that do not suffer from this problem.

40

Joint Causal Inference from Multiple Contexts

Proof We refer the reader to Appendix B, and in particular to Theorems 35, 37, and 38.
In practice, one often does not have access to the joint distribution PM (X, C), but only toa finite sample of it. In that case we have:
Corollary 21 The FCI variants mentioned in Theorem 20 are also asymptotically consistent under the assumptions stated in Theorem 20 if the conditional independence test(including the choice of the threshold to decide between independence and dependence) is
consistent.
Proof Direct application of Lemma 36 in Appendix B.

4.3. Related WorkIn this section we provide a more detailed comparison with related work. Since the pioneering work by Fisher (1935), many different causal discovery methods that can deal withdata from different contexts have been proposed. Table 4 provides an overview of some ofthese methods and the features they offer. Note that JCI offers most features of all methods. By implementing the JCI framework using sophisticated causal discovery methods forobservational data (plus background knowledge) one obtains versatile and powerful causal
discovery algorithms for multiple contexts. We will now discuss in detail some of the aspectsof the related work.
4.3.1. Latent Confounders
Most score-based methods that combine multiple contexts (like the ones by Cooper andYoo, 1999; Tian and Pearl, 2001; Sachs et al., 2005; Eaton and Murphy, 2007; Hauser andBühlmann, 2012; Mooij and Heskes, 2013; Oates et al., 2016a) and some constraint-based
methods (Zhang et al., 2017; Yang et al., 2018) assume causal sufficiency, i.e., that nolatent confounders are present. This simplifies the causal discovery problem considerably,
but the assumption is likely violated in practice and may lead to wrong conclusions. Thisis well-known for causal discovery from a single observational data set, but also applies to
the JCI setting.
4.3.2. Cycles
As we have seen in Proposition 10, the method by Fisher (1935) can handle cycles. Less
well-known is that also LCD (Cooper, 1997) and Trigger (Chen et al., 2007) can handle
cycles (see Proposition 16). Hyttinen et al. (2012) provide an algorithm for linear SCMs
with cycles and confounders that deals with perfect interventions. The methods by Hyttinen
et al. (2014) and Mooij and Heskes (2013) can deal with cycles in a linear (or approximately
linear) setting. The method by Hyttinen et al. (2014) relies on d-separation, which only
applies in certain settings (see Theorem 9). The method can be modified to use σ-separation
instead (Forré and Mooij, 2018). The way Mooij and Heskes (2013) handle cycles is not
as straightforward. Generally, their method could handle nonlinear cyclic models, but for41

Nonlinear mechanisms

CyclesPerfect interventions

Mechanism changes

Activity interventions

Other context changes

Unknown intervention/context targets

Learns intervention/context targets

Global causal discovery

Different variables in each context

Combination strategy(Fisher, 1935)
LCD (Cooper, 1997)(Cooper and Yoo, 1999)(Tian and Pearl, 2001)(Sachs et al., 2005)(Eaton and Murphy, 2007)
Trigger (Chen et al., 2007)(Claassen and Heskes, 2010)(Tillman and Spirtes, 2011)(Hauser and Bühlmann, 2012)(Hyttinen et al., 2012)(Mooij and Heskes, 2013)(Hyttinen et al., 2014)(Triantafillou and Tsamardinos, 2015)(Rothenhäusler et al., 2015)(Oates et al., 2016a)
ICP (Peters et al., 2016)(Zhang et al., 2017)(Yang et al., 2018)(Forré and Mooij, 2018)
Joint Causal Inference (this work)
FCI-JCI (this work)
ASD-JCI (this work)

Latent confounders

Mooij, Magliacane and Claassen

+
+
+
+
+
+
+
+
+
+
+
+
+
+

+
+
+
+
+
+
+
+
+
+
±
+
+
+
+
+
+
+
+
+

+
+
+
+
±
±
±
+
+
+
?
+

+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

+
+
+
+
+
+
+
+
+
+
+
+
+
+

+
+
+
+
+
+
+
+
+
+
+
+

+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

+
+
+
+
+
+
+
+
+
+
+
+

+
+
+
+
+
+
+

+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+

+
+
+
+
+
-

b
b
b
b
b
b
b
a
a
b
a
b
a
a
a
b
b
b
a/b
a
b
b
b

Table 4: Overview of causal discovery methods that can combine data from multiple contexts. Features offered by the original implementations of these methods are indicated. When a feature is offered only under additional restrictive assumptions,
it is indicated with a ± sign. Combination strategies (right-most column) are:(a) obtain statistics or constraints from each context separately and then construct a single causal graph based on the combined statistics, (b) pool all dataand construct a single causal graph directly from the pooled data.

42

Joint Causal Inference from Multiple Contexts

computational reasons, their implementation linearizes the SCMs around each (contextdependent) equilibrium, thereby basically assuming that d-separation holds within eachcontext. The method by Rothenhäusler et al. (2015) assumes linearity and can deal withcycles in that case, under a certain condition that suffices to prove identifiability of themethod. The method by Peters et al. (2016) can handle cycles, as our Corollary 17 shows.
The JCI framework in general allows for cycles, but requires its implementation to support
this.
4.3.3. Selection Bias
The only causal discovery method for multiple data sets that is explicitly claimed to be able
to deal with selection bias (i.e., conditioning on a latent variable that is a common effectof one or more of the observed variables), at least to some extent, is the IOD algorithm(Tillman and Spirtes, 2011). It allows for different sets of observed (system) variables in eachcontext and for different distributions in each context, while assuming that each context
can be described by a MAG that is the marginal of a common MAG defined on the union of
all system variables. It performs conditional independence tests in each data set separately,
and merges the p-values of the test results using Fisher’s method. It then constructs the
PAG that represents simultaneously all contexts. Since it does not assume invariance ofthe distribution across contexts, it can deal with a single (latent) context variable thatmodels mechanism changes or other “soft” interventions that do not change the conditional
independences in the distribution. It can also deal with perfect interventions since Fisher’s
method is used to test for independence in all contexts (see also Figure 7).
4.3.4. Imperfect Interventions and Other Context Changes
Cooper and Yoo (1999) provided the first score-based causal discovery algorithm that could
deal with data from multiple contexts, focusing on perfect interventions with known targets.
They describe in detail how to handle perfect interventions and introduced the idea ofadding explicit context variables to deal with mechanism changes, which was later refinedby Eaton and Murphy (2007), who provide an algorithm that can handle (stochastic) perfectinterventions with unknown targets, soft interventions, and mechanism changes. Also Sachs
et al. (2005) use a score-based causal discovery algorithm based on the ideas of Cooper andYoo (1999) that uses a greedy search strategy through the space of DAGs.
Another recent approach to constraint-based causal discovery in a JCI setting is the one
by Yang et al. (2018); these authors propose the algorithm IGSP that can be seen as animplementation of JCI for causally sufficient, acyclic models with a diagonal experimental
design under JCI Assumptions 1, 2, for mechanism changes with intervention targets assumed to be known. An advantage of IGSP over our JCI approach is that essentially no
assumptions on the context distribution need to be made (apart from positivity) since it
relies on a weaker faithfulness assumption.
Tian and Pearl (2001) were the first to consider mechanism changes. They deal withsequences of mechanism changes, exploiting changes in the distribution to infer descendantsof the changed mechanism. This is followed by a constraint-based approach from observational data that also takes into account the background knowledge on the causal orderingof the system variables inferred from analysing the interventional data. A similar approach
43

Mooij, Magliacane and Claassen

using the differences between data from experimental conditions and an observational baseline as background knowledge for a constraint-based approach was applied by Magliacane
et al. (2016b) on the data of Sachs et al. (2005).
Claassen and Heskes (2010) handle certain environment changes: direct causal relations between system variables are assumed to be invariant across contexts, but latent
confounding (and more generally, the exogenous distribution) may differ between contexts.
Rothenhäusler et al. (2015) assume stochastic shift interventions in which the mean of atarget variable is shifted by an (independent) random amount. Various multi-task “structure learning” (i.e., Bayesian network learning) approaches that put a prior on the similarityof the DAGs in multiple contexts which encourages them to be similar have been proposed(e.g., Oates et al., 2014, 2016b).
Some methods which allow for a single context variable have been applied in settings ontime-series data, by using time as the context variable (Friedman et al., 2000; Zhang et al.,2017). This extends the more usual approach of treating time-series data by assuming
invariance of the causal structure across time as in dynamic Bayesian networks (DBNs)(Murphy, 2002), methods based on Granger causality (Granger, 1969), or constraint-based
approaches (Entner and Hoyer, 2010).
JCI allows one to handle all interventions and context changes discussed above in a
unified way.
4.3.5. Multiple Context Variables
Some causal discovery methods for combining data from different contexts that explicitly
consider a context variable, allow for a single context variable only, for example, LCD,
ICP, and the method by Zhang et al. (2017). There is an important advantage to allowing
multiple context variables, as JCI does generally. One might argue that the case of multiple
context variables can always be reduced to a case with a single context variable, by simply
combining all context variables {Ck }k∈K into a single tuple C = (Ck )k∈K . However, thisreduction to a single context variable typically loses information. This is illustrated in
Figure 13. When using only a single context variable in that case, the DMG cannot be
identified from conditional independences in the data. On the other hand, when using all
three context variables with JCI, the complete DMG can be identified, even when the causal
relations between context and system variables are unknown.
4.3.6. Dependent Context Variables
If one allows for multiple context variables and considers the joint distribution on context
and system variables, as we do in JCI, one should account for possible dependencies between
the context variables. Indeed, incorrectly assuming the context variables to be independenta priori may lead to wrong conclusions. An example is provided in Figure 14. In thatexample, incorrectly assuming the contexts to be independent leads to the wrong conclusion that context variable Cβ causes the system variable X0 , at least for causal discovery
algorithms that are tolerant to faithfulness violations.
This issue was recognized and addressed in recent work (Oates et al., 2016a) by introducing a novel graphical modeling framework, Conditional DAGs (CDAGs), which bears
some similarity with our approach. However, a disadvantage of the CDAG framework is
44

Joint Causal Inference from Multiple Contexts

C
X1

X2

X3

Cα

Cβ

Cγ

X1

X2

X3

Figure 13: Example that shows that allowing multiple context variables (right) has advantages over considering a single context variable only (left). The causal graph on
the right can be identified by JCI (with JCI Assumptions 1, 2, 3) from conditional independences in the pooled data, whereas the causal graph on the left is
not identifiable.

Cα

Cβ

Cα

X1

Cβ

X1

Figure 14: Example that shows that incorrectly assuming independent context variables canlead to wrong conclusions when using causal discovery algorithms that are tolerant to faithfulness violations. Left: true causal graph, with dependent context
variables, which is identifiable by JCI. Right: causal graph that reproduces allconditional dependencies (except Cα ⊥6 ⊥ Cβ ) and minimizes the number of faithfulness violations (faithfulness here implies Cβ ⊥6 ⊥ X1 | Cα ), when (incorrectly)
assuming that context variables are independent.

that existing causal discovery methods cannot be directly applied to learn a CDAG from
data, and the wealth of results on causal modeling with SCMs cannot be used directly. Oneof the key advantages of the JCI framework is that it utilizes existing theory and methods,
as it reduces a causal discovery problem from multiple contexts to a purely observational
one with background knowledge. This is one of the reasons why JCI offers many morefeatures than the approach by Oates et al. (2016a). We also note that CDAGs can be dealt
with as a special case of the JCI framework.
4.3.7. Partially Overlapping Sets of Variables
A particular case of missing data that has been addressed by some of the methods is when theset of observed variables differ between data sets, while still having some overlap. The first
one to address this using constraint-based causal discovery was Tillman (2009), and several
other methods have been proposed over the years (Claassen and Heskes, 2010; Tillman andSpirtes, 2011; Hyttinen et al., 2014; Triantafillou and Tsamardinos, 2015). JCI can only
deal with this when strengthening its faithfulness assumption: one would need to assume
that the context variables are discrete, and that every conditional distribution P(X | C = c)
for c ∈ C with P(C = c) > 0 is faithful with respect to the marginalization of the same45

Mooij, Magliacane and Claassen

DMG GI on the system variables that were observed in that context (or, in case of perfectinterventions with known targets, the corresponding marginal intervened DMG).
4.3.8. Influence Diagrams
Our representation of a system within a context imposed by its environment bears strong
similarities with influence diagrams (Dawid, 2002). A formal difference is that we consider
the context variables to be random variables that reflect the empirical distribution of theexperimental design, whereas in influence diagrams they are interpreted as non-random
decision variables. The advantage of treating context variables as random variables isthat this allows one to apply standard causal discovery techniques (designed for randomvariables) jointly on system and context variables. In particular, the standard notion of
statistical conditional independence (Dawid, 1979) suffices. If one would like to treat thecontext variables as decision (i.e., non-random) variables, extended notions of conditional
independence would be necessary (Forré and Mooij, 2019). Since we can always view thecontext variables as random variables in the empirical distribution of the experimental
design (see also Footnote 17), this allows us to make use of the standard notion of conditional
independence for the purposes of causal discovery.
4.3.9. Selection Diagrams
Our representation of a system within a context imposed by its environment also bears somesimilarities with selection diagrams (Bareinboim and Pearl, 2013). Selection diagrams havealso been used for causal modeling in different contexts, but one crucial difference is that
we are modeling the joint distribution on the intervention and system variables, whereas aselection diagram represents the conditional distribution of the system variables given theintervention (“selection”) variables. Because we are modeling the joint distribution and not
only the conditional one, we can apply standard causal discovery techniques directly onpooled data, something that would not be as trivial when using selection diagrams instead.
Bareinboim and Pearl (2013) define a selection diagram as follows:
Definition 22 Let M = hI, J , H, X , E, f , PE i and M∗ = hI, J , H, X , E, f ∗ , P∗E i be two
acyclic SCMs corresponding to two different contexts, that only differ with respect to theircausal mechanisms and exogenous distributions. In particular, they share the same augmented graph H and hence also their graphs are identical, G(M) = G(M∗ ). The selection
˙ wherediagram S induced by hM, M∗ i is the acyclic directed mixed graph with nodes I ∪I,
I := {ī : i ∈ I} is a copy of I of selection variable indices, such that
1. the induced subgraph of S on I equals the common graph of M and M∗ , i.e., SI =G(M) = G(M∗ ), and
2. for each i ∈ I such that fi 6= fi∗ or PE pa

G(M) (i)6= P∗E pa

G(M) (i)there is an edge ī → i in

S.
From the definition, it is apparent that a selection diagram essentially models two contexts, and that the selection variables in the selection diagram correspond to the childrenof the context variable in our representation. Indeed, consider a JCI model of the form(5) (p. 21) with a single binary context variable C. The joint SCM can be split into two
46

Joint Causal Inference from Multiple Contexts

context-specific SCMs, M0 and M1 , and the induced selection diagram D can be obtained
from the causal graph G(M) as follows: (i) each edge i1 → i2 or i1 ↔ i2 in G(M) between
system variables i1 , i2 ∈ I is also in D; (ii) if C → i in G(M) for i ∈ I then ī → i is in D.
Since the JCI framework can be used to learn (features of the) causal graph G(M) from
data, this means that we can thereby learn (features of) the selection diagram from data.
It is not clear how a selection diagram could be used to represent the same information that an SCM with multiple context variables can represent. Indeed, even though theselection diagram has multiple selection variables, it is still modeling only two contexts,
corresponding with just a single binary context variable in the JCI framework.

5. Experiments
In this section we report on the experiments we performed with JCI, comparing various implementations of the framework with several baselines and state-of-the-art causal discovery
methods. We experimented both with simulated data with perfectly known ground truth
and with real-world data where the ground truth is only known approximately. The sourcecode that we used for producing the results and plots in this section is provided under afree and open source license as Online Appendix 1.
5.1. Methods and Baselines
In our experiments we study different implementations of JCI, based on two existing causal
discovery algorithms: ASD (Hyttinen et al., 2014; Magliacane et al., 2016b; Forré and Mooij,
2018) and FCI (Spirtes et al., 1999; Zhang, 2008a). The ASD algorithm is accurate but
slow, while FCI is faster but less accurate due to its “greedy” approach. Another difference
between both methods is that ASD can deal with partial inputs, while for FCI it is necessary to provide all independence test results it asks for. Although FCI (and with somesmall extensions, also ASD) can deal with selection bias, we ignore this additional complication here and use simplified implementations that assume that there is no selection bias.
For the cyclic case, we used the adaptation of the ASD algorithm proposed by Forré and
Mooij (2018) that replaces d-separation with its general cyclic generalization, σ-separation.
Adapting FCI to the cyclic case seems less straightforward and is beyond the scope of thispaper.21
Table 5 provides an overview of all implementations that we have studied here. The
methods will be discussed in more detail in the next few subsections. The “CI Tests” column
in the table describes what conditional independence test are performed, and how, and canhave the following values:
A : use all variables, including context variables; the conditional independence tests performed are of the form X̃a ⊥⊥ X̃b | X̃S with {a} ∪ {b} ∪ S ⊆ I ∪ K and {a}, {b}, S
mutually disjoint.
S : use only system variables; the conditional independence tests performed are of the form
Xa ⊥⊥ Xb | XS with {a} ∪ {b} ∪ S ⊆ I and {a}, {b}, S mutually disjoint.
21. The d-separation case has recently been addressed by Strobl (2018).

47

Mooij, Magliacane and Claassen

SS : system variables only, separately for each context; the conditional independence tests
performed are of the form Xa ⊥⊥ Xb | XS , C = c for {a} ∪ {b} ∪ S ⊆ I and {a}, {b}, S
mutually disjoint, for all values c ∈ C with P(c) > 0. Note that this method assumes
that the context domain C is discrete.
SF : system variables only, separately for each context, using Fisher’s method;22 theconditional independence tests performed are of the form Xa ⊥⊥ Xb | XS , C = c for{a} ∪ {b} ∪ S ⊆ I and {a}, {b}, S mutually disjoint, for all values c ∈ C with P(c) > 0.
Note that this method assumes that the context domain C is discrete.
NC : test all variables, except for conditional independences between context variables; theconditional independence tests performed are of the form Xa ⊥⊥ X̃b | X̃S with a ∈ I,{b} ∪ S ⊆ I ∪ K and {a}, {b}, S mutually disjoint.
PF : test all pairs of context and system variables, conditioning on the remaining context variables; the conditional independence tests that are performed are of the form
Xi ⊥⊥ Ck | CK\{k} for i ∈ I, k ∈ K.
5.1.1. ASD Variants
For ASD, we implemented different variants as described in Table 5. The variants ASD-obs,
ASD-pooled and ASD-pikt use the original implementation of Hyttinen et al. (2014) inthe acyclic case and the σ-separation adaptation of Forré and Mooij (2018) in the cyclic
case, with the weights and query method of Magliacane et al. (2016b). ASD-obs onlyuses the observational context and ignores data from the other contexts, ASD-pooled pools
data from all contexts but does not add context variables, and ASD-pikt uses data fromall contexts and assumes that the contexts correspond to perfect interventions with known
targets. Inspired by the approach of Tillman (2009) and Tillman and Spirtes (2011), we also
implemented a variant ASD-meta that uses Fisher’s method as a “meta-analysis” method
to combine the p-values of conditional independence tests performed with data from eachcontext separately into a single overall p-value, which was then used as input for ASD. We
use these implementations as state-of-the-art causal discovery methods for comparison.
For the JCI approach, we study several variants that differ in terms of which JCI assumptions they make, whether they also test conditional independences between context
variables, whether all context variables are used or whether they are first merged into a single
context variable, and whether the intervention targets of the context variables are considered to be known or not. The different implementations are described in detail in Table 5.
ASD-JCI123sc and ASD-JCI1sc both use a single merged context variable C = (Ck )k∈K ,whereas the other ASD-JCI variants use all context variables {Ck }k∈K as separate variables. The background knowledge for all ASD-JCI variants is a subset of the three JCI
Assumptions 1, 2 and 3. The only ASD-JCI variant that uses background knowledge onintervention targets is ASD-JCI123kt (but it does not make any assumptions on the type of
the intervention).
22. Fisher’s method (Fisher, P
1925) aggregates N independent p-values {pi }N
i=1 by computing the p-value
2
for the statistic F := −2 N
i=1 log pi , which has a χ distribution with 2N degrees of freedom if all N
p-values pi are independent.

48

Joint Causal Inference from Multiple Contexts

5.1.2. FCI Variants
We implemented different variants of the FCI algorithm by adapting the implementation in
the R package pcalg (Kalisch et al., 2012). We used the default configuration, i.e., an orderindependent (“stable”) skeleton phase, and no conservative or majority rule modifications(Colombo and Maathuis, 2014). For simplicity, we assumed that no selection bias is present,
which means that the rules R5–R7 in Zhang (2008a) can be ignored in the FCI algorithm,
and only PAGs without (possibly) undirected edges need to be considered. We consider two
variants of FCI as the current state-of-the-art: FCI-obs, which uses only the observational
context, and FCI-pooled, which uses pooled data from all contexts but does not add
context variables. We also implemented a “meta-analysis” approach FCI-meta that uses
Fisher’s method to combine the p-values from separate contexts into overall p-values thatare used as input for the FCI algorithm. Finally, we have three variants (FCI-JCI123,
FCI-JCI1 and FCI-JCI0) referring to the JCI adaptation of the FCI algorithm as describedin Section 4.2.4, for three different combinations of JCI Assumptions (we have not yet
implemented and evaluated the speedup for FCI-JCI123r but leave that for future work).
The output of FCI is a PAG. Here we will not evaluate the PAG itself, since estimating
the PAG is often not the ultimate task in causal discovery, but instead we will evaluatepresence and absence of ancestral relations that can be identified from the estimated PAG.
The procedure we use for this task is explained in Appendix B.4. For the special case of
FCI-JCI123, we also read off the direct intervention targets (and non-targets) from theestimated PAG, as explained in Appendix B.5. We encode identified presence of a featureby +1, identified absence by −1, and an unidentifiable feature by 0. These predictions canthen also easily be bootstrapped (or more precisely, bagged). The bagged feature scores
can then be used as a score for the confidence that the feature is present.
5.1.3. LCD Variants
The LCD implementation simply iterates over all context variables and ordered pairs ofsystem variables and tests for the LCD pattern. As conditional independence test we
test whether the partial correlation vanishes. As confidence measure for an LCD pattern hC, X, Y i (see also Figure 12, where hX1 , X2 , X3 i corresponds with hC, X, Y i), we use− log pC ⊥⊥ Y . Note that LCD predicts the presence of an ancestral relation X ∈ an(Y ), theabsence of a confounder between X and Y , and the absence of a direct causal effect of C
on Y .23
5.1.4. ICP Variants
We also compare with the ICP function in the R package InvariantCausalPrediction(Peters et al., 2016).24 Under the assumption of causal sufficiency, ICP returns directcausal relations. However, in our setting, in which causal sufficiency cannot be assumed,
ICP will generally return ancestral causal relations, as shown in Corollary 17. Note that
ICP assumes a single context variable, whereas one typically may have multiple context
23. The absence of a bidirected edge X ↔ Y in the marginalization of the graph G on {C, X, Y } implies also
the absence of the bidirected edge X ↔ Y in the full graph G. Furthermore, the absence of the directedge C → Y in this marginalization implies the absence of the direct edge C → Y in the full graph G.
24. We used the default arguments, except that we set stopIfEmpty to TRUE.

49

Mooij, Magliacane and Claassen

variables in the data. One way to deal with this is to merge all context variables before thepooled data is input to ICP, which we do in ICP-sc. Another way is to run ICP on eachcontext variable separately, hiding the other context variables when feeding the pooled data
to ICP, and finally merging all predictions. This is done in ICP-mc.
The conditional independence tests that ICP performs are of the form C ⊥⊥ Xi | XS
with S ⊆ I \ {i}, for i ∈ I. The default conditional independence test used in ICP firstlinearly regresses Xi on XS for each context C = c individually, and once globally. It thentests whether there exists a context c in which the mean or the variance of the regression
residuals is different from the global mean or variance of the residuals. All p-values are then
combined using a Bonferroni correction. Note that this test can detect more conditional
dependencies than a simple partial correlation test, as it also considers the variation between
the variances across contexts.
As a confidence score for the ancestral causal relation i ∈ an(j), we use − log pi99Kj ,where pi99Kj is the p-value returned by ICP for system variable i being ancestor of systemvariable j.
5.1.5. Fisher’s Test for Causality
This is a very simple and immensely popular baseline in which we simply go through allpairs (i, k) of a system variable i ∈ I and a context variable k ∈ K, perform the conditional
independence test Xi ⊥⊥ Ck | CK\{k} on the pooled data resulting in p-value pik . It is limited
to discovery of ancestral causal relations from context to system variables. As confidence
value for the ancestral causal relation k ∈ an(i) we report − log pik + log α, where α is thethreshold for the independence test.
5.1.6. Bootstrapping
A simple way to improve the stability of causal discovery algorithms is bootstrapping. Fora method that outputs a confidence measure for a certain prediction we simply average theconfidence measures over bootstrap samples. For FCI, as confidence measure we simply takea {−1, 0, 1}-valued variable encoding the identifiable absence/unidentifiability/identifiable
presence of an ancestral relation (or direct intervention target, for FCI-JCI123). For LCD,
we average − log pC ⊥⊥ Y over bootstrap samples. For ICP, we similarly average the negative
logarithm of the p-values for the discovered ancestral causal relations. We do not bootstrap
ASD variants because of the high computational complexity. In our experiments, we use
100 bootstrap samples. Bootstrapped methods are indicated with a suffix “-bs”.
5.1.7. Conditional Independence Tests
Using an appropriate conditional independence test is important to obtain good causal
discovery results. In this work we will use two different conditional independence tests,
both relying on the assumption that the context variables are discrete and the systemvariables have a multivariate Gaussian distribution given the context. However, the JCI
framework imposes no principled restrictions on the conditional independence tests used,
so one could also use non-parametric tests instead, for example.
The default conditional independence test that we used is the following. For testingX̃a ⊥⊥ X̃b | X̃S , we distinguish two cases:
50

Joint Causal Inference from Multiple Contexts

Name

Data

Baselines
ASD-obs
observational
ASD-pooled
pooledASD-metaall
ASD-piktall
FCI-obs
observational
FCI-pooled
pooledFCI-metaall
JCI Implementations
ASD-JCI0
pooledASD-JCI1
pooledASD-JCI12
pooledASD-JCI123
pooledASD-JCI123kt
pooledFCI-JCI123
pooledFCI-JCI1
pooledFCI-JCI0
pooledLCD-scpooledICP-scpooledASD-JCI1-sc
pooledASD-JCI123-sc pooled
LCD-mcpooledICP-mcpooledFisherpooledContext variables

JCI Assumptions
1 2 3Interventions

CI Tests

none
none
none
none
none
none
none

−
−
−
−
−
−
−

−
−
−
−
−
−
−

−
−
−
−
−
−
−

none
any
any
perfect, KT
none
any
any

S
S
SF
SS
S
S
SF

all
all
all
all
all
all
all
all
single (merged)
single (merged)
single (merged)
single (merged)
all (one-by-one)all (one-by-one)all (one-by-one)

−
+
+
+
+
+
+
−
+
+
+
+
+
+
+

−
−
+
+
+
+
−
−
−
−
−
+
−
−
+

−
−
−
+
+
+
−
−
−
−
−
+
−
−
−

any
any
any
any
any, KT
any
any
any
any
any
any
any
any
any
any

A
A
A
NC
NC
NC
A
A
A∗
A∗
A∗
A∗
A
A∗
PF

Table 5: Variants of implemented JCI algorithms and baselines used in our experiments.
JCI Assumption 0 is always assumed. “KT” is an abbreviation of “known targets”.
The meaning of “CI Tests” is (more detailed explanation in main text): A: use
all variables, including context variables; S: use only system variables; SS: systemvariables only, separately for each context; SF: system variables only, separately
for each context, using Fisher’s method to combine them into a single p-value; NC:
test all variables, except for conditional independences between context variables;
PF: test all pairs of context and system variables, conditioning on the remaining
context variables. The meaning of the superscript ∗ is explained in Section 5.1.7.
Bootstrapped versions of methods will be indicated with a suffix “-bs” (and havebeen omitted from this table for clarity).

51

Mooij, Magliacane and Claassen• S ∩ K = ∅: the test then reduces to a standard partial correlations test.
• Otherwise, we go through all observed values cS∩K of CS∩K , and use a standard
partial correlations test to calculate a p-value for X̃a ⊥⊥ X̃b | X̃S\K , CS∩K = cS∩K . We
then aggregate the p-values corresponding to observed values of CS∩K into one overall
p-value for X̃a ⊥⊥ X̃b | X̃S using Fisher’s method for aggregating p-values.25
Note that for the case with zero context variables (i.e., a single context), this test reduces
to a standard partial correlations test.
For the ICP implementations, we make use of the implementation in the R package
InvariantCausalPrediction. This by default makes use of a conditional independence test
that uses more than just partial correlations. We extended this conditional independence
test to allow for conditioning on context variables, and make use of this extended test in
all the algorithms that assume a single merged context, i.e., the ones marked with a ∗in Table 5. It assumes that there is a single context variable, i.e., |K| = 1. For testingX̃a ⊥⊥ X̃b | X̃S , it distinguishes the following cases:• If ({a} ∪ {b} ∪ S) ∩ K = ∅, it reduces to a standard partial correlations test.
• If X̃a = C is the context variable, it uses linear regression to fit X̃b as a linear functionof XS , using data pooled over all contexts, and calculates the corresponding residuals.
It then goes through all observed values c of C, and tests whether the residuals incontext c have a different distribution than the residuals in the other contexts (i.e.,with C 6= c). This two-sample test is performed by comparing the means by a ttest, and the variances by an F -test. The two resulting p-values are combined witha Bonferroni correction. The resulting p-values, one for each context, are then also
combined with a Bonferroni correction.26• If X̃b = C is the context variable, proceed similarly as in the previous case.
• If the context variable C is part of X̃S , we go through all observed values c of C, anduse a standard partial correlations test to calculate a p-value for Xa ⊥⊥ Xb | XS\K , C =c. We then aggregate the p-values corresponding to observed values of C into one
overall p-value for Xa ⊥⊥ Xb | X̃S using Fisher’s method for aggregating p-values.
As a final note regarding the conditional independence tests, we state that ASD-pikt
uses a standard partial correlation test to calculate a p-value for each context separately. It
subsequently combines all these p-values by adding the log p-values, but taking into account
how the graph structure is changed through perfect interventions.
The choice of the p-value threshold α for rejecting the null hypothesis of independence
is an important one. To obtain consistent results, one should let α decrease to 0 withincreasing sample size. In our experiments, we used fixed sample size and simply used a
global threshold α = 0.01.
25. Not to be confused with Fisher’s test for causality that we described in Section 5.1.5.
26. This is the default test for continuous data in the ICP function of theInvariantCausalPrediction.

52

R

packageJoint Causal Inference from Multiple Contexts

5.2. Simulations
We simulated random linear-Gaussian SCMs with p system variables and q context variables. We considered both the acyclic setting and the cyclic one. We simulated stochastic
interventions of two different intervention types: mechanism changes, and perfect interventions.
Random causal graphs were simulated by drawing directed edges independently between
system variables with probability . For the acyclic models, we only allowed directed edges
i1 → i2 for i1 < i2 with i1 , i2 ∈ I. For cyclic models, we allowed directed edges i1 → i2
for i1 6= i2 with i1 , i2 ∈ I, and subsequently selected only the graphs in which at least one
cycle exists. We drew bidirected edges independently between all unordered pairs of systemvariables with probability η, and associated each bidirected edge with a separate latent
confounding variable. For each context variable, we randomly selected a single systemvariable as its target, while ensuring that each system variable has at most one context
variable as its direct cause. We sampled all linear coefficients between system variables,
context variables and confounders from the uniform distribution on [−1.5, −0.5] ∪ [0.5, 1.5].
The exogenous variables (“error terms”) were sampled independently from the standardnormal distribution. To ensure that system variables have comparable scales, we rescaled
the weight matrix such that each system variable would have variance 1 if all its directcauses would be i.i.d. standard-normal.
We used binary context variables in a “diagonal” design. This means that for eachrandom SCM, we simulated q + 1 contexts, with the first context being purely observational(i.e., Ck = 0 for all k ∈ {1, . . . , q}), and the other q contexts corresponding with one of thecontext variables turned on (say Ck0 = 1 for some k 0 ∈ {1, . . . , q}) and the others turnedoff (Ck = 0 for the other k ∈ {1, . . . , q} \ {k 0 }). We either took all interventions to bemechanism changes, or all interventions to be perfect. For mechanism changes, we simply
add the value of the parent context variable to the structural equation (i.e., this correspondswith adding a constant offset of 1 to the intervention target variable when the intervention isturned on). For perfect interventions, we additionally set the linear coefficients of incomingedges on the intervention target to zero. Finally, we sampled N observed values of systemvariables from each context and combined all samples into one pooled data set. This wasdone for each random SCM separately.
5.3. Evaluation
In evaluating the results, we consider different prediction tasks: establishing the absence orpresence of ancestral causal relations between system variables, the absence or presence of
direct causal relations between system variables, and the absence or presence of confounders
between system variables. In addition, we consider predicting the absence or presence of
indirect intervention targets (i.e., whether or not some context variable is ancestor of somesystem variable) and of direct intervention targets (i.e., whether or not some context variable
is parent of some system variable).
Each method outputs a confidence score for each feature of interest, where positivescores mean that it is more likely that a feature is present in the causal graph G, whereas
negative scores mean that it is more likely that a feature is absent. The higher the absolutevalue of the score, the more likely its presence or absence is. The predictions are pooled
53

Mooij, Magliacane and Claassen

both within model instances (e.g., all possible ancestral relations i ∈ anG (j) for all orderedpairs of system variables i, j ∈ I) and across model instances to gather more statistics.
The scores are then ranked and turned into ROC curves and PR curves (one PR curve forthe presence, and one for the absence of the features) by comparing with the true features.
In the ROC curves, we use solid lines for positive scores (feature present) and negativescores (feature absent), and dotted lines for vanishing scores (feature presence/absence isunknown).
5.4. Results: Small Simulated Models
We first present results for small models with p = 4 system variables and 0 ≤ q ≤ 4 (as adefault, q = 2) context variables. We used  = 0.5, η = 0.5, and sampled Nc = 500 samples
for each context, i.e., N = 500(q + 1) samples in total.
5.4.1. ASD-JCI vs. Baselines (Causal Mechanism Changes)
We start by showing off the advantage that JCI can offer over existing methods. We
first consider only ASD variants because this most clearly shows the impact of how one
merges data from different contexts and how one treats the context variables, since theother aspects of the causal discovery algorithm are the same for all ASD variants. In
Figure 15 we present results for several ASD variants for acyclic models with causal mechanism changes. We compare the JCI variants ASD-JCI123 (unknown intervention targets)
and ASD-JCI123kt (known intervention targets) with the available baselines, ASD-obs (observational data only), ASD-pooled (pooled data from all context treated as if they were allobservational, context variables not included), ASD-meta (using Fisher’s method to combine
p-values from conditional independence tests in separate contexts), and ASD-pikt (whichassumes that interventions are perfect and uses knowledge of intervention targets). The
tasks of predicting ancestral causal relations and direct causal relations show relatively
similar ROC and PR curves. Predicting the absence or presence of confounders is a morechallenging task.
The three baselines ASD-obs, ASD-pooled and ASD-meta show very similar performance
behaviors. In particular, for the tasks of predicting the presence of the features, these baselines perform poorly, not much better than random guessing. This is partially due to thesmall sample size, but also to the fact that many relationships are simply not identifiablefrom purely observational data alone. ASD-pikt even performs poorly on nearly all prediction tasks in this simulation setting because it incorrectly assumes that the interventions
are perfect. The two JCI variants, on the other hand, strongly outperform the baselines
and obtain very high precisions. In particular, even without knowing the intervention targets, ASD-JCI123 manages to predict the presence of (direct and indirect) causal relations
at maximum precision for low recall. Exploiting knowledge of the intervention targets,
ASD-JCI123-kt obtains an even more impressive precision for predicting ancestral causal
relations for a large range of recalls. This illustrates the significant improvement in precisionthat JCI can yield.
Figure 16 shows a similar picture in the cyclic setting, where all ASD variants make useof σ-separation (Forré and Mooij, 2018). The task of predicting the presence of ancestral
relations is easier than in the acyclic setting, because for most pairs of system variables, one
54

Joint Causal Inference from Multiple Contexts

Figure 15: Results of some ASD variants (acyclic, causal mechanism changes) forsmall models. The two JCI variants (with unknown/known intervention targets)strongly outperform the baselines in this setting. From top to bottom: ancestral
relations, direct causes, confounders. From left to right: ROC curves, PR curves
for presence of feature, PR curves for absence of feature.

is ancestor of the other due to the cycles. The task of predicting the absence of ancestral
relations, on the other hand, is more challenging. Detecting the presence or absence ofconfounders in this setting has become nearly impossible, for any method. For the othertasks, the JCI approach again shows substantially improved precisions compared to thebaselines.
5.4.2. ASD-JCI vs. Baselines (Perfect Interventions)Figures 17 and 18 show results for respectively the acyclic and cyclic settings, but now for
perfect interventions with known targets rather than causal mechanism changes.
In these perfect intervention scenarios, the JCI variants again obtain a much higher
precision than any of the baselines, with the sole exception of ASD-pikt. In this setting,the latter method successfully exploits the assumed perfect nature of the interventions,55

Mooij, Magliacane and Claassen

Figure 16: Results of some ASD variants (cyclic, causal mechanism changes) forsmall models. The two JCI variants substantially outperform the baselines inthis setting.

thereby outperforming the JCI variants that do not make any assumption about the natureof the intervention. However, a significant disadvantage of ASD-pikt in practice is that its
assumption of perfect interventions with known targets may not be valid. As we alreadysaw in Section 5.4.1, ASD-pikt then breaks down, in contrast to the JCI variants.
For the cyclic setting with perfect interventions (Figure 18), we observe that predicting
the presence of confounders still seems impossible for all methods, but predicting theirabsence seems at least feasible in principle (although it seems a very challenging task).
ASD-pikt again obtains the highest precisions, followed by the JCI variants.
5.4.3. Influence of JCI Assumptions
We now investigate in more detail which JCI assumptions are responsible for the excellent
performance of the JCI variants of ASD. Figure 19 shows that, as expected, the more priorknowledge about the context variables is used, the better the predictions become.
56

Joint Causal Inference from Multiple Contexts

Figure 17: Results of some ASD variants (acyclic, perfect interventions) for small
models. ASD-pikt, which takes into account the perfect nature of the interventions, is the best performing method in this setting. The JCI variants do notassume perfect interventions, but still yield a vast improvement of the precisionof the predicted features with respect to the other baselines.

However, surprisingly, the largest boost in precision with respect to the observational
baseline is due to simply pooling the data and adding the context variables: ASD-JCI0
already strongly improves over ASD-obs. Adding more background knowledge regarding
the nature of the context variables helps to improve the results further. JCI Assumption 1
yields a marginal improvement. JCI Assumptions 2 (and 3) do not lead to any further
improvements for discovering the causal relations between system variables in this setting,though. Exploiting knowledge of the intervention targets, on the other hand, turns out
to be very helpful for getting highly accurate predictions for ancestral relations between
system variables, and also significantly improves the precision of predicting direct causal
relations and confounders. We have shown here only the results for the acyclic setting withcausal mechanism changes because we obtained qualitatively similar results for the othersimulation settings.
57

Mooij, Magliacane and Claassen

Figure 18: Results of some ASD variants (cyclic, perfect interventions) for small
models. We see a similar picture as in the acyclic case in Figure 17.

We also investigated variants of ASD-JCI0, ASD-JCI1 and ASD-JCI12 where we did not
perform any independence tests on the context variables, i.e., using conditional independence testing scheme “NC” rather than “A”. This is possible because ASD is capable ofhandling incomplete inputs. We obtained almost identical results (in all simulation settings
considered) to the standard variants of those methods in which we do perform conditional
independence tests on the context variables themselves (not shown here).
5.4.4. Multiple Context Variables vs. Single (Merged) Context Variable
Figure 20 shows that for ASD-JCI variants, exploiting multiple context variables leads tobetter results than using a single (merged) context variable, as expected. Similar results
hold for the cyclic settings and for causal mechanism changes (not shown).
5.4.5. FCI Variants
Figure 21 shows the results for the various FCI variants. FCI is seen to be somewhat lessaccurate than ASD, but bootstrapping helps to boost precision for lower recalls. Similarly
58

Joint Causal Inference from Multiple Contexts

Figure 19: Influence of JCI assumptions for ASD-JCI (acyclic, causal mechanism
changes) for small models. Exploiting more prior knowledge is seen to lead tobetter results. For reference, the non-JCI baseline ASD-obs is shown which onlyuses observational data.

to ASD, we conclude that the JCI variants of FCI substantially outperform the non-JCI
variants. FCI-JCI1 and FCI-JCI123 seem to yield identical results in this setting.
Results for causal mechanism changes are very similar to those for perfect interventions,and therefore are not shown here. The results for the cyclic setting are also not shown,because FCI was designed for the acyclic setting.
5.4.6. LCD and ICP
Figure 22 shows the results of the LCD and ICP variants for the task of predicting ancestral
relations. We only show the results for perfect interventions, as the results for causal mechanism changes are similar. LCD and ICP both only can predict the presence of ancestral
relations, not their absence. LCD and ICP apparently benefit from merging the context
variables into a single one. A possible explanation of this phenomenon could be that a
combination of conditional independence tests of the form Ck ⊥⊥ Xi0 | Xi with each Ck bi59

Mooij, Magliacane and Claassen

Figure 20: Multiple context variables vs. single (merged) context variable for
ASD-JCI (acyclic, perfect interventions) for small models. Exploiting
multiple context variables leads to better results than using a single (merged)
context variable.

nary (and both Xi and Xi0 real-valued) might be less reliable than a single test C ⊥⊥ Xi0 | Xi
where C is categorical with  2 states. Another observation we made is that bootstrapping
does lead to only marginal improvements for these methods (not shown).
5.4.7. Varying the Number of Context Variables
As we have seen, discovery of causal relations between system variables can benefit strongly
from observing the system in multiple contexts. As Figure 23 shows, the more context variables are taken into account, the better the predictions for ASD-JCI123 become. Although
not shown here, the same conclusion holds as well for the other JCI variants ASD-JCI0,
ASD-JCI1 and ASD-JCI123-kt. It does not hold for ASD baselines in general, but it does
for ASD-pikt if interventions are perfect. For the JCI variants of FCI (FCI-JCI0, FCI-JCI1
and FCI-JCI123) we also observed that the more context variables are available, the moreaccurate the predictions become. This is in line with our expectation that jointly analyz60

Joint Causal Inference from Multiple Contexts

Figure 21: Results of FCI variants (acyclic, perfect interventions) for small models.

Figure 22: Multiple context variables vs. single (merged) one for LCD and ICP(perfect interventions) for small models. LCD and ICP both benefit frommerging the context variables into a single one for the task of predicting thepresence of ancestral relations. Top: acyclic; bottom: cyclic.

61

Mooij, Magliacane and Claassen

Figure 23: Results for different numbers of context variables for ASD-JCI123(acyclic, causal mechanism changes) for small models. Taking into accountmore context variables leads to better predictions.

ing data from multiple experiments makes it easier to estimate the causal structure of thesystem.
Interestingly, the same conclusion does not hold for ASD-JCI1-sc and ASD-JCI123-sc.
This suggests that having multiple contexts is mostly beneficial if each context variabletargets only a small subset of system variables, and only for methods that can explicitly
take into account multiple context variables. For LCD and ICP, precision also does not
improve monotonically with the number of context variables. Although LCD-mc and ICP-mc
in principle allow for multiple context variables, they suffer from a drop in recall because
they focus on detecting a certain causal pattern that becomes increasingly rare with more
context variables. Indeed, for the extreme case q = p, each system variable is targeted by asingle context variable in our simulation setting, and hence one would expect LCD-mc and
ICP-mc to make no predictions at all. Any predictions they make must therefore be false
positives, resulting in low precision.
62

Joint Causal Inference from Multiple Contexts

Figure 24: Discovering indirect intervention targets (acyclic, causal mechanism
changes) on small models. Top: ASD variants; Bottom: FCI variants.

Figure 25: Discovering indirect intervention targets (cyclic, causal mechanism
changes) on small models. Top: ASD variants; Bottom: FCI variants.

63

Mooij, Magliacane and Claassen

5.4.8. Discovering Indirect Intervention Targets
Figure 24 shows for the acyclic setting with causal mechanism changes how accurately
indirect intervention targets (i.e., which system variables are descendants of each context
variable) can be discovered by various methods. Baselines ASD-obs, ASD-pooled, ASD-meta
and ASD-pikt cannot learn intervention targets (neither direct ones nor indirect ones), since
they do not represent context variables explicitly, and are therefore excluded.27 LCD and
ICP also cannot address this task.28 Although ASD-JCI123-kt makes use of known directintervention targets (i.e., which system variables are children of each context variable?),this means that there is still a non-trivial task of learning the indirect ones.
The task of deciding that a system variable is targeted is an easier one than deciding thata system variable is not targeted by an intervention. Although Fisher’s test is generally
hard to beat when it comes to predicting indirect intervention targets, ASD-JCI123-kt
outperforms it in this setting by exploiting the knowledge about direct intervention targets.
While JCI Assumption 2 turned out to be unimportant for learning the causal relations
between system variables, it is seen to be very useful for this task of learning causal relations
between context and system variables.
Figure 25 shows a largely similar picture for the cyclic setting with causal mechanism
changes. Surprisingly, FCI variants are also performing quite well in the cyclic setting. We
do not show the results for perfect interventions here, as we observed that this task is easier,and the results are generally better, but otherwise mostly similar conclusions are obtained.
The only exception is that for perfect interventions, the best method is Fisher’s approach.
5.4.9. Discovering Direct Intervention Targets
Fisher’s test is not able to predict direct intervention targets (i.e., which system variables
are children of each intervention variable?), but the ASD-JCI variants can, as well as
FCI-JCI123. Figure 26 shows the results for these algorithms in the four different simulation settings. The task is considerably easier in the acyclic setting than in the cyclic
setting. Having perfect interventions makes it slightly easier than with causal mechanism
changes. We again notice that exploiting JCI Assumption 2 considerably improves performance on this task. Surprisingly, FCI-JCI123 obtains almost perfect precision in all
scenarios, outperforming ASD-JCI123 notably in the cyclic cases. We do not understand
why this is the case.
5.4.10. Computation Time
So far we have only considered the accuracy of the predictions. Another interesting aspect isthe computation time that various methods need. Figure 27 shows total computation time(for all prediction tasks together) for all methods considered thus far. Note the logarithmic
scale on the x-axis. We only show runtimes for causal mechanism changes since those for
27. However, it would be trivial to extend ASD-pikt such that it can predict indirect intervention targets, by
combining the known direct intervention targets of a certain intervention variable with the descendantsof those assumed targets as predicted by the method as a postprocessing step.
28. However, when assuming also JCI Assumption 2, both LCD and ICP could be used to learn indirectintervention targets.

64

Joint Causal Inference from Multiple Contexts

Figure 26: Discovering direct intervention targets on small models. From top tobottom: acyclic, causal mechanism changes; acyclic, perfect interventions; cyclic,
causal mechanism changes; cyclic, perfect interventions.

perfect interventions are nearly identical. On the other hand, we do see that the cyclic
setting is more computationally demanding in general than the acyclic one.
Already for small models of p+q = 4+2 = 6 variables, the ASD algorithms become slow
because they are performing an optimization over a large discrete space. The availability ofmore background knowledge makes the search space considerably smaller, and hence leads
to reduced computation time for the ASD variants. Also, the search space is considerablylarger in the cyclic setting than in the acyclic one. By design, FCI variants are much faster,65

Mooij, Magliacane and Claassen

Figure 27: Runtimes for various methods on small models. Shown are runtimes forcausal mechanism changes; for perfect interventions, runtimes are similar. Left:acyclic; Right: cyclic.

Figure 28: Runtimes for three different algorithms as a function of the numberof context variables on small models. Note the considerably different rangesof the (logarithmic) x-axis. Results are omitted if the computation took too
long to finish.

but bootstrapping also takes its toll. The fastest methods are Fisher’s test, LCD and ICP.
Figure 28 shows how computation time scales with the number of context variables, forthree JCI implementations (ASD-JCI123, ASD-JCI1 and FCI-JCI123).
5.5. Results: Larger Simulated Models
We now present results for larger models, with p = 10 system variables and q = 10 context
variables (the meaning of the simulation parameters is explained in Section 5.2). We only
consider causal mechanism changes here, but we do distinguish the acyclic and cyclic settings. We again used 500 samples per context. For the acyclic setting, we used  = η = 0.25,
while for the cyclic setting we used  = η = 0.15 to get more or less similarly dense graphs
66

Joint Causal Inference from Multiple Contexts

Figure 29: FCI results for discovering ancestral causal relations between systemvariables in larger models. Top: acyclic; bottom: cyclic.

in both scenarios. The motivation for these parameter choices is that they are somewhat
comparable to the setting of the real-world data set that we will study in Section 5.8.
For these larger models, computation time for bootstrapped FCI methods became prohibitive for the default conditional independence test (described in Section 5.1.7). Although
we can speed up the implementation of this test that we were using considerably by implementing it more efficiently, we here simply replaced it by a standard partial correlation test.
This led to a speedup of about one order of magnitude at no apparent loss of accuracy.
5.5.1. FCI Variants
Figure 29 shows the accuracy for the task of predicting ancestral causal relations between
system variables for various FCI-JCI variants and for various FCI baselines, in both acyclic
and cyclic settings. The conclusions are in line with what we already observed for smaller
models. Again, bootstrapping FCI helps considerably to boost the accuracy of its predictions. As before, FCI-obs (which uses only observational data) performs worst. The
two baselines FCI-pooled and FCI-meta (that make use of all data) lead to a moderate
improvement. The JCI variants (FCI-JCI0, FCI-JCI1 and FCI-JCI123) perform the best,delivering almost maximum precision for a considerable recall range on the task of predicting the presence of an ancestral relation. JCI Assumption 1 does not seem to help much,as the results for FCI-JCI0 and FCI-JCI1 seem to be identical. Assuming in addition JCI
Assumption 2 (and 3) does help to obtain slightly higher precision. The good performanceof FCI variants in the cyclic setting is surprising, since FCI was designed for the acyclic
case.
67

Mooij, Magliacane and Claassen

Figure 30: Bootstrapped LCD and ICP results for discovering ancestral causal
relations between system variables for larger models. Top: acyclic; bottom:
cyclic.

5.5.2. LCD and ICP
In Figure 30, we show the accuracy of bootstrapped LCD and ICP for the task of predicting ancestral causal relations between system variables, in both the acyclic and the cyclic
setting. As for FCI, bootstrapping improves the accuracy of LCD and ICP results, and we
decided to only show the bootstrapped results here. Contrary to what we observed for small
models, the “multiple context” (“-mc”) versions of both algorithms now clearly outperformthe versions that use only a single (merged) context (“-sc”) in these settings. Interestingly,
the accuracy of LCD is quite similar to that of ICP. The additional complexity of ICP
apparently does not lead to substantially better results than the LCD algorithm alreadyoffers in these settings. Also, the precision of LCD-mc is comparable to that of FCI-JCI123,
the most accurate of the JCI variants of FCI (cf. Figure 29), on the task of predicting thepresence of ancestral relations.
5.5.3. Discovering Intervention Targets
Figure 31 shows the performance of FCI-JCI variants (and as a baseline, Fisher’s test) on
the task of discovering indirect intervention targets, for both the acyclic and cyclic setting.
Interestingly, JCI Assumption 2 seems necessary to obtain good results on this task. Still,
FCI-JCI123-bs is outperformed by Fisher’s test.
On the other hand, Fisher’s test cannot identify direct intervention targets, whereas
FCI-JCI123 can. Figure 32 shows that FCI-JCI123 can identify the direct intervention
targets as well as the non-targets with high precision. Surprisingly, this works also in thecyclic case.
68

Joint Causal Inference from Multiple Contexts

Figure 31: FCI results for discovering indirect intervention targets in larger models.
Top: acyclic; bottom: cyclic.

Figure 32: FCI results for discovering direct intervention targets in larger models.
Top: acyclic; bottom: cyclic.

69

Mooij, Magliacane and Claassen

Figure 33: Runtimes for various methods on larger models. Left: acyclic; Right: cyclic.

5.5.4. Computation Time
Figure 33 shows total runtimes of the methods that we ran on the larger simulated models.
First, we observe no big differences between the runtimes for the cyclic setting with respect
to the acyclic one. However, we do observe huge differences in runtime between various
methods. LCD variants and Fisher’s test are by far the fastest. ICP variants come second.
Bootstrapping puts a large toll on computation time for FCI variants. JCI variants of
FCI are much slower than non-JCI variants. This seems to be mostly due to an exponential
increase in the number of conditional independence tests. Indeed, we observed that FCI-JCI
variants are conditioning on a substantial fraction of all 210 subsets of all context variablesin the skeleton search phase. Nevertheless, JCI variants of FCI are still computationallyfeasible in this setting, even with bootstrapping.
5.6. Results: Large Simulated Models
We now present results for large simulated models, with p = 100 system variables and q = 10
context variables. We only consider causal mechanism changes with unknown targets andonly the acyclic setting. We used  = η = 0.02, which yields rather sparse graphs, in order toavoid that the computations would take too long (the meaning of the simulation parameters
is explained in Section 5.2). We used only 100 samples per context, because the tasks would
become too easy otherwise due to the sparsity of the graphs. We again used the standard
partial correlation test for FCI variants in this setting instead of the default conditional
independence test described in Section 5.1.7 for computational efficiency reasons.29
29. It is still possible to use the standard test, and the results are slightly better, but the small gain in
precision does not seem to justify the large increase in computation time. Implementing FCI-JCI123r
as proposed in Section 4.2.5 would probably yield a significant reduction in computation time. Also,
alternatives for the skeleton search phase such as FCI+ (Claassen et al., 2013) could be employed togain further speedups. Last but not least, a more efficient implementation of the standard conditional
independence test would help considerably.

70

Joint Causal Inference from Multiple Contexts

Figure 34: Discovering ancestral causal relations between system variables in large
models. Top: FCI variants; Bottom: Bootstrapped LCD and ICP variants. Notethat we zoomed in on the PR curves.

5.6.1. Ancestral Causal Relations between System Variables
Figure 34 shows the accuracy for the task of discovering ancestral causal relations between
system variables for various methods and baselines. Overall, we see that JCI variants
outperform non-JCI baselines. On a detailed level, the conclusions are somewhat different
than what we saw for smaller models.
We start by discussing the results for the task of predicting the presence of causal
relations. Like before, we see that bootstrapping FCI helps considerably to increase theprecision of the predictions for the lower recall range. On the other hand, it reduces recall,possibly because only half of the available data is used and the independence threshold
was not adjusted. As before, FCI-obs (which uses only observational data) performs worst.
However, FCI-obs outperforms random guessing by a large margin in this setting. This isespecially noteworthy given that it is only using 100 observational samples. Interestingly,
the main improvement in this setting is obtained by pooling the data; whether one includes
the context variables (FCI-JCI0) or not (FCI-pooled) does not seem to make much of adifference. Also, using more JCI background knowledge yields only small improvements:
the differences between FCI-JCI0, FCI-JCI1 and FCI-JCI123 are small.
We observe that LCD-mc-bs obtains a higher precision at low recall than the FCI-JCI
variants. On the other hand, the FCI-JCI variants maintain a decent precision over a larger
recall range, contrary to the LCD and ICP variants that only look for very specific patterns,
which may explain that their precision drops off at lower recall than it does for FCI-JCI.
The “multiple context” (“-mc”) versions of LCD and ICP outperform the versions that use
71

Mooij, Magliacane and Claassen

Figure 35: Results for discovering intervention targets in large models. Top: indirectintervention targets; Bottom: direct intervention targets. Note that we zoomed
in on the PR curves.

only a single (merged) context (“-sc”) in these settings. Interestingly, both variants of
LCD outperform the corresponding variant of the more complicated ICP algorithm.
For predicting the absence of causal relations, the random guessing baseline alreadyobtains a high precision, because of the sparsity of the graphs. FCI variants improve on
this, roughly halving the error for the most confident predictions when using the bootstrapped versions. FCI-JCI variants again obtain the highest precision on this task, but
don’t significantly outperform FCI-pooled.
5.6.2. Discovering Intervention Targets
Figure 35 shows the performance of FCI-JCI variants (and as a baseline, Fisher’s test) on
the task of discovering intervention targets. For discovering indirect intervention targets,
Fisher’s test is now slightly outperformed by FCI-JCI123-bs. Interestingly, JCI Assumption 2 seems necessary to obtain any results on that particular task. Investigating the
PAGs shows that the edges between context and system variables are mostly bidirected for
FCI-JCI1 and FCI-JCI0, which explains why these two algorithms yield no predictions atall. We do not have a good explanation for this behavior, but speculate that the sparsesetting with many nodes makes certain empirical violations of faithfulness quite likely. Oneof the features of FCI-JCI123 is that it can discover direct intervention targets, somethingthat Fisher’s test cannot. We find that FCI-JCI123 identifies with high precision directintervention targets as well as non-targets in this simulation setting.
72

Joint Causal Inference from Multiple Contexts

Figure 36: Runtimes for various methods on large models.
5.6.3. Computation Time
Figure 36 shows total runtimes of the methods that we ran on the large simulated models. For most methods, the largest part of the total running time is spent on performing
independence tests. In particular, the tests that subdivide data according to context are
relatively slow since we have not seriously optimized their implementation.
5.7. Summary of Results on Simulated Data
We have seen in our experiments with simulated data that JCI methods typically outperform
non-JCI methods, in some settings by a large margin. For certain tasks, our newly proposed
FCI-JCI algorithms provide the new state-of-the-art. Interestingly, Fisher’s baseline turned
out to be hard to beat on the task of discovering indirect intervention targets. However,there are other tasks for which Fisher’s baseline cannot be applied but for which our newly
proposed methods do apply, such as the task of discovering direct intervention targets.
As expected, LCD, ICP and ASD variants work in both the acyclic and cyclic setting.
While FCI variants were expected to work only in the acyclic setting, we were surprisedby how well they perform in the cyclic setting.30 Often, but not always, adding more
context variables leads to better results. Having multiple contexts was seen to be mostly
beneficial if each context variable targets only a small subset of system variables, andthen only for methods that can explicitly take into account multiple context variables. Aninteresting exception to this are LCD and ICP, which due to their sensitivity to very specificpatterns actually degrade in performance when too many system variables become directly
targeted by context variables. Exploiting more JCI background knowledge typically led tobetter results, but it depends on the task and simulation setting how large the benefits are.
Interestingly, the largest boost in accuracy for discovery of causal relations between systemvariables comes already from JCI Assumption 0 (i.e., from pooling the data and adding thecontext variables).
30. These empirical observations led us to conjecture that FCI does not need to be adapted for the σseparation setting. It was shown very recently that this is indeed the case: FCI is also sound andcomplete in that setting (Mooij and Claassen, 2020).

73

Mooij, Magliacane and Claassen

As to the relative merits of the various JCI methods that we compared, it is difficult
to state this concisely. Different methods behave differently on different tasks in different
simulation settings, both in terms of precision and recall as well as in terms of computation
time, and for many methods there is a combination of a task and a simulation setting inwhich they do relatively well. Generally speaking, LCD and ICP behave quite similarly, arerelatively fast and obtain high precision but have low recall; ASD-JCI variants are among
the most accurate and have highest recall, but computation time explodes for more thana handful of variables; performance of FCI-JCI methods turns out to be somewhere in
between, both in terms of accuracy and in terms of computation time. When implemented
properly, their scalability is comparable to that of standard FCI, where the total numberof variables |I| + |K| and the sparsity of the underlying MAGs are important factors in thecomputation time.
Some aspects that should be kept in mind when interpreting these results is that all
simulations have been done under JCI Assumptions 1, 2, 3, and there was no model misspecification. From that perspective, it was to be expected that JCI methods would work
well. Also, we make no claims as to how our conclusions would generalize to different settings, for example, with non-Gaussian noise distributions, discrete variables, or continuous
variables with non-linear interactions. For those settings, the choice of the conditional independence tests could have a large influence on the results, for example. A detailed studyof that is beyond the scope of this paper.
5.8. Results: Real-world (Flow Cytometry) Data
In this subsection, we present an application of the Joint Causal Inference framework on
real-world data: the flow cytometry data of Sachs et al. (2005). The data consists of acollection of data sets, where each data set corresponds with a different experimental condition in which a system was perturbed and subsequently measured. The system consistsof an individual cell, drawn randomly from a collection of primary human immune systemcells. The system variables measure the abundances of several phosphorylated protein andphospholipid components in an individual cell using a measurement technology known as
flow cytometry. Performing the measurement destroys the cell, and hence, it is not possible
to obtain multiple measurements over time from the same cell. Instead, snap-shot measurements of thousands of individual cells are available, obtained in different experimental
conditions in which the cells were perturbed with molecular interventions, performed byadministering certain reagents to the cells. Most of these interventions are not perfect,
but rather change the activity of some component by an unknown amount. There is priorknowledge about the targets of these interventions (see Table 3), but it is not clear whether
interventions are as specific as claimed. Many existing causal discovery approaches assume
that the true causal graph is acyclic and that the system variables are causally sufficient.
However, it is known that these cellular signaling networks contain strong feedback loops,
and it is quite likely that some of the variables may be subject to latent confounding. Thus,
this type of experimental data constitutes a compelling motivation for the Joint Causal
Inference framework.
Over the years, this particular flow cytometry data set has become a “benchmark” incausal discovery (see e.g. Ramsey and Andrews (2018) for some references). Many causal
74

Joint Causal Inference from Multiple Contexts

Cγ
0
1
0
0
0
0
0
0

Cδ
0
0
1
0
0
0
0
0

C
0
0
0
1
0
0
0
0

Cζ
0
0
0
0
1
0
0
0

Cη
0
0
0
0
0
1
0
0

(Cα , Cθ , Cι )(1,0,0)(1,0,0)(1,0,0)(1,0,0)(1,0,0)(1,0,0)(0,1,0)(0,0,1)

NC
853
911
723
810
799
848
913
707

Reagents added
α-CD3, α-CD28
α-CD3, α-CD28, AKT inhibitor
α-CD3, α-CD28, G0076
α-CD3, α-CD28, Psitectorigenin
α-CD3, α-CD28, U0126
α-CD3, α-CD28, LY294002
PMA
β2CAMP

Table 6: Experimental design for part of the Sachs et al. (2005) flow cytometry data usedin our experiments.

discovery methods have been applied to this data, and in many cases, the “consensus network” in Sachs et al. (2005), visualized in Figure 38(a), was used as a ground truth toevaluate the results of the causal discovery procedure. However, we would like to pointout that there are good reasons to be skeptical about the assumption that the “consensus
network” represents the true causal graph of the system (as acknowledged by several domain experts we spoke). Indeed, by inspecting the data one can find many examples where
the data is incompatible with the hypothesis that the “consensus network” is a realistic
and complete description of the underlying system. For example, according to the “consensus network”, an intervention that inhibits the activity of Mek should have no effecton Raf (because Raf is a direct cause of Mek and there is no feedback loop from Mek toRaf). However, in the data we see an increase of more than one order of magnitude in theabundance of Raf when U0126 (a reagent assumed to inhibit Mek activity) is added to thecells (see Figure 37). So either U0126 also directly targets Raf, or Mek must be a cause ofRaf, in both cases contradicting the “consensus network”. In the literature regarding thissignaling pathway (which has been studied in great detail since it plays an important role inmany human cancers) it is often suggested that there should be a feedback loop back from
Erk to Raf (whose molecular mechanism is still unknown). This would be in line with ourobservations from the data in Figure 37, and also imply that the “consensus network” is incomplete. This is just one example illustrating that the data is not entirely compatible withthe “consensus network”, and it is easy to find more of these examples via visual inspectionof the data. Therefore, we will not make use of the “consensus network” as a ground truth
to compare with when evaluating the output of the causal discovery algorithms.
Sachs et al. (2005) use an MCMC method to estimate the structure of a causal Bayesian
network from the combined observational and interventional data, making use of the modified BDe score proposed by Cooper and Yoo (1999). Eaton and Murphy (2007) later useda dynamic programming algorithm to solve the estimation problem exactly. Like the original analysis by Sachs et al. (2005), many causal discovery methods that have been appliedon this data rely on the background knowledge about the intervention types and targets,for which we provide (our interpretation) in Table 3. A notable exception is Eaton andMurphy (2007), who were the first to estimate the intervention targets directly from thedata. Exploiting the background knowledge on intervention targets and types simplifies thecausal discovery problem considerably. However, the accuracy of this background knowl75

10

10

8

8

6

6

ln Erkln MekMooij, Magliacane and Claassen

4
2

4
2

0

0
0

2

4

6

8

10

0

ln Raf

2

4

6

8

10

ln MekFigure 37: Log-abundances of Mek vs. Raf (left) and of Erk vs. Mek (right). Blue: observational baseline (Cα = 1, C\α = 0); Red: reagent U0126 added (Cζ = 1). We
observe: (i) the measurement noise is quite small; (ii) Raf and Mek are highly
correlated (“consensus network”: Raf is a direct cause of Mek); (iii) strong evidence for feedback (intervening on Mek changes Raf abundance) if we assume
that U0126 directly targets Mek but not Raf; (iv) the Mek inhibitor U0126 increases Mek abundance (so modeling this as a perfect intervention would not
be realistic); (v) Mek and Erk are independent in both contexts (even though
Mek is a direct cause of Erk according to the “consensus network”), an apparent
violation of faithfulness.

edge is not universally accepted. In particular, many biologists that we spoke with were
skeptical about the assumed specificity of the interventions (i.e., the interventions may haveadditional direct effects that are not listed in the table).
We ran various FCI-JCI variants on a subset of the flow cytometry data.31 The experimental design of the original data is described in Table 2 (left), p. 34. In order toavoid deterministic relations between context variables, we merged context variables Cα(α-CD3/CD28), Cθ (PMA) and Cι (β2CAMP), as discussed in Section 4.1, leading to theexperimental design in Table 2 (right). This means that we must interpret the merged context variable as referring to the addition of PMA or β2CAMP, combined with the omissionof α-CD3/CD28. However, note that there are still approximate conditional independencesin this experimental design of the form Cβ ⊥⊥ Ck | (Cα , Cθ , Cι ) for k ∈ {γ, δ, , ζ, η}. This
could lead to problems with JCI Assumption 3. For that reason, but also in order to enable
comparisons with other results reported in the literature, we only used the 8 (out of 14)
experimental conditions in the data set in which no ICAM.2 had been administered (i.e.,with Cβ = 0), and ignored the others.
Similarly to Eaton and Murphy (2007), we do not use the background knowledge regarding intervention types or targets. We only assume that the experimental setting is capturedby the JCI framework. JCI Assumption 1 should be true because the intervention is performed some time (approximately 20 minutes) before the measurements are done. We havealready discussed the validity of JCI Assumption 2 for this particular experimental setting
31. As preprocessing, we simply took the logarithm of the raw values.

76

Joint Causal Inference from Multiple Contexts

in Section 3.4. Assuming that the context variables provide a complete causal descriptionof the context (in particular, that there are no unintended batch effects), JCI Assumption 2
applies. JCI Assumption 3 then also applies since there are no conditional independencesin the context distribution (after merging Cα , Cθ and Cι and leaving out all contexts withCβ = 1). When using FCI-JCI123, we can then learn the intervention targets from the data
itself, without making use of the background knowledge on intervention types and targets.
For comparison, we also ran FCI-obs, i.e., standard FCI using only the observational
data set (i.e., the one in which only global activators α-CD3 and α-CD28 have been administered). We also ran FCI-meta, which uses Fisher’s method to combine p-values ofconditional independence tests in the 8 separate experimental conditions, which are thenused as input for standard FCI. Finally, we ran FCI-pooled, i.e., standard FCI on the 8experimental conditions pooled together (but excluding context variables). In all those FCI
variants, we assumed that no selection bias would be present. We additionally compare withother JCI implementations, in particular, multiple variants of LCD and ICP. Computation
times of various implementations are reported in Figure 42.
The “consensus network” and the PAGs obtained by the FCI baselines are shown in
Figure 38. Figure 39 shows PAGs obtained by the FCI-JCI variants. Note that we show the
PAGs obtained without bootstrapping, although these are not necessarily stable. Therefore,
we show in Figure 40 also the bootstrapped results for the learned ancestral causal relations
between system variables, and the learned intervention targets. One can see here that most,but not all, of these features are stably predicted by the FCI-JCI variants. In particular,
the causal relations32 Mek99KRaf, PLCg99KPIP2, Akt99KErk, P3899KPKC are predictedby both FCI-JCI123 and FCI-JCI1 with high confidence.
Regarding the learned indirect intervention targets, FCI-JCI123 has an advantage over
FCI-JCI1 because it can exclude bidirected edges between context and system variables.
Nonetheless, FCI-JCI1 predicts (in accordance with FCI-JCI123) that G0076 affects PLCg,
PIP2, PKC and P38. For discovering indirect intervention targets, Fisher’s test for causality is simple and powerful. However, it is not able to learn the direct intervention targets,like FCI-JCI123 can. Although the direct intervention targets learned by FCI-JCI123 do
not all correspond with the “consensus network”, they do agree for example on Psitectorigenin directly targeting PIP2. Most direct intervention targets learned by FCI-JCI123 werealso found by Eaton and Murphy (2007). It is interesting to see that FCI-JCI123 considers
both Mek and Erk as direct intervention targets of U0126, while Raf is identified as anindirect target.
In the absence of a reliable ground truth, we draw the following conclusions from these results. First, the reference methods that exploit knowledge on intervention types and targets
obtain rather consistent results. Second, the JCI methods that do not assume knowledge
on intervention types and targets (and the approach by Eaton and Murphy (2007)) showless consistent results. This could indicate that the data contain not enough signal in order to solve this more ambitious task reliably. In that case, some model misspecification(for example, strong non-linearities or deviations from Gaussianity, which makes a simple
partial correlation based conditional independence test inadequate) could lead to inconsistencies between methods. Nonetheless, the performance of the FCI-JCI variants appears
32. We write i 99K j if i is a cause of j.

77

Mooij, Magliacane and Claassen

LY294002

AKT.inh

PIP3

PLCg

Psitectorigenin

PLCg

PIP2

Raf

beta2CAMP + noAlphaCD3/28

PKA

PKC

G0076

Mek

Akt

Raf

JNK

Erk

PIP3

p38
PKC

U0126

PIP2

PMA + noAlphaCD3/28

Akt

PKA

Mek
p38

Erk

JNK

(a) “Consensus Network”
Raf

Mek

PLCg

PIP2

Erk

PIP3

Akt

(b) FCI-pooled
Raf

PKC

p38

Mek

JNK

PLCg

PIP2

Akt

PIP3

PKA

(c) FCI-obs

Erk

PKA

PKC

p38

JNK

(d) FCI-meta

Figure 38: PAGs resulting from various FCI baselines on the flow cytometry dataof Sachs et al. (2005). Intervention variables are denoted with rectangles,
system variables with ellipses. From top to bottom: (a) “Consensus network”
according to Sachs et al. (2005); (b) FCI on pooled data (without adding thecontext variables); (c) FCI on the first (“observational”) data set in which only
global activators α-CD3 and α-CD28 have been administered; (d) FCI with as
input the result of Fisher’s method for combining conditional independence test
results from all data sets.

to be a considerable improvement over the simple FCI baselines FCI-obs, FCI-meta and
FCI-pooled. Remarkably, FCI-JCI123 manages to orient most of the edges. The output
also resembles the consensus network, although some of the edges seem to be reversed.
Considering that we have not taken into account the available background knowledge on
the intervention types and targets (Table 3), and that we have not used any tuning of theparameters we consider this still an impressive and encouraging result that illustrates thepotential that JCI has for analyzing complex scientific data sets.
78

Joint Causal Inference from Multiple Contexts

Raf

PMA/beta2CAMP + noAlphaCD3/28

Mek

U0126

Erk

Akt

AKT.inh

LY294002

G0076

PLCg

Psitectorigenin

PKC

PKA

p38

PIP2

JNK

PIP3

(a) FCI-JCI123
PLCg

PIP2

PIP3

Raf

Psitectorigenin

PKC

Mek

AKT.inh

Erk

U0126

Akt

LY294002

p38

PKA

JNK

G0076

PMA/beta2CAMP + noAlphaCD3/28(b) FCI-JCI1
PLCg

PIP2

PIP3

Psitectorigenin

AKT.inh

Raf

PKC

Mek

Erk

U0126

Akt

LY294002

p38

PKA

JNK

G0076

PMA/beta2CAMP + noAlphaCD3/28(c) FCI-JCI0

Figure 39: PAGs resulting from various FCI-JCI variants on the flow cytometry
data of Sachs et al. (2005). These causal discovery methods do not make useof the biological prior knowledge regarding intervention types and targets, butlearn the intervention targets from the data. Intervention variables are denotedwith rectangles, system variables with ellipses. From top to bottom, less JCI
Assumptions are made. Note that these are individual PAGs that have notbeen bootstrapped. To get an idea of the robustness, Figure 40 shows also thecorresponding bootstrap estimates for certain features of the PAGs.

79

Mooij, Magliacane and Claassen

Direct causal relations

1.00

0.75

Raf->Mek
Raf->PLCg
Raf->PIP2
Raf->PIP3
Raf->Erk
Raf->Akt
Raf->PKA
Raf->PKC
Raf->p38
Raf->JNK
Mek->Raf
Mek->PLCg
Mek->PIP2
Mek->PIP3
Mek->Erk
Mek->Akt
Mek->PKA
Mek->PKC
Mek->p38
Mek->JNK
PLCg->Raf
PLCg->Mek
PLCg->PIP2
PLCg->PIP3
PLCg->Erk
PLCg->Akt
PLCg->PKA
PLCg->PKC
PLCg->p38
PLCg->JNK
PIP2->Raf
PIP2->Mek
PIP2->PLCg
PIP2->PIP3
PIP2->Erk
PIP2->Akt
PIP2->PKA
PIP2->PKC
PIP2->p38
PIP2->JNK
PIP3->Raf
PIP3->Mek
PIP3->PLCg
PIP3->PIP2
PIP3->Erk
PIP3->Akt
PIP3->PKA
PIP3->PKC
PIP3->p38
PIP3->JNK
Erk->Raf
Erk->Mek
Erk->PLCg
Erk->PIP2
Erk->PIP3
Erk->Akt
Erk->PKA
Erk->PKC
Erk->p38
Erk->JNK
Akt->Raf
Akt->Mek
Akt->PLCg
Akt->PIP2
Akt->PIP3
Akt->Erk
Akt->PKA
Akt->PKC
Akt->p38
Akt->JNK
PKA->Raf
PKA->Mek
PKA->PLCg
PKA->PIP2
PKA->PIP3
PKA->Erk
PKA->Akt
PKA->PKC
PKA->p38
PKA->JNK
PKC->Raf
PKC->Mek
PKC->PLCg
PKC->PIP2
PKC->PIP3
PKC->Erk
PKC->Akt
PKC->PKA
PKC->p38
PKC->JNK
p38->Raf
p38->Mek
p38->PLCg
p38->PIP2
p38->PIP3
p38->Erk
p38->Aktp38->PKA
p38->PKC
p38->JNKJNK->Raf
JNK->Mek
JNK->PLCg
JNK->PIP2
JNK->PIP3
JNK->Erk
JNK->Akt
JNK->PKA
JNK->PKC
JNK->p38

Ancestral relationsCon
Es sens
tim usate (Sd ac
Ea (Sa hs et c
A
hid IC CI Moon & hs et al.,de P (M oij Mu t al. 20
nIC (M ag & rp , 2 05e
P ( in liac He hy 00 )Me sh an ske (2 5)ins aus e e s ( 007ha en t a 20 )us et l., 2 13
en al. 0 )et , 2 16)
al. 01, 6
FC FCI-JC2016)
I-JC I1 )I12 23
3
FC FCI-JC-bs
I-JC I1
I1
FC FCI-JC-bs
I-JC I0
I0
FC FCI-o-bs
I-o b
FC FCI-p bs-bs
I-p oo s
oo lele d
FC FCI-md-bs
I-m e
eta taLC LCD -bs
D-m -mc- c
ICPICP-mbs-m ccLC LCD bs
D-s -scc-b
ICPICP-ss
-sc c
-bs

Co
n
Es sens
tim usate (Sd ac
hid IC Eat (Sachs ede P M on hs t al.
nIC (M o & et , 2
P ( ein oij Mu al. 00
Me sh & H rph , 2 5)
ins aus es y ( 005
ha en ke 20 )us et s (2 07
en al. 0 )et , 2 13)
al. 01, 2 6)
01
6)

Raf->Mek
Raf->PLCg
Raf->PIP2
Raf->PIP3
Raf->Erk
Raf->Akt
Raf->PKA
Raf->PKC
Raf->p38
Raf->JNK
Mek->Raf
Mek->PLCg
Mek->PIP2
Mek->PIP3
Mek->Erk
Mek->Akt
Mek->PKA
Mek->PKC
Mek->p38
Mek->JNK
PLCg->Raf
PLCg->Mek
PLCg->PIP2
PLCg->PIP3
PLCg->Erk
PLCg->Akt
PLCg->PKA
PLCg->PKC
PLCg->p38
PLCg->JNK
PIP2->Raf
PIP2->Mek
PIP2->PLCg
PIP2->PIP3
PIP2->Erk
PIP2->Akt
PIP2->PKA
PIP2->PKC
PIP2->p38
PIP2->JNK
PIP3->Raf
PIP3->Mek
PIP3->PLCg
PIP3->PIP2
PIP3->Erk
PIP3->Akt
PIP3->PKA
PIP3->PKC
PIP3->p38
PIP3->JNK
Erk->Raf
Erk->Mek
Erk->PLCg
Erk->PIP2
Erk->PIP3
Erk->Akt
Erk->PKA
Erk->PKC
Erk->p38
Erk->JNK
Akt->Raf
Akt->Mek
Akt->PLCg
Akt->PIP2
Akt->PIP3
Akt->Erk
Akt->PKA
Akt->PKC
Akt->p38
Akt->JNK
PKA->Raf
PKA->Mek
PKA->PLCg
PKA->PIP2
PKA->PIP3
PKA->Erk
PKA->Akt
PKA->PKC
PKA->p38
PKA->JNK
PKC->Raf
PKC->Mek
PKC->PLCg
PKC->PIP2
PKC->PIP3
PKC->Erk
PKC->Akt
PKC->PKA
PKC->p38
PKC->JNK
p38->Raf
p38->Mek
p38->PLCg
p38->PIP2
p38->PIP3
p38->Erk
p38->Aktp38->PKA
p38->PKC
p38->JNKJNK->Raf
JNK->Mek
JNK->PLCg
JNK->PIP2
JNK->PIP3
JNK->Erk
JNK->Akt
JNK->PKA
JNK->PKC
JNK->p38

0.50

0.25

0.00
0.25
Confidence Score

0.50

0.75

1.00

Figure 40: Causal relationships between the biochemical agents in the flow cytometry data of Sachs et al. (2005), according to different causal discovery
methods and the “consensus network” according to Sachs et al. (2005) (which we
do not consider as a reliable complete ground truth). We also included resultsof causal discovery methods reported in other works.

80

Joint Causal Inference from Multiple Contexts

Direct intervention targets

Indirect intervention targets

1.00

sen
Co
n

Co
n

sen

su
s
Ea (Saton chs& M eturp al.,hy 200(20 5)07F
FC CI-JFishe)
I-JC CI1 r
I12 23
3
FC FCI-J -bs
I-JC CI1
I
1
FC FCI-J -bs
I-JC CI0
I0bs

AKT.inh->Raf
AKT.inh->Mek
AKT.inh->PLCg
AKT.inh->PIP2
AKT.inh->PIP3
AKT.inh->Erk
AKT.inh->Akt
AKT.inh->PKA
AKT.inh->PKC
AKT.inh->p38
AKT.inh->JNK
G0076->Raf
G0076->Mek
G0076->PLCg
G0076->PIP2
G0076->PIP3
G0076->Erk
G0076->Akt
G0076->PKA
G0076->PKC
G0076->p38
G0076->JNK
Psitectorigenin->Raf
Psitectorigenin->Mek
Psitectorigenin->PLCg
Psitectorigenin->PIP2
Psitectorigenin->PIP3
Psitectorigenin->Erk
Psitectorigenin->Akt
Psitectorigenin->PKA
Psitectorigenin->PKC
Psitectorigenin->p38
Psitectorigenin->JNK
U0126->Raf
U0126->Mek
U0126->PLCg
U0126->PIP2
U0126->PIP3
U0126->Erk
U0126->Akt
U0126->PKA
U0126->PKC
U0126->p38
U0126->JNK
LY294002->Raf
LY294002->Mek
LY294002->PLCg
LY294002->PIP2
LY294002->PIP3
LY294002->Erk
LY294002->Akt
LY294002->PKA
LY294002->PKC
LY294002->p38
LY294002->JNK
PMA/beta2CAMP + noAlphaCD3/28->Raf
PMA/beta2CAMP + noAlphaCD3/28->Mek
PMA/beta2CAMP + noAlphaCD3/28->PLCg
PMA/beta2CAMP + noAlphaCD3/28->PIP2
PMA/beta2CAMP + noAlphaCD3/28->PIP3
PMA/beta2CAMP + noAlphaCD3/28->Erk
PMA/beta2CAMP + noAlphaCD3/28->Akt
PMA/beta2CAMP + noAlphaCD3/28->PKA
PMA/beta2CAMP + noAlphaCD3/28->PKC
PMA/beta2CAMP + noAlphaCD3/28->p38
PMA/beta2CAMP + noAlphaCD3/28->JNK

sus
Ea (Saton chs& M eturp al.,hy 200( 5
FC FCI-J2007)
I-JC CI1 )I12 233-bs

AKT.inh->Raf
AKT.inh->Mek
AKT.inh->PLCg
AKT.inh->PIP2
AKT.inh->PIP3
AKT.inh->Erk
AKT.inh->Akt
AKT.inh->PKA
AKT.inh->PKC
AKT.inh->p38
AKT.inh->JNK
G0076->Raf
G0076->Mek
G0076->PLCg
G0076->PIP2
G0076->PIP3
G0076->Erk
G0076->Akt
G0076->PKA
G0076->PKC
G0076->p38
G0076->JNK
Psitectorigenin->Raf
Psitectorigenin->Mek
Psitectorigenin->PLCg
Psitectorigenin->PIP2
Psitectorigenin->PIP3
Psitectorigenin->Erk
Psitectorigenin->Akt
Psitectorigenin->PKA
Psitectorigenin->PKC
Psitectorigenin->p38
Psitectorigenin->JNK
U0126->Raf
U0126->Mek
U0126->PLCg
U0126->PIP2
U0126->PIP3
U0126->Erk
U0126->Akt
U0126->PKA
U0126->PKC
U0126->p38
U0126->JNK
LY294002->Raf
LY294002->Mek
LY294002->PLCg
LY294002->PIP2
LY294002->PIP3
LY294002->Erk
LY294002->Akt
LY294002->PKA
LY294002->PKC
LY294002->p38
LY294002->JNK
PMA/beta2CAMP + noAlphaCD3/28->Raf
PMA/beta2CAMP + noAlphaCD3/28->Mek
PMA/beta2CAMP + noAlphaCD3/28->PLCg
PMA/beta2CAMP + noAlphaCD3/28->PIP2
PMA/beta2CAMP + noAlphaCD3/28->PIP3
PMA/beta2CAMP + noAlphaCD3/28->Erk
PMA/beta2CAMP + noAlphaCD3/28->Akt
PMA/beta2CAMP + noAlphaCD3/28->PKA
PMA/beta2CAMP + noAlphaCD3/28->PKC
PMA/beta2CAMP + noAlphaCD3/28->p38
PMA/beta2CAMP + noAlphaCD3/28->JNK

0.75

0.50

0.25

0.00
0.25
Confidence Score

0.50

0.75

1.00

Figure 41: Intervention effects on biochemical agents in the flow cytometry dataof Sachs et al. (2005), according to different causal discovery methods and the“consensus network” according to Sachs et al. (2005) (which we do not consider
as a reliable complete ground truth). We also included results of causal discovery
methods reported in other works.

81

Mooij, Magliacane and Claassen

Figure 42: Runtimes for various methods on the flow cytometry data of Sachs
et al. (2005).

6. Conclusions and Discussion
In this work, we proposed Joint Causal Inference (JCI), a powerful and elegant framework
for causal discovery from data sets from multiple contexts. JCI generalizes the ideas ofcausal discovery based on experimentation (as in randomized controlled trials and A/Btesting) to multiple context and system variables. Seen from another perspective, it also
generalizes the ideas of causal discovery from purely observational data to the setting of data
sets from multiple contexts—for example, different interventional regimes—by reducing thelatter to a special case of the former, with additional background knowledge on the causal
relationships involving the context variables. We proposed different flavours of JCI thatdiffer in the amount of background knowledge that is assumed, some being more conservative
than others. JCI can be implemented with any causal discovery method that can take into
account the background knowledge. Surprisingly, we saw that one can even apply an offthe-shelf causal discovery algorithm for purely observational data on the pooled data (withcontext variables included), completely ignoring the background knowledge, and thereby
already obtain significant improvements in the accuracy of the discovered causal relations.
We have seen how JCI deals with different types of interventions in a unified fashion,
how it reduces learning intervention targets to learning the causal relations between context
and system variables, and that it allows one to fully exploit all the information in the joint
distribution on system and context variables. JCI was partially inspired by the approachby Eaton and Murphy (2007), but is much more generally applicable, as it allows for latent
confounders and cycles, which are both important in many application domains. Especially
noteworthy is that more conservative flavours of JCI allow for confounders between systemand context variables, which cannot always be excluded, for example because the relevant
aspects of the system’s context were only partially observed.
We have investigated various implementations of JCI, amongst which some existingalgorithms (LCD, ICP, and standard estimators for the presence of a causal effect in arandomized controlled trial), and also proposed novel implementations that are adaptations
82

Joint Causal Inference from Multiple Contextsof algorithms for causal discovery from purely observational data to the JCI setting. In
particular, we proposed ASD-JCI, an adaptation of the method of Hyttinen et al. (2014)
combined with ideas from Magliacane et al. (2016b), which is very flexible and accurate.
By replacing d-separation with σ-separation (Forré and Mooij, 2018), ASD-JCI can also
be used in general nonlinear cyclic settings. A major disadvantage of ASD-JCI is that it
becomes computationally extremely expensive already for as few as about 7 variables. We
also proposed FCI-JCI, an adaptation of the FCI algorithm that enables it to exploit theapplicable JCI background knowledge. This algorithm is less accurate than ASD-JCI, butmuch faster.
We evaluated different implementations of the JCI approach on synthetic data. We sawthat JCI implementations outperform other state-of-the-art causal discovery algorithms inmost settings. In some cases, the gains were quite extreme; for example, while purely
observational causal discovery methods did not perform better than random guessing on
small models, JCI variants were able to discover with almost perfect precision ancestralcausal relations between system variables. The only case in which all JCI implementations
were outperformed by another causal discovery algorithm that combines data from different
contexts, was the setting in which the contexts correspond with perfect interventions withknown targets. The reason is that none of the JCI implementations exploited the perfectnature of the interventions. However, we also saw that if interventions are not perfect (forexample, in the case of causal mechanism changes), JCI implementations still perform very
well, while algorithms relying on the perfect nature of interventions may suffer from model
misspecification. Another interesting observation we made in the experiments on synthetic
data is that for the task of discovering indirect (ancestral) causal relations, the classic (andvery simple and fast) LCD algorithm can be competitive with more sophisticated algorithms,
like ICP and bootstrapped FCI-JCI.
We further illustrated the use of JCI by analyzing flow cytometry protein expression data
(Sachs et al., 2005), a famous “benchmark” in the field of causal discovery. Unfortunately,
applying ASD-JCI on the 11 system and 6 context variables would take excessive amountsof computation time, so we had to resort to FCI-JCI instead for causal discovery on a
global scale. We compared with LCD and ICP variants that do causal discovery locally.
The results of various methods differ considerably, but show also some consistent patterns.
This suggests that there is indeed a strong causal signal in the data, but it seems hard toconclude which of the various methods is best equipped to extract this signal most reliably,
because the ground truth is only partially known. In future work, we plan to analyze more
recent cytometry data sets that will allow for a more principled validation. Because often
the true causal structure is not known, while interventional data is available, this requires
to extend JCI-based causal discovery with causal prediction techniques, enabling one topredict the results of a particular intervention (Magliacane et al., 2018).
JCI offers increased flexibility when it comes to designing experiments for the purposeof causal discovery, as the JCI framework facilitates analysis of data from almost arbitraryexperimental designs. This allows researchers to trade off the number and complexity ofexperiments to be done with the reliability of the analysis of the data for the purpose of
causal discovery. Compared with existing methods, the framework offered by JCI is themost generally applicable, handling various intervention types and other context changes ina unified and non-parametric way, allowing for latent variables and cycles, and also applies
83Mooij, Magliacane and Claassen

when intervention types and targets are unknown, a common situation in causal discovery
for complex systems.
As future work, we plan to (i) weaken the faithfulness assumption of JCI with respect
to the context variables to allow for even more general experimental designs, (ii) address
the problem of learning from data sets with non-identical (but overlapping) sets of observed
variables, (iii) address selection bias, (iv) develop algorithms that need less computation
time for delivering reliable results, (v) work on more applications on real-world data.

Acknowledgments
We thank Thijs van Ommen for useful discussions and the reviewers and editor for theirconstructive comments. SM, JMM and TC were supported by NWO, the Netherlands
Organization for Scientific Research (VIDI grant 639.072.410). SM was also supported by
the Dutch programme COMMIT/ under the Data2Semantics project. TC was supported
by EU-FP7 grant agreement n.603016 (MATRICS).

Appendix A. Proofs
In this appendix we provide the proofs that were omitted from the main text.
A.1. JCI Foundations
Theorem 13 Assume that JCI Assumptions 0, 1 and 2 hold for SCM M:k ∈ K,
Ck = fk (CpaH (k)∩K , EpaH (k)∩J ),M : Xi = fi (XpaH (i)∩I , CpaH (i)∩K , EpaH (i)∩J ),i ∈ I,
QP(E) = j∈J P(Ej ),For any other SCM M̃ satisfying JCI Assumptions 0, 1 and 2 that is the same as M exceptthat it models the context differently, i.e., of the form˜k ∈ K,
Ck = fk (CpaH̃ (k)∩K , EpaH̃ (k)∩J˜ ),M̃ : Xi = fi (XpaH (i)∩I , CpaH (i)∩K , EpaH (i)∩J ),i ∈ I,
QP(E) = j∈J˜ P(Ej ),with J ⊆ J˜ and paH (i) = paH̃ (i) for all i ∈ I, we have that(i) the conditional system graphs coincide: G(M)do(K) = G(M̃)do(K) ;(ii) if M̃ and M induce the same context distribution, i.e., PM (C) = PM̃ (C), then forany perfect intervention on the system variables do(I, ξI ) with I ⊆ I (including thenon-intervention I = ∅), M̃do(I,ξI ) is observationally equivalent to Mdo(I,ξI ) .
(iii) if the context graphs G(M̃)K and G(M)K induce the same separations, then also G(M̃)
and G(M) induce the same separations (where “separations” can refer to either dseparations or σ-separations).
84

Joint Causal Inference from Multiple Contexts

Proof Let M be an SCM of the form (5). Under JCI Assumption 1, the structural
equations for the context variables do not depend on the system variables:
Ck = fk (CpaH (k)∩K , EpaH (k)∩J ),k ∈ K.

Because of JCI Assumption 2, paH (K) ∩ paH (I) ∩ J = ∅, i.e., the context variables do notshare any exogenous variable with the system variables. This means that in G(M), anyedge between a context variable and a system variable must be a directed edge pointing
from context to system variable, i.e., of the form k → i with k ∈ K, i ∈ I.
Since the structural equations for the system variables of M̃ coincide with those of M,
their solutions (in terms of the context and exogenous variables) also coincide, even after
any perfect intervention on a subset of the system variables. Since C is independent of
EpaH (I) (both for M as well as for M̃), and since PE = PẼJ by assumption, this implies
that the interventional distributions of M and M̃ coincide for any perfect intervention ona subset of system variables if PM (C) = PM̃ (C).
Assume now that G(M̃)K and G(M)K induce the same separations. In the remainderof this proof, “open” can be read either consistently as “σ-open” or as “d-open”. Note thatby assumption, G(M)do(K) = G(M̃)do(K) , and that the edges in G(M)do(K) are a subset of
those in G(M), and of those in G(M̃). We will prove that G(M̃) and G(M) induce thesame separations by first showing that for any two context nodes connected by a path πin G(M)K such that π is A-open in G(M) for some A ⊆ I ∪ K, we can find a path π 0in G(M)K between the two nodes that is A0 -open in G(M) where A0 = A ∩ K ∪ B withB ⊆ K ∩ anG(M)do(K) (A \ K). For π to be A-open in G(M), any collider on π that is nota G(M)-ancestor of A ∩ K must be a G(M)-ancestor of A \ K. Since the latter does notnecessarily imply that the collider must also be G(M̃)-ancestor of A \ K, the idea will be toreplace the variables from A\K in the conditioning set by variables in K∩anG(M)do(K) (A\K)(i.e., context nodes that are guaranteed to be both G(M)-ancestors and G(M̃)-ancestors of
A \ K) that are G(M)-descendants of those colliders that are not already G(M)-ancestorsof A ∩ K. It will turn out that this is not always possible to achieve for π, but that we canconstruct another path π 0 for which this can be done.
Consider a path π in G(M)K between k0 ∈ K and kn ∈ K that is A-open in G(M)for some A ⊆ I ∪ K. We will iteratively construct a walk in G(M)K between the sametwo nodes k0 and kn that is both A-open in G(M) and (A ∩ K) ∪ B-open in G(M), whereB ⊆ K ∩ anG(M)do(K) (A \ K). We will proceed by induction. Suppose a walk πm in G(M)K
between k0 and kn is A-open in G(M). Then it is A ∪ Bm -open in G(M) where Bm =(K ∩ anG(M)do(K) (A \ K)) \ ncol(πm ). Consider the “problematic” colliders col(πm ) \anG(M) (A ∩ K ∪ Bm ) on πm , i.e., the ones that are not ancestors of A ∩ K ∪ Bm . If there
are any, choose one such problematic collider c ∈ K on πm . Since c is not G(M)-ancestor ofA ∩ K ∪ Bm , but πm is A-open, it has to be G(M)-ancestor of A \ K. This means that there
is a directed path in G(M) that starts at c, passes through zero or more context nodes,
none of which lie in A ∩ K ∪ Bm by assumption, and then through zero or more systemnodes, until it ends at a system node in A \ K. Let kc ∈ K be the last context node on thisdirected path before the path crosses the context-system boundary. By assumption, kc must
exist as a non-collider on πm (otherwise it would be in Bm and c would be G(M)-ancestorof Bm ), hence we can make a shortcut by replacing the subwalk of πm between c and kc
85

Mooij, Magliacane and Claassen

by a directed path c → · · · → kc in G(M)K , which necessarily entirely consists of context
nodes that are not in A. If kc occurs more than once on this new walk, remove the entire
subwalk between the two outermost occurrences of kc , such that kc only occurs once. This
new walk πm+1 must be A-open: c (if still present) is now a non-collider that is not in A,
none of the (non-collider) nodes on the directed path (if still present) between c and kc arein A, and kc itself is not in A and is a G(M)-ancestor of A, so it does not matter whether
it is a collider or non-collider. The number of problematic colliders on πm+1 is at least one
less than on πm : c is no longer a collider, and if kc became a collider on πm+1 , it won’t be
problematic (as it is itself in (K ∩ anG(M)do(K) (A \ K)) and cannot also occur as non-collider
on πm+1 ), We repeat this procedure until no problematic colliders are present anymore.
This yields a walk πM that is both A-open and A0 -open, with A0 = (A ∩ K) ∪ B where
B = BM = (K ∩ anG(M)do(K) (A \ K)) \ ncol(πM ). We now shorten this A0 -open walk πM
in G(M)K into an A0 -open path π 0 in G(M)K .
This implies that there must be an A0 -open path π̃ 0 in G(M̃)K connecting k0 and kn ,by assumption. Every collider on π̃ 0 is a G(M̃)K -ancestor of A ∩ K ∪ B, and hence G(M̃)ancestor of A ∩ K ∪ B, and hence G(M̃)-ancestor of A. Therefore, π̃ 0 is also A0 -open inG(M̃). But then it must also be A-open, as we can add A\K to the conditioning set without
blocking any non-collider on π̃ 0 , and then remove B \ (A \ K) from the conditioning set as
all colliders are still kept open due to either being G(M̃)-ancestor of A ∩ K or of A \ K.
Consider now any path in G(M) that is A-open, for A ⊆ I ∪ K. Any edge on thepath between a system node and a context node must be of the form i ← k (with i ∈ I,
k ∈ K) or k → i, where i is in another strongly-connected component than k and k cannot
be in A (because the path was assumed to be A-open). Replacing each longest subpath
consisting entirely of context nodes k0 . . . kn (with all k0 , . . . , kn ∈ K) by a corresponding
A-open path in G(M̃)K between k0 and kn gives a walk in G(M̃) that by construction is also
A-open in G(M̃). Any system collider on this walk must be a collider on the original path,and therefore G(M)-ancestor of A, and therefore also G(M̃)-ancestor of A. Any system
non-collider on this walk is also a system non-collider on the original path and therefore
not in A or, in case of σ-separation, pointing only to nodes in the same strongly-connectedcomponent of G(M), and hence of G(M̃). Any context non-collider on this walk cannot be
in A, or, in case of σ-separation, points to the same strongly-connected component in G(M̃),
since the replacing path in G(M̃)K was A-open by construction. Any context collider on
this walk that is a G(M)K -ancestor of (A ∩ K) ∪ B, and therefore must be G(M̃)-ancestorof A. The walk can be shortened into an A-open path in G(M̃).
Similarly, one can show that any path in G(M̃) that is A-open, there must be a corresponding path in G(M) that is A-open.

Corollary 14 Assume that JCI Assumptions 0, 1 and 2 hold for SCM M. Then there
exists an SCM M̃ that satisfies JCI Assumptions 0, 1 and 2 and 3, such that(i) the conditional system graphs coincide: G(M)do(K) = G(M̃)do(K) ;(ii) for any perfect intervention on the system variables do(I, ξI ) with I ⊆ I (including
the non-intervention I = ∅), M̃do(I,ξI ) is observationally equivalent to Mdo(I,ξI ) ;86

Joint Causal Inference from Multiple Contexts(iii) if the context distribution PM (C) contains no conditional or marginal independences,then the same σ-separations hold in G(M̃) as in G(M); if in addition, the Directed
Global Markov Property holds for M, then also the same d-separations hold in G(M̃)
as in G(M).
Proof Let M be an SCM of the form (5). Under JCI Assumption 1, the structural
equations for the context variables do not depend on the system variables:
Ck = fk (CpaH (k)∩K , EpaH (k)∩J ),k ∈ K.

Because of JCI Assumption 2, paH (K) ∩ paH (I) ∩ J = ∅, i.e., the context variables do notshare any exogenous variable with the system variables.
Consider now the modified SCM M̃ of the form:k∈KCk = gk (EC ),M̃ : Xi = fi (XpaH (i)∩I , CpaH (i)∩K , EpaH (i)∩J ),i ∈ I,
QP(E) = j∈J˜ P(Ej ),Q
where J˜ = J ∪ {C} contains an additional exogenous variable EC ∈ k∈K Ek with components (EC )k ∈ Ck with distribution P(EC ) = PM (C) and gk the projection on the k thcomponent gk : EC 7→ (EC )k . By construction, this SCM M̃ satisfies JCI Assumptions 1
and 2. The only aspect that requires some work is to prove that M̃ as constructed above
is simple (Definition 4).
Take O ⊆ I and consider the solution function for O according to M:
gO : X (paH (O)\O)∩I × C (paH (O)\O)∩K × E paH (O)∩J → X I\O .
This solves the structural equations for O \ I, and since these are the same for M̃ as for M,
the same solution function works also for M̃. Now take Q ⊆ K and consider the solutionfunction gQ : EC → C Q with components EC → Ck : eC 7→ gk (eC ), k ∈ Q. Any othersolution function can be obtained by composition. We conclude that M̃ is simple. M̃ also
induces the same context distribution PM̃ (E) = PM (E) and satisfies JCI Assumption 3 byconstruction. The other statements now follow by applying Theorem 13, where the only
thing left to show is that G(M̃)K and G(M)K induce the same σ-separations if the context
distribution PM (C) contains no conditional independences, and the same d-separations if
in addition the Directed Global Markov Property holds for M.
Marginalizing out the system variables (both in M as well as in M̃) yields M\I andM̃\I , with graphs G(M\I ) = G(M)K and G(M̃\I ) = G(M̃)K , respectively. By the Generalized Directed Global Markov property, since PM (C) has no conditional independences,
there must be a K-σ-open path in G(M)K between any two context nodes k 6= k 0 ∈ K,
for any K ⊆ K with {k, k 0 } ∩ K = ∅. If the Directed Global Markov property holds forM, then it holds for G(M\I ), and hence there must even be a K-d-open path in G(M)K
between any two context nodes k 6= k 0 ∈ K, for any K ⊆ K with {k, k 0 } ∩ K = ∅. Since byconstruction G(M̃)K contains all bidirected edges k ↔ k 0 , there is a K-d-open path in M̃\Ibetween any two context nodes k 6= k 0 ∈ K, for any K ⊆ K with {k, k 0 } ∩ K = ∅.

87

Mooij, Magliacane and Claassen

The following Lemma and Corollary extend these fundamental results further, whichenables one to state a precise relationship between our JCI approach of jointly modelingsystem and context with alternative approaches based on modeling the system conditionalon its context (e.g., Yang et al. (2018)).
Lemma 23 Let M be an SCM that satisfies JCI Assumptions 0, 1, 2. Then the same
restricted separations hold in G(M) as in the conditional system graph G(M)do(K) , i.e.,X⊥

G(M)do(K)Y | Z ⇐⇒ X ⊥ Y | Z
G(M)whenever X, Y, Z ⊆ I ∪ K with X ∩ K = ∅ and K ⊆ Y ∪ Z (where “separations” can refer
to either d-separations or σ-separations).
Proof Let X, Y, Z ⊆ I ∪ K be such that X ∩ K = ∅ and K ⊆ Y ∪ Z. Let G1 , G2 be twographs in {G(M), G(M)do(K) }. Let π be a path in G1 between a node in X and a nodein Y that is open in G1 and that does not contain any non-endpoint nodes in X ∪ Y . It
cannot have non-endpoint nodes in K, because those would be either in Y (a contradiction),or in Z (and since they would be non-colliders with an outgoing directed edge pointing toanother strongly-connected component, they would block the path, another contradiction).
But then the same path π must be present in G2 as well. It is easy to see that it must also
be open in G2 , since for each i ∈ I, deG1 (i) = deG2 (i) and scG1 (i) = scG2 (i).
We can now formulate the following slightly adapted version of Corollary 14:
Corollary 24 Assume that JCI Assumptions 0, 1 and 2 hold for SCM M. Then there
exists an SCM M̃ that satisfies JCI Assumptions 0, 1 and 2 and 3, such that(i) the conditional system graphs coincide: G(M)do(K) = G(M̃)do(K) ;(ii) as a consequence, the same restricted separations hold in G(M̃) as in G(M) and intheir corresponding conditional system graphs, i.e.,
X ⊥ Y | Z ⇐⇒ X
G(M̃)

⊥

Y | Z ⇐⇒ X

G(M̃)do(K)⊥G(M)do(K)Y | Z ⇐⇒ X ⊥ Y | Z
G(M)whenever X, Y, Z ⊆ I ∪ K with X ∩ K = ∅ and K ⊆ Y ∪ Z (where “separations” canrefer to either d-separations or σ-separations);(iii) for any perfect intervention on the system variables do(I, ξI ) with I ⊆ I (including
the non-intervention I = ∅), and any perfect intervention on all context variablesdo(K, c):PM̃ (X | do(K, c), do(I, ξI )) = PM (X | do(K, c), do(I, ξI ));(iv) as a consequence, PM (X | C) = PM̃ (X | C), and in particular, the same restrictedconditional independences hold, i.e.,
X ⊥⊥ Y | Z ⇐⇒ X
G(M̃)

⊥
⊥

Y | Z ⇐⇒ X

G(M̃)do(K)⊥⊥G(M)do(K)Y | Z ⇐⇒ X ⊥⊥ Y |Z

whenever X, Y, Z ⊆ I ∪ K with X ∩ K = ∅ and K ⊆ Y ∪ Z;
88

G(M)

Joint Causal Inference from Multiple Contexts(v) the context distribution PM̃ (C) contains no conditional or marginal independences.
Proof The same SCM M̃ as constructed in the proof of Corollary 14, but with a generic
distribution of P(EC ) that contains no conditional or marginal independences, is easily seen
to fulfill all requirements.

A.2. Minimal Conditional (In)Dependencies
In this section we generalize two useful Lemmas from Claassen and Heskes (2011) to thecyclic setting.
Definition 25 Let X, Y, Z, S ⊆ V be sets of nodes in a DMG G = hV, E, Fi. Let ⊥ denotea DMG-separation property, e.g., d-separation (⊥d ) or σ-separation (⊥σ ). We say that theminimal separation
X ⊥ Y | S ∪ [Z]G

holds if and only if
X ⊥Y |S ∪ Z
G∧∀Q ( Z : X 6⊥ Y | S ∪ Q.
G

In words: all nodes in Z are required (in the context of the nodes in S) to separate X from
Y . Similarly: we say that the minimal connection
X 6⊥ Y | S ∪ [Z]G

holds if and only if
X 6⊥ Y | S ∪ Z
G∧∀Q ( Z : X ⊥ Y | S ∪ Q.
G

In words: all nodes in Z are required (in the context of the nodes in S) to connect X withY.
Note that despite the notation, a minimal connection is not the logical negation of a minimal
separation.
Minimal connections imply the absence of certain ancestral relations:Lemma 26 Let {X}, {Y }, S, {Z} ⊆ V be mutually disjoint sets of nodes in a DMG G =hV, E, Fi. For both d-separation (⊥d ) and σ-separation (⊥σ ), we have that:
X 6⊥ Y | S ∪ [{Z}] =⇒ Z ∈/ anG ({X, Y } ∪ S).
G

Proof The minimal connection means that all paths between X and Y are closed whenconditioning on S and there exists at least one path between X and Y that is open whenconditioning on S ∪ {Z}. For d-separation, this means that such a path (i) contains acollider not in anG (S), (ii) every collider is in anG (S ∪ {Z}), (iii) every non-collider is not in
S ∪{Z}. For σ-separation, this means that such a path (i) contains a collider not in anG (S),89

Mooij, Magliacane and Claassen(ii) every collider is in anG (S ∪ {Z}), (iii) every non-collider is either not in S ∪ {Z}, or if
it is, it points to neighboring nodes in the same strongly-connected component only.
Thus there exists a path between X and Y that contains a collider in anG ({Z}) thatis not in anG (S). If Z ∈ anG (S) this would be a contradiction. If Z ∈ anG (X), thenwe can consider the walk between X and Y obtained from composing the subpath of theoriginal path between Y and the first collider (starting from Y ) in anG ({Z}) \ anG (S) witha directed path to Z and then on to X, without passing through nodes in S. This walk
between X and Y must be open when conditioning on S, and hence there exists a path
between X and Y that is open when conditioning on S, a contradiction. Similarly we obtaina contradiction if Z ∈ anG (Y ).
On the other hand, minimal separations imply the presence of certain ancestral relations:Lemma 27 Let {X}, {Y }, S, Z ⊆ V be mutually disjoint sets of nodes in a DMG G =hV, E, Fi. For both d-separation (⊥d ) and σ-separation (⊥σ ), we have that:X ⊥ Y | S ∪ [Z] =⇒ Z ⊆ anG ({X, Y } ∪ S).
G

Proof Let Q ( Z. Consider a path between X and Y that is open when conditioning on
S ∪ Q, but becomes blocked when conditioning on S ∪ Z. For d-separation, this means that(i) every collider on the path is in anG (S ∪ Q), (ii) every non-collider is not in S ∪ Q, and(iii) it contains a non-collider in S ∪ Z. For σ-separation, this means that (i) every collider
on the path is in anG (S ∪ Q), (ii) every non-collider is either not in S ∪ Q or if it is, it points
to a neighboring node on the path in another strongly-connected component, and (iii) it
contains a non-collider in S ∪ Z that points to a neighboring node on the path in anotherstrongly-connected component. In both cases, we have that (i) every collider on the path
is in anG (S ∪ Q) and (ii) it contains a non-collider in Z \ Q. Consider a maximal directedsubpath of the path starting at a non-collider U in Z \ Q and stopping at a collider or atan end node. Then U ∈ anG ({X, Y } ∪ S ∪ Q).
So, for each Q ( Z, there exists a U ∈ Z \ Q with U ∈ anG ({X, Y } ∪ S ∪ Q). Thus
for every Zi ∈ Z, we either obtain (taking Q = Z \ {Zi }) an ancestral relation of the form
Zi ∈ anG ({X, Y } ∪ S), or, otherwise, at least Zi ∈ anG (Zj ) for some Zj ∈ Z \ {Zi }. Define a
directed graph A with nodes Z ∪{ω} (where ω represents {X, Y }∪S) and add an edge Zi →ω whenever our construction yields an ancestral relation of the form Zi ∈ anG ({X, Y } ∪ S),
or otherwise, an edge Zi → Zj if our construction yields Zi ∈ anG (Zj ).
Then, taking the transitive closure of the constructed directed graph A and using transitivity of ancestral relations, for any Zi ∈ Z we either obtain Zi ∈ anG ({X, Y } ∪ S), or
Zi is in some strongly-connected component C ⊆ Z in A. In the latter case, we can apply
the reasoning above (taking now Q = Z \ C) to conclude that there exists a Zj ∈ C with
Zj ∈ anG ({X, Y } ∪ S) or Zj ∈ anG (C 0 ) where C 0 ⊆ Z is another strongly-connected component. Since the strongly-connected components of Z form an acyclic structure, repeating
this reasoning a finite number of times, we ultimately conclude that Zi ∈ anG ({X, Y }∪S).
An implication of this is that the intersection of all sets that separate a node X from anode Y can only consist of ancestors of X or Y :90Joint Causal Inference from Multiple Contexts

Proposition 28 Let X, Y ∈ V be different nodes
T in a DMG G = hV, E, Fi. For d-separation(⊥d ) or σ-separation (⊥σ ), consider Z ∗ := {Z ⊆ V : X ∈/ Z, Y ∈/ Z, X ⊥G Y | Z}. Then∗Z ⊆ anG ({X, Y }).
T
Proof First, note that Z ∗ = {Z ⊆ V : X ∈/ Z, Y ∈/ Z, X ⊥G Y | [Z]}. From Lemma 27,
X ⊥G Y | [Z] implies Z ⊆ anG ({X, Y }). Hence Z ∗ ⊆ anG ({X, Y }).

Appendix B. Soundness, Consistency and Completeness Properties of
FCI-JCIIn this appendix we will formulate and prove various results concerning soundness andcompleteness of FCI-JCI variants.
B.1. Preliminaries on MAGs and PAGs
We start by summarizing the basic definitions and results from the theory of maximal
ancestral graphs and partial ancestral graphs (Spirtes et al., 2000; Richardson and Spirtes,
2002; Zhang, 2006, 2008a,b) that we will need.
B.1.1. Directed Maximal Ancestral Graphs
Richardson and Spirtes (2002) introduced a class of graphs known as maximal ancestralgraphs (MAGs). The general formulation of MAGs allows for undirected edges which areuseful when modeling selection bias, but here we will only use directed maximal ancestralgraphs without undirected edges (sometimes abbreviated as DMAGs in the literature) as we
assume for simplicity that there is no selection bias. In order to define a directed maximal
ancestral graph, we need the notion of inducing path.
Definition 29 Let G = hV, E, Fi be an acyclic directed mixed graph (ADMG). An inducing
path between two nodes u, v ∈ V is a path in G between u and v on which every node (except
for the end nodes) is a collider on the path and an ancestor in G of an end node of the path.
We can now state:
Definition 30 A directed mixed graph G = hV, E, Fi is called a directed maximal ancestralgraph (DMAG) if all of the following conditions hold:1. Between any two different nodes there is at most one edge, and there are no self-cycles;
2. The graph contains no directed or almost directed cycles (“ancestral”);
3. There is no inducing path between any two non-adjacent nodes (“maximal”).
Given an ADMG, we can define a corresponding DMAG (Richardson and Spirtes, 2002):
Definition 31 Let G = hV, E, Fi be an ADMG. The directed maximal ancestral graphinduced by G is denoted DMAG(G) and is defined as DMAG(G) = hṼ, Ẽ, F̃i such thatṼ = V and for each pair u, v ∈ V with u 6= v, there is an edge in DMAG(G) between u andv if and only if there is an inducing path between u and v in G, and in that case the edge
91

Mooij, Magliacane and Claassen

in DMAG(G) connecting u and v is:u → v if u ∈ anG (v),u ← v if v ∈ anG (u),u ↔ v if u 6∈ anG (v) and v 6∈ anG (u).
An important property of the induced DMAG is that it preserves all ancestral and nonancestral relations. More precisely, for two nodes u, v in ADMG G: u ∈ anG (v) if and only ifu ∈ anDMAG(G) (v). Another important property of the induced DMAG is that it preserves
all d-separations. Indeed, A ⊥dDMAG(G) B | C ⇐⇒ A ⊥dG B | C for all A, B, C ⊆ V. We
sometimes identify a DMAG H with the set of ADMGs G such that DMAG(G) = H. For
an acyclic SCM M, we will define its induced DMAG as DMAG(M) := DMAG(G(M)).
For a directed maximal ancestral graph H, define its independence model to bedIM(H) := {hA, B, Ci : A, B, C ⊆ V, A ⊥ B | C},
Hi.e., the set of all d-separations entailed by the DMAG. For a simple SCM M with endogenous index set I and distribution PM (X), we define its independence model to be
IM(M) := {hA, B, Ci : A, B, C ⊆ I, XA ⊥⊥ XB | XC },PM

i.e., the set of all (conditional) independences that hold in its (observational) distribution.
If M is acyclic, then by the Markov property, IM(M) ⊇ IM(DMAG(M)); the faithfulness
assumption then means that IM(M) ⊆ IM(DMAG(M)).
B.1.2. Directed Partial Ancestral Graphs
Since in many cases, the true DMAG is unknown, it is often convenient when performing
causal reasoning to be able to represent a set of hypothetical DMAGs in a compact way. For
this purpose, partial ancestral graphs (PAGs) have been introduced (Zhang, 2006). Again,
since we are assuming no selection bias for simplicity, we will only discuss directed PAGs
(that is, PAGs without undirected or circle-tail edges, i.e., edges of the form {−−, −−◦ , ◦−−}).
Definition 32 We call a mixed graph G = hV, Ei with nodes V and edges E of the types{→, ←, ←◦ , ↔, ◦−◦ , ◦→} a directed partial ancestral graph (DPAG) if:
1. Between any two different nodes there is at most one edge, and there are no self-cycles;
2. The graph contains no directed or almost directed cycles (“ancestral”);
3. There is no inducing path between any two non-adjacent nodes (“maximal”).
Given a DMAG or DPAG, its induced skeleton is an undirected graph with the same nodes
and with an edge between any pair of nodes if and only if the two nodes are adjacent in the
DMAG or DPAG. We often identify a DPAG with the set of all DMAGs that have the same
skeleton as the DPAG, have an arrowhead (tail) on each edge mark for which the DPAG
has an arrowhead (tail) at that corresponding edge mark, and for each circle in the DPAG,
have either an arrowhead or a tail at the corresponding edge mark. Hence, the circles in a
DPAG can be thought of as to represent either an arrowhead or a tail.
92

Joint Causal Inference from Multiple Contexts

We extend the definitions of (directed) walks, (directed) paths and colliders for directed
mixed graphs to apply also to DPAGs. Edges of the form i ← j, i ←◦ j, i ↔ j are called into
i, and similarly, edges of the form i → j, i ◦→ j, i ↔ j are called into j. Edges of the formi → j and j ← i are called out of i. In addition, we define:
Definition 33 A path v0 , e1 , v1 , . . . , vn between nodes v0 and vn in a DPAG G = hV, Ei iscalled a possibly directed path from v0 to vn if for each i = 1, . . . , n, the edge ei between
vi−1 and vi is not into vi−1 (i.e., is of the form vi−1 ◦−◦ vi , vi−1 ◦→ vi , or vi−1 → vi ). The
path is called uncovered if every subsequent triple is unshielded, i.e., vi and vi−2 are notadjacent in G for i = 2, . . . , n.
If H1 and H2 are DMAGs, then we call them Markov equivalent if IM(H1 ) = IM(H2 ).
One can show that this implies that H1 and H2 must have the same skeleton and the sameunshielded colliders. The FCI algorithm maps the independence model IM(H) of a DMAG
H to a DPAG P. Zhang (2008a) showed that FCI is sound and complete, which means that• P has the same skeleton as H;• As a set of DMAGs, P contains H and all Markov equivalent DMAGs;• For every circle edge mark in P, there exists a DMAG in P Markov equivalent to
H that has a tail at the corresponding edge mark, and there exists a DMAG in P
Markov equivalent to H that has an arrowhead at the corresponding edge mark.
We will denote the completely oriented directed partial ancestral graph that contains theMarkov equivalence class of a DMAG H by CDPAG(H). For an acyclic SCM M we will denote its corresponding CDPAG representation as CDPAG(M) := CDPAG(DMAG(G(M))).
We will make use of the notion of (in)visible edges in a DMAG (Zhang, 2008b):
Definition 34 A directed edge i → j in a DMAG is said to be visible if there is a node k
not adjacent to j, such that either there is an edge between k and i that is into i, or there
is a collider path between k and i that is into i and every collider on the path is a parent ofj. Otherwise i → j is said to be invisible.
We will use the same notion in a DPAG, but call it definitely visible (and its negation
possibly invisible). If a directed edge in a DPAG is definitely visible, it must be visible inall DMAGs in the DPAG.
B.2. Soundness and Consistency of FCI-JCI
We are now equipped to prove the soundness (and for some cases, completeness) of FCI-JCI,
the adaptation of FCI that incorporates the JCI background knowledge that we introducedin Section 4.2.4.
First, the soundness of FCI-JCI is easy to prove by checking that the soundness of the
FCI orientation rules is not invalidated by the JCI background knowledge.
Theorem 35 Let M be an acyclic SCM that satisfies JCI Assumption 0 and a subset of
JCI Assumptions 1, 2, 3. Suppose that its distribution PM (X, C) is faithful w.r.t. the graph
G(M). With input IM(M), and with the right JCI Assumptions, FCI-JCI outputs a DPAG
that contains DMAG(M).
93

Mooij, Magliacane and Claassen

Proof First note that the skeleton obtained by FCI-JCI must coincide with that ofDMAG(M) (as it would for standard FCI). Indeed, if JCI Assumption 3 is made, thecontext nodes are all adjacent in DMAG(M) by assumption. For all other edges, and also
for edges between context nodes if JCI Assumption 3 is not made: the edge is in the skeleton
found by FCI-JCI if and only if it is in the skeleton of DMAG(M), for the same reason as
for standard FCI.
Furthermore, one can easily see that FCI rule R0 is still sound and will not conflict withthe application of the background knowledge stemming from the JCI Assumptions.
The extra orientation rules to incorporate the JCI background knowledge are easily
seen to be sound themselves, as they just impose additional features on the DPAG that aresatisfied by DMAG(M) by assumption.
By checking the soundness proofs of each of the standard FCI orientation rules R1-R4
and R8-R10 in Zhang (2006), it is obvious that all these rules are sound when applied
on any DPAG as long as (i) it contains the true DMAG and (ii) FCI rule R0 has been
completely applied, i.e., all unshielded colliders have been oriented as such. Rules R5-R7
are not needed since we assumed no selection bias.
Hence all the rules applied by FCI-JCI are sound (as long as they are applied in theprescribed ordering), and hence the final DPAG must contain the true DMAG(M).
In general, soundness of a constraint-based causal discovery algorithm implies consistency of the algorithm when using appropriate conditional independence tests.

Lemma 36 If a conditional independence test, including the choice of the sample-size dependent threshold (to decide between the null and alternative hypothesis), is consistent, thenany sound constraint-based causal discovery algorithm based on the test is asymptoticallyconsistent.
Proof Consistency of the conditional independence test means that for any distribution,
the probability of a Type I or Type II error converges to 0 as sample size N → ∞. Since
the number of tests is finite for a fixed number of variables, and the number of possible
predictions made by the algorithm is finite, the probability of any error then converges to
0. Any constraint-based algorithm that is sound (i.e., that would return correct answers
when using an independence oracle, including the possible answer “unknown”) is therefore
asymptotically consistent if it makes use of that conditional independence test.
As an example of a consistent test, Kalisch and Bühlmann (2007) provide a choice of thethreshold for the standard partial correlation test that ensures asymptotic consistency of
the test under the assumption that the distribution is multivariate Gaussian. Anotherexample of a distribution-free and even strongly-consistent conditional independence test is
proposed by Györfi and Walk (2012). In general, when using the p-value as a test statistic,
one should choose the sample-size dependent threshold αN (where the p-value of the test
result is used to decide “dependence” if p ≤ αN and “independence” otherwise) in such a
way that αN → 0 as N → ∞ at a suitable rate.
94

Joint Causal Inference from Multiple Contexts

B.3. Completeness of FCI-JCI
Regarding completeness, we currently only know how to prove the completeness of thevariants FCI-JCI0 and FCI-JCI123. In particular, we do not know whether FCI-JCI1 iscomplete. The completeness of FCI-JCI0 is obvious, because FCI-JCI0 reduces to thestandard FCI algorithm without additional background knowledge.
Theorem 37 Let M be an acyclic SCM that satisfies JCI Assumption 0. Suppose that its
distribution PM (X, C) is faithful w.r.t. the graph G(M). With input IM(M), the output of
FCI-JCI0 is a CDPAG in which all edge marks that can possibly be identified from IM(M)have been oriented.
Proof Follows immediately from the completeness of FCI (Zhang, 2008a) under the additional assumption of no selection bias.
Proving the completeness of FCI-JCI123 is more work. The proof strategy is to introduce
additional variables that mimic the JCI background knowledge. We can then apply thecompleteness results for the standard FCI algorithm (Zhang, 2008a).
Theorem 38 Let M be an acyclic SCM that satisfies JCI Assumptions 0, 1, 2, 3. Supposethat its distribution PM (X, C) is faithful w.r.t. the graph G(M). With input IM(M), theoutput of FCI-JCI123 is a DPAG in which all edge marks that can possibly be identified
from IM(M) and the JCI background knowledge have been oriented.
Proof The adjacency phase (skeleton search) of FCI-JCI123 and the orientation of unshielded triples by applying FCI rule R0 are both sound, as we have seen in Theorem 35.
Furthermore, the skeleton and unshielded colliders found by FCI-JCI123 will be the same
as found by standard FCI (in particular, note that FCI would not orient any unshielded
colliders on a context node, since the true DMAG(M) does not have these).33
Before continuing with the FCI orientation rules R1–R4 and R8–R10, FCI-JCI123 now
uses the JCI background knowledge to orient the following edges:• k ↔ k 0 for all k 6= k 0 ∈ K;
• k → i for k ∈ K, i ∈ I if k and i are adjacent.
After this background orientation step, the only edges in the skeleton that remain to be(further) oriented are the ones connecting two system variables. Denote the DPAG identified
by FCI-JCI123 so far by P.
Each DMAG H with IM(H) = IM(M) and that satisfies the JCI Assumptions 0, 1, 2, 3
must be contained in P. Consider any such DMAG H. We can extend it to a DMAG H∗ ,˙ ∪{r}˙
defined over an extended set of variables I ∪K
∪˙ K̄ where K̄ := {k̄ : k ∈ K} is a copy ofK, by adding edges r → k for all k ∈ K, adding edges k̄ → k for all k ∈ K, and removing
all bidirected edges k ↔ k 0 for all k 6= k 0 ∈ K (see also Figure 43). By construction, themarginal DMAG of H∗ on I ∪ K is H.34
If we run FCI on IM(H∗ ) then the first stages of the algorithm yield a DPAG P ∗ with:
33. If we would use the results of statistical conditional independence tests on a finite data sample, thenthere could be differences between the DPAGs constructed by FCI and FCI-JCI123 after these stages.
34. For the notion of marginal DMAG, see Richardson and Spirtes (2002).

95

Mooij, Magliacane and Claassen• the skeleton of H∗ , which equals the skeleton of P together with the additionallyconstructed edges (r ∗−∗ k for all k ∈ K, and k̄ ∗−∗ k for all k ∈ K) but without theedges between the context variables (k ∗−∗ k 0 for k 6= k 0 ∈ K);• the unshielded colliders in P plus the additionally constructed unshielded colliders(r ∗→ k ←∗ k̄ for all k ∈ K), which are all identified by rule R0;• the arrowheads identified by rule R1 that could also be found in P, plus the edge
orientations k → i for k ∈ K, i ∈ I ∩ chG(M) (k) that are obtained from rule R1;
This means that the subgraph of DPAG P (obtained by FCI-JCI123 so far) induced on thesystem nodes I is identical to the subgraph of DPAG P ∗ (obtained by FCI so far) inducedon the system nodes I. In addition, all pairs of a system node i ∈ I and a context nodek ∈ K are identically connected in the two DPAGs.
By examining rules R1-R4 and R8-R10 of FCI (which are used without modificationsin FCI-JCI123) in detail, one can check that they will perform exactly the same edge markorientations on P as on P ∗ . For rules R2, R3 and R8-R10 this is obvious because the only
subsets of nodes that play a role in those rules necessarily must be in I. For rules R1 andR4 the situation is only slightly more complicated: a single node appearing in those rules
can be in I ∪ K, while all others must be in I. Hence each of these rules is applicable tosome tuple of nodes in P if and only if it is applicable to the same tuple of nodes in P ∗ .
Hence, the final DPAG obtained by FCI-JCI123 from IM(M) and the final DPAG obtainedby FCI from IM(H∗ ) induce identical subgraphs on the system nodes I.
Thus, if all DMAGs H in IM(M) that satisfy JCI Assumptions 0, 1, 2, 3 have a certain
invariant edge mark on an edge between two system variables, then all extended DMAGs H∗
must have the same invariant edge mark. All these extended DMAGs H∗ must be Markov
equivalent, since FCI arrives at the same CDPAG for all IM(H∗ ). Now suppose FCI-JCI123
left an edge mark on an edge between two system variables unoriented. Then FCI must also
leave the corresponding edge mark unoriented when it is run on IM(DMAG(M)∗ ). Thismeans that there must exist DMAGs that are Markov equivalent to DMAG(M)∗ that havean arrowhead at that spot, but also DMAGs that are Markov equivalent to DMAG(M)∗that have a tail at that spot. Marginalizing those DMAGs down to I ∪ K gives DMAGs
that are Markov equivalent to DMAG(M), satisfy JCI Assumptions 0, 1, 2, 3 and havean arrowhead respectively a tail at that spot. This means that all edge marks between
system variables that could possibly be oriented, have been oriented by FCI-JCI123. This
completes the proof of arrowhead and tail completeness of FCI-JCI123.

B.4. Reading off Definite (Non-)Ancestors From a DPAG
Zhang (2006) conjectured the soundness and completeness of a criterion to read off definite
ancestral relations from a CDPAG. Roumpelaki et al. (2016) proved soundness of thiscriterion.35 We will need a slightly stronger result (with a similar proof) for DPAGs:
35. Roumpelaki et al. (2016) also claim to have proved completeness, but their proof is flawed: the last partof the proof that aims to prove that u, v are non-adjacent appears to be incomplete.

96

Joint Causal Inference from Multiple Contexts

(a)

(b)
Cr
Cα

Cβ

Cγ

Cα

X0

X1

X2

X0

X3

X4

Cᾱ

Cβ

Cβ̄

X1
X3

Cγ

Cγ̄

X2
X4

Figure 43: (a) Example DMAG satisfying JCI Assumptions 1, 2, 3 and (b) corresponding
extended DMAG with additional variables as used in the proof of Theorem 38.

Proposition 39 Let M be an acyclic SCM. Let P be a DPAG that contains DMAG(M),and in which all unshielded colliders in DMAG(M) have been oriented. For two nodesi, j ∈ P: If• there is a directed path from i to j in P, or• there exist uncovered possibly directed paths from i to j in P of the form i, u, . . . , j
and i, v, . . . , j such that u, v are non-adjacent nodes in P,
then i causes j according to M, i.e., i ∈ anG(M) (j).
Proof First, if there is a directed path from i to j in P, it must be in any DMAG in
P, hence there must be a directed path from i to j in DMAG(M) as well. Thereforei ∈ anG(M) (j).
Second, assume that there exist uncovered possibly directed paths from i to j in Pof the form i, u, . . . , j and i, v, . . . , j such that u, v are non-adjacent in P. If DMAG(M)has i → u, the path i, u, . . . , j must actually correspond to a directed path in DMAG(M)because otherwise it would contain unshielded colliders that were not oriented, contradicting the assumptions. If DMAG(M) has i ←∗ u instead, it must have i → v to avoid an
unshielded collider u ∗→ i ←∗ v that was not oriented, and hence must have a directed pathi, v, . . . , j. In both cases, DMAG(M) must have a directed path from i to j, and hencei ∈ anG(M) (j).
Zhang (2006, p. 137) provides a sound and complete criterion to read off definite nonancestors from a CDPAG. It is easy to prove the soundness of the criterion also for (arbitrary) DPAGs:
Proposition 40 Let M be an acyclic SCM. Let P be a DPAG that contains DMAG(M).
For two nodes i, j ∈ P: if there is no possibly-directed path from i to j in P then i ∈/anG(M) (j).
97

Mooij, Magliacane and Claassen

(a)

(b)
C
X1

(c)
C

X2

C

X1

X2

X1

X2

Figure 44: Example to illustrate that directed edges in the DPAG obtained by FCI-JCI123
do not necessarily correspond with a direct cause. (a) Graph G(M), satisfying JCI Assumptions 1, 2, 3; (b) corresponding DMAG(M); (c) corresponding
DPAG P output by FCI-JCI123. The directed edge C → X2 in P identified by
FCI-JCI123 does not correspond with a direct causal effect of C on X2 (notethat there is no directed edge C → X2 in G(M)).

Proof If i ∈ anG(M) (j) then there is a directed path from i to j in DMAG(M). Since P
contains DMAG(M), this must correspond with a possibly-directed path from i to j in P.
These two propositions allow us to read off (a subset of the) ancestral and non-ancestral
relations that are identifiable from the conditional independences in the joint distribution
and the JCI background knowledge (if applicable) from the DPAGs output by the various
FCI variants (FCI-JCI123, FCI-JCI1, FCI-JCI0, FCI).
B.5. Discovering Direct Intervention Targets with FCI-JCI123
One of the features of FCI-JCI123 is that it allows one to read off direct intervention targets
from the DPAG that it outputs. Naı̈vely interpreting a directed edge k → i from context
node k ∈ K to system node i ∈ I in the DPAG output by FCI-JCI123 as meaning thatk directly targets i is incorrect, as can be seen from the example in Figure 44. Here we
propose a provably correct (but possibly incomplete) procedure.
We will first consider how to read direct intervention targets from DMAGs before applying this to DPAGs.
Lemma 41 Let M be an acyclic SCM that satisfies JCI Assumptions 0, 1, 2, 3. For k ∈ K,i ∈ I: if• k → i in DMAG(M), and• for all nodes j in DMAG(M) s.t. k → j → i in DMAG(M), j → i is visible,
then k is a direct cause of i according to M, i.e., k → i ∈ G(M).
Proof Because the edge k → i is present in DMAG(M), there must be an inducing path
between k and i in G(M). This must be a collider path into i where each collider is ancestorof i. First suppose that the path consists of more than a single edge. Denote its first collider(the one adjacent to k) by j. If the first edge on the path would be into k, then it must
be k ↔ j and j must be a context node (because of JCI Assumptions 1, 2). Similarly, allsubsequent nodes on the inducing path (except for the final node i) must be collider nodes
98

Joint Causal Inference from Multiple Contexts

and hence in K. But then the final edge is between a context node and system node i andinto the context node, contradicting JCI Assumption 1 or 2. Hence the first edge on the
inducing path must be k → j. The same edge k → j must then occur in DMAG(M). The
remainder of the inducing path is actually an inducing path between j and i that is into j.
By Zhang (2008b, Lemma 9), j → i is in DMAG(M) and it is invisible. This contradicts
the assumption. Therefore the inducing path in G(M) between k and i must consist of asingle edge. This must be out of k because of JCI Assumptions 1, 2, and is thus necessarilyof the form k → i. Hence k → i is in G(M).
The following result enables us to read off direct intervention targets from the DPAG output
by FCI-JCI123.
Proposition 42 Let M be an acyclic SCM that satisfies JCI Assumptions 0, 1, 2, 3.
Suppose that its distribution PM (X, C) is faithful w.r.t. the graph G(M). Let P be the
DPAG output by FCI-JCI123 with input IM(M). Let k ∈ K, i ∈ I.
• If k is not adjacent to i in P, k is not a direct cause of i according to M, i.e.,k → i 6∈ G(M).
• If:
1. k → i in P, and
2. for all system nodes j ∈ I s.t. k → j in P and j ◦−◦ i or j ◦→ i or j → i in P,
the edge j → i is definitely visible in the DPAG obtained from P by replacing theedge between j and i by j → i,then k is a direct cause of i according to M, i.e., k → i ∈ G(M).
Proof Because FCI-JCI123 is sound (Theorem 38), P contains DMAG(M).
For the first statement: if k is not adjacent to i in P, then the two nodes are not adjacentin any DMAG in P, and in particular, in DMAG(M). This means that k → i 6∈ G(M),because otherwise, k → i would be in DMAG(M), a contradiction.
The second statement follows from Lemma 41 and from the JCI Assumptions 1, 2, 3.
If the context variables represent interventions, then this allows us to learn (a subset of)the direct targets and non-targets of each intervention. While it is easy to see that thiscriterion is sound, we do not know whether it is complete.

Appendix C. ASD: Accounting for Strong Dependencies
The causal discovery and reasoning algorithm that we refer to as ASD (Accounting forStrong Dependencies) used in this work is based on an algorithm proposed by Hyttinen
et al. (2014) and extensions proposed by Magliacane et al. (2016b) and Forré and Mooij(2018). To make this paper more self-contained, we give here a short description of thisalgorithm, referring the reader for more details to the original publications.
Hyttinen et al. (2014) formulate causal discovery as an optimization problem where aloss function is minimized over possible causal graphs. Intuitively, the loss function can bethought of as measuring the amount of evidence against the hypothesis that the data was
generated by an SCM with a particular graph. The loss function depends on the hypothetical causal graph and on a list of input statements. For the purely observational case,99

Mooij, Magliacane and Claassen

nthe input consists of a list S = (aj , bj , Zj , λj ) j=1 of weighted conditional independence
statements. Here, the weighted statement (aj , bj , Zj , λj ) with {aj }, {bj }, Zj disjoint sets ofendogenous variable indices and λj ∈ R̄ := R ∪ {−∞, +∞} encodes that the conditional
independence Xaj ⊥⊥ Xbj | XZj holds with “confidence” λj , where a finite value of λj gives a“soft constraint” and a value of λj = ±∞ imposes a “hard constraint”. Positive weights encode that we have empirical support in favor of the independence, whereas negative weights
encode empirical support against the independence (in other words, in favor of dependence).
The loss function simply sums the absolute weights of all the input statements that would
be violated if the true causal graph would consist of the hypothetical one:X
λj (1λj >0 − 1aj ⊥G bj | Zj ),L(G, S) :=(aj ,bj ,Zj ,λj )∈S

where 1 is the indicator function. While the original implementation by Hyttinen et al.
(2014) makes use of d-separation, Forré and Mooij (2018) show how this can be modified
for σ-separation. Causal discovery can now be formulated as the optimization problem:G ∗ = arg min L(G, S),(7)

G∈G(I)

where G(I) denotes the set of all possible causal graphs with nodes I (ADMGs in theacyclic case, and DMGs in the cyclic case).
The optimization problem (7) may have multiple minima, for example because theunderlying causal graph is not identifiable from the inputs. Nonetheless, some of the featuresof the causal graph (e.g., the presence or absence of a certain directed edge) may still beidentifiable. Let f : G(I) → {0, 1} be a feature, i.e., a Boolean function of the causal graph
G. We employ the method proposed by Magliacane et al. (2016b) for scoring the confidence
that feature f is present in the causal graph by calculating the difference between the
optimal losses under the additional hard constraints that the feature f is present vs. thatthe feature f is absent in the causal graph:C(f, S) :=min

G∈G(I):¬f (G)L(G, S) −min

G∈G(I):f (G)L(G, S).

(8)This confidence is positive if there is less evidence against its presence than against its
absence, negative if there is less evidence against its absence than its presence, and vanishes
if there is as much evidence against its presence as there is against its absence. As features,
we can consider for example the presence of a direct causal relation, the presence of an(ancestral) causal relation, and the presence of a latent confounder.
Magliacane et al. (2016b) showed that this scoring method is sound for oracle inputs.
Theorem 43 For any feature f : G(I) → {0, 1}, the ASD confidence score C(f, S) of (8)is sound and complete for oracle inputs with infinite weights. In other words, C(f, S) = ∞if f is identifiable from the inputs, C(f, S) = −∞ if ¬f is identifiable from the inputs, andC(f, S) = 0 if f is unidentifiable from the inputs.
Additionally, Magliacane et al. (2016b) showed that the scoring method is asymptoticallyconsistent under a consistency condition on the weights that encode the confidence of conditional (in)dependence.
100

Joint Causal Inference from Multiple Contexts

Theorem 44 Assume that the weights are asymptotically consistent, meaning that(−∞ Xaj ⊥6 ⊥ Xbj | XZjPλj →+∞ Xaj ⊥⊥ Xbj | XZj ,(9)P(where → means convergence in probability) as the number of samples N → ∞. Then for
any feature f : G(I) → {0, 1}, the ASD confidence score C(f, S) of (8) is asymptoticallyPPconsistent, i.e., C(f, S) → ∞ if f is identifiably true, C(f, S) → −∞ if f is identifiablyP

false, and C(f, S) → 0 otherwise.
In our experiments, we used the weights proposed in Magliacane et al. (2016b): λj =log pj − log α, where pj is the p-value of a statistical test with independence as null hypothesis, and α is a significance level (e.g., 1%). These weights have the desirable property thatindependences typically get a smaller absolute weight than strong dependencies. This leads
to the strong dependencies dominating the loss function, which explains the acronym ASD(Accounting for Strong Dependencies) that we use here to describe this method. By choosing a sample-size dependent threshold αN such that αN → 0 as N → ∞ at a suitable rate,
these weights may become asymptotically consistent. Kalisch and Bühlmann (2007) provide
a choice of αN for partial correlation tests that ensures asymptotic consistency under theassumption that the distribution is multivariate Gaussian. Another possibility for obtaining
consistent weights would be to base them on the distribution-free and strongly-consistentconditional independence test proposed by Györfi and Walk (2012).

References
R. Ayesha Ali, Thomas S. Richardson, and Peter Spirtes. Markov equivalence for ancestralgraphs. The Annals of Statistics, 37(5B):2808–2837, 2009.
Elias Bareinboim and Judea Pearl. A general algorithm for deciding transportability ofexperimental results. Journal of Causal Inference, 1:107–134, 2013.
Tineke Blom, Anna Klimovskaia, Sara Magliacane, and Joris M. Mooij. An upper boundfor random measurement error in causal discovery. In Proceedings of the 34th Annual
Conference on Uncertainty in Artificial Intelligence (UAI 2018), 2018.
Stephan Bongers, Patrick Forré, Jonas Peters, Bernhard Schölkopf, and Joris M. Mooij.
Foundations of structural causal models with cycles and latent variables. arXiv.org
preprint, arXiv:1611.06221v3 [stat.ME], May 2020. URL https://arxiv.org/abs/1611.
06221v3.
Lin S. Chen, Frank Emmert-Streib, and John D. Storey. Harnessing naturally randomized
transcription to infer regulatory relationships among genes. Genome Biology, 8(10):R219,
2007.
David M. Chickering. Optimal structure identification with greedy search. Journal of
Machine Learning Research, 3:507–554, 2002.
101

Mooij, Magliacane and Claassen

Tom Claassen and Tom Heskes. Causal discovery in multiple models from different experiments. In Advances in Neural Information Processing Systems 23 (NIPS 2010), pages
415–423, Vancouver, British Columbia, Canada, 2010.
Tom Claassen and Tom Heskes. A logical characterization of constraint-based causal discovery. In Proceedings of the 27th Annual Conference on Uncertainty in Artificial Intelligence(UAI 2011), pages 135–144, 2011.
Tom Claassen, Joris M. Mooij, and Tom Heskes. Learning sparse causal models is not
NP-hard. In Ann Nicholson and Padhraic Smyth, editors, Proceedings of the 29th Annual
Conference on Uncertainty in Artificial Intelligence (UAI 2013), pages 172–181. AUAI
Press, 2013.
Diego Colombo and Marloes H. Maathuis. Order-independent constraint-based causal structure learning. The Journal of Machine Learning Research, 15(1):3741–3782, 2014.
Gregory F. Cooper. A simple constraint-based algorithm for efficiently mining observational
databases for causal relationships. Data Mining and Knowledge Discovery, 1(2):203–224,
1997.
Gregory F. Cooper and Edward Herskovits. A Bayesian method for the induction of probabilistic networks from data. Machine Learning, 9:309–347, 1992.
Gregory F. Cooper and Changwon Yoo. Causal discovery from a mixture of experimental
and observational data. In Proceedings of the Fifteenth Conference on Uncertainty in
Artificial Intelligence (UAI 1999), pages 116–125, Stockholm, Sweden, 1999.
A. Philip Dawid. Conditional independence in statistical theory. Journal of the Royal
Statistical Society B, 41(1):1–31, 1979.
A. Philip Dawid. Influence diagrams for causal modelling and inference. International
Statistical Review, 70(2):161–189, 2002.
Vanessa Didelez, A. Philip Dawid, and Sara Geneletti. Direct and indirect effects of sequential treatments. In Proceedings of the 22nd Conference on Uncertainty in Artificial
Intelligence (UAI 2006), pages 138–146. AUAI Press, 2006.
Daniel Eaton and Kevin Murphy. Exact Bayesian structure learning from uncertain interventions. In Proceedings of the Eleventh International Conference on Artificial Intelligenceand Statistics (AISTATS 2007), San Juan, Puerto Rico, 2007.
Doris Entner and Patrik O. Hoyer. On causal discovery from time series data using FCI.
In Proceedings of the Fifth European Workshop on Probabilistic Graphical Models (PGM2010), pages 121–128, 2010.
Ronald A. Fisher. Statistical Methods for Research Workers. Oliver and Boyd, Edinburgh,
1925.
Ronald A. Fisher. The Design of Experiments. Hafner, 1935.
102

Joint Causal Inference from Multiple Contexts

Patrick Forré and Joris M. Mooij. Markov properties for graphical models with cycles andlatent variables. arXiv.org preprint, arXiv:1710.08775 [math.ST], October 2017. URL
https://arxiv.org/abs/1710.08775.
Patrick Forré and Joris M. Mooij. Constraint-based causal discovery for non-linear structural causal models with cycles and latent confounders. In Proceedings of the 34th Annual
Conference on Uncertainty in Artificial Intelligence (UAI 2018), 2018.
Patrick Forré and Joris M. Mooij. Causal calculus in the presence of cycles, latent confounders and selection bias. In Proceedings of the 35th Annual Conference on Uncertainty
in Artificial Intelligence (UAI 2019), 2019.
Nir Friedman, Michal Linial, Iftach Nachman, and Dana Pe’er. Using Bayesian networks
to analyze expression data. Journal of Computational Biology, 7(3-4):601–620, 2000.
Dan Geiger, Thomas Verma, and Judea Pearl. Identifying independence in Bayesian networks. Networks, 20(5):507–534, 1990.
Clive W. J. Granger. Investigating causal relations by econometric models and cross-spectral
methods. Econometrica, 37(3):424–438, 1969.
László Györfi and Harro Walk. Strongly consistent nonparametric tests of conditional
independence. Statistics & Probability Letters, 82:1145–1150, June 2012.
Alain Hauser and Peter Bühlmann. Characterization and greedy learning of interventional
Markov equivalence classes of directed acyclic graphs. Journal of Machine Learning
Research, 13:2409–2464, 2012.
David Heckerman, Dan Geiger, and David M. Chickering. Learning Bayesian networks:
The combination of knowledge and statistical data. Machine Learning, 20:197–243, 1995.
Christina Heinze-Deml, Jonas Peters, and Nicolai Meinshausen. Invariant causal prediction
for nonlinear models. Journal of Causal Inference, 6(2):20170016, 2018.
Antti Hyttinen, Frederick Eberhardt, and Patrick O. Hoyer. Learning linear cyclic causal
models with latent variables. Journal of Machine Learning Research, 13:3387–3439, 2012.
Antti Hyttinen, Frederick Eberhardt, and Matti Järvisalo. Constraint-based causal discovery: Conflict resolution with answer set programming. In Proceedings of the 13th
Conference on Uncertainty in Artificial Intelligence (UAI 2014), pages 340–349, Quebec
City, Quebec, Canada, 2014.
Markus Kalisch and Peter Bühlmann. Estimating high-dimensional directed acyclic graphswith the PC-algorithm. Journal of Machine Learning Research, 8:613–636, March 2007.
Markus Kalisch, Martin Mächler, Diego Colombo, Marloes H. Maathuis, and Peter
Bühlmann. Causal inference using graphical models with the R package pcalg. Journal of Statistical Software, 47(11):1–26, 2012.
103

Mooij, Magliacane and Claassen

Yutaka Kano and Shohei Shimizu. Causal inference using nonnormality. In Proceedingsof the International Symposium on Science of Modeling, the 30th Anniversary of theInformation Criterion, pages 261–270, 2003.
Patrick Kemmeren, Katrin Sameith, Loes A.L. van de Pasch, Joris J. Benschop, Tineke L.
Lenstra, Thanasis Margaritis, Eoghan ODuibhir, Eva Apweiler, Sake van Wageningen, Cheuk W. Ko, Sebastiaan van Heesch, Mehdi M. Kashani, Giannis AmpatziadisMichailidis, Mariel O. Brok, Nathalie A.C.H. Brabers, Anthony J. Miles, Diane
Bouwmeester, Sander R. van Hooff, Harm van Bakel, Erik Sluiters, Linda V. Bakker,
Berend Snel, Philip Lijnzaad, Dik van Leenen, Marian J.A. Groot Koerkamp, andFrank C.P. Holstege. Large-scale genetic perturbations reveal regulatory networks andan abundance of gene-specific repressors. Cell, 157(3):740–752, 2014.
Mikko Koivisto and Kismat Sood. Exact Bayesian structure discovery in Bayesian networks.
Journal of Machine Learning Research, 5:549–573, 2004.
Jan T. A. Koster. On the validity of the Markov interpretation of path diagrams of Gaussian
structural equations systems with correlated errors. Scandinavian Journal of Statistics,
26:413–431, 1999.
Jan T.A. Koster. Markov properties of nonrecursive causal models. Annals of Statistics, 24(5):2148–2177, 1996.
Sara Magliacane, Tom Claassen, and Joris M. Mooij. Joint causal inference from observational and interventional datasets. arXiv.org preprint, arXiv:1611.10351v1 [cs.LG],
November 2016a. URL http://arxiv.org/abs/1611.10351v1.
Sara Magliacane, Tom Claassen, and Joris M. Mooij. Ancestral causal inference. In Advancesin Neural Information Processing Systems 29 (NIPS 2016), pages 4466–4474, Barcelona,
Spain, 2016b.
Sara Magliacane, Thijs van Ommen, Tom Claassen, Stephan Bongers, Philip Versteeg,and Joris M Mooij. Domain adaptation by using causal inference to predict invariant
conditional distributions. In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. CesaBianchi, and R. Garnett, editors, Advances in Neural Information Processing Systems 31(NeurIPS 2018), pages 10869–10879. Curran Associates, Inc., 2018.
Subramani Mani. A Bayesian Local Causal Discovery Framework. PhD thesis, Universityof Pittsburg, March 2006. URL http://d-scholarship.pitt.edu/10181/.
Florian Markowetz, Steffen Grossmann, and Rainer Spang. Probabilistic soft interventions
in conditional Gaussian networks. In Proceedings of the Tenth International Workshop
on Artificial Intelligence and Statistics (AISTATS 2005), Bridgetown, Barbados, 2005.
Christopher Meek. Strong completeness and faithfulness in Bayesian networks. In Proceedings of the 11th Conference on Uncertainty in Artificial Intelligence (UAI 1995), pages
411–419, 1995.
104

Joint Causal Inference from Multiple Contexts

Nicolai Meinshausen, Alain Hauser, Joris M. Mooij, Jonas Peters, Philip Versteeg, andPeter Bühlmann. Methods for causal inference from gene perturbation experiments andvalidation. Proceedings of the National Academy of Sciences of the United States ofAmerica, 113(27):7361–7368, 2016.
Joris M. Mooij and Tom Claassen. Constraint-based causal discovery in the presence ofcycles. arXiv.org preprint, arXiv:2005.00610 [math.ST], May 2020. URL https://arxiv.
org/abs/2005.00610.
Joris M. Mooij and Tom Heskes. Cyclic causal discovery from continuous equilibrium data.
In Proceedings of the 29th Annual Conference on Uncertainty in Artificial Intelligence(UAI 2013), pages 431–439, 2013.
Joris M. Mooij, Jonas Peters, Dominik Janzing, Jakob Zscheischler, and Bernhard
Schölkopf. Distinguishing cause from effect using observational data: Methods and benchmarks. Journal of Machine Learning Research, 17(32):1–102, 2016.
Kevin Murphy. Dynamic Bayesian Networks: Representation, Inference and Learning.
PhD thesis, University of Pittsburg, July 2002. URL http://www.cs.ubc.ca/~murphyk/
Thesis/thesis.pdf.
Radford M. Neal. On deducing conditional independence from d-separation in causal graphswith feedback. Journal of Artificial Intelligence Research, 12:87–91, 2000.
Chris J. Oates, Jim Korkola, Joe W. Gray, and Sach Mukherjee. Joint estimation of multiple
related biological networks. Annals of Applied Statistics, 8(3):1892–1919, 2014.
Chris J. Oates, Jim Q. Smith, and Sach Mukherjee. Estimating causal structure using
conditional DAG models. Journal of Machine Learning Research, 17(1):1880–1903, 2016a.
Chris J. Oates, Jim Q. Smith, Sach Mukherjee, and James Cussens. Exact estimation ofmultiple directed acyclic graphs. Statistics and Computing, 26(4):797–811, 2016b.
Judea Pearl. A constraint propagation approach to probabilistic reasoning. In Proceedings ofthe First Conference on Uncertainty in Artificial Intelligence (UAI 1985), pages 357–370,
1986.
Judea Pearl. Comment: Graphical models, causality, and intervention. Statistical Science,
8:266–269, 1993.
Judea Pearl. Causality: Models, Reasoning and Inference. Cambridge University Press,
2009.
Judea Pearl and Rina Dechter. Identifying independence in causal graphs with feedback. In
Proceedings of the 12th Annual Conference on Uncertainty in Artificial Intelligence (UAI1996), pages 420–426, 1996.
Jonas Peters, Joris M. Mooij, Dominik Janzing, and Bernhard Schölkopf. Causal discovery
with continuous additive noise models. Journal of Machine Learning Research, 15:2009–
2053, 2014.
105

Mooij, Magliacane and Claassen

Jonas Peters, Peter Bühlmann, and Nicolai Meinshausen. Causal inference using invariant
prediction: identification and confidence intervals. Journal of the Royal Statistical Society,
Series B, 78(5):947–1012, 2016.
Jonas Peters, Dominik Janzing, and Bernhard Schölkopf. Elements of Causal Inference:
Foundations and Learning Algorithms. MIT Press, Cambridge, MA, USA, 2017.
Joseph Ramsey and Bryan Andrews. FASK with interventional knowledge recovers edges
from the Sachs model. arXiv.org preprint, arXiv:1805.03108 [q-bio.MN], 2018. URL
https://arxiv.org/abs/1805.03108.
Thomas S. Richardson and Peter Spirtes. Ancestral graph Markov models. The Annals of
Statistics, 30(4):962–1030, August 2002.
Dominik Rothenhäusler, Christina Heinze, Jonas Peters, and Nicolai Meinshausen. BACKSHIFT: Learning causal cyclic graphs from unknown shift interventions. In Advancesin Neural Information Processing Systems 28 (NIPS 2015), pages 1513–1521. Curran
Associates, Inc., 2015.
Anna Roumpelaki, Giorgos Borboudakis, Sofia Triantafillou, and Ioannis Tsamardinos.
Marginal causal consistency in constraint-based causal learning. In Frederick Eberhardt, Elias Bareinboim, Marloes Maathuis, Joris Mooij, and Ricardo Silva, editors,
Proceedings of the UAI 2016 Workshop on Causation: Foundation to Application, number 1792 in CEUR Workshop Proceedings, pages 39–47, Aachen, 2016. URL http://ceur-ws.org/Vol-1792/paper5.pdf.
Karen Sachs, Omar Perez, Dana Pe’er, Douglas A. Lauffenburger, and Garry P. Nolan.
Causal protein-signaling networks derived from multiparameter single-cell data. Science,
308, 2005.
Ilya Shpitser, Robin J. Evans, Thomas S. Richardson, and James M. Robins. Introduction
to nested Markov models. Behaviormetrika, 41(1):3–39, 2014.
Peter Spirtes. Conditional independence in directed cyclic graphical models for feedback.
Technical Report CMU-PHIL-54, Carnegie Mellon University, 1994.
Peter Spirtes. Directed cyclic graphical representations of feedback models. In Philippe
Besnard and Steve Hanks, editors, Proceedings of the Eleventh Conference on Uncertainty
in Artificial Intelligence (UAI 1995), pages 499–506, San Francisco, CA, USA, 1995.
Morgan Kaufmann.
Peter Spirtes, Thomas Richardson, Cristopher Meek, Richard Scheines, and Clark Glymour.
Using path diagrams as a structural equation modelling tool. Sociological Methods &Research, 27:182–225, 1998.
Peter Spirtes, Christopher Meek, and Thomas S. Richardson. An algorithm for causal
inference in the presence of latent variables and selection bias. In Clark Glymour andGregory F. Cooper, editors, Computation, Causation and Discovery, chapter 6, pages
211–252. The MIT Press, 1999.
106

Joint Causal Inference from Multiple Contexts

Peter Spirtes, Clark Glymour, and Richard Scheines. Causation, Prediction, and Search.
MIT press, 2nd edition, 2000.
Eric V. Strobl. A constraint-based algorithm for causal discovery with cycles, latent variables and selection bias. International Journal of Data Science and Analytics, 2018.
Jin Tian and Judea Pearl. Causal discovery from changes. In Proceedings of the 17th
Conference in Uncertainty in Artificial Intelligence (UAI 2001), Seattle, Washington,
USA, 2001.
Robert E. Tillman. Structure learning with independent non-identically distributed data.
In Proceedings of the 26th Annual International Conference on Machine Learning (ICML2009), pages 1041–1048, 2009.
Robert E. Tillman and Peter Spirtes. Learning equivalence classes of acyclic models with latent and selection variables from multiple datasets with overlapping variables. In Proceedings of the 14th International Conference on Artificial Intelligence and Statistics (AISTATS 2011), 2011.
Sofia Triantafillou and Ioannis Tsamardinos. Constraint-based causal discovery from multiple interventions over overlapping variable sets. Journal of Machine Learning Research,
16:2147–2205, 2015.
Sofia Triantafillou, Vincenzo Lagani, Christina Heinze-Deml, Angelika Schmidt, Jesper Tegner, and Ioannis Tsamardinos. Predicting causal relationships from biological data: Applying automated causal discovery on mass cytometry data of human immune cells. Scientific Reports, 7:12724, 2017.
Thijs van Ommen and Joris M. Mooij. Algebraic equivalence of linear structural equation models. In Proceedings of the 33rd Annual Conference on Uncertainty in Artificial
Intelligence (UAI 2017), 2017.
Philip J.J.P. Versteeg and Joris M. Mooij. Boosting local causal discovery in highdimensional expression data. arXiv.org preprint, arXiv:1910.02505v2 [stat.ML], November 2019. URL https://arxiv.org/abs/1910.02505v2. Accepted for publication in
BIBM 2019.
Larry Wasserman. All of Statistics. Springer Texts in Statistics. Springer, 2004.
Sewall Wright. Correlation and causation. Journal of Agricultural Research, 20:557–585,
1921.
Karren D. Yang, Abigail Katcoff, and Caroline Uhler. Characterizing and learning equivalence classes of causal DAGs under interventions. In Proceedings of Machine Learning
Research 80 (ICML 2018), pages 5537–5546, 2018.
Jiji Zhang. Causal Inference and Reasoning in Causally Insufficient Systems. PhD thesis,
Carnegie Mellon University, July 2006. URL http://citeseerx.ist.psu.edu/viewdoc/
download?doi=10.1.1.466.7206&rep=rep1&type=pdf.
107

Mooij, Magliacane and Claassen

Jiji Zhang. On the completeness of orientation rules for causal discovery in the presenceof latent confounders and selection bias. Artificial Intelligence, 172(16-17):1873–1896,
2008a.
Jiji Zhang. Causal reasoning with ancestral graphs. Journal of Machine Learning Research,
9:1437–1474, 2008b.
Kun Zhang, Biwei Huang, Jiji Zhang, Clark Glymour, and Bernhard Schölkopf. Causal
discovery from nonstationary/heterogeneous data: Skeleton estimation and orientation
determination. In Proceedings of the Twenty-Sixth International Joint Conference on
Artificial Intelligence (IJCAI 2017), pages 1347–1353, 2017.

108