Distributional Equivalence and Structure Learning for Bow-free Acyclic Path Diagrams

We consider the problem of structure learning for bow-free acyclic path diagrams (BAPs). BAPs can be viewed as a generalization of linear Gaussian DAG models that allow for certain hidden variables. We present a first method for this problem using a greedy score-based search algorithm. We also prove some necessary and some sufficient conditions for distributional equivalence of BAPs which are used in an algorithmic ap- proach to compute (nearly) equivalent model structures. This allows us to infer lower bounds of causal effects. We also present applications to real and simulated datasets using our publicly available R-package.


Introduction
We consider learning the causal structure among a set of variables from observational data. In general, the data can be modelled with a structural equation model (SEM) over the observed and unobserved variables, which expresses each variable as a function of its direct causes and a noise term, where the noise terms are assumed to be mutually independent. The structure of the SEM can be visualized as a directed graph, with vertices representing variables and edges representing direct causal relationships. We assume the structure to be recursive (acyclic), which results in a directed acyclic graph (DAG). DAGs can be understood as models of conditional independence, and many structure learning algorithms use this to find all DAGs which are compatible with the observed conditional independencies (Spirtes et al., 1993). Often, however, not all relevant variables are observed. The resulting marginal distribution over the observed variables might still satisfy some conditional independencies, but in general these will not have a DAG representation (Richardson and Spirtes, 2002). Also, there generally are additional constraints resulting from the marginalization of some of the variables (Evans, 2016;Shpitser et al., 2014).
In this paper we consider a model class which can accommodate certain hidden variables. Specifically, we assume that the graph over the observed variables is a bow-free acyclic path diagram (BAP). This means it can have directed as well as bidirected edges (with the directed part being acyclic), where the directed edges represent direct causal effects, and the bidirected edges represent hidden confounders. The bow-freeness condition means there cannot be both a directed and a bidirected edge between the same pair of variables. The BAP can be obtained from the underlying DAG over all (hidden and observed) variables via a latent projection operation (Pearl, 2000) (if the bow-freeness condition admits this). We furthermore assume a parametrization with linear structural equations and Gaussian noise, where two noise terms are correlated only if there is a bidirected edge between the two respective nodes. In certain situations, it is beneficial to consider this restricted class of hidden variable models, as it forms a middle ground between DAG models that don't allow any hidden variables and maximal ancestral graph (MAG) models (Richardson and Spirtes, 2002) that allow arbitrarily many and general hidden variables. Such a restricted model class, if not heavily misspecified, results in a smaller distributional equivalence class, and estimation is expected to be more accurate than for more general hidden variable methods like FCI (Spirtes et al., 1993), RFCI (Colombo et al., 2012), or FCI+ (Claassen et al., 2013).
The goal of this paper is structure learning with BAPs, that is, finding the set of BAPs that best explains some observational data. Just like in other models, there is typically an equivalence class of BAPs that are statistically indistinguishable, so a meaningful structure search result should represent this equivalence class. We propose a penalized likelihood score that is greedily optimized and a heuristic algorithm (supported by some theoretical results) for finding equivalent models once an optimum is found. This method is the first of its kind for BAP models.
(c) Figure 1: (a) DAG with hidden variables H 1 , H 2 , H 3 , (b) resulting BAP over the observed variables X 1 , . . . , X 4 with annotated edge weights, and (c) resulting graph if X 3 is also not observed, which is not a BAP.

