Identifying Causal Effects with the R Package causaleffect

Do-calculus is concerned with estimating the interventional distribution of an action from the observed joint probability distribution of the variables in a given causal structure. All identifiable causal effects can be derived using the rules of do-calculus, but the rules themselves do not give any direct indication whether the effect in question is identifiable or not. Shpitser and Pearl constructed an algorithm for identifying joint interventional distributions in causal models, which contain unobserved variables and induce directed acyclic graphs. This algorithm can be seen as a repeated application of the rules of do-calculus and known properties of probabilities, and it ultimately either derives an expression for the causal distribution, or fails to identify the effect, in which case the effect is non-identifiable. In this paper, the R package causaleffect is presented, which provides an implementation of this algorithm. Functionality of causaleffect is also demonstrated through examples.


Introduction
When discussing causality, one often means the relationships between events, where a set of events directly or indirectly causes another set of events. The aim of causal inference is to draw conclusions from these relationships by using available data and prior knowledge. Causal inference can also be applied when determining the effects of actions on some variables of interest. These types of actions are often called interventions and the results of the interventions are referred to as causal effects.
The causal inference can be divided into three sub-areas: discovering the causal model from the data, identifying the causal effect when the causal structure is known and estimating an identifiable causal effect from the data. Our contribution belongs to the second category, arXiv:1806.07161v1 [stat.ME] 19 Jun 2018 identification of causal effects. As a starting point, we assume that the causal relationships between the variables are known in a non-parametric form and formally presented as a probabilistic causal model (Pearl 1995). Part of the variables may be latent. The causal structure, i.e., the non-parametric causal relationships, can be described using a directed acyclic graph (DAG). A causal effect is called identifiable if it can be uniquely determined from the causal structure on basis of the observations only.
Do-calculus (Pearl 1995) consist of a set of inference rules, which can be used to express the interventional probability distribution using only observational distributions. The rules of docalculus do not themselves indicate the order in which they should be applied. This problem is solved in the algorithm developed by Tian and Pearl (2003) and Shpitser and Pearl (2006b). The algorithm is proved to determine the interventional distribution of an identifiable causal effect. When faced with an unidentifiable effect, the algorithm provides a problematic graph structure called a hedge, which can be thought of as the cause of unidentifiability.
Other R packages for causal inference are summarized in Table 1. It can be seen that in addition to causaleffect, only pcalg (Kalisch et al. 2012) supports the identification of causal effects. pcalg supports the generalized back-door criterion but does not support the frontdoor criterion. Thus, according to our knowledge, causaleffect is the only R package that implements a complete algorithm for the identification of causal effects.
An algorithm equivalent to the one developed by (Shpitser and Pearl 2006b) has been implemented earlier by Lexin Liu in the CIBN software using JavaBayes, which is a graphical software interface written in Java by Fabio Gagliardi Cozman. In addition to causal effect identification CIBN also provides tools for creating and editing graphical models. CIBN is freely available from http://web.cs.iastate.edu/~jtian/Software/CIBN.htm. DAGitty (Textor et al. 2011) provides another free interface for causal inference and causal modeling. One of the main features of DAGitty is finding sufficient adjustment sets for the minimization of bias in causal effect estimation. DAGitty can also be used to determine instrumental variables, which is a feature currently not provided by causaleffect. However, DAGitty does not provide a complete criterion for identifiability.
Familiarity of Pearl's causal model, do-calculus and basic graph theory is assumed throughout the paper. These concepts are briefly reviewed in Appendix A. A more detailed description can be found in (Pearl 2009) and (Koller and Friedman 2009). Notation similar to that of (Shpitser and Pearl 2006b) is also utilized repeatedly in this paper. Capital letters denote variables and small letters denote their values. Bold letters denote sets which are formed of the previous two. The abbreviations P a(Y) G , An(Y) G , and De(Y) G denote the set of observable parents, ancestors and descendants of the node set Y while also containing Y itself. It should also be noted that the shorthand notation of bidirected edges is used to represent the direct effects of an unobserved confounding variable on the two variables at the endpoints of the bidirected edge.
A motivating example is presented in Section 2. The identification algorithm is presented in Section 3 and the details of its R implementation are described in Section 4. Section 5 showcases the usage of causaleffect in R with some simple examples, and describes some curious special cases arising from the nature of the algorithm itself. Section 6 concludes this paper by providing some examples of similar algorithms, where the work of this paper could be applicable.

Packages for specific applications ASPBay
Bayesian inference on causal genetic variants using affected sib-pairs data (Dandine-Roulland 2015) cin Causal inference for neuroscience (Luo et al. 2011) mwa Causal inference in spatiotemporal event data (Schutte and Donnay 2015) qtlnet Causal inference of QTL networks (Neto and Yandell 2014)

Packages for estimation of causal effects from data CausalGAM
Estimation of causal effects with generalized additive models (Glynn and Quinn 2010) InvariantCausalPrediction Invariant causal prediction (Meinshausen 2016) iWeigReg Improved methods for causal inference and missing data problems (Tan and Shu 2013)

