Identifying Causal Effects with the R Package causaleffect

This is an electronic reprint of the original article. This reprint may differ from the original in pagination and typographic detail. Author(s): Title: Year: Version: Please cite the original version: All material supplied via JYX is protected by copyright and other intellectual property rights, and duplication or sale of all or part of any of the repository collections is not permitted, except that material may be duplicated by you for your research use or educational purposes in electronic or print form. You must obtain permission for any other use. Electronic or print copies may not be offered, whether for sale or otherwise to anyone who is not an authorised user. Identifying Causal Effects with the R Package causaleffect Tikka, Santtu; Karvanen, Juha Tikka, S., & Karvanen, J. (2017). Identifying Causal Effects with the R Package causaleffect. Journal of Statistical Software, 76(12), 1-30. https://doi.org/10.18637/jss.v076.i12 2017 JSS Journal of Statistical Software February 2017, Volume 76, Issue 12. doi: 10.18637/jss.v076.i12 Identifying Causal Effects with the RPackage causaleffect Santtu Tikka University of Jyvaskyla Juha Karvanen University of Jyvaskyla Abstract 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 (2006b) 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 Rpackage causaleffect is presented, which provides an implementation of this algorithm. Functionality of causaleffect is also demonstrated through examples. Keywords: DAG, do-calculus, causality, causal model, identifiability, graph, C-component, hedge, d-separation. 1. 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, 2Identifying Causal Effects with the RPackage causaleffect 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 Rpackages 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 Rpackage 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)Gdenote the set of observable parents, ancestors and descendants of the node set Ywhile 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 3and the details of its Rimplementation are described in Section 4. Section 5 showcases the usage of causaleffect in Rwith some simple examples, and describes some curious special cases arising from the nature of the algorithm itself. Section 6concludes this paper by providing some examples of similar algorithms, where the work of this paper could be applicable. Journal of Statistical Software 3 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) pcalg Methods for graphical models and causal inference SVMMatch Causal effect estimation and diagnostics with support vector machines (Ratkovic 2015) wfe Weighted linear fixed effects regression models for causal inference (Kim and Imai 2014) Packages for sensitivity analysis and other specific problems in causal inference causalsens Selection bias approach to sensitivity analysis for causal effects (Blackwell 2015) cit Causal inference test (Millstein 2016) ImpactIV Identifying causal effect for multi-component intervention using instrumental variable method (Ding 2012) inferference Methods for causal inference with interference (Saul 2015) MatchingFrontier Computation of the balance – sample size frontier in matching methods for causal inference (King et al. 2015) mediation Causal mediation analysis (Tingley et al. 2014) qualCI Causal inference with qualitative and ordinal information on outcomes(Kashin et al. 2014) SimpleTable Bayesian inference and sensitivity analysis for causal effects from 2 ×2 and 2 ×2×Ktables in the presence of unmeasured confounding (Quinn 2012) treatSens Sensitivity analysis for causal inference (Carnegie et al. 2016) Packages for causal discovery CAM Causal additive model (CAM) (Peters and Ernest 2015) D2C Predicting causal direction from dependency features (Bontempi et al. 2015) pcalg Methods for graphical models and causal inference Packages for identification of causal effects causaleffect Deriving expressions of joint interventional distributions in causal models pcalg Methods for graphical models and causal inference Table 1: Rpackages for causal inference. 4Identifying Causal Effects with the RPackage causaleffect 2. Example on do-calculus Consider identification of causal effect Px(y)in the graph Gof 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 Px(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,Yand Zbe pairwise disjoint sets of nodes in the graph Ginduced by a causal model M. Here GX,Zmeans the graph that is obtained from Gby removing all incoming edges of Xand all outgoing edges of Z. Let Pbe 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: Px(y|z,w) = Px(y|w),if (Y |= Z|X,W)GX. 2. Exchanging actions and observations: Px,z(y|w) = Px(y|z,w),if (Y |= Z|X,W)GX,Z. 3. Insertion and deletion of actions: Px,z(y|w) = Px(y|w),if (Y |= Z|X,W)GX,Z(W), where Z(W) = Z\An(W)GX. 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 Px(y)in the causal model of Figure 1, we begin with the factorization Px(y) = X w,z Px(y|w, z)Px(z|w)Px(w).(1) Let us start by focusing on the first term in the sum. Because (Y |= Z|X, W)GX,Z rule 2 implies that Px(y|w, z) = Px,z(y|w) and by noting that (Y |= X|Z, W)GZ,X rule 3 allows us to write Px,z(y|w) = Pz(y|w). By expanding the previous expression we get Pz(y|w) = X x Pz(y|w, x)Pz(x|w).(2) Journal of Statistical Software 5 Figure 1: Graph Gfor the illustrative example. Rule 2 and the fact that (Y |= Z|X, W)GZtogether imply Pz(y|w, x) = P(y|w, x, z).(3) The condition (X |= Z|W)GZand rule 3 allow us to write Pz(x|w) = P(x|w).(4) Inserting (3) and (4) into (2) yields Pz(y|w) = X x P(y|w, x, z)P(x|w).(5) Focusing now on the second term of (1) we see that because (Z |= X|W)GXrule 2 implies that Px(z|w) = P(z|x, w).(6) Similarly, the third term simplifies by using rule 3 and the condition (W |= X)GXrule 3. Px(w) = P(w).(7) Finally, we combine the results above by inserting (5), (6) and (7) into (1) which yields the expression for the causal effect. Px(y) = X w,z X x P(y|w, x, z)P(x|w)!P(z|x, w)P(w) 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 Px(y) = X s X x P(y|x,s)P(x)!P(s|x), whenever the set Sblocks all directed paths from Xto Y, there are no unblocked back-door paths from Xto Sand Xblocks all back-door paths from Sto 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. 3. Identifiability algorithm 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 6Identifying Causal Effects with the RPackage causaleffect (a) Graph G. (b) A subgraph of Ginduced by the set {X, Z1, Z2}. Figure 2: An example illustrating the definition of an induced subgraph. 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). 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. 3.1. 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=hW,Fiand G=hV,Eibe graphs such that W⊂V. If every pair of nodes X, Y ∈Wis connected by an edge in graph Hprecisely when they are connected by an edge of the same direction in graph G, then His an induced subgraph induced by the set Wand 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, Z1 and Z2from Gin 2(a). Perhaps the most important definition is C-component (confounded component). Definition 2 (C-component, (Shpitser and Pearl 2006b) 3).Let G=hV,Eibe a graph. If there exists a set Bsuch that B⊂Eand Bcontains only bidirected edges, and the graph hV,Biis connected, then Gis a C-component. Both graphs in Figure 2are 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. Definition 3 (Maximal C-component).Let Gbe a graph and C=hV,Eia C-component such that C⊂G. C-component Cis maximal (with respect to graph G) if H⊂Cfor every bidirected path Hof graph Gwhich contains at least one node of the set V. Journal of Statistical Software 7 Figure 3: Path H. Tian (2002) proved, that the joint probability distribution P(V)of the observed variables of graph Gcan 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 Gis 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 Gif 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 Gcompletely 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 PP a(Y)(Y). Definition 4 (C-tree, (Shpitser and Pearl 2006b) 4).Let Gbe a C-component such that every observed node has at most one child. If there is a node Ysuch that G[An(Y)G] = G, then Gis 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 Gbe a graph and Yits root set. If Gis a C-component, and every observed node has at most one child, then Gis Y-rooted C-forest. Both C-components in Figure 2are 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 Px(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 Gcontains a hedge for Px(y), then the effect is not identifiable. Definition 6 (Hedge, (Shpitser and Pearl 2006b) 6).Let G=hV,Eibe a graph, and X,Y⊂ Vdisjoint subsets. If there are two R-rooted C-forests F=hVF,EFiand F0=hVF0,EF0i such that VF∩X6=∅,VF0∩X=∅, F0⊂F, and R⊂An(Y)GX, then Fand F0form hedge for Px(y)in G. 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 8Identifying Causal Effects with the RPackage causaleffect form Px(y)in (Tian and Pearl 2002). The result states that Px(y)is identifiable if and only if there are no bidirected paths between Xand any of its children in G[An(Y)G]. Consider the graph H=hV,Eiin Figure 3containing the nodes Xand Yand a bidirected path connecting them formed by the intermediary nodes {Z1, . . . , Zk}. One can observe, that the C-forests Hand H[V\ {X}]form a hedge for Px(Y, Z1, . . . , Zk). 3.2. 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 1for Px(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 Gcontains 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=hV,Eiis an ordering of its nodes, where either X > Y or Y > X for all pairs of nodes X, Y ∈V, X 6=Y in G. In addition, no node can be greater than its descendants in π. In other words, if Xis an ancestor of Yin 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 1is 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 Yin Gare eliminated. This is possible due to the fact that the input of the algorithm assumes that Gis an I-map of Gand 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)GX,W. It is possible to index the nodes of Gand the nodes of any subgraph of Gusing the topological ordering. This property is utilized on lines four, six and seven. The notation V(i−1) πrefers to all nodes in Gthat are smaller than Viin π. Any topological ordering of Gis also a topological ordering for any subgraph of G. This means, that it is unnecessary to determine Journal of Statistical Software 15 object that represents the denominator and fraction is set to TRUE when it is necessary to represent the expression as a fraction. 4.3. Maximal C-components In Section 3.1 it was shown, that for every causal diagram Gthere 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 Gis utilized to find the bidirected edges of G. Definition 8 (adjacency matrix).An adjacency matrix of a graph G=hV,Eiis a n×n matrix A= [aij], where nis the number of nodes of G,V={V1, V2, . . . , Vn}and aij is the number of edges from Vito Vj. Because Gis 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 Vi and Vjare connected by at least one bidirected edge if and only if aij ≥1and aji ≥1. Thus all bidirected edges can be determined by comparing Ato its transpose A>, and by choosing only those edges which correspond to indices with aij ≥1and aji ≥1. The subgraph of Gcontaining 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. 4.4. 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 1consists of the sets xand yincluding the graph G, and returns a probability object, which is a list structure that describes the expression of the causal distribution Px(y)in terms of P(V). The returned object can be further parsed into a character representation. The Rfunction of Algorithm 1is called id. This function takes five parameters as input: a string vector y, a string vector x, a distribution object P, an igraph graph Gand a string vector to. The first four parameters correspond to their mathematical counterparts, namely the vectors x,y,Pand 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 Ras the functions intersect,setdiff and union. The observed portion of Gis saved as G.obs. This graph contains all the observed nodes of Gand the edges between them. In addition, the observed nodes are saved into vector v, and the ancestors of yare saved into vector anc. The implementation of each line of Algorithm 1 is presented next. 16 Identifying Causal Effects with the RPackage causaleffect 1: if x=∅,then return Pv∈v\yP(v). 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 vand ywith the sumset of Pand returns P. 2: if V6=An(Y)G,then return ID(y,x∩An(Y)G, P(An(Y)G), G[An(Y)G)]. 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. 3: let W= (V\X)\An(Y)GX. if W6=∅,then return ID(y,x∪w, P, G). To construct a vector wwhich represents the node set W, one must first construct the subgraph GX. To accomplish this, all incoming edges of Xhave 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 xwith the command E(G) [1:length(E(G)) %->% x], where Eis a function that returns all edges of G. When the subgraph has been constructed, wcan 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). 4: if C(G[V\X]) = {G[S1], . . . , G[Sk]},then return Pv∈v\(y∪x)Qk i=1 ID(si,v\si, P, G). The set C(G[V\X]) can be found with the function c.components. This function determines the node set of every maximal C-component of the input graph, and returns them as a list s. If the length of this list is larger than one, then id returns a new distribution object with sumset = setdiff(v, union(y, x)), recursive = TRUE, children = productlist, where every object in productlist is determined by a new recursive call for every C-component G[Si], i = 1, . . . , k that was found. These components are constructed by calling id with the arguments id(s[[i]], setdiff(v, s[[i]]), P, G, to),i= 1, . . . , k. If the algorithm did not proceed to any of the previous lines, then the additional condition C(G[V\X]) = {G[S]}must be true. The node set of the single C-component G[S]is now saved in the vector s, which was previously a list. This means that sis replaced by s[[1]]. 5: if C(G) = {G},then throw FAIL(G, G[S]). 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 Gitself, then line five is triggered. This is checked by comparing sand v. If they are equal, then the computation is interrupted by the stop function and an error message is produced. The error message Journal of Statistical Software 17 describes the C-forests which form the problematic hedge structure for the causal effect of the current recursion stage. 6: if G[S]∈C(G),then return Pv∈s\yQVi∈SP(vi|v(i−1) π). 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 Viin S. The conditioning nodes are the ones that precede Viin the topological ordering to. 7: if (∃S0)S⊂S0such that G[S0]∈C(G),then return ID(y,x∩s0,QVi∈S0P(Vi|V(i−1) π∩S0, v(i−1) π\s0), G[S0]). 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[S0]. Vector sis replaced by a vector corresponding to the nodes of S0, since the nodes of Sare 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[S0]and every distribution object in productlist is constructed by setting var <- s[i] and cond <- v[0:(ind[i]-1)] for every node Viin S0. Algorithm 2is also implemented in causaleffect as the function idc. This function iterates through the nodes zwhich 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. 5. 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 Rwhile retaining any attributes that have been set for the nodes or edges of the graph. 5.1. 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,xand Gare identical to those of id. The parameter zis optional and it is used to represent the conditioning variables of idc. The initial probability object Pwhich 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,xand z actually contain nodes that are present in G. If Gis not a DAG then causal.effect will also terminate. The last parameter expr is a logical variable. If assigned to TRUE,causal.effect 18 Identifying Causal Effects with the RPackage causaleffect will return the expression in L A T EX 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 EX 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 Ritself. This is demonstrated by replicating some of the graphs of Section 3.3. The graph of Figure 1is created as follows. R> library("igraph") R> fig1 <- graph.formula(W -+ X, W -+ Z, X -+ Z, Z -+ Y, X -+ Y, Y -+ X, + simplify = FALSE) R> fig1 <- set.edge.attribute(graph = fig1, name = "description", + index = c(5,6), value = "U") R> ce1 <- causal.effect(y = "Y", x = "X", z = NULL, G = fig1, expr = TRUE) R> ce1 [1] "\\left(\\sum_{W,Z}P(W)P(Z|W,X)\\left(\\sum_{X}P(Y|W,X,Z)P(X|W)\\right)\\right)" Here X -+ Z denotes a directed edge from Xto Z. The argument simplify = FALSE allows the insertion of duplicate edges for the purposes of forming bidirected arcs. Recalling the internal notation from Section 4.1 we must denote the unidirected edges that correspond to a bidirected edge with a special description parameter, and assign its value to "U". This can be done with the set.edge.attribute function of the igraph package. Finally, the expression for the interventional distribution is obtained by using the causal.effect function. Usually one needs to apply the standard Rfunction cat to obtain the expression with only singular slash symbols. R> cat(ce1) \left(\sum_{W,Z}P(W)P(Z|W,X)\left(\sum_{X}P(Y|W,X,Z)P(X|W)\right)\right) To observe unidentifiability, the graph of Figure 5(a) is also constructed and an attempt is made to identify Px(y). R> fig5 <- graph.formula(Z_1 -+ X, X -+ Z_2, Z_2 -+ Y, Z_1 -+ X, X -+ Z_1, + Z_1 -+ Z_2, Z_2 -+ Z_1, Z_1 -+ Y, Y -+ Z_1, X -+ Y, Y -+ X, + simplify = FALSE) R> fig5 <- set.edge.attribute(graph = fig5, name = "description", + index = 4:11, value = "U") R> causal.effect(y = "Y", x = "X", z = NULL, G = fig5, expr = TRUE) Error: Graph contains a hedge formed by C-forests of nodes: {Z_1,X,Z_2} and {Z_2}. Journal of Statistical Software 19 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 1in 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 Px(y). We omit these warnings from the code for clarity. R> gml1 <- parse.graphml("g1.graphml", format = "standard") R> ce2 <- causal.effect(y = "Y", x = "X", z = NULL, G = gml1, expr = TRUE) R> cat(ce2) \left(\sum_{W,Z}P(W)P(Z|W,X)\left(\sum_{X}P(Y|W,X,Z)P(X|W)\right)\right) We see that the result agrees with the one derived from the manually constructed graph. For conditional causal effects, we simply utilize the parameter zof the causal.effect function. For example, we can obtain the formula for Px(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 P(z|w, x) Pz[P(z|w, x)]. 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). 5.2. A complex expression The conditional distributions P(vi|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 20 Identifying Causal Effects with the RPackage causaleffect Figure 7: An example of a graph, where an identifiable causal effect results in a complex expression. of the simplification rules which were described in the previous section, and their inability to handle every situation. The graph Gof Figure 7serves to demonstrate this phenomenon. An attempt is made to identify Px(z1, z2, z3, y)in this graph. Tian (2002) proved this effect to be identifiable, and showed that its expression is Px(z1, z2, z3, y) = P(z1|x, z2)X x P(y, z3|x, z1, z2)P(x, z2). When applying Algorithm 1to this causal effect, it is necessary to compute a conditional distribution P∗(Y|Z2, Z3),where P∗(y, z2, z3) = X x P(y|z2, x, z3, z1)P(z3|z2, x)P(x|z2)P(z2) and Pis the joint distribution of the observed variables of G. Now, the function causal.effect is applied as follows. R> fig7 <- graph.formula(X -+ Z_1, Z_1 -+ Y, Z_3 -+ Y, Z_2 -+ X, + Z_2 -+ Z_1, Z_2 -+ Z_3, X -+ Y, Y -+ X, X -+ Z_3, Z_3 -+ X, + X -+ Z_2, Z_2 -+ X, Y -+ Z_2, Z_2 -+ Y, simplify = FALSE) R> fig7 <- set.edge.attribute(graph = fig7, name = "description", + index = 7:14, value = "U") R> ce3 <- causal.effect(y = c("Z_1", "Z_2", "Z_3", "Y"), x = "X", + z = NULL, G = fig7, expr = TRUE) R> cat(ce3) This results in the expression P(z1|z2, x)(PxP(y|z2, x, z3, z1)P(z3|z2, x)P(x|z2)P(z2)) Px,y P(y|z2, x, z3, z1)P(z3|z2, x)P(x|z2)P(z2) × X x,z3,y P(y|z2, x, z3, z1)P(z3|z2, x)P(x|z2)P(z2)!P(z3|z2) This result is clearly more cumbersome than the one determined by Tian. However, it can be shown that this expression is correct by using do-calculus. Because the set {X, Z2}d-separates Journal of Statistical Software 21 all paths from Z1to Z3, it follows that (Z3 |= Z1|X, Z2)G,so P(z1|z2, x)X x P(y|z2, x, z3, z1)P(z3|z2, x)P(x|z2)P(z2) =P(z1|z2, x)X x P(y|z2, x, z3, z1)P(z3|z2, x, z1)P(x, z2) =P(z1|z2, x)X x P(y, z3|z2, x, z1)P(x, z2), where the second equality is due to the conditional independence of Z1and Z3given Xand Z2. The last line is equivalent with Tian’s expression up to the ordering of terms. It can be shown, that the remaining terms are subtracted from the expression. P(z3|z2)Px,z3,y P(y|z2, x, z3, z1)P(z3|z2, x)P(x|z2)P(z2) Px,y P(y|z2, x, z3, z1)P(z3|z2, x)P(x|z2)P(z2) =P(z3|z2)P(z2) Px,y P(y|z2, x, z3, z1)P(z3|z2, x)P(x, z2). By applying the same logic to the denominator, it follows that P(z3|z2)P(z2) Px,y P(y|z2, x, z3, z1)P(z3|z2, x)P(x, z2)=P(z3|z2)P(z2) Px,y P(y, z3|z2, x, z1)P(x, z2). By using the conditional independence of Z1and Z3given Xand Z2one gets P(z3, z2) PxP(z3|z2, x, z1)P(x, z2)=P(z3, z2) PxP(z3|z2, x)P(x, z2) =P(z3, z2) PxP(z3|z2, x)P(x, z2)=P(z3, z2) PxP(z3, z2, x)=P(z3, z2) P(z3, z2)= 1. The expression produced by causal.effect is correct despite its complexity. 5.3. d-separation Algorithm 1does 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 Viin the topological ordering π, even though at least some nodes on paths preceding Viare often d-separated by some sets of nodes. In these cases, the nodes that are d-separated with Vi could be excluded from the expression, because they are conditionally independent from Vi in G. This situation is demonstrated by determining the expression of Px,w(y)in the graph Gof 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 1would 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). 22 Identifying Causal Effects with the RPackage causaleffect Figure 8: An example of a graph with additional conditional independences. 6. Discussion We have introduced Rpackage 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 IDZ(Bareinboim and Pearl 2012). z-identifiability deals with a situation, where it is possible to utilize a set Zthat is disjoint from Xto 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. •Counterfactual and conditional counterfactual identifiability algorithms ID* and IDC* (Shpitser and Pearl 2007). The work presented in this paper could be utilized to implement these algorithms. Journal of Statistical Software 23 References Bareinboim E, Pearl J (2012). “Causal Inference by Surrogate Experiments: z-Identifiability.” In N de Freitas, K Murphy (eds.), Proceedings of the Twenty-Eight Conference on Uncertainty in Artificial Intelligence, pp. 113–120. AUAI Press. Bareinboim E, Pearl J (2013a). “A General Algorithm for Deciding Transportability of Experimental Results.” Journal of Causal Inference,1, 107–134. doi:10.1515/jci-2012-0004. Bareinboim E, Pearl J (2013b). “Meta-Transportability of Causal Effects: A Formal Approach.” In Proceedings of the 16th International Conference on Artificial Intelligence and Statistics (AISTATS), pp. 135–143. Blackwell M (2015). causalsens: Selection Bias Approach to Sensitivity Analysis for Causal Effects.Rpackage version 0.1.1, URL https://CRAN.R-project.org/package=causalsens. Bontempi G, Olsen C, Flauder M (2015). D2C: Predicting Causal Direction from Dependency Features.Rpackage version 1.2.1, URL https://CRAN.R-project.org/package=D2C. Brandes U, Eiglsperger M, Herman I, Himsolt M, Marshall MS (2002). “GraphML Progress Report Structural Layer Proposal.” In Graph Drawing, volume 2265 of Lecture Notes in Computer Science, pp. 501–512. Springer-Verlag. doi:10.1007/3-540-45848-4_59. Carnegie NB, Harada M, Hill J (2016). treatSens: A Package to Assess Sensitivity of Causal Analyses to Unmeasured Confounding. New York. Rpackage version 2.0.1, URL https: //CRAN.R-project.org/package=treatSens. Csardi G, Nepusz T (2006). “The igraph Software Package For Complex Network Research.” InterJournal,Complex Systems, 1695. doi:10.1142/s0219525914500064. Dandine-Roulland C (2015). ASPBay: Bayesian Inference on Causal Genetic Variants Using Affected Sib-Pairs Data.Rpackage version 1.2, URL https://CRAN.R-project.org/ package=ASPBay. Dawid AP (1979). “Conditional Independence in Statistical Theory.” Journal of the Royal Statistical Society B,41, 1–31. Ding P (2012). ImpactIV: Identifying Causal Effect for Multi-Component Intervention Using Instrumental Variable Method.Rpackage version 1.0, URL https://CRAN.R-project. org/package=ImpactIV. Glynn A, Quinn K (2010). CausalGAM: Estimation of Causal Effects with Generalized Additive Models.Rpackage version 0.1-3, URL https://CRAN.R-project.org/package= CausalGAM. Huang Y, Valtorta M (2006). “Pearl’s Calculus of Intervention Is Complete.” In Proceedings of the Twenty-Second Conference on Uncertainty in Artificial Intelligence, pp. 217–224. AUAI Press. Hyttinen A, Eberhardt F, Järvisalo M (2015). “Do-Calculus When the True Graph Is Unknown.” In Proceedings of the 31st Conference on Uncertainty in Artificial Intelligence, pp. 395–404. AUAI Press. 24 Identifying Causal Effects with the RPackage causaleffect Kalisch M, Mächler M, Colombo D, Maathuis MH, Bühlmann P (2012). “Causal Inference Using Graphical Models with the RPackage pcalg.” Journal of Statistical Software,47(11), 1–26. doi:10.18637/jss.v47.i011. Karvanen J (2015). “Study Design in Causal Models.” Scandinavian Journal of Statistics, 42(2), 361–377. doi:10.1111/sjos.12110. Kashin K, Glynn A, Ichino N (2014). qualCI: Causal Inference with Qualitative and Ordinal Information on Outcomes.Rpackage version 0.1, URL https://CRAN.R-project.org/ package=qualCI. Kim IS, Imai K (2014). wfe: Weighted Linear Fixed Effects Regression Models for Causal Inference.Rpackage version 1.3, URL https://CRAN.R-project.org/package=wfe. King G, Lucas C, Nielsen R (2015). MatchingFrontier:RPackage for Computing the Matching Frontier.Rpackage version 1.0.0, URL https://CRAN.R-project.org/package= MatchingFrontier. Koller D, Friedman N (2009). Probabilistic Graphical Models: Principles and Techniques. The MIT Press. Luo X, Small D, shan Li C, Rosenbaum P (2011). cin: Causal Inference for Neuroscience. Rpackage version 0.1, URL https://CRAN.R-project.org/package=cin. Maler E, Paoli J, Sperberg-McQueen CM, Yergeau F, Bray T (2004). “Extensible Markup Language (XML) 1.0 (Third Edition).” Technical report, W3C. URL http://www.w3.org/ TR/2004/REC-xml-20040204. Marchetti GM, Drton M, Sadeghi K (2015). ggm: Functions for Graphical Markov Models. Rpackage version 2.3, URL https://CRAN.R-project.org/package=ggm. Meinshausen N (2016). InvariantCausalPrediction: Invariant Causal Prediction.Rpackage version 0.6-0, URL https://CRAN.R-project.org/package= InvariantCausalPrediction. Millstein J (2016). cit: Causal Inference Test.Rpackage version 2.1, URL https://CRAN. R-project.org/package=cit. Neto EC, Yandell BS (2014). qtlnet: Causal Inference of QTL Networks.Rpackage version 1.3.6, URL https://CRAN.R-project.org/package=qtlnet. Pearl J (1995). “Causal Diagrams for Empirical Research.” Biometrika,82, 669–688. doi: 10.1093/biomet/82.4.669. Pearl J (2009). Causality: Models, Reasoning and Inference. 2nd edition. Cambridge University Press, New York. Peters J, Ernest J (2015). CAM: Causal Additive Model (CAM).Rpackage version 1.0, URL https://CRAN.R-project.org/package=CAM. Quinn KM (2012). SimpleTable: Bayesian Inference and Sensitivity Analysis for Causal Effects from 2×2and 2×2×KTables in the Presence of Unmeasured Confounding. Rpackage version 0.1-2, URL https://CRAN.R-project.org/package=SimpleTable.