Bayesian Probabilities for Constraint-Based Causal Discovery

1 Introduction: Robust Causal Discovery

In real-world systems interactions between a set of variables V are often modeled in the form of a causal DAG (directed ayclic graph) GC. A directed path from A to B in such a graph GC indicates a causal relation A ⇒ B in the system, whereas an edge A→B in GC indicates a direct causal link.

The causal Markov and faithfulness assumptions link the structure of the graph GC to observed probabilistic in/dependencies through d-separation [Pearl, 2000], which forms the basis behind most existing causal discovery procedures. Together, they imply that the causal DAG GC is also minimal, in the sense that no proper subgraph can satisfy both assumptions and produce the same probability distribution [Zhang and Spirtes, 2008].

If some of the variables in the causal DAG are hidden then the independence relations between the observed variables may be represented in the form of an ancestral graph (AG) [Richardson and Spirtes, 2002]; intuitively similar to a DAG

∗The paper on which this extended abstract is based was the recipient of the best paper award of the 2012 Conference on Uncertainty in Artificial Intelligence (UAI) [Claassen and Heskes, 2012].

except that it can also contain bi-directed arcs X ↔ Y (unobserved common cause between X and Y ) and undirected edges X--Y (selection effects, ignored in this article).

Different graphs can produce the same set of independencies: the equivalence class [G] is the set of all graphs that are indistinguishable from G in terms of implied independencies. The invariant features, common to all members of [G], can be represented in the form of a partial ancestral graph (PAG) P: an ancestral graph where circle marks '◦' on edges indicate 'unknown tail or arrowhead'. A complete PAG encodes all identifiable, present or absent causal relations; see e.g. [Zhang, 2008] for details on how to read PAGs.

With this in mind, the task of a causal discovery algorithm is to find as many invariant features of the equivalence class over observed variables from a given data set as possible.

Causal discovery paradigms

The so-called constraint-based algorithms search for conditional independencies X⊥⊥Y |Z in order to eliminate direct causal links X - Y from the causal PAG P. An efficient search strategy can uncover the entire skeleton, after which a number of orientation rules are executed to find invariant tails '-' and arrowheads '>' on edges. Members of this group include the IC-algorithm [Pearl and Verma, 1991], PC/FCI [Spirtes et al., 2000], and many others; see e.g. [Glymour et al., 2004; Kalisch et al., 2011]. Of these, the FCI algorithm in conjunction with the orientation rules in [Zhang, 2008] is sound and complete in the large-sample limit when hidden common causes and/or selection bias may be present.

They tend to output a crisp and clear causal model. The downside is that for realistic, finite data sets they give little indication of which parts of the network are stable (reliable), and which are not: if unchecked, even one erroneous borderline independence decision may lead to multiple incorrect orientations [Spirtes, 2010]. To tackle this lack of robustness, [Ramsey et al., 2006] proposed a conservative approach involving explicit validation of certain orientation rules. The resulting Conservative PC algorithm is indeed more robust, but also significantly less informative than vanilla PC.

The score-based algorithms build on the implied minimality of the causal graph. They define a scoring criterion S(G,D), often corresponding to a (Bayesian) likelihood, that measures how well a Bayesian network with structure G fits the observed data D, while preferring simpler networks, with fewer free parameters, over more complex ones. For DAG structures with either discrete or Gaussian variables closed form solutions exist that can be computed efficiently from data. These are employed in algorithms such as K2 [Cooper and Herskovits, 1992] and the Greedy Equivalence Search (GES) [Chickering, 2002] to search for an optimal structure.

Score-based procedures can output a range of high-scoring models. Multiple alternatives are arguably less straightforward to interpret, but it does allow for a measured interpretation of the reliability of inferred causal relations, and is less susceptible to incorrect categorical decisions [Heckerman et al., 1999]. The main drawback is the need to rely on the causal sufficiency assumption (no latent confounders).

2 The Best of BothWorlds

The strength of constraint-based algorithms lies in the ability to handle data from arbitrary causal DAGs and turn it into clear, unambiguous causal output. The strength of Bayesian score-based approaches lies in the robustness and implicit confidence measure that a likelihood-weighted combination of multiple models can bring.

▶ Our idea is to improve on both methods by using a Bayesian approach to estimate the reliability of different constraints, and use this to decide if, when, and how that information should be used.

Instead of classifying pieces of information as 'true' or not, we want to rank and process independence constraints according to a principled confidence measure, and build up a global causal model starting from the most reliable information down to a discretionary minimum confidence level.

For that we use a recently developed method to break up the causal inference process into a series of modular steps that can be executed in arbitrary order. It works by translating observed in/dependence constraints into basic logical causal statements L via:

where square brackets indicate a minimal set of nodes.

The logical causal statement L : (Z ⇒ X) ∨ (Z ⇒ Y ) states that there is either a causal relation from Z to X, or from Z to Y , or both. Similarly, the logical causal statement L : (Z ⇏ X) states that there is no directed path from Z to X in the underlying causal DAG GC; see [Claassen and Heskes, 2011] for details. Subsequent statements follow from straightforward deduction on the causal properties transitivity and acyclicity.