Example on do-calculus
Consider identification of causal effect P x (y) in the graph G of Figure 1. We show how this causal effect can be identified by applying do-calculus (Pearl 2009) manually. Later the same example is reconsidered using the identification algorithm.
First, the rules of do-calculus are shortly reviewed. The purpose of do-calculus is to represent the interventional distribution P x (y) by using only observational probabilities. A causal effect is identifiable, if such an expression can be found by applying the rules of do-calculus repeatedly. This result follow directly from the definition of identifiability due to the fact that all observational distributions are assumed identical for the causal models that induce G.
Let X, Y and Z be pairwise disjoint sets of nodes in the graph G induced by a causal model M .
Here G X,Z means the graph that is obtained from G by removing all incoming edges of X and all outgoing edges of Z. Let P be the joint distribution of all observed and unobserved variables of M . Now, the following three rules hold (Pearl 1995): 1. Insertion and deletion of observations: 2. Exchanging actions and observations: 3. Insertion and deletion of actions: The rules of do-calculus can be shown to be true by using d-separation and the definition of the do(·)-operator. Pearl presented proofs for these three rules (Pearl 1995). Do-calculus has also been shown to be complete, meaning that the expressions of all identifiable causal effects can be derived by using the three rules (Shpitser and Pearl 2006b;Huang and Valtorta 2006).
To identify P x (y) in the causal model of Figure 1, we begin with the factorization Let us start by focusing on the first term in the sum.
and by noting that (Y |= X|Z, W ) G Z,X rule 3 allows us to write P x,z (y|w) = P z (y|w).
By expanding the previous expression we get Rule 2 and the fact that (Y |= Z|X, W ) G Z together imply The condition (X |= Z|W ) G Z and rule 3 allow us to write Inserting (3) and (4) into (2) yields Focusing now on the second term of (1) we see that because (Z |= X|W ) G X rule 2 implies that P x (z|w) = P (z|x, w).
Similarly, the third term simplifies by using rule 3 and the condition (W |= X) G X rule 3.
Finally, we combine the results above by inserting (5), (6) and (7) into (1) which yields the expression for the causal effect.
In Section 3.3 we will see how the causal effect can be identified by applying the algorithm of (Shpitser and Pearl 2006b). The previous result highly resembles the front-door criterion, which states that whenever the set S blocks all directed paths from X to Y, there are no unblocked back-door paths from X to S and X blocks all back-door paths from S to Y. However, neither W , Z, or {W, Z} satisfy the role of the set S. The criterion would certainly hold if we removed W from the graph. Even if a causal effect is identifiable, the rules of do-calculus themselves do not guarantee that they could be used to form an expression for the interventional distribution, and that it would contain only observed quantities. It is also not self-evident in which order the rules of do-calculus should be applied to reach the desired expression from the joint distribution of the observed variables P (V).

Identifiability algorithm
To overcome these limitations an identifiability algorithm has been developed by Shpitser and Pearl (2006b). This algorithm can be used to determine the identifiability of any causal effect, in addition of generating the expression for the interventional distribution in the case of an identifiable effect.

Definitions
Some graph theoretic definitions are necessary in order to present the algorithm. The notation mostly follows that of (Shpitser and Pearl 2006b) with some slight alterations for the benefit of the reader.
Definition 1 (Induced Subgraph). Let H = W, F and G = V, E be graphs such that W ⊂ V. If every pair of nodes X, Y ∈ W is connected by an edge in graph H precisely when they are connected by an edge of the same direction in graph G, then H is an induced subgraph induced by the set W and H = G [W].
Defining new graphs using only a set of nodes can easily be achieved using induced subgraphs. For example, the graph in Figure 2(b) is an induced subgraph induced by the nodes X, Z 1 and Z 2 from G in 2(a).
Perhaps the most important definition is C-component (confounded component).
Definition 2 (C-component, (Shpitser and Pearl 2006b) 3). Let G = V, E be a graph. If there exists a set B such that B ⊂ E and B contains only bidirected edges, and the graph V, B is connected, then G is a C-component. Figure 2 are examples of C-components. Even if a graph is not a C-component, at least one of its subgraphs is guaranteed to be a C-component because every subgraph induced by a single node is always a C-component. It is often of greater interest to determine how a given graph can be partitioned in C-components that contain as many nodes as possible.