Example of a BAP
Consider the DAG in Figure 1a, where we observe variables X 1 , . . . , X 4 , but do not observe H 1 , H 2 , H 3 . The only (conditional) independency over the observed variables is X 1 ⊥ ⊥ X 3 | X 2 , which is also represented in the corresponding BAP in Figure 1b. The parametrization of this BAP would be with (  Hence the model parameters in this case are B 21 , B 32 , B 43 , Ω 11 , Ω 22 , Ω 33 , Ω 44 , and Ω 24 . An example of an acyclic path diagram that is not bow-free is shown in Figure 1c.

Challenges
The main challenge, like with all structure search problems in graphical modelling, is the vastness of the model space. The number of BAPs grows super-exponentially. Hence, as is the case for DAGs, exhaustively scoring all BAPs and finding the global score optimum is very challenging. For DAGs, Silander and Myllymäki (2006) proposed a surprisingly simple algorithm whose runtime is exponential in the number of nodes and which is feasible for problems with up to about 30 nodes. However, extending their idea to BAPs is not straightforward, and we aim to deal with settings where the number of nodes can be significantly larger.
Another major challenge, specifically for our setting, is the fact that a graphical characterization of the (distributional) equivalence classes for BAP models is not yet known. In the DAG case, for example, it is known that models are equivalent if and only if they share the same skeleton and v-structures (Verma and Pearl, 1991). A similar result is not known for BAPs (or the more general acyclic directed mixed graphs). This makes it hard to traverse the search space efficiently, since one cannot search over the equivalence classes (like the greedy equivalence search for DAGs, see Chickering (2002)). It also makes it difficult to evaluate simulation results, since the graphs corresponding to the ground truth and the optimal solution may be distinct and yet still represent the same model.

Contributions
We provide the first structure learning algorithm for BAPs. It is a scorebased algorithm and uses greedy hill climbing to optimize a penalized likelihood score. We are able to achieve a significant computational speedup by decomposing the score over the bidirected connected components of the graph and caching the score of each component. To mitigate the problem of local optima, we perform many random restarts of the greedy search.
We propose to approximate the distributional equivalence class of a BAP by using a greedy strategy for likelihood scoring. If two BAPs are similar with respect to their penalized likelihoods within a tolerance, they should be treated as statistically indistinguishable and hence as belonging to the same class of (nearly) equivalent BAPs. Based on such greedily computed (near) equivalence classes, we can then infer bounds of total causal effects, in the spirit of Maathuis et al. (2009Maathuis et al. ( , 2010. We present some theoretical results towards equivalence properties in BAP models, some of which generalize to acyclic path diagrams. In particular, we prove some necessary and some sufficient conditions for BAP equivalence. Furthermore, we present a Markov Chain Monte Carlo method for uniformly sampling BAPs based on ideas from Kuipers and Moffa (2015).
We obtain promising results on simulated data sets despite the challenges listed above. Comparing the highest-scoring BAPs and DAGs on real datasets exemplifies the more conservative nature of BAP models.

Related Work
There are two main research communities that intersect at this topic. On the one side there are the path diagram models, going back to Wright (1934) and then being mainly developed in the behavioral sciences (Jöreskog, 1970;Duncan, 1975;Glymour and Scheines, 1986;Jöreskog, 2001). In this setting a model for the edge functions is assumed, usually a parametric model with linear edge functions and Gaussian noise. In a very general formulation, the graph over the observed variables is assumed to be an acyclic directed mixed graph (ADMG), which can have bows. While in general the parameters for these models are not identified, Drton et al. (2011) give necessary and sufficient conditions for global identifiability. Complete necessary and sufficient conditions for the more useful almost everwhere identifiability remain unknown (however, see Foygel et al. (2012) for some necessary and some sufficient conditions). BAP models are a useful subclass, since they are almost everywhere identified (Brito and Pearl, 2002). Drton et al. (2009) provided an algorithm, called residual iterative conditional fitting (RICF), for maximum likelihood estimation of the parameters for a given BAP.
On the other side there are the non-parametric hidden variable models, which are defined as marginalized DAG models (Pearl, 2000) 1 . The marginalized distributions are constrained by conditional independencies, as well as additional equality and inequality constraints (Evans, 2016). When just modelling the conditional independence constraints, the class of maximal ancestral graphs (MAGs) is sufficient (Richardson and Spirtes, 2002). Shpitser et al. (2014) have proposed the nested Markov model using ADMGs to also include the additional equality constraints. Finally, mDAGs were introduced to model all resulting constraints (Evans, 2016). In general BAPs induce independence constraints and also Verma constraints (Richardson and Spirtes, 2002, Sections 7.3 and 8), as well as other restrictions that do not apply in the non-parametric case. The BAP in Figure 1b, for example, implies a Verma constraint. Gaussian BAPs are also 'maximal', in the sense that every missing edge induces a constraint. In the non-parametric case, with each additional layer of constraints learning the graphical structure from data becomes more complicated, but at the same time more available information is utilized and a possibly more detailed structure can be learned. In the Gaussian case, however, all models are parameteric, and fitting BAPs that do not correspond to conditional independence models is essentially no different to fitting those that do. At the graphical level the search is perhaps easier, since we do not have to place the restriction of ancestrality on the structure of the graph. However, unlike for MAGs, the equivalence class of BAPs is not known, which means that one may end up fitting the same model multiple times in the form of different graphs. Furthermore, BAPs are easier to interpret as hidden variable models. This can be seen when comparing the BAP in Figure 1b with the corresponding MAG. The latter would have an additional edge between X 1 and X 4 since there is no (con-ditional) independency of these two variables. As can be verified, the BAP and the MAG in this example are not distributionally equivalent, since the former encodes additional non-independence constraints.
Structure search for MAGs can be done with the FCI (Spirtes et al., 1993), RFCI (Maathuis et al., 2009), or FCI+ (Claassen et al., 2013 algorithms. Silva and Ghahramani (2006) propose a fully Bayesian method for structure search in linear Gaussian ADMGs, sampling from the posterior distribution using an MCMC approach. Shpitser et al. (2012) employ a greedy approach to optimize a penalized likelihood over ADMGs for discrete parametrizations.

Outline of this Paper
In Section 2 we give an in-depth overview of the model and its estimation from data, as well as some distributional equivalence properties. In Section 3 we present the details of our greedy algorithm with various computational speedups. In Section 4 we present empirical results on simulated and real datasets. All proofs as well as further theoretical results and justifications can be found in the Appendix.

Graph Terminology
Let X 1 , . . . , X d be a set of random variables and V = {1, . . . , d} be their index set. The elements of V are also called nodes or vertices. A mixed graph or path diagram G on V is an ordered tuple we say there is a directed edge from i to j and write i → j ∈ G. If (i, j) ∈ E B , we must also have (j, i) ∈ E B , and we say there is a bidirected edge between i and j and write i ↔ j ∈ G. The set pa G (i) := {j | j → i ∈ G} is called the parents of i. This definition extends to sets of nodes S in the obvious way: pa G (S) := i∈S pa G (i). The in-degree of i is the number of arrowheads at i.
is called a subgraph of G, and we write G ⊆ G. The induced subgraph G W for some vertex set W ⊂ V is the restriction of V to vertices W . If G ⊆ G but G = G, we call G a strict subgraph of G and write G ⊂ G. The skeleton of G is the undirected graph over the same node set V and with edges i − j if and only if i → j ∈ G or i ↔ j ∈ G (or both).
A path π between i and j is an ordered tuple of (not necessarily distinct) nodes π = (v 0 = i, . . . , v l = j) such that there is an edge between v k and v k+1 for all k = 0, . . . , l − 1. If the nodes are distinct, the path is called non-overlapping (note that in the literature a path is mostly defined as nonoverlapping). The length of π is the number of edges λ(π) = l. If π consists only of directed edges pointing in the direction of j, it is called a directed path from i to j. A node j on a non-overlapping path π is called a collider if π contains a non-overlapping subpath (i, j, k) with two arrowheads into j 2 . Otherwise j is called a non-collider on the path. If j is a collider on a non-overlapping path (i, j, k), we call (i, j, k) a collider triple. Moreover, if (i, j, k) is a collider triple and i and k are not adjacent in the graph, then (i, j, k) is called a v-structure. A path without colliders is called a trek.
Let A, B, C ⊂ V be three disjoint sets of nodes. The set an(C) := C ∪ {i ∈ V | there exists a directed path from i to c for some c ∈ C} is called the ancestors of C. A non-overlapping path π from a ∈ A to b ∈ B is said to be m-connecting given C if every non-collider on π is not in C and every collider on π is in an(C). If there are no such paths, A and B are m-separated given C, and we write A ⊥ ⊥ m B | C. We use a similar notation for denoting conditional independence of the corresponding set of variables X A and X B given X C : A graph G is called cyclic if there are at least two distinct nodes i and j such that there are directed paths both from i to j and from j to i. Otherwise G is called acyclic or recursive. An acyclic path diagram is also called an acyclic directed mixed graph (ADMG). An acyclic path diagram having at most one edge between each pair of nodes is called a bow-free 3 acyclic path diagram (BAP). An ADMG without any bidirected edges is called a directed acyclic graph (DAG).

The Model
A linear structural equation model (SEM) M is a set of linear equations involving the variables X = (X 1 , . . . , X d ) T and some error terms = ( 1 , . . . , d ) T : where B is a real matrix, cov( ) = Ω is a positive semi-definite matrix, and we assume that all variables X have been normalized to mean zero. M has an associated graph G that reflects the structure of B and Ω. For every non-zero entry B ij there is a directed edge from j to i, and for every nonzero entry Ω ij there is a bidirected edge between i and j. Thus we can also write (1) as: with cov( i , j ) = Ω ij for all i, j ∈ V .
Our model is a special type of SEM, which we refer to as Gaussian BAP model 4 . In particular, we make the following assumptions: (A1) The errors follow a multivariate Normal distribution N (0, Ω).
(A2) The associated graph G is a BAP.
Assumption (A1) is not strictly needed for our equivalence results, but we rely on it for fitting the models in practice using the RICF method of Drton et al. (2009).
Often M is specified via its graph G, and we are interested to find parameters θ G compatible with G. We thus define the parameter spaces for the edge weight matrices B (directed edges) and Ω (bidirected edges) for a given BAP G as Ω is symmetric positive semi-definite} and the combined parameter space as The covariance matrix for X is given by: where φ : Θ G → S G maps parameters to covariance matrices, and S G := φ(Θ G ) is the set of covariance matrices compatible with G. Note that φ(θ) exists since G is acyclic by (A2) and therefore I − B is invertible. We assume that the variables are normalized to have variance 1, that is, we are interested in the subsetS G ⊂ S G , whereS G = {Σ ∈ S G | Σ ii = 1 for all i = 1, . . . , d}, and its preimage under φ, One of the main motivations of working with BAP models is parameter identifiability. This is defined below: Brito and Pearl (2002) show that for any BAP G, the set of normalized non-identifiable parameters has measure zero.
The causal interpretation of BAPs is the following. A directed edge from X to Y represents a direct causal effect of X on Y . A bidirected edge between X and Y represents a hidden variable which is a cause of both X and Y , see also Figure 1. In practice, one is often interested in predicting the effect of an intervention at X j on another variable X i . This is called the total causal effect of X j on X i and can be defined as Pearl, 2000). For linear Gaussian path diagrams such as in (1) or (2), this is a constant quantity given by

Penalized Maximum Likelihood
Consider a BAP G. A first objective is to estimate the parameters θ G from n i.i.d. samples of model (2), denoted by {x (s) i } (i = 1, . . . , d and s = 1, . . . , n). This can be done by maximum likelihood estimation using the RICF method of Drton et al. (2009). Given the Gaussianity assumption (A1) and the covariance formula (3), one can express the log-likelihood for some given parameters θ G and the sample covariance matrix S as: where Σ G = φ(θ G ) is the covariance matrix implied by parameters θ G , see for example Mardia et al. (1979, (4.1.9)). However, due to the structural constraints on B and Ω it is not straightforward to maximize this for θ G . RICF is an iterative method to do so, yielding the maximum likelihood estimate:θ We now extend this to the scenario where the graph G is also unknown, using a regularized likelihood score with a BIC-like penalty term that increases with the number of edges. Concretely, we use the following score for a given BAP G: We have scaled the log-likelihood and penalty with 1/n so that the score is expected to be O(1) as n increases. Compared with the usual BIC penalty, we chose our penalty to be twice as large, since this led to better performance in simulations studies 5 . The number of nodes is typically fixed, so does not matter for comparing graphs over the same vertex set. We included it to make explicit the penalization of the model parameters (which correspond to nodes and edges).
In our search for the true causal graph G, we assume that if Σ ∈S H for any other graph H, thenS H ⊇S G , that is H represents a strict supermodel of G. This rules out the possibility that 'by chance' we land on a distribution contained in a submodel, and is a minimal requirement for causal learning. The set of matrices that violate the requirement has measure zero within the modelS G (assuming entries in B and Ω are generated according to a positive joint density with respect to the Lebesgue measure) 6 . This requirement is analogous to the faithfulness assumption of Spirtes et al. (1993), though faithfulness applies separately to individual conditional independence constraints rather than to the entire model.

Equivalence Properties
There is an important issue when doing structure learning with graphical models: typically the maximizers of (7) will not be unique. This is a fundamental problem for most model classes and a consequence of the model being underdetermined. In general, there are sets of graphs that are statistically indistinguishable (in the sense that they can all parametrize the same joint distributions over the variables). These graphs are called distributionally equivalent. For nonparametric DAG models (without nonlinearity or non-Gaussianity constraints), for example, the distributional equivalence classes are characterized by conditional independencies and are called Markov equivalence classes. For BAPs, distributional equivalence is not completely characterized yet (see Spirtes et al. (1998) or Williams (2012) for a discussion of the linear Gaussian ADMG case), but we present some necessary and some sufficient conditions, that can be used to simplify structure search in practice. Let us first make precise the different notions of model equivalence.
Definition 2. Two BAPs G 1 , G 2 over a set of nodes V are Markov equivalent if they imply the same m-separation relationships.
This essentially means they imply the same conditional independencies, and the definition coincides with the classical notion of Markov equivalence when G 1 and G 2 are both DAGs. The following notion of distributional equivalence is stronger.
We now present some sufficient and some necessary conditions for distributional equivalence in BAP models. Note that the Gaussianity assumption (A1) is not required for these to hold. Spirtes et al. (1998) showed the following global Markov property for general linear path diagrams: if there are nodes a, b ∈ V and a possibly empty set C ⊂ V such that a ⊥ ⊥ m b | C, then the partial correlation of X a and X b given X C is zero. In addition, if such an m-separation does not hold then the partial correlation is non-zero for almost all distributions. As a direct consequence, we get the following first result:

Necessary Conditions
Lemma 1. If two BAPs G 1 , G 2 do not share the same m-separations, they are not distributionally equivalent.
Unlike for DAGs, the converse is not true, as the counterexample in Figure 2 shows. For DAGs it is trivial to show that having the same skeleton is necessary for Markov equivalence, since a missing edge between two nodes means they can be d-separated, and thus a conditional independency would have to be present in the corresponding distribution. For BAPs a missing edge does not necessarily result in an m-separation, as the counterexample in Figure 2 shows. However, the following result will allow us to improve the necessary condition of same m-separations for BAPs to the same as for DAGs.
Theorem 1. Let G and G be distributionally equivalent BAPs on vertices V . Then, for any subset W ⊆ V , the induced subgraphs G W and G W are also distributionally equivalent.
If we in particular look at the induced subgraphs of size two and three we obtain the following necessary conditions for distributional equivalence.
Corollary 1. Let G and G be distributionally equivalent BAPs. Then they have the same skeleton and v-structures.
Since m-separations are not fully determined by the skeleton and the vstructures of a graph, it is also worthwhile to look at larger subgraphs. This leads, for example, to the following result: if two graphs are distributionally equivalent and a particular path is a so-called discriminating path in both graphs, then the discriminated triple will be a collider in both or in neither (see Ali et al., 2009, Section 3.4).
The criteria given above are not complete, in the sense that there exist BAPs that are not distributionally equivalent and yet this cannot be proven by applying these results. For example, the BAPs in Figure 3 are not distributionally equivalent, which can be shown using the results of Shpitser et al. (2014). However, they both have no m-separations. A complete characterization remains an open problem.  (a) and (c) share the same m-separations (none) but are not distributionally equivalent since they have different skeletons (using Corollary 1). Figure 3: The two BAPs in (a) and (b) differ only in the direction of the X 1 , X 3 edge; both have no m-separations, and every induced subgraph leads to models which are distributionally equivalent. However, by using the results of Shpitser et al. (2014) one can show that these models are not distributionally equivalent.

Sufficient Conditions
To prove sufficient conditions, we first give a characterization of the equivalence class in terms of treks (collider-free paths) using Wright's path tracing formula (Wright, 1960). Wright's formula expresses the covariance between any two variables in a path diagram as the sum-product over the edge labels of the treks between those variables, as long as all variables are normalized to variance 1. A precise statement as well as a proof of a more general version of Wright's formula can be found in the Appendix (Theorems 3 and 4). As an example, consider the BAP in Figure 1b. There are two treks between X 2 and X 4 : X 2 → X 3 → X 4 and X 2 ↔ X 4 . Hence cov(X 2 , X 4 ) = B 32 B 43 + Ω 24 , assuming normalized parameters. Similarly we have cov(X 1 , X 4 ) = B 21 B 32 B 43 . As a consequence of Wright's formula, we can show that having the same skeleton and collider triples is sufficient for two acyclic path diagrams to be distributionally equivalent: Theorem 2. Let G 1 , G 2 be two acyclic path diagrams that have the same skeleton and collider triples. Then G 1 and G 2 are distributionally equivalent.
Considering Figure 3b, for example, this result shows that if we replace the X 1 ↔ X 2 edge with X 1 → X 2 , the resulting graph is distributionally equivalent to the original.
For DAGs, it is known that the weaker condition of having the same skeleton and v-structures is sufficient for being Markov equivalent. For BAPs this is not true, as the counterexample in Figure 2 (together with Lemma 1) shows.
We therefore have that the distributional equivalence class of a BAP G: • is contained in the set of BAPs with the same skeleton and v-structures as G and • contains the set of BAPs with the same skeleton and collider triples as G.
We know that the first relation is strict by the counterexample mentioned above and have strong evidence that the second relation is strict as well (Nowzohour, 2015, Appendix B) 7 .

Greedy Search
We aim to find the maximizer of (7) over all graphs over the node set V = {1, . . . , d}. Since exhaustive search is infeasible, we use greedy hillclimbing. Starting from some graph G 0 , this method obtains increasingly better estimates by exploring the local neighborhood of the current graph. At the end of each exploration, the highest-scoring graph is selected as the next estimate. This approach is also called greedy search and is often used for combinatorial optimization problems. Greedy search converges to a local optimum, although typically not the global one. To alleviate this we repeat it multiple times with different (random) starting points. We use the following neighborhood relation. A BAP G is in the local neighborhood of G if it differs by exactly one edge, that is, the number of edges differs by at most one, and one of the following holds: 1. G ⊂ G (edge addition), 2. G ⊂ G (edge deletion), or 3. G and G have the same skeleton (edge change).
If we only admit the first condition, the procedure is called forward search, and it is usually started with the empty graph. Instead of searching through the complete local neighborhood at each step (which can become prohibitive for large graphs), we can also select a random subset of neighbors and only consider those.
In Sections 3.1 and 3.2 we describe some adaptations of this general scheme, that are specific to the problem of BAP learning. In Section 3.3 we describe our greedy equivalence class algorithm.

Score Decomposition
Greedy search becomes much more efficient when the score separates over the nodes or parts of the nodes. For DAGs, for example, the log-likelihood can be written as a sum of components, each of which only depends on one node and its parents. Hence, when considering a neighbor of some given DAG, one only needs to update the components affected by the respective edge change. A similar property holds for BAPs. Here, however, the components are not the nodes themselves, but rather the connected components of the bidirected part of the graph (that is, the partition of V into sets of vertices that are reachable from each other by only traversing the bidirected edges). For example, in Figure 1b the bidirected connected components (sometimes also called districts) are {X 1 }, {X 2 , X 4 }, {X 3 }. This decomposition property is known (Tian, 2005;Richardson, 2009), but for completeness we give a derivation in appendix A.2. We write out the special case of the Gaussian parametrization below.
Let us write p X G for the joint density of X under the model (2), and p G for the corresponding joint density of . Let C 1 , . . . , C K be the connected components of the bidirected part of G. We separate the model G into submodels G 1 , . . . , G K of the full SEM (2), where each G k consists only of nodes in V k = C k ∪ pa(C k ) and without any edges between nodes in pa(C k ) \ C k . Then, as we show in appendix A.2, the log-likelihood of the model with joint density p X G given data D = {x (s) i } (with 1 ≤ i ≤ d and 1 ≤ s ≤ n) can be written as: ..,n ) refers to the likelihood of the X j -marginal of p X G k .
For our Gaussian parametrization, using (5), this becomes where S G k is the restriction of S to the rows and columns corresponding to C k , and σ 2 kj is the diagonal entry of Σ G k corresponding to parent node j. Note that now the log-likelihood depends on {x (s) i } and p X G only via S and Σ G 1 , . . . , Σ G K . Furthermore, the log-likelihood is now a sum of contributions from the submodels G k . This means we only need to re-compute the likelihood of the submodels that are affected by an edge change when scoring the local neighborhood. In practice, we also cache the submodel scores, that is, we assign each encountered submodel a unique hash and store the respective scores, so they can be re-used.

Uniformly Random Restarts
To restart the greedy search we need random starting points (BAPs), and it seems desirable to sample them uniformly at random 8 . Just like for DAGs, it is not straightforward to achieve this. What is often done in practice is uniform sampling of triangular (adjacency) matrices and subsequent uniform permutation of the nodes. However, this does not result in uniformly distributed graphs, since for some triangular matrices many permutations yield the same graph (the empty graph is an extreme example). The consequence is a shift of weight to more symmetric graphs, that are invariant under some permutations of their adjacency matrices. A simple example with BAPs for d = 3 is shown in Figure 4. One way around this is to design a random process with graphs as states and a uniform limiting distribution. A corresponding Markov Chain Monte Carlo (MCMC) approach is described for example in Melançon et al. (2001) for the case of DAGs. See also Kuipers and Moffa (2015) for an overview of different sampling schemes.
We adapted the MCMC algorithm for BAPs as described below.
be a position sampled uniformly at random and let σ ∈ Bernoulli(0.5) be a single Bernoulli draw. We then form G k+1 = (V, E D , E B ) by applying the following rules.

If there is an edge at
It is easy to check that the resulting transition matrix is irreducible and symmetric (see Appendix A.3) and hence the Markov chain has a (unique) uniform stationary distribution. Thus, starting from any graph, after an initial burn-in period, the distribution of the visited states will be approximately uniform over the set of all BAPs. In practice, we start the process from the empty graph and sample after taking O(d 4 ) steps (c.f. Kuipers and Moffa (2015)).
It is straightforward to adapt this sampling scheme to a number of constraints, for example uniform sampling over all BAPs with a given maximal in-degree.

Greedy Equivalence Class Construction
We propose the following recursive algorithm to greedily estimate the distributional equivalence class EC(G) of a given BAP G with score ζ. We start by populating the empirical equivalence class EC(G) with graphs that have the same skeleton and collider triples as G, since these are guaranteed to be equivalent by Theorem 2. This is a significant computational shortcut, since these graphs do not have to be found greedily anymore. Then, starting once from each of the graphs in EC(G) found above, at each recursion level we search all edge-change neighbors of the current BAP for BAPs that have a score within of ζ (edge additions or deletions would result in non-equivalent graphs by Corollary 1). For each such BAP, we start a new recursive search until a maximum depth of d(d − 1)/2 (corresponding to the maximum number of possible edges) is reached, and always comparing against the original score ζ. Already visited states are stored and ignored. Finally, all found graphs are added to EC(G). The main tuning parameter here is , essentially specifying the threshold for numerical error, as well as statistically indistinguishable graphs. This results in conservative estimates for total causal effects using the methods discussed in Section 4.1 by also including neighboring equivalence classes, that are statistically indistiguishable from the given data.

Implementation
Our implementation is done in the statistical computing language R (R Core Team, 2015), and the code is available as an R package on github (Nowzohour, 2017). We make heavy use of the RICF implementation fitAncestralGraph 9 in the ggm package (Marchetti et al., 2015). We noted that there are sometimes convergence issues, so we adapted the implementation of RICF to include a maximal iteration limit (which we set to 10 by default).

Empirical Results
In this section we present some empirical results to show the effectiveness of our method. First, we consider a simulation setting where we can compare against the ground truth. Then we turn to a well known genomic data set, where the ground truth is unknown, but the likelihood of the fitted models can be compared against other methods.

Causal Effects Discovery on Simulated Data
To validate our method, we randomly generate ground truths, simulate data from them, and try to recover the true total causal effects from the gener-ated datasets. This procedure is repeated N = 100 times and the results are averaged. We now discuss each step in more detail.
Randomly generate a BAP G. We do this uniformly at random (for a fixed model size d = 10 and maximal in-degree α = 2). The sampling procedure is described in Section 3.2.
Randomly generate parameters θ G . We sample the directed edge labels in B independently from a standard Normal distribution. We do the same for the bidirected edge labels in Ω, and set the error variances (diagonal entries of Ω) to the respective row-sums of absolute values plus an independently sampled χ 2 (1) value 10 .

Simulate data {x
This is straightforward, since we just need to sample from a multivariate Normal distribution with mean 0 and covariance φ(θ G ). We use the function rmvnorm() from the package mvtnorm (Genz et al., 2014).

Find an estimateĜ from {x
We use greedy search with R = 100 uniformly random restarts (as outlined in Section 3), as well as one greedy forward search starting from the empty model.

Compare G andĜ.
A direct comparison of the graphs does not make sense since they could be different but in the same equivalence class. We therefore estimate the equivalence classes of both G andĜ using the greedy approach described in Section 3.3 with = 10 −10 to get EC(G) and EC(Ĝ).
Since the estimated equivalence classes are empirical, it is not straightforward to compare them. For one, they might be intersecting, but not equal (if the recursion level was set too low and they were started from different graphs for example). More relevantly, they might be entirely different, but still agree in large areas of the graph. We therefore chose to evaluate not the graph structure but the identifiability of causal effects. Often this is also more relevant in practice. Maathuis et al. (2009) developed a method (which they called IDA) to find identifiable causal effects in a multiset of DAGs. We apply the same idea in our setting. Specifically, this means we estimate the causal effects matrixÊ for each graph G ∈ EC(G) (using the estimated parametersθ G = (B ,Ω ) and (4)). We then take absolute values and take the entry-wise minima over allÊ to obtainÊ min G , the minimal absolute causal effects matrix (if an entry E ij is nonzero, there is a nonzero causal effect from X i to X j for every graph in the equivalence class). We do the same forĜ to getÊ min G . What is left is to compare the minimal absolute causal effects matrix E min G of the ground truth to the minimal absolute causal effects matrix E min G of the estimate. Thus, our target set consists of all pairs (i, j), such that (Ê min G ) ij > 0. We score the pairs according to our estimatedÊ min G values, and we report the area under the ROC curve (AUC, see Hanley and McNeil (1983)). The AUC ranges from 0 to 1, with 1 meaning perfect classification and 0.5 being equivalent to random guessing 11 . In our case, we have a separate ROC curve for each graph. The points on the curve correspond to the thresholding on the estimated absolute value of the causal effects; the k-th point shows the situation when we classify the largest k − 1 values as causal, and the rest as non-causal.
The results for 100 simulations can be seen in Figure 5; the average AUC is 0.75. While this suggests that perfect graph discovery is usually not achieved, causal effects can be identified to some extent. We also note that our simulation setting is challenging, in the sense that non-zero edge weights can be arbitrarily close to zero. The computations took 2.5 hours on an AMD Opteron 6174 processor using 20 cores.

Genomic Data
We also applied our method to a well-known genomics data set (Sachs et al., 2005), where the expression of 11 proteins in human T-cells was measured under 14 different experimental conditions (the sample size varies between 707 and 927). There are likely hidden confounders, which makes this setting suitable for hidden variable models. However, it is questionable whether the bow-freeness, linearity, and Gaussianity assumptions hold to a reasonable approximation (in fact the data seem not to be multivariate normal). Furthermore, there does not exist a ground truth network (although some individual links between pairs of proteins are reported as known in the original paper). So we abstain from comparing a "best" network with reported links in literature, but instead use this as an example for comparing highestscoring BAPs and DAGs.
To do this, we first log-transform all variables since they are heavily skewed. We then run two sets of greedy searches for each of the 14 data sets: one with BAPs and one with DAGs. We use 100 random restarts in both cases. The results can be seen in Figures 6 and 7. The computations 11 Some care has to be taken because of the fact that the cases (Ê min G )ij > 0 and (Ê min G )ji > 0 exclude each other, but we took this into account when computing the false positive rate. took 4 hours for the BAP models and 1.5 hours for the DAG models on an AMD Opteron 6174 processor using 20 cores. Note that while the BAPs and DAGs look very similar in many cases, the BAPs are more conservative in identifying causal effects. Eg for dataset 4 there is a v-structure at pip3 (with pip2 and plcg) in both the highestscoring BAP and DAG. However, by Theorem 2, this part of the BAP is equivalent to versions with different edge directions (as long as the collider is preserved). This is not the case for the DAG. Hence, in the DAG model these edges are identifiable, but this identifiability disappears in the presence of potential hidden confounders in BAPs. This exemplifies the more conservative nature of BAP models. Another example is the v-structure at pakts473 (with pka and pkc) in dataset 8.

Conclusions
We have presented a structure learning method for BAPs, which can be viewed as a generalization of Gaussian linear DAG models that allow for certain latent variables. Our method is computationally feasible and the first of its kind. The results on simulated data are promising, keeping in mind that structure learning and inferring causal effects are difficult, even for the easier case with DAGs. The main sources of errors (given the model assumptions are fulfilled) are sampling variability, finding a local optimum  Figure 7). For simplicity only one highest-scoring graph is shown per example while equivalent and equally high-scoring graphs are omitted. Note that the equivalence classes in the corresponding BAPs and DAGs are similar but some v-structures lead to identifiablity in DAGs but not in BAPs. Dataset Names: 1: cd3cd28 2: cd3cd28icam2+aktinhib 3: cd3cd28icam2+g0076 4: cd3cd28icam2+psit 5: cd3cd28icam2+u0126 6: cd3cd28icam2+ly 7: cd3cd28icam2 8: cd3cd28+aktinhib 9: cd3cd28+g0076 10: cd3cd28+psitect 11: cd3cd28+u0126 12: cd3cd28+ly 13: pma 14: b2camp only, and not knowing the equivalence classes. Local optima are a general weakness of many structure learning methods in Bayesian networks since this problem is NP-hard in general (Chickering, 1996). In our simulations, overestimating the equivalence class leads to too few causal effects, while the opposite happens if we underestimate it. On the other hand, our approach of greedily approximating the empirical equivalence class is building on the idea that some models are statistically indistinguishable, due to limited sample size and estimation error. Therefore, our approach has the advantage that it can include neighboring equivalence classes, that score almost as well, which is desirable from a statistical point of view. Our theoretical results about model equivalence go some way towards characterizing the distributional equivalence classes in BAP models and allow us to efficiently approximate them empirically.
In many applications, not all relevant variables are observed, calling for hidden variable models. While there have been structure learning methods for general hidden variable models for many years (FCI, RFCI, FCI+, see Spirtes et al. (1993); Colombo et al. (2012); Claassen et al. (2013)), causal inference based on these models is very conservative (Malinsky and Spirtes, 2017). BAP models are restricted hidden variable models, where the restriction comes from the bow-freeness constaint. As such, they form an interesting middle ground between general hidden variable models and models that do not allow any hidden variables. In particular, the bowfreeness constraint leads to improved identifiability of causal effects when compared to general hidden variable models, while being more conservative than models without hidden variables. This makes our structure learning algorithm for BAPs a useful addition to existing structure learning methods. Structure learning for a different type of restricted hidden variable models is considered in Frot et al. (2017), and it will be interesting to compare our results with this method.

A.1.1 Necessary Conditions
The following Lemma shows that the point (B, Ω) = (0, I) is non-singular for the map φ : Θ G → S G for any BAP G. The result also appears in Brito and Pearl (2002) and Drton et al. (2011).
Proof. We proceed by induction on d, the number of vertices in G. If d = 1 then the result is trivial since B = 0.
Otherwise assume without loss of generality that the last vertex d has no children. The result holds for the subgraph of the remaining vertices by the induction hypothesis and the fact that the distribution of X is defined recursively. We know that Σ is of the form and we may deduce by the induction hypothesis that the first d − 1 rows of B are zero, and the upper (d − 1) × (d − 1)-sub matrix of Ω is the identity matrix. Hence where ω T = (ω 1d , . . . , ω d−1,d ), and but note that β T ω = 0 by the bow-free assumption, so we get and hence β + ω = 0. Now note that for each j, either β dj = 0 or ω jd = 0 by the bow-free assumption; hence β + ω = 0 implies that β = ω = 0, leaving Σ = I d−1 0 0 T ω dd and hence ω dd = 1. This completes the result.
Corollary 2. Let G be a BAP. For some neighborhood U of the set of covariance matrices containing I, if Σ ∈ S G ∩ U with σ ij = 0 (for i = j), then this implies that ω ij = β ij = 0.
Note that Lemma 2 allows a direct proof of the fact that having the same skeleton is necessary for BAPs to be distributionally equivalent by looking at the tangent spaces of the models at Σ = I and showing that they are determined by the skeletons of the graphs.
In the proof of Theorem 1 we make use of the language of polynomial varieties (see Cox et al. (2007) for an overview). A variety is a set defined by the zeros of some collection of polynomials (in our case polynomials in the entries of Σ), and all SEM models are varieties.
Let G be a BAP with vertices V = W∪W where W ∩ W = ∅. Let B W G be the set of matrices B ∈ B G such that only entries corresponding to directed edges in G between vertices in W have non-zero coefficients. Similarly, let O W G be the set of Ω ∈ O G such that entries corresponding to edges outside W are zero and diagonal entries outside W are 1. Define a modelS W G as the image of the map φ applied to (B W G , O W G ). So in other words, we only manipulate parameters in G that correspond to vertices and edges in G W . The resulting model is canonically isomorphic to S G W via a simple projection, since this is the same setup as for the BAP G W , but with the matrices extended to include independent vertices in W ≡ V \ W .
Let T W be the set of covariance matrices Σ on V such that Σ W W = I and Σ W W = 0 (i.e. so that vertices outside W are completely independent).
We will show that looking at the set of covariance matrices in S G that are also in T W is essentially the same as the setS W G . Since the first set is a property of the full model, and the second set is determined by the subgraph G W , this will be enough to prove Theorem 1.
Proof of Theorem 1. What we need to prove is that S G = S G implies S G W = S G W . Consider again the variety T W defined above.
Clearly S G = S G implies T W ∩ S G = T W ∩ S G . We will show that the irreducible component of T W ∩ S G which contains Σ = I is the same asS W G ; since this last quantity is isomorphic to S G W we will prove the result.
First, note thatS W G ⊆ T W andS W G ⊆ S G , so clearlyS W G ⊆ S G ∩ T W . In addition, note that by Corollary 2, in a neighborhood of Σ = I every element of T W ∩ S G is also contained inS W G . It follows that the entire irreducible component of T W ∩ S G containing Σ = I is contained withinS W G , and therefore that the irreducible component of T W ∩ S G containing Σ = I isS W G .
We can now prove the Theorem giving necessary conditions for BAP equivalence.
Proof of Corollary 1. Let us first consider vertex pairs, i.e. W = {i, j}. By Theorem 1 we have G W being distributionally equivalent to G W . If G W = G W we would have i ⊥ ⊥ m j in one of the graphs but not the other, and using Lemma 1 this would lead to a contradiction. Hence G W = G W for any vertex pair, and hence G and G must have the same skeleton.
Let us now consider vertex triplets W = {i, j, k}, such that (without loss of generality) there is a v-structure at j in G W . Then i ⊥ ⊥ m k in G W and by the same argument as above we must have i ⊥ ⊥ m k also in G W . This is only possible if there is a v-structure at j in G W . Hence G and G must have the same v-structures.

A.1.2 Sufficient Conditions
We first make precise the definition of an important class of paths: treks. These are paths that do not contain colliders. We adopt the notation of Foygel et al. (2012). A trek τ from i to j can have one of the following forms: where v L l = i and v R r = j and in the second case v 0 = v L 0 = v R 0 . Accordingly, we define the left-hand side of τ as Left(τ ) = v L l ← · · · ← v L 0 and the righthand side of τ as Right(τ ) = v R 0 → · · · → v R r . Note that there is nothing inherently directional about a trek other that the (arbitrary) definition which end node is on the left. That is, every trek from i to j is also a trek from j to i just with the left and right sides switched. We denote the lengths of the left-and right-hand sides of a trek τ by λ L (τ ) and λ R (τ ) respectively. If τ does not contain a bidirected edge, we define its head to be H τ = v 0 . If the left-and right-hand sides of τ do not intersect (except possibly at H τ ), we call τ simple 12 . We define the following sets that will be useful later: We will usually drop the subscript if it is clear from context which is the reference graph. We now show some intermediate results that are well-known, but we prove them here for completeness nevertheless. All of these apply more generally to path diagrams (possibly cyclic and with bows). Lemma 4. Let G be a path diagram over d nodes and B ∈ B G . Then Proof. By induction on l. For l = 1 the claim follows from the definition of B G . Using the inductive hypothesis we get and the claim follows, since every directed path from j to i of length l can be decomposed into a directed path π of length l − 1 from j to some node k and the edge k → i. Proof. Since G is acyclic, there is an ordering of the nodes, such that B is strictly lower triangular and hence all its eigenvalues are zero. Furthermore, the longest directed path in G has length d − 1. Therefore the result follows from Lemma 3 and Lemma 4.
The following theorem is a version of Wright's theorem that applies to non-standardized variables. It does not require a proper parametrization (in the sense that Ω needs to be positive definite). This result is probably known to experts, but we could not find a proof in the literature.
Theorem 3. Let G be a (possibly cyclic) path diagram over d nodes, B ∈ B G , and Ω ∈ R d×d such that Ω is symmetric (but not necessarily positive definite) and Ω ij = 0 if i ↔ j is not an edge in G. Then the entries of the matrix φ = (I − B) −1 Ω(I − B) −T are given by as a shorthand for the edge contribution 13 of a trek τ given parameter matrices B and Ω. We write c(τ ; B, Ω) = c e (τ ; B, Ω) · Ω Hτ Hτ for the total contribution of τ (where we define Ω Hτ Hτ to be 1 if τ contains a bidirected edge and therefore H τ = ∅).
Using Lemma 3, we can expand φ as φ = ∞ k=0 ∞ l=0 B k Ω(B l ) T . We now first show the following intermediate result, which interprets the entries of these matrices as contributions of certain treks: for integers k ≥ 0, l ≥ 0. To see this, we expand the double matrix product and use Lemma 4 to get and (8) follows since each bracketed expression corresponds to one side of the trek from i to j via a and b (and the diagonal entries of Ω do not correspond to a trek, so they are separate). Now summing over k and l gives the following which gives the result for the diagonal entries φ ii .
For the off-diagonal entries φ ij , we can get a simpler expression involving only simple treks and the diagonal entries φ ii . Note that every trek τ can be uniquely decomposed into a simple part ξ(τ ) and a (possibly empty) non-simple part ρ(τ ) (we just split at the point, where the right-and the left-hand sides of τ first intersect). Since (dropping the parameter matrices B and Ω in our notation), we can factor out the contributions of the simple parts. Note that if the simple part ξ(τ ) contains a bidirected edge, then ρ(τ ) must be empty and Ω H ξ (τ )H ξ (τ ) = 1. Hence (9) becomes and the result follows.
The following version for standardized parameters is often quoted as Wright's theorem.
Theorem 4. Let G be a (not necessarily acyclic) path diagram over d nodes, B ∈ B G , and Ω ∈ R d×d such that Ω is symmetric (but not necessarily positive definite) and Ω ij = 0 if i ↔ j is not an edge in G. Furthermore assume that we have standardized parameters B, Ω such that (φ(B, Ω)) ii = 1 for all i. Then the off-diagonal entries of φ(B, Ω) are given by Proof. This is a direct consequence of Theorem 3.
We can now prove Theorem 2, which is a consequence of Wright's formula.
Proof of Theorem 2. Let θ G 1 ∈Θ G 1 and choose θ G 2 = (B 2 , Ω 2 ) such that their edge labels agree, that is, This is possible since G 1 and G 2 have the same skeleton: we just assign the edge labels of G 1 to G 2 , irrespective of the edge type. The diagonal entries of Ω 2 are still free-we now show that they can be used to enforce (φ(B 2 , Ω 2 )) ii = 1 for all i, which defines a linear system for the diagonal entries of Ω 2 . Let d = diag(Ω 2 ) be the vector consisting of the diagonal elements of Ω 2 , and write (10) as M d + c = 1, where M is the coefficient matrix of the linear system, and c is constant. To show that (10) always has a solution, we need to show that det(M ) = 0. Without loss of generality, assume that the nodes are topologically ordered according to G 2 (this is possible since G 2 is assumed to be acyclic), that is, there is no directed path from i to j if i > j.
Then we have H τ < i (or H τ = ∅) for all τ ∈ T ii , and using the expression for φ ii in Theorem 3 we see that M must be lower triangular with diagonal equal to 1. Thus det(M ) = 1, and we can enforce (10).
Since G 1 and G 2 share the same collider triples, the sets of simple treks between any two nodes are the same in both graphs: S ij G 1 = S ij G 2 ∀i, j. Together with Theorem 4 and the fact that the edge labels agree this shows that φ(θ G 1 ) = φ(θ G 2 ).
What is left to show is that Ω 2 is a valid covariance matrix, that is, it is positive semi-definite. By (3) and (11) we have that where Σ 1 = φ(θ G 1 ). Since Σ 1 is positive semidefinite, so is Ω 2 .
The joint density for the errors can be factorized according to the independence structure implied by Ω. Let us adopt the notation {X} I := {X i } i∈I and { } I := { i } i∈I for some index set I. Then we have { } C k ⊥ ⊥ { } C l Writing D for the observed data {x s i } (with 1 ≤ i ≤ d and 1 ≤ s ≤ n), the log-likelihood can then be written as where l(p X G k ; {x (s) j } s=1,...,n ) refers to the likelihood of the X j -marginal of p X G k .

A.3 Symmetry and Irreducibility of Markov Chain
We show that the transition matrix of the Markov Chain described in Algorithm 1 is symmetric and irreducible. For two BAPs G, G , let P (G, G ) be the probability of a single step transition from G to G .
2. Irreducibility: ∃G 1 , . . . , G n such that P (G, G 1 ) Proof. Let p be the probability of sampling one position (i, j), i.e. p = 1/(d(d − 1)). Let us first consider the case where G and G only differ by one edge addition, i.e. WLOG either In both cases we get P (G, G ) = p/2 = P (G , G). In the first case, by multiplying the probabilities along the branches of 1a and 2a in Algorithm 1 respectively, and since G has no cycles. In the second case, by multiplying the probabilities along the branches of 1a and 2b respectively. Hence symmetry holds, and irreducibility is trivially true in this case. For the general case, note that the transitions described in Algorithm 1 involve either edge additions or deletions, so if G, G do not differ by only one edge addition, we have P (G, G ) = P (G , G) = 0. Furthermore, we can always find a collection of graphs G 1 , . . . , G n , such that irreducibility holds, e.g. by successively removing edges from G until the graph is empty and then successively adding edges until we arrive at G . Then we have P (G i , G i+1 ) = p/2 > 0 for all 1 ≤ i < n by the case considered above, and the claim follows.