Crucial in this connection is that these logical causal statements are implied directly by the structure of the underlying causal graph GC, and so in turn by the induced ancestral graph over the observed variables. As a result, we can obtain probability estimates on logical causal statements L ∈ L by summing the normalized posterior likelihoods of all structures G that entail L through d-separation:

For the likelihood estimates p(D|G) on possible DAG structures we employ the well-known Bayesian Dirichlet (BD) metric for discrete variables [Heckerman et al., 1995].

with n the number of variables, ri the multiplicity of variable Xi, qi the number of possible instantiations of the parents of Xi in G, Nijk the number of cases in data set D in which variable Xi has the value ri(k) while its parents are instantiated as qi(j), and with N ′ ij = ∑ri k=1 N ′ ijk pseudocounts for a Dirichlet prior over the free parameters.

For the prior over structures p(G) we can opt for a uniform distribution over all possible graphs, or even include additional background information. To obtain proper probability estimates we still need to normalize eq.(1) by dividing through the sum of all contributions from all graphs. It is well-known that the number of possible graphs increases super-exponentially with the number of nodes. However, eq.(1) equally applies to graphs over arbitrary subsets X ⊂ V of the observed variables. For sparse graphs we can limit the search to (small) subsets of max. size K ≪ |V| in order to keep the computations feasible, without losing much information on the possible logical causal statements that can be inferred. It implies that we can employ an efficient search strategy over increasing subsets X of variables, similar to PC/FCI. It does mean that there are now multiple ways (different subsets of variables) to obtain a given logical causal statement L, but it can be shown that it is sufficient to only keep track of the maximum probability estimates obtained so far. The resulting p(L|D(X)) form a conservative estimate of the theoretical optimal probability that could be obtained from the entire data set D.

When considering graphs over subsets of variables, we still need to account for the fact that the minimal DAG over an arbitrary subset X ⊂ V may be unfaithful (a 'uDAG') to the underlying causal structure. For example in Figure 1(a) the corresponding ancestral graph contains invariant bi-directed edges W ↔ Z and V ↔ T from hidden common causes that cannot be accommodated in a DAG. As a result, apparent direct causal links such as X → Z and V → T appear in the minimal uDAG in (b).

To avoid drawing incorrect conclusions we need to rely on a modified d-separation inference rule:

Lemma 1. Let G be a uDAG for some distribution p(X). Let GX‖Y be the graph obtained by eliminating the edge X - Y from G (if present). Then, if X ⊥⊥GX‖Y Y | Z then:

In words: independence from d-separation remains valid, but the identifiable dependencies are restricted. The rule can be extended to indirect dependencies in a uDAG G by showing that if π = 〈X, .., Y 〉 is the only unblocked path from X to Y given Z in G, then X ⊥⊥⧸p Y |Z; see [Bouckaert, 1995; Claassen and Heskes, 2012] for details. From this we build a (pre-computed) mapping G → L from each possibly unfaithful uDAG G to all valid logical causal statements L.

Finally, as we now consider graphs over subsets of different sizes, it becomes essential to ensure that we also have a consistent prior p(G) for structures of different sizes. Perhaps surprisingly, this is not obtained by applying the same strategy at different levels: a uniform distribution over DAGs over {X, Y,Z} implies p("X ⊥⊥ Y ") = 6/25, whereas a uniform distribution over two-node DAGs implies p("X ⊥⊥ Y ") = 1/3. We obtain a consistent multi-level prior by starting from a preselected level K, and then extend to different sized structures through marginalization.

3 Implementation and results

We can now turn the results from the previous section into a working algorithm, called the Bayesian Constraint-based Causal Discovery (BCCD) algorithm, depicted below.

It starts from a data set D and available background information I, and outputs a matrix of identified causal relations MC, indicating explicitly for each pair of variables whether there is a causal relation X ⇒ Y , absence of causal relation X ⇏ Y , or 'unknown'. It also produces a graphical causal model in the form of a PAG P.

A crucial step in the algorithm is the mapping G × L from (possibly unfaithful) DAG structures to logical causal statements in line 9, which converts posterior likelihoods for structures into probability estimates for causal relations. This mapping is the same for each run, so it can be precomputed once from the rules such as in Lemma 1, and stored for use afterwards (l.1) The uDAGs G are represented as adjacency matrices. For speed and efficiency purposes, we choose to limit the structures to size K ≤ 5, which gives a list of 29,281 uDAGs at the highest level. For details about representation and rules, see [Claassen and Heskes, 2012].

The adjacency search (l.3-15), loops over subsets from neigbouring nodes, looking for identifiable causal information, while keeping track of adjacencies that can be eliminated (l.11). As the set W = {X, Y } ∪ Z can be encountered in different ways, line (7) checks if the test on that set has been performed already. A list of probability estimates p(L|D) for each logical causal statement is built up (l.10), until no more information is found.

The inference stage (l.16-21) then processes the list L in decreasing order of reliability, until the threshold is reached.

Algorithm 1 Bayesian Constraint-based Causal Discovery