Both graphs in
if H ⊂ C for every bidirected path H of graph G which contains at least one node of the set V. Tian (2002) proved, that the joint probability distribution P (V) of the observed variables of graph G can always be factorized in such a way, that each term of the resulting product corresponds to a maximal C-component. This property is in a fundamental role in the algorithm, since it can be used to recursively divide the expression of the interventional distribution into simpler expressions.
If a given graph G is not a C-component, it can still be divided into a unique set C(G) of subgraphs, each a maximal C-component of G. This follows from the fact, that there exists a bidirected path between two nodes in G if and only if they belong in the same maximal C-component, which in turn follows from the definition of a maximal C-component. This means, that the bidirected paths of graph G completely define its maximal C-components.
C-trees are a special case of C-components. They are closely related to direct effects, which are causal effects of the form P P a(Y ) (Y ).
Definition 4 (C-tree, (Shpitser and Pearl 2006b) 4). Let G be a C-component such that every observed node has at most one child. If there is a node Y such that G[An(Y ) G ] = G, then G is a Y -rooted C-tree.
Using only C-trees and C-components it is already possible to characterize identifiability of effects on a single variable. C-forest is the multivariate generalization of a C-tree in such a way that the root set, which is the set of nodes {X ∈ G | De(X) G \ {X} = ∅}, contains one or more nodes.
Definition 5 (C-forest, (Shpitser and Pearl 2006b) 5). Let G be a graph and Y its root set. If G is a C-component, and every observed node has at most one child, then G is Y-rooted C-forest.
Both C-components in Figure 2 are also C-forests, because every observed node has at most one child in both graphs. In addition, their root sets consist only of a single node. There exists a connection between C-forests and general causal effects of the form P x (Y). A graph structure formed by a pair of C-trees is used to determine such effects. Shpitser and Pearl (2006b) proved, that if a graph G contains a hedge for P x (y), then the effect is not identifiable.
Definition 6 (Hedge, (Shpitser and Pearl 2006b) 6). Let G = V, E be a graph, and X, Y ⊂ V disjoint subsets. If there are two R-rooted C-forests Hedges are a remarkable structure, since they generalize certain results regarding identifiability. One example of such a result is the condition for identification of a causal effect of the form P x (y) in (Tian and Pearl 2002). The result states that P x (y) is identifiable if and only if there are no bidirected paths between X and any of its children in G[An(Y) G ]. Consider the graph H = V, E in Figure 3 containing the nodes X and Y and a bidirected path connecting them formed by the intermediary nodes {Z 1 , . . . , Z k }. One can observe, that the

Algorithm
Using the previously presented definitions it is now possible to define Algorithm 1, which completely characterizes the identifiability problem of general causal effects. Shpitser and Pearl (2006b) showed, that the expression returned by Algorithm 1 for P x (y) is always correct if the effect in question is identifiable. They also showed, that if the algorithm is interrupted on line five, then the original graph G contains a hedge, preventing the identifiability of the effect. The existence of a hedge is therefore equivalent with unidentifiability. This result also shows the completeness of do-calculus, because the algorithm only applies standard rules of probability manipulations and the three rules of do-calculus. All variables are assumed to be discrete, but the algorithm can also be applied in a continuous case, when the respective sums are replaced with integrals.
The algorithm is required to be able to iteratively process the nodes of the graph, which means that the nodes have to be ordered in some meaningful fashion. This ordering must be able to take the directions of the edges into account, and at least one such ordering must always exist for any given graph. Topological ordering has all of these prerequisite properties.
Definition 7 (Topological Ordering). Topological ordering π of a DAG G = V, E is an ordering of its nodes, where either X > Y or Y > X for all pairs of nodes X, Y ∈ V, X = Y in G. In addition, no node can be greater than its descendants in π. In other words, if X is an ancestor of Y in G, then X < Y .
There exists at least one topological ordering for any DAG, but in some cases there can be multiple orderings. One way to always construct an ordering for a given graph is to begin by determining all nodes without parents, and ordering them arbitrarily. Next, all nodes without parents excluding the nodes found in previous step are determined and again ordered arbitrarily. It is also assigned, that the largest node in the previous step is smaller than the smallest node in the current step. This process is iterated, until all nodes have been ordered.
Algorithm 1 is simple in a sense that at each recursion stage the computation proceeds to exactly one line only. This is easy to see from the fact that after a condition regarding any of the line has been checked, either a return or a FAIL command will be executed. If x = ∅ on line one, then the marginal distribution P (y) is computed instead of a causal effect. This can be achieved by marginalizing over the joint distribution P (V). On line two, all non-ancestors of Y in G are eliminated. This is possible due to the fact that the input of the algorithm assumes that G is an I-map of G and thus all necessary conditional independences hold. On line three, interventions are added to the original causal effect, which is feasible due to the third rule of do-calculus, because (Y |= W|X) G X,W .
It is possible to index the nodes of G and the nodes of any subgraph of G using the topological ordering. This property is utilized on lines four, six and seven. The notation V (i−1) π refers INPUT: Value assignments x and y, joint distribution P (v) and a DAG G = V, E . G is an I-map of P . OUTPUT: Expression for P x (y) in terms of P (v) or FAIL(F, F ).
to all nodes in G that are smaller than V i in π. Any topological ordering of G is also a topological ordering for any subgraph of G. This means, that it is unnecessary to determine a new ordering for each subgraph of G. Instead, one can fix the ordering before applying the algorithm.
The maximal C-components of G[V \ X] are determined on line four and their factorization property is utilized. If more than one C-components were found, it is now necessary to calculate a new causal effect for every C-component. The algorithm proceeds to either line five, six or seven in the case if only one C-component was found.
If Algorithm 1 throws FAIL, then the original graph G contains a hedge formed by graph G and G[S] of the current recursion stage, due to which the original effect is not identifiable and computation terminates. If the algorithm continues, then it is necessary to determine whether G[S] is a maximal C-component of G. If this is the case, then the condition of line six has been satisfied. In the other case, the computation of the intervention can be limited to the intersection of sets X and S on line seven.
Identifiability of conditional interventional distributions is characterized by Algorithm 2. This algorithm is a generalization of Algorithm 1 and in fact it utilizes the function ID in the computation. It was constructed by Shpitser and Pearl (2006a) for identifying conditional causal effects i.e., causal effects of the form P x (y|z). They showed, that this algorithm is also sound and complete for identifying all such effects.
The primary focus of this paper however, is the implementation of Algorithm 1. The imple-

