We compare score-based and constraint-based learning in the presence of latent confounders. We use a greedy search strategy to identify the best fitting maximal ancestral graph (MAG) from continuous data, under the assumption of multivariate normality. Scoring maximal ancestral graphs is based on (a) residual iterative conditional fitting [Drton et al., 2009] for obtaining maximum likelihood estimates for the parameters of a given MAG and (b) factorization and score decomposition results for mixed causal graphs [Richardson, 2009, Nowzohour et al., 2015]. We compare the score-based approach in simulated settings with two standard constraintbased algorithms: FCI and conservative FCI. Results show a promising performance of the greedy search algorithm.
Causal graphs can capture the probabilistic and causal properties of multivariate distributions. Under the assumptions of causal Markov condition and faithfulness, the graph induces a factorization for the joint probability distribution, and a graphical criterion (d-separation) can be used to identify all and only the conditional independencies that hold in the joint probability distribution.
The simplest case of a causal graph is a directed acyclic
graph (DAG). A causal DAG G and faithful probability
distribution P constitute a causal Bayesian network (CBN)
[
There exist two main approaches for learning causal graphs
from data. Constraint-based approaches infer the
conditional independencies imprinted in the data and search
for a DAG/MAG that entails all (and only) of these
independences according to d/m-separation. Score-based
approaches try to find the graph G that maximizes the
likelihood of the data given G (or the posterior), according
to the factorization imposed by G. In general, a class of
causal graphs, that are called Markov equivalent, fit the
data equally well. Constraint-based approaches are more
efficient and output a single graph with clear semantics,
but give no indication the relative confidence in the model.
Moreover, they have been shown to be sensitive to error
Numerous constraint-based and score-based algorithms exist that learn causal DAGs (classes of Markov equivalent DAGs) from data. Learning MAGs on the other hand is typically done with constraint-based algorithms. A score-based method for mixed causal graphs (not necessarily MAGs) has recently been proposed [Nowzohour et al., 2015] based on relative factorization results [Tian and Pearl, 2003, Richardson, 2009].
Using these decomposition results, we implemented a
simple greedy search for learning MAGs from data. We
compare the results of this approach with FCI [S
The rest of the paper is organized as follows: Section 2 briefly reviews causal graphs with and without causal sufficiency. Section 3 gives an overview of constraint-based and score-based methods for DAGs and MAGs. Section 4 describes a greedy search algorithm for learning MAGs. Related work is discussed in Section 5. Section 6 compares the performance of the algorithm against FCI and CFCI. Conclusions and future work are presented in 7. 2
We begin with some graphical notation: A mixed graph (MG) is a collection of nodes (interchangeably variables) V, along with a collection of edges E. Edges can be directed (X ! Y ) or bi-directed (X $ Y ). A path is a sequence of adjacent edges (without repetition). The first and last node of a path are called endpoints of the path. A bi-directed path is a path where every edge is bi-directed. A directed path is a path where every edge is directed and oriented in the same direction. We use X 99K Y to symbolize a directed path from X to Y . A directed cycle occurs when there exists a directed path X 99K X. An almost directed cycle is occurs when X $ Y and X 99K Y . A triplet hX; Y; Zi on consequent nodes on a path are form a collider if X ! Y Z. If X and Z are not adjacent, the triplet is an unshielded collider.
A mixed graph is called ancestral if it has no directed and
almost directed cycles. An ancestral graph without
bidirected edges is a DAG. X is a parent of Y in a MG G
if X ! Y in G. We use the notation P aG (X), AnG (X) to
denote the set of parents and ancestors of X in G.
Under causal sufficiency, DAGs can be used to model
causal relationships: For a graph G over a set of variables
V, X ! Y in G if X causes Y directly (no variables in
V mediate this relationship). Under the causal Markov
condition and faithfulness
P (V) = Y P (V jP aG (V ))
V 2V (1) Thus, the parameters of the joint probability distribution describe the probability density function of each variable given its parents in the graph. An interesting property of CBNs, that constitutes the basis of constraint-based learning, is the following: Every missing edge in a DAG of a CBN corresponds to a conditional independence. Hence, if X is independent from Y given Z (symb. X ?? Y jZ) in P, then X and Y are not adjacent in G.
In general, a class of Markov equivalent DAGs fit the data equally well. DAGs in a Markov equivalent class share the same skeleton and unshielded colliders. A Pattern DAG (PDAG) can be used to represent the Markov equivalent class of DAGs: It has the same edges as every DAG in the Markov equivalence class, and the orientations that are shared by all DAGs in the Markov equivalence class. Confounded relationships cannot be represented in DAGs, and mixed causal graphs were introduced to tackle this problem. The most straightforward approach is with semiMarkov causal models (SMCMs) Tian and Pearl [2003]. The graphs of semi-Markov causal models are acyclic directed mixed graph (ADMGs). Bi-directed edges are used to denote confounded variables, and directed edges denote direct causation. A pair of variables can share up to two edges (one directed, one bi-directed). The conditional independencies that hold in a faithful distribution are represented through the criterion of m-separation:
graph G = (E; V), a path p between A and B is mconnecting given (conditioned on) a set of nodes Z , Z V n fA; Bg if 1. Every non-collider on p is not a member of Z. 2. Every collider on the path is an ancestor of some member of Z.
A and B are said to be m-separated by Z if there is no m-connecting path between A and B relative to Z. Otherwise, they are said to be m-connected given Z.For graphs without bi-directed edges, m-separation is reduced to the d-separation criterion.
Markov equivalence classes of semi-Markov causal
models do not have a simple characterization, because Markov
equivalent SMCMs do not necessarily share the same
edges: Absence of an edge in a SMCM does not necessarily
correspond to an m-separation. Figure 1 shows an example
of two SMCMs that encode the same m-separations but do
not have the same edges
Maximal ancestral graphs are also used to model causality and conditional independencies in causally insufficient systems. MAGs are mixed ancestral graphs, which means that they can have no directed or almost directed cycles. Every pair of variables X; Y in an ancestral graph is joined by at most one edge. The orientation of this edge represents (non) causal ancestry. A directed edge X ! Y denotes that X is an ancestor of Y , but the relation is not necessarily direct in the context of modeled variables (see for example edge A ! D in MAG M1 of Figure 1). Moreover, X and Y may also be confounded (e.g. edge B ! D in MAG M1 of Figure 1). A bi-directed edge X $ Y denotes that X and Y are confounded.
Like SMCMs, ancestral graphs encode the conditional
independencies of a faithful distribution according to the
criterion of m-separation. Maximal ancestral graphs are
graphs in which every missing edge (non-adjacency)
corresponds to a conditional independence. Every ancestral
graph can be extended to a maximal ancestral graph by
adding some bi-directed edges [Ri
Figure 1 illustrates some differences in SMCMs and MAGs that represent the same marginal of a DAG. For example, A is a causal ancestor of D in DAG G1, but not a direct cause (in the context of observed variables). Therefore, the two are not adjacent in the corresponding SMCM S1 over fA; B; C; Dg. However, the two cannot be rendered independent given any subset of fB; Cg, and therefore A ! D is in the respective MAG M1.
On the same DAG, B is another causal ancestor (but not a direct cause) of D. The two variables share the common cause L. Thus, in the corresponding SMCM S1 over fA; B; C; Dg B $ D is present. However, a bi-directed edge between B and D is not allowed in MAG M1, since it would create an almost directed cycle. Thus, B ! D is in M1.
Overall, a SMCM has a subset of the adjacencies of its
MAG counterpart. These extra adjacencies in MAGs
correspond to pairs of variables that cannot be m-separated given
any subset of observed variables, but neither directly causes
the other, and the two are not confounded. These
adjacencies can be checked in a SMCM using a special type of
path, called inducing path [Ri
A M1: A P1: A
B B B B C C C C
D D
D D A S2:
A M2: A P2: A
B B B B
C C C C
D D
D D
As mentioned above, there are two main approaches for learning causal networks from data, constraint-based and score-based. Constraint-based methods estimate from the data which conditional independencies hold, using appropriate tests of conditional independence. Each conditional independence corresponds to an m(d)-separation constraint. Constraint-based algorithms try to eliminate all graphs that are inconsistent with the observed constraints, and ultimately return only the statistically equivalent graphs consistent with all the tests.
Notice that the number of possible conditional
independencies are exponential to the number of variables. For graphs
that are maximal, i.e. every missing edge corresponds to a
conditional independence (DAGs and MAGs but not
SMCMs), there exist efficient procedures that can return the
skeleton and invariant orientations of the Markov
equivalence class of graphs that are consistent with a data set,
using only a subset of conditional independencies.
The
Given faithfulness, the subset of conditional
independencies that have been identified during the skeleton
identification stage are sufficient to make all invariant orientation and
return the PDAG or PAG that represents the Markov
equivalence class of causal graphs that are consistent with all
and only the cached conditional independencies. In
practice, the orientation stage is sensitive to error
Score-based methods on the other hand search over the
space of possible graphs trying to maximize a score that
reflects how well the graph fits the data. This score is
typically related to the likelihood of the graph given the
data, P (DjG). For multinomial and Gaussian
parametrizations, respectively, B
The number of possible graphs is super-exponential to the
number of variables. For DAGs, efficient search-and-score
learning is based on the factorization in Equation 1. Thus,
the likelihood of the graph can be decomposed into a
product of individual likelihoods of each node given its parents.
This can make greedy search very efficient: In each step
of the search, all the graphs that occur with single changes
of the current graph are considered. Using the score
decomposition, one only needs to recompute the scores of
nodes that are affected by the change (i.e. the set of
parents changes). Unfortunately, the factorization presented
in Equation 1 does not apply to probability distributions
that are faithful to mixed graphs. This happens because
variables connected with bi-directed edges (confounded)
are no longer independent given their parents in the graph.
Thus, SMCMs and MAGs do not admit a simple
factorization, where each node has a single contribution to the
likelihood. However, the joint probability of a set of variables
V according to an ADMG G can be factorized based on
sets of variables, called the c-components [Tian and Pearl,
2003], or distric
Parametrizations of the set of distributions obeying the
conditional independence relations given by an ADMG
are available for multivariate discrete [Richardson, 2009]
and multivariate Gaussian distributions [Ri
Let V = fVi : i = 1; : : : ; V g be a random vector of V variables following a multivariate normal distribution N (O; ) with positive definite covariance matrix . Let G be a MAG. Then graph G defines a system of linear equations: ij Vj + i; i 2 f1; : : : ; V g (2) Vi =
X j2P aG(Vi) Let B(G) be the collection of all real V V matrices B = ( ij ) such that (i) ij = 0 when j ! i is not in G, and (ii) (I B) is invertible. Let (G) be all the V V matrices = (!ij ) such that (i) is positive definite and (ii) !ij = 0 if j $ i is not in G.
Then the system of linear equations (2) can be written as V = BV + , and for B 2 B(G); Cov( ) = 2 (G) it has a unique solution that is a multivariate normal vector with covariance matrix = (I B) 1 (I B) T , where the superscript T denotes the transpose inverse. The family of distributions with covariance matrix in the set f = (I B) 1 (I B) T g is called the normal linear model associated with G (symb N(G)). For MAGs, the normal linear model is everywhere identifiable. Let D be a V N matrix of observations for the variables V . Then the empirical covariance matrix is defined as 1 n S =
DDT : For N V + 1, S is almost surely positive definite. For a MAG G, the log likelihood of the model is
N
1 N lG (B; jS) = ln(j2
j N 2 tr[(I
B)T
where lG ( ^ jG) is the likelihood of the graph G with the
MLE parameters B^; ^ . BIC is an asymptotically
correct criterion for selecting among Gaussian ancestral graph
models [Ri
A simple greedy strategy starts from a MAG G with a score sc and then checks the local neighborhood (i.e. the graphs that stem from the current graph after making a single edge change) for the lowest-scoring network. The algorithm continues this “hill-climbing” until no single edge change reduces the score.
Algorithm 1 begins with the empty graph, where each node is a component. At every subsequent step, every possible edge change is considered: For absent edges, the possible actions are addLeft, addRight, addBidirected. For directed edges, the possible actions are reverse, orientBidirected, remove. For bi-directed edges the possible actions are oriend end end end until no action reduces curScore; sc curScore ;
Maximum likelihood estimates B^; ^ that maximize (3)
can be found using the residual iterative conditional
fitting (RICF) algorith
Based on the factorization of
N2 X jCkjln(2 ) + ln Q j Gk j
2 j2P aGk kj tr[ Gk1SGk
jP aG (Ck) n fCkgj] ; where the Gk is the graph consisting only of nodes in Ck [ P aG (Ck) without any edges among variables in P aG (Ck)n Ck, and the subscript Gk denotes the restriction of a matrix to the rows and columns participating in Gk. k2j denotes the diagonal entry of Gk corresponding to parent node k. The log likelihood is now a sum of c-component scores. The scoring function is typically the negative log likelihood entLeft, orientRight, remove.
Score decomposition described in Equation 4 is used to avoid re-fitting the entire MAG. Instead, only the likelihood of the c-components affected by the change need to be reestimated. Algorithm 1 describes a simple greedy search strategy for learning MAG structure.
Only actions that do not create directed or almost directed cycles are attempted. To efficiently check for cycle creation, a matrix of ancestral relationships1 of the current MAG is cached. Edge removals can never create directed cycles. Using the cached ancestral matrix, it is straightforward to check whether the addition of a directed edge will create a directed cycle, or if the addition of a bi-directed edge will create an almost directed cycle. Almost directed cycles can also be created when adding directed edges: For each edge X $ Y , adding edge J ! I will create a semidirected cycle if I is an ancestor if X (Y ) and Y (X ) is an ancestor of J .
At the end of each iteration, the matrix of ancestral
relationships is updated. If a previously missing edge is added,
the update takes O(V 2) time. If an edge is removed, the
matrix is recomputed using Warshall’ s algorithm for
tran
The c-components are only updated in case a bi-directed edge is added or altered in any way. When adding a bidirected edge, the corresponding c-components of the endpoints are merged if separate. When an existing bi-directed edge is removed (completely or becomes directed), the corresponding c-component Ck is divided in the new connected components. The scores of the affected components are recomputed using new RICF estimates. In any other case, the c-components remain the same, and the score of the c-component whose corresponding graph Gk is affected by the change is recomputed. This procedure is described in Algorithm 2.
When no single-edge change improves the current score, the algorithm terminates and the current network is returned. Greedy hill-climbing procedures can be stuck in local optima (minima). To tackle this problem, they are often augmented with meta-heuristics such as random restarts, TABU lists or simulated annealing. For the scope of this work we use no such heuristic. In preliminary experiments, however, we found that augmenting Algorithm 1 with a TABU heuristic did not significantly improve performance.
1Changing edge orientation is equivalent to a removing the edge and then adding it re-oriented. To test for possible cycles efficiently, a matrix of all the non-trivial ancestral relationships (more than one variable in the path) is also cached. Reversing an edge X ! Y creates a directed cycle only if X is a non-trivial ancestor of Y .
Several constraint-based algorithms exist for learning a
Markov equivalence class of MAGs from an observational
data set: FCI [S
Nowzohour et al. [2015] propose a greedy search with random restarts for learning “bow-free” ADMGs from data and introduce the score decomposition showed in Equation 4. Since Markov equivalence for ADMGs that are not MAGs has not yet been characterized, they use a greedy strategy for obtaining the empirical Markov equivalence class, based on score similarity. The authors use the estimated ADMGs to compute causal effects and show promising results. However, since they do not necessarily find maximal ancestral graphs, they do not compare against constraint-based methods or evaluate the accuracy of the learnt graphs.
Marginalizing out variables from causal DAGs results in
some additional equality constraints that are not
conditional independencies. Nested Markov models Shpi
We compared the performance of Algorithm 1 against FCI and CFCI in simulated data. We simulated 100 random DAGs over 10, 20 and 50 variables. To control the sparseness of the DAGs, we set the maximum parents of each node. We present results for sparse networks, where each variable was allowed to have up to 3 parents, and denser networks where each variables was allowed to have up to 5 parents. For each DAG, 10% of the variables were marginalized (1, 2 and 5 variables respectively). The ground truth PAG PGT was then created for each marginal DAG.
Data sets with 100, 1000 and 5000 samples were simulated for each DAG and random parameters with absolute values in f0:1; 0:9g. The corresponding marginal data sets were input in Algorithm 1, FCI and CFCI. FCI and CFCI were run with a significance threshold of 0.05 and a maximum conditioning size 5. Algorithm 1 outputs a MAG. To compare the outputs of the algorithms, the corresponding PAG was created for each MAG output. We use PF CI ; PCF CI and PGS to denote the outputs of FCI, CFCI and Algorithm 1, respectively.
Summarizing PAG differences is not trivial, and many different approaches are used in the literature. As a general metric of how different two PAGs are, we use the structural hamming distance (shd) for PAGs, defined as follows: Let P^ be the output PAG and P be the ground truth PAG. For each change (edge addition, edge deletion, change arrowhead, change tail) required to transform P^ into P, shd is increased by 1. We also use precision and recall, as described in Tillman and Spirtes [2011]: Precision is defined as the number of edges in the output PAG with the correct orientations, divided by the number of edges in the output PAG. Recall is defined as the number of edges in the output PAG with correct orientations, divided by the number of edges in the ground truth PAG. These metrics are very conservative, since they penalize even small differences. For example, an edge that is ! in the ground truth but in the output PAG will be classified as a false positive. Greedy search has larger structural hamming distances than FCI and CFCI. More specifically, out of 900 cases (over all variable and sample sizes), FCI outperforms GSMAG in 657 cases for sparse networks and in 715 cases for dense networks, while CFCI outperforms GSMAG in 682 cases for sparse networks and in 710 cases for dense networks. In terms of precision, CFCI is again the best of the three (outperforms GSMAG in 601 and 659 out of 900 cases for sparse and dense networks, respectively). FCI has the poorest precision: it outperforms GSMAG in 318 and 221 cases for sparse and dense networks, respectively). Finally, GSMAG has the best recall out of all algorithms, with CFCI being second. Specifically, terms of recall, FCI outperforms GSMAG in 139 and 83 cases, while CFCI outperforms GSMAG in 357 and 249 cases for sparse networks and dense networks, respectively. Naturally, greedy search is much slower than both conservative and plain FCI. Intriguingly, GSMAG’s performance declines for the largest attempted sample size (5000 samples), particularly for larger networks. This happens because greedy search tends to include many false positive edges. It is possible that this is related to the naive greedy search, and could be improved by augmenting some kind of heuristic for escaping local minima, or by adjusting the scoring criterion. Figure 6 shows the score of the output MAG for Algorithm 1 and the ground truth MAG. To compare also with FCI, we used the method presented in Zhang [2008] to obtain a MAG from PF CI . Notice that this cannot be applied to the output of CFCI, since it is not a complete PAG (due to unfaithful triplets). Greedy search typically scores closer to the ground truth, particularly for denser networks. 7
We present an implementation of a greedy search algorithm for learning MAGs from observations, and compare it to FCI and CFCI. To the best of our knowledge, this is the first comparison of score-based and constraint-based search in the presence of confounders.
The algorithm uses the decomposition presented in Nowzohour et al. [2015] for bow-free SMCMs. Compared to SMCMs, MAGs are less expressive in terms of causal statements. However, since they have no almost directed cycles, fitting procedures for obtaining maximum likelihood estimates always converge. Semi-Markov causal models that are Markov equivalent to the output MAG could be identified as a post-processing step.
Heuristic procedures for escaping local minima could also be explored, to improve the performance of GSMAG. Algorithm efficiency could also possibly be improved by updating without recomputing inverse matrices and where applicable.
Other interesting directions include taking weighted averages for specific PAG features, or using both constraintbased and score-based techniques for hybrid learning. Greedy search in the space of PAGs instead of MAGs could also be explored, since a transformational characterization for Markov equivalent MAGs exists [Zhang and Spirtes, 2005].
This work was funded by European Research Council (ERC) and is part of the CAUSALPATH - Next Generation Causal Analysis project, No 617393.