In : database D over variables V, backgr.info I Out: causal relations matrixMC, causal PAG P Stage 0 - Mapping 1: G ×L←Get uDAG Mapping(V,Kmax = 5) 2: p(G) ←Get Prior(I) Stage 1 - Search 3: fully connected P, empty list L, K = 0, θ = 0.5 4: while K ≤Kmax do 5: for all X ∈ V, Y ∈ Adj(X) in P do 6: for all Z ⊆ Adj(X)\Y , |Z| = K do 7: W ← Check Unprocessed(X, Y, Z) 8: ∀G ∈ GW : compute p(G|DW) 9: ∀L : p(LW|DW) ← ∑G→LW p(G|DW) 10: ∀L : p(L) ← max(p(L), p(LW|DW)) 11: P ← p("Wi -⧹⧸-Wj "|DW) > θ 12: end for 13: end for 14: K =K +1 15: end while Stage 2 - Inference 16: LC = empty 3D-matrix size |V|3 , i = 1 17: Ł←Sort Descending (Ł, p(L)) 18: while p(Li) > θ do 19: LC ←Run Causal Logic(LC, Li) 20: i ← i + 1 21: end while 22: MC ←Get Causal Matrix(LC) 23: P ←Map To PAG(P,MC)

Statements in L are added one-by-one to the matrix of logical causal statements LC (encoding identical to L, see [Claassen and Heskes, 2011]), with additional information inferred from the causal logic rules. Basic conflict resolution is achieved by not overriding information already derived from more reliable statements. The final step (l.22,23) retrieves all explicit causal relations in the form of a causal matrix MC, and maps this onto the skeleton P obtained from Stage 1 to return a graphical PAG representation.

4 Experimental Evaluation

We have tested various aspects of the BCCD algorithm in many different circumstances, and against various other methods. The principal aim in this paper is to verify the viability of the Bayesian approach. We compare our results and that of other methods from data against known ground-truth causal models. For that, we generate random causal graphs with certain predefined properties (adapted from [Melancon et al., 2000; Chung and Lu, 2002]), generate random data from this model, and marginalize out one or more hidden confounders. We looked at the impact of the number of data points, size of the models, sparseness, choices for parameter settings etc. on the performance to get a good feel for expected strengths and weaknesses in real-world situations.

It is well-known that the relative performance of different causal discovery methods can depend strongly on the performance metric and/or specific test problems used in the evaluation. Therefore, we will not claim that our method is inherently better than others based on the experimental results below, but simply note that the fact that in nearly all test cases the BCCD algorithm performed as good or better than other methods, is a clear indication of its viability and potential.

First we implemented the BCCD approach as a simple independence test through a modified mapping G → I of structures G to implied in/dependence statements I. Figure 2 shows a typical example in the form of ROC-curves for different sized data sets, compared against a chi-squared test and a Bayesian log-odds test from [Margaritis and Bromberg, 2009], with the prior on independence as the tuning parameter for BCCD. For 'regular' conditional independence there was no significant difference (BCCD slightly ahead, more as conditioning set increases). But for minimal independencies other methods reject for both high and low decision thresholds, resulting in the looped curves in (b); BCCD has no such problem.

Figure 3 shows typical results for the BCCD algorithm itself: for a data set of 1000 records the PAG accuracy for both FCI and conservative FCI peaks around a threshold α ≈ 0.05 - lower for more records, higher for less - with conservative FCI consistently outperforming standard FCI. The BCCD algorithm peaks at a cut-off value θ ∈ [0.5, 0.7] with an accuracy that is slightly higher than the maximum for conservative FCI. The PAG accuracy tends not to vary much over this interval, making the default choice θ = 0.5 fairly safe, even though the number of invariant edge marks does increase significantly (more decisions).

Figure 4 depicts the causal accuracy as a function of the tuning parameter for the three methods. The BCCD dependency is set against (1-θ) so that going from 0 → 1 matches processing the list of statements in decreasing order of reliability. As hoped/expected: changing the decision parameter θ allows to access a range of accuracies, from a few very reliable causal relations to more but less certain indications. In contrast, the accuracy of the two FCI algorithms cannot be tuned effectively through the decision parameter α. The reason behind this is apparent from Figure 2(b): changing the decision threshold in an independence test shifts the balance between dependence and independence decisions, but it cannot identify or alter the balance in favor of more reliable decisions. We consider the fact that the BCCD can do exactly that as the most promising aspect of the Bayesian approach.

5 Discussion

The experimental results confirm that the Bayesian approach is both viable and promising: even in a basic implementation the BCCD algorithm already outperforms other state-of-theart causal discovery algorithms. It yields slightly better accuracy, and comes with an easy tuning parameter that can be used to vary from making just a few but very reliable causal statements to many possibly less certain decisions.

An interesting question is how far off from the theoretical optimumwe are: at the moment it is not clear whether we are fighting for the last few percent or if sizeable gains can still be made. Obvious improvements include scoring equivalence classes (instead of all DAGs), handle larger substructures through sampling, and to score ADMGs directly [Silva and Ghahramani, 2009; Evans and Richardson, 2010]. Finally we want to extend the method to handle continuous/mixed data without the need for discretization.