INPUT:
Value assignments x, y and z, joint distribution P (v) and a DAG G = V, E . G is an I-map of P . OUTPUT: Expression for P x (y|z) in terms of P (v) or FAIL(F, F ). mentation of Algorithm 2 follows seamlessly from this implementation, because at the bottom of any recursive stack of IDC the function ID is ultimately called, which determines if the original conditional effect is identifiable. The only additional task is to determine whether a suitable node for the d-separation condition exists on line 1.

Application in practice
We return to the example presented in Section 2. The graph of the example along with some subgraphs are shown here in Figure 4. Let G = V, E be a graph such as in Figure 4(a) and a causal effect of interest P x (y), which is to be identified from the joint distribution P (X, Y, Z, W ). Only a single topological ordering exists for the nodes of G, and it is W < X < Z < Y . Clearly x = ∅, V = An(Y ) G and W = ∅, so the first three lines are ignored and line four is triggered, since it is now necessary to identify three new causal effects in the following expression: Consider the first term of the product. Because V = An(W ) G , line two is triggered, and non-ancestors of W are ignored. This results in the first term simplifying to P (w) because Line two is also triggered when computing the second term, and in a subgraph induced by ancestors of Z as in Figure 4(b). Observing that the algorithm proceeds to line 6 and the second term simplifies again The last term P w,x,z (y) triggers line four, because So it is mandatory to compute P x (y) from P (X|w)P (Y |X, w, z) in the graph corresponding to Figure 4(c). It should be noted, that this causal effect differs from the original effect P x (y), because the joint distribution P (V) of observed variables of G is not the same as the distribution P (X|w)P (Y |X, w, z) of the subgraph of the current recursion stage.
Line two is triggered next, and since Y has no observed ancestors in the graph corresponding to 4(c), it follows that P x (y) = x P (x|w)P (y|x, w, z).
An expression for the original causal effect is obtained by combining the previous results The result agrees with the result derived in Section 2.
Algorithm 1 can also be used to detect unidentifiability. Let F = V, E be a graph of Figure  5(a) and a causal effect of interest P x (y), which is to be identified from P (X, Y, Z 1 , Z 2 ). Let the topological ordering of the nodes of F be The computation starts from line three Z 1 is added to the original intervention, so P z 1 ,x (y) has to be identified. Line four is triggered next, because Since v \ ({y} ∪ {z 1 , x}) = {z 2 }, two new causal effects have to be identified following expression: Consider the first term of the product. Clearly V = An(Z 2 ) F , so the algorithm proceeds to line two, which means that in a subgraph formed by ancestors of Z 2 as in Figure 5 Thus the original effect P x (y) is not identifiable.

Implementation using R
The programming language R (R Core Team 2016) was chosen for the implementation of Algorithm 1. The R packages XML (Temple Lang 2016), igraph (Csardi and Nepusz 2006) and ggm (Marchetti et al. 2015) are utilized repeatedly throughout the implementation.

Graph files
A graph G induced by the causal model is a crucial argument of Algorithm 1. Many file formats for visualizing graphs are available, each with their own strengths and weaknesses. Some of these formats are very simple, and do not differentiate directed and undirected graphs. Some formats offer excessive features for describing causal models, or they might require handling complex syntax, which can be time consuming.
GraphML (Brandes et al. 2002) is a user-friendly file format for graphs. Its features include support for directed graphs and visualizations. GraphML is based on the extensible markup language XML (Maler et al. 2004), which makes processing of graphs files almost effortless. One can also include the names of the nodes within the GraphML file itself, so the user is not limited to having to input the node names themselves inside the R environment. Graphical editors for creating GraphML files are freely available for the user. A special function called parse.graphml has been developed for processing GraphML files. However, the implementation of Algorithm 1 is not limited to GraphML files alone. Any file format supported by the igraph package can be used, as long as the graph follows one of the following notations for bidirected edges.  Bidirected edges can be separated from unidirected edges by using graphical parameters. For this purpose, three distinct notations have been selected to describe bidirected edges, which correspond to unobserved nodes.
The available notations for bidirected edges are shown in Figures 6(a), 6(b) and 6(c). It should be noted, that notations 1 and 2 are almost identical. Because of their similarity, both notations 1 and 2 are referred to as standard notation. Notation 3, as shown in Figure 6(c) differs from the previous two. It is apparent, that this notation cannot be used as such, because it induces loops in the graph which is not allowed in the context of DAGs. However, GraphML format enables the assignment of parameters for the edges, which in turn allows one to separate these edges from their unidirected counterparts. When using notation 3, one must define a parameter called description for the two unidirected edges corresponding to the bidirected edge, and assign its value to "U" (Unobserved). Notation 3 is used in the implementation itself, which is why it is referred to as the internal notation.
The process of importing GraphML files created by a graphical editor is handled by using the R package XML. This package contains the function xmlParse, which is utilized to import graph files into R objects. It should be noted, that these objects only reflect their internal C objects and are thus different from ordinary R objects. This means that the memory reserved by the XML objects has to be freed after the files have been imported. Normally R does this automatically.
Algorithm 1 requires only a small portion of the XML content, and the unnecessary content is removed in the process of searching for the important items. Items of importance are those that contain data about the node names, node count, edge count and the values of the description parameters of the edges. If notation 1 or 2 of Figures 6(a) and 6(b) was used for the bidirected arcs, it is converted to match the internal format of Figure 6(c). The XML search is implemented using the function getNodeSet of the XML package. This function uses XPath, which is a processing language for XML content search (Simpson 2002).
When the crucial information has been extracted, an igraph graph is formed from the remaining content. igraph is a tool for visualizing and processing graphs, and it can handle graphs which may contain millions of nodes due to its implementation in C. This package also offers many useful functions related to Algorithm 1, such as determining the ancestors of a node, constructing a topological ordering and generating induced subgraphs from a set of edges or nodes. One of the main goals of igraph is the effortless implementation of graph algorithms.

Distribution objects
An important question regarding the Algorithm 1 of Section 3.2, is how the probability distribution which changes at each recursive stage should be implemented. An intuitive solution is to construct a distribution object, which maintains the terms currently present in the expression. Distribution objects are recursive by construction as is the algorithm itself.
In practice this means that when any of the lines four, six or seven is triggered, sub objects are formed, which correspond to the product terms of the expression. These sub objects can further branch into sub objects of their own and so forth. causaleffect implements an R class called probability to represent the distribution objects.
Multiple attributes have to be set for the distribution objects in order to present the probability distribution precisely. The string vectors var and cond are one of the most common attributes, because they enable the definition of a simple conditional distribution. A distribution is formed by the variables described in var conditioned on those of cond. For example, let p be a distribution object, and let the values of its attributes be var = "Y" and cond = "X". Therefore object p represents the conditional distribution P (Y |X).
When the distribution is a product, the individual terms are defined in a list of distribution objects called children and a logical variable recursive is set to TRUE to differentiate this object from those containing only a single term. For example, for a distribution which represents the distribution P * = P (Z|X)P (X|Y )P (Y ) one has to set children = list(a,b,c), where the objects a, b and c represent the distributions P (Z|X), P (X|Y ) and P (Y ) respectively.
For marginal distributions a string vector sumset has been defined. The contents of this vector correspond to the variables which the distribution is to be summed over in the discrete case, or integrated over in the continuous case. In simple situations this parameter is not needed, but often with more complex graphs one encounters instances, where the computation of conditionals is no longer straightorward. Suppose one had to compute the marginal distribution P * (X) of X from the joint distribution P * (X, Y, Z) of the previous example.
To achieve this, one has to set sumset = c("Y","Z") for the matching distribution object, because P * (X) = Y,Z P (Z|X)P (X|Y )P (Y ).
The level of complexity increases further when computing conditionals from distributions which consist of multiple product terms. The previously presented attributes are often insufficient to form an expression for the corresponding distribution object. Consider once more the joint distribution P * . Computing the marginal conditional distribution P * (X|Y ) results in The implementation is able to handle similar situations, where the expression can easily be simplified using the following procedure. Any term which does not depend on the summation index, will be placed outside of the sum. Next, it is checked whether any expressions can be simplified by changing the order of summation. Corresponding terms are subtracted if possible.
These simplification rules are not sufficient to handle every situation. For example, the expression X P (Y |X)P (X) cannot be simplified using the procedure above. One cannot remove any terms from within the sum and the summation order is clearly fixed. In situations, where the denominator is necessary in order to correctly form the expression, one needs to include additional attributes called divisor and fraction. These attributes are similar to the attributes children and recursive in a sense that divisor contains the distribution object that represents the denominator and fraction is set to TRUE when it is necessary to represent the expression as a fraction.

Maximal C-components
In Section 3.1 it was shown, that for every causal diagram G there exists a unique set C(G) of maximal C-components of G. To construct this set, one has to begin by determining all bidirected edges of G. Afterwards, a subgraph containing only bidirected edges is formed. This subgraph will contain one or more components, which are connected subgraphs of G. Because these components are disjoint and every pair of nodes within a component is connected by a bidirected path, it follows that they must be the maximal C-components of G. The adjacency matrix of G is utilized to find the bidirected edges of G.
Definition 8 (adjacency matrix). An adjacency matrix of a graph G = V, E is a n × n matrix A = [a ij ], where n is the number of nodes of G, V = {V 1 , V 2 , . . . , V n } and a ij is the number of edges from V i to V j .
Because G is directed, its adjacency matrix is not necessarily symmetric. When notation 3 of Figure 6(c) is used to describe the bidirected edges, it is easy to confirm that two nodes V i and V j are connected by at least one bidirected edge if and only if a ij ≥ 1 and a ji ≥ 1. Thus all bidirected edges can be determined by comparing A to its transpose A , and by choosing only those edges which correspond to indices with a ij ≥ 1 and a ji ≥ 1.
The subgraph of G containing only bidirected edges is constructed by using the function subgraph.edges of the igraph package. This function retains all nodes of the input graph, but removes all the edges that were not given as input. The subgraph returned by this function is further divided into components by using the function decompose.graph which is also provided by igraph.

Implementation
All necessary preparations have been presented to implement Algorithm 1. Any probability distribution can be represented with a corresponding distribution object, and the adjacency matrix provides a method to determine the maximal C-components of G. Other important methods are provided by the igraph package, such as constructing subgraphs and determining the ancestors of a given set of nodes. In this implementation, the input of Algorithm 1 consists of the sets x and y including the graph G, and returns a probability object, which is a list structure that describes the expression of the causal distribution P x (y) in terms of P (V). The returned object can be further parsed into a character representation.
The R function of Algorithm 1 is called id. This function takes five parameters as input: a string vector y, a string vector x, a distribution object P, an igraph graph G and a string vector to. The first four parameters correspond to their mathematical counterparts, namely the vectors x, y, P and G. The last parameter to is a string vector representing some topological ordering of the nodes of G. All required set theoretic operations are included in R as the functions intersect, setdiff and union.
The observed portion of G is saved as G.obs. This graph contains all the observed nodes of G and the edges between them. In addition, the observed nodes are saved into vector v, and the ancestors of y are saved into vector anc. The implementation of each line of Algorithm 1 is presented next.
The truth value of the expression x = ∅ is determined on line 1. This is done by computing the length of x. If the length is zero, then id combines the difference of the sets v and y with the sumset of P and returns P.
The truth value of the condition on line 2 is determined by computing the length of the vector setdiff(v, anc). If the length is not zero, then id is called with the arguments id(y, intersect(x, anc), P, anc.graph, to), where anc.graph is the induced subgraph G[An(Y) G ], which is constructed by using the induced.subgraph function of the igraph package. This function takes a set of nodes and a graph as input, and constructs a subgraph, which retains all of the nodes given as input, and all of the edges between them in the original graph.
To construct a vector w which represents the node set W, one must first construct the subgraph G X . To accomplish this, all incoming edges of X have to be determined. A useful operator is provided by the igraph package to accomplish this. The operator %->% can be used to find incoming or outgoing edges of a node. In this case, one finds the incoming nodes of x with the command E(G) [1:length(E(G)) %->% x], where E is a function that returns all edges of G. When the subgraph has been constructed, w can also be constructed. If the length of w is not zero, then id is called with the arguments id(y, union(x, w), P, G, to).  The function c.components is utilized again in order to find the maximal C-components of G.
If in addition to having only a single C-component this C-component is G itself, then line five is triggered. This is checked by comparing s and v. If they are equal, then the computation is interrupted by the stop function and an error message is produced. The error message describes the C-forests which form the problematic hedge structure for the causal effect of the current recursion stage.
If the single C-component found on line four is one of the maximal C-components of G, then the function id returns a new distribution object. The sumset of this object is set to setdiff(s, y). The distribution is a product so it must also be set, that recursive = TRUE for this new object. The objects in the list children are determined by new recursive calls for every node V i in S. The conditioning nodes are the ones that precede V i in the topological ordering to.
If the single C-component found on line four is not one of the maximal C-components of G, then it must be a subgraph of some maximal C-component G [S ]. Vector s is replaced by a vector corresponding to the nodes of S , since the nodes of S are no longer required. The function id is called with the following attributes id(y, intersect(x, s), probability(recursive = TRUE, children = productlist), s.graph, to), where s.graph is the induced subgraph G[S ] and every distribution object in productlist is constructed by setting var <- Algorithm 2 is also implemented in causaleffect as the function idc. This function iterates through the nodes z which it receives as input in addition to the parameters that were previously defined for the id function. The d-separation condition on line 1 is checked by using the function dSep from the ggm package.

Package causaleffect
The primary goal of the causaleffect package is to provide the implementation described in Section 4. The package also provides a means of importing GraphML files into R while retaining any attributes that have been set for the nodes or edges of the graph.

Using causaleffect in R
The primary function which serves as a wrapper for the functions id and idc is called causal.effect. This function can be called as causal.effect(y, x, z = NULL, G, expr = TRUE) where the parameters y, x and G are identical to those of id. The parameter z is optional and it is used to represent the conditioning variables of idc. The initial probability object P which is a parameter of id does not have to be specified by the user. In essence, causal.effect starts from an empty distribution object, and gradually builds the final expression if possible. Also, the topological ordering to of the function id is automatically generated by the topological.sort function of the igraph package. It is verified, that the vectors y, x and z actually contain nodes that are present in G. If G is not a DAG then causal.effect will also terminate. The last parameter expr is a logical variable. If assigned to TRUE, causal.effect will return the expression in L A T E X syntax. Otherwise, the probability object used internally by id is returned, which can be manually parsed by the user to gain the desired output. The function get.expression is also provided to get a string representation of a probability object. This function currently supports L A T E X syntax only.
First, causaleffect is loaded to demonstrate the usage of the package.

R> library("causaleffect")
The causal.effect function can be utilized without first importing a graph file. One can utilize the igraph package to construct graphs within R itself. This is demonstrated by replicating some of the graphs of Section 3.3. The graph of Figure 1 is created as follows.
The identification fails in this case due to a hedge present in the graph.
Another function provided by causaleffect is parse.graphml which can be called as parse.graphml(file, format = c("standard", "internal"), nodes = c(), use.names = TRUE) Parameter file is the path to the GraphML file the user wishes to convert into an igraph graph. Parameter format should match the notation that is used to denote bidirected edges in the input graph. The vector nodes can be used to give names to the nodes of the graph if they have not been specified in the file itself or alternatively, to replace them. Finally, use.names is a logical vector indicating whether the names of the nodes should be read from the file or not.
We provide an example GraphML file in the replication materials to demonstrate the use of the parse.graphml function. The file g1.graphml contains the graph of Figure 1 in standard notation. This means that we do not have to provide names for the nodes or set the unidentified edges manually. First, we read the file into R. This produces several warnings which can be ignored because they are related to the visual attributes created by the graphical editor that was used to produce g1.graphml. These attributes play no role in the identification of P x (y). We omit these warnings from the code for clarity.
For conditional causal effects, we simply utilize the parameter z of the causal.effect function. For example, we can obtain the formula for P x (z|w) in the graph of Figure 1.
R> cond1 <-causal.effect(y = "Z", x = "X", z = "W", G = gml1, expr = TRUE) R> cat(cond1) \frac{P(Z|W,X)}{\left(\sum_{Z}P(Z|W,X)\right)} In mathematical notation the result reads This is a typical case where the resulting expression is slightly awkward due to the incompleteness of the simplification rules. However, in this case it is easy to see that the expression can be simplified into P (z|w, x).

A complex expression
The conditional distributions P (v i |v (i−1) π ) that are computed on line 6 can sometimes produce difficult expressions when causal effects are determined from complex graphs. This is a result of the simplification rules which were described in the previous section, and their inability to handle every situation. The graph G of Figure 7 serves to demonstrate this phenomenon. An attempt is made to identify P x (z 1 , z 2 , z 3 , y) in this graph. Tian (2002) proved this effect to be identifiable, and showed that its expression is When applying Algorithm 1 to this causal effect, it is necessary to compute a conditional distribution P * (Y |Z 2 , Z 3 ), where P * (y, z 2 , z 3 ) = x P (y|z 2 , x, z 3 , z 1 )P (z 3 |z 2 , x)P (x|z 2 )P (z 2 ) and P is the joint distribution of the observed variables of G. Now, the function causal.effect is applied as follows.
By applying the same logic to the denominator, it follows that x,y P (y, z 3 |z 2 , x, z 1 )P (x, z 2 ) . By using the conditional independence of Z 1 and Z 3 given X and Z 2 one gets The expression produced by causal.effect is correct despite its complexity.

d-separation
Algorithm 1 does not utilize every possible independence property of a given graph G. For example, the conditional distribution of line six is conditioned on all nodes preceding V i in the topological ordering π, even though at least some nodes on paths preceding V i are often d-separated by some sets of nodes. In these cases, the nodes that are d-separated with V i could be excluded from the expression, because they are conditionally independent from V i in G. This situation is demonstrated by determining the expression of P x,w (y) in the graph G of Figure 8. The function causal.effect is utilized R> fig8 <-graph.formula(z -+ x, z -+ w, x -+ y, w -+ y) R> ce3 <-causal.effect(y = "y", x = c("x", "w"), z = NULL, G = fig8, + expr = TRUE) R> cat(ce3) The function returns P (y|x, w) even though Algorithm 1 would return P (y|z,x,w). This is possible because (Y |= Z|X, W ) G . This means that our implementation is able to simplify the expression into P (y|x, w).

Discussion
We have introduced R package causaleffect for deriving expressions of joint interventional distributions in causal models. The task is a specific but important part of causal inference. We believe that our implementation has two practical use cases. First, causaleffect can be simply used to derive expressions of interventional distributions for complex causal models or to check manual derivations. This is an important step in the estimation of causal effects in complicated settings (Karvanen 2015). Second, causaleffect can be used as a building block in simulation studies and automated systems where identifiability needs to be checked for a large number of causal models. An example of this kind usage is already given by Hyttinen et al. (2015). The efficiency of the presented implementation causaleffect could be analyzed further for example by simulation studies. However, an attempt to maximize performance was made by utilizing the most efficient packages available for the processing of graph files and for the objects corresponding to them. The existing simplification rules of the expressions could also be further improved, but it should be noted that sometimes the more complex expression can prove useful.
There have been many recent developments in the field of causality resulting in graph theoretic algorithms similar to ID and IDC. These include for example: • Causal effect z-identifiability algorithm ID Z (Bareinboim and Pearl 2012). z-identifiability deals with a situation, where it is possible to utilize a set Z that is disjoint from X to achieve identifiability.
• Causal effect transportability algorithm sID (Bareinboim and Pearl 2013a). Transportability means, that results obtained from experimental data can be generalized into a larger population, where only observational studies are applicable.
• Causal effect meta-transportability algorithm µsID (Bareinboim and Pearl 2013b). Metatransportability is an extension of the concept of transportability, where the results are to be generalized from multiple experimental studies simultaneously.
The work presented in this paper could be utilized to implement these algorithms.

A.1. Graphs
The definitions that are presented here follow those of (Koller and Friedman 2009). Graph is an ordered pair G = V, E , where V and E are sets such that The elements of V are the nodes of G, and the elements of E are the edges of G. A graph where V 1 is a parent of V 2 and V 2 is a child of V 1 . This can also be notated as V 2 ← V 1 .
In the first case, H is a path from V 1 to V n . In the second case H is a cycle. If n = 1, then H = {V 1 }, ∅ is also a path. A path H is a directed path if all of its edges are directed and point to the same direction, which means that either A node V 2 is a descendant of V 1 in G, if there exists a directed path H from V 1 to V 2 and H ⊂ G. Respectively, V 2 is an ancestor of V 1 in G, if there exists a directed path H from V 2 to V 1 and H ⊂ G. If a graph G does not contain any cycles, it is acyclic. A graph G = V, E is connected if there exists a path H ⊂ G between every pair of nodes V i , V j ∈ V. Examples of paths are cycles are presented in Figure 9.
If a graph is directed it is also possible to consider its subgraphs as undirected graphs, when all of the edges of the graph are regarded as undirected edges. For example, a directed graph contain paths, even if it does not contain any directed paths. The directed graph in Figure  10 contains a path connecting the nodes X and Y , even though they are not connected by a directed path.
Let G = V, E be a graph and Y ⊂ V. Assume that the nodes of Y correspond to some observed variables, and that the set V can also contain nodes, which in turn correspond to some unobserved variables. Then the abbreviations P a(Y) G , An(Y) G , and De(Y) G denote the set of observable parents, ancestors and descendants of the node set Y while also containing Y.

A.2. Causal model
Causal model can be used to describe the functional relationships between variables of interest. In addition, the model enables the formal treatment of actions or interventions on the   variables of the model. Judea Pearl defined the deterministic causal model and its probabilistic counterpart (Pearl 2009, p. 203-205), which are presented in this section.  where pa i is any realization of the unique minimal set of variables PA V i in V \ V i (connoting parents) sufficient for representing f i . Likewise, U V i ⊆ U stand for the unique minimal set of variables in U sufficient for representing f i .
For each causal model M there is a corresponding graph G = W, E . The node set W contains a node for each observed and unobserved variable of M . The edge set E is determined by the functional relationships between the variables of V and U in the causal model M . The set E contains an edge from X to Y if X ∈ PA Y , which means that there is an edge coming into V i from every node required to uniquely define f V i . Likewise, the set E contains an edge from U to every node V i such that U ∈ U V i .
The definition of causal model does not set any limitations for the unobserved variables. Thus any unobserved node can be a parent of an arbitrary number of observed nodes. If every unobserved node is a parent of exactly two observed nodes, then the causal model is a semi-Markovian causal model. Verma (1993) showed, that for any causal model with unobserved variables one can construct a semi-Markovian causal model that encodes the same set of conditional independences. This is why only semi-Markovian models are considered in this paper.
The edges coming from unobserved variables are sometimes denoted as in Figure 11. However, it is common not to include the unobserved nodes in the visual representation of the graph, which serves to simplify the notation. Instead, it is said that there exists a bidirected edge between X and Y , which corresponds to the effect of the unobserved variable. Thus the notation of Figure 12 is utilized instead of the one in Figure 11.
This notation is used in (Huang and Valtorta 2006;Shpitser and Pearl 2006b;Tian 2002). It should be noted, that a bidirected edge is not the same as two directed edges between two nodes, as this would induce a cycle in the graph which is not allowed. Next, the definition of the causal model is expanded by defining a probability distribution for the unobserved variables.
Definition 10 (Probabilistic Causal Model, (Pearl 2009)  where M D is a (deterministic) causal model and P (U) is the joint distribution of the variables in U.
Henceforth in the paper, the term causal model refers to a probabilistic semi-Markovian causal model without exception. Similarly, any graphs discussed will also refer to the graphs induced by these causal models. A graph G induced by a causal model is strongly related to the joint distribution P of all variables in the model, where P = n i=1 P (v i |pa * (V i ) G ) k j=1 P (u j ), and P a * (.) G also contains all unobserved parents. If this relationship holds, then G is an I-map (independence map) of P . Independence properties of G and P are closely related through the following definition Disjoint sets X and Y are said to be d-separated by Z in G if every path from X to Y is d-separated by Z in G.
If X and Y are d-separated by Z in G, then X is independent of Y given Z in every P for which G is an I-map of P . The notation of (Dawid 1979) is used to denote this statement as (X |= Y | Z) G .

A.3. Causal effects
Interventions on a causal model alter the functional relationships between its variables. Any intervention do(X = x) on a causal model M produces a new model M x = U, V, F x , P (U) , where F x is obtained by replacing f X ∈ F for each X ∈ X with a constant function, where the constants are defined as the x values of do(X = x). It is now feasible to formalize the notion of causal effects as follows.
Definition 12 (Causal Effect, (Shpitser and Pearl 2006b)). Let M = U, V, F, P (U) be a causal model and Y, X ⊂ V. The causal effect of do(X = x) on the set Y in M is the marginal distribution of Y in M x , which is noted by P (Y|do(X = x)) = P x (Y).
For every action do(X = x) it is required that P (x|P a(X) G \ X) > 0. This limitation ensures that P x (V) and its marginals are well defined. The restriction stems from the fact that it is unfeasible to force X to attain values which cannot be observed. No inference can be made from the distribution of such an intervention using observational data.
Definition 13 (Causal Effect Identifiability, (Shpitser and Pearl 2006b) 2). Let G = V, E be a graph and Y, X ⊂ V. The causal effect of do(X = x) on the set Y, where Y ∩ X = ∅, is identifiable in G if P 1 x (Y) = P 2 x (Y) for every pair of causal models M 1 and M 2 such that P 1 (V) = P 2 (V) and P 1 (x|P a(X) G \ X) > 0.
It is often impossible to show that a causal effect is identifiable by using solely the definition, because one would have to compare every causal model that agree on the distribution of the