§ INTRODUCTION

In many fields of science, a primary interest is determining causal effects, that is, distributions $P(\+Y \given \doo(\+X),\+Z)$, where variables $\+ Y$ are observed, variables $\+ X$ are intervened upon (forced to values irrespective of their natural causes) and variables $\+ Z$ are conditioned on [Pearl, 2009]. In this paper, instead of placing various parametric restrictions based on background knowledge, we are interested in the question of identifiability: can the causal effect be uniquely determined from the distributions (data) we have and a graph representing our structural knowledge on the generating causal system.

In the most basic setting we are identifying causal effects from a single observational input distribution, corresponding to passively observed data. To solve such problems more generally than what is possible with the back-door adjustment [Spirtes et al., 1993, Pearl, 2009, Greenland et al., 1999], Pearl, 1995 introduced do-calculus, a set of three rules that together with probability theory enable the manipulation of interventional distributions. Shpitser and Pearl, 2006 and Huang and Valtorta, 2006 showed that do-calculus is complete by presenting polynomial-time algorithms whose each step can be seen as a rule of do-calculus or as an operation based on basic probability theory. The algorithms have a high practical value because the rules of do-calculus do not by themselves provide an indication on the order in which the rules should be applied. The algorithms save us from manual application of do-calculus, which is a tedious task in all but the simplest problems.

Since then many extensions of the basic identifiability problem have appeared. In identifiability using surrogate experiments [Bareinboim and Pearl, 2012], or $z$-identifiability, an experimental distribution is available in addition to the observed probability distribution. For data observed in the presence of selection bias, both algorithmic and graphical identifiability results have been derived [Bareinboim and Tian, 2015, Correa et al., 2018]. More generally, the presence of missing data necessitates the representation of the missingness mechanism, which poses additional challenges [Mohan et al., 2013, Shpitser et al., 2015, Bhattacharya et al., 2019]. Another dimension of complexity is the number of available data sources. Identification from a mixture of observational and interventional distributions that originate from multiple conceptual domains is known as transportability for which complete solutions exist in a specific setting [Bareinboim and Pearl, 2014].

While completeness has been accomplished for a number of basic identifiability problems, there are still many challenging but important extensions to the identifiability problem that have not been studied so far. Table <ref> recaps the current state of the art identifiability results; it also describes generalizations that we aim to investigate in this paper. To find solutions to the more complicated identifiability problems, we present a unified approach to the identification of observational and interventional causal queries by constructing a search algorithm that directly applies the rules of do-calculus.  We impose no restrictions on the number or type of known input distributions: we thus provide a solution to problems for which no other algorithmic solutions exist (row 8 in Table <ref>). We also extend to identifiability under missing data together with mechanisms related to selection bias and transportability (row 11 in Table <ref>).

The following introductory example does not fall under any of the previously solved special cases of Table <ref>, thus necessitating our approach. Consider human resource management analyzing the remuneration policy in a company. The graph of Figure <ref> shows the key variables. The salary ($Y$) of an employee consist of a base salary (not modeled explicitly) and a bonus ($B$). The base salary depends on education ($E$) and performance ($X$) of the employee. The level of performance is evaluated by the supervisor of the employee and it is one of the factors that affects the level of bonus. The performance depends on the education as well as the attitude ($A$) of the employee and the team. The attitude may also have a direct effect on the bonus.

We are interested in estimating the total causal effect of performance on the salary, i.e., quantifying how changes in performance affect the salary. The salary ($Y$), bonus ($B$), education ($E$) and performance ($X$) of the employees are available from the registry of the human resource department. Data on attitude ($A$), bonus ($B$) and performance ($X$) have been collected in an anonymized survey that cannot be linked to the registry on the personal level. The question of interest is the identifiability of $P(Y \given \doo(X))$ from the data sources $P(Y,B,E,X)$ and $P(A,B,X)$. Using the machinery developed in this paper, we can identify the causal effect with the formula

\[
P(Y \given \doo(X)) = \sum_{B,A}P(A)P(B \given X,A)\sum_{E} P(E)P(Y|X,B,E),
\]

where the conditional distributions can be directly determined from the available data sources (see Section <ref> for further details).

[dot5 = 00XleftX : Performance] at (0,0) ;
[dot5 = 00EbelowE : Education] at (1,-0.7) ;
[dot5 = 00YrightY : Salary] at (2,0) ;
[dot5 = 00Babove rightB : Bonus] at (1,0.7) ;
[dot5 = 00AaboveA : Attitude] at (0,0.7) ;

[->] (X) – (Y);
[->] (X) – (B);
[->] (E) – (Y);
[->] (E) – (X);
[->] (B) – (Y);
[->] (A) – (X);
[->] (A) – (B);

Graph for the example on human resource management in a company.

To combat the inherent computational complexity of the search-based approach, we derive rules and techniques that avoid unnecessary computational steps. We are able to detect trivial queries where non-identifiability can be determined directly from the inputs. We also present a search heuristic that considerably speeds up the search in cases where the effect is indeed identifiable. The approach, called , is provably sound and it retains the completeness in the cases previously proven to be solved by the rules of do-calculus. We can easily scale up to the problem sizes commonly reported in the literature. The R package dosearch [R Core Team, 2019, Tikka et al., 2020] provides an implementation of  and is available on CRAN at <https://CRAN.R-project.org/package=dosearch>.

Problem         Input     data

(Reference)      Target     (assumptions)     pattern     (complete)

1      Causal effect identifiability     $P(\+Y \given \doo(\+ X))$     $P(\+V)$     None     ID

[Shpitser and Pearl, 2006]                 (Yes)
2      Causal effect identifiability      $P(\+Y \given \doo(\+ X), \+ Z)$     $P(\+V)$     None     IDC

[Shpitser and Pearl, 2006]                 (Yes)
3      $z$-identifiability     $P(\+Y \given \doo(\+ X), \+ Z)$        $P(\+V)$, $P(\+V \setminus \+B \given \doo(\+B))$     None     zID

[Bareinboim and Pearl, 2012]         (NE, ED)         (Yes)
4     g-identifiability     $P(\+Y \given \doo(\+ X))$         $\{P(\+V \setminus \+B_i \given \doo(\+ B_i) \}$      None     gID

[Lee et al., 2019]         (ED)         (Yes)
5      Surrogate outcome      $P(\+Y \given \doo(\+ X), \+ Z)$         $\{P(\+ A_i \given \doo(\+ B_i), \+ C_i) \}$      None     TRSO

identifiability         (NE, SO)         (No)

[Tikka and Karvanen, 2019]
6      $mz$-transportability     $P(\+Y \given \doo(\+ X), \+ Z)$     $\{P(\+V \setminus (\+ B_i \cup \+ T_i) \given \doo(\+B_i), \+ T_i)\}$      None     TR$^{\textrm{mz}}$

[Bareinboim and Pearl, 2014]         (NEDD, ED)        (Yes)
7      Selection bias recoverability     $P(\+Y \given \doo(\+ X), \+ Z)$     $P(\+V \setminus S \given S)$     Selection     RC

[Bareinboim and Tian, 2015]                 (Unknown)
8     Generalized identifiability      $P(\+Y \given \doo(\+ X), \+ Z)$         $\{ P(\+ A_i \given \doo(\+ B_i), \+ C_i) \}$      None     do-search

9      Missing data recoverability     $P(\+V)$     $P(\+ V^*)$     Restricted     –

[Mohan et al., 2013]                 (Yes)

10      Missing data recoverability     $P(\+V)$     $P(\+ V^*)$     Arbitrary     –
[Bhattacharya et al., 2019]                 (Unknown)
11     Generalized identifiability      $P(\+ Y \given \doo(\+ X), \+ Z)$         $\{ P(\+ A_i^* \given \doo(\+ B_i), \+ C_i^*) \}$          Arbitrary      do-search

with missing data                   (No)

Solved and unsolved problems in causal effect identification. Bold-italic denotes the previously unsolved problems for which  can now be used. Input $P(\+ V)$ stands for passively observed joint distribution of all variables. Input $P(\+V^*)$ is the joint distribution with missing data (see Section <ref>). The variable sets present in the same distribution are disjoint.  Input $P(\+V \setminus \+ B \given \doo(\+B))$ stands for an experiment where all variables are measured and input $P(\+ A \given \doo(\+ B))$ stands for an experiment where only a subset $\+ A \subset \+V$ of the variables is measured. Notation $\{\cdot\}$ denotes a set of inputs enumerated by the index $i$. The assumptions of nested experiments (NE), entire distributions (ED) and nested experiments in different domains (NEDD) are explained in Section <ref>. Assumptions related to surrogate outcomes (SO) can be found in [Tikka and Karvanen, 2019]. Input $P(\+V \given S)$ means the joint distribution under selection bias. The last column gives the name of an algorithm that can be used to solve the problem if one exists and whether it (or a theorem when no algorithm is provided) provides a complete solution to the problem, or whether the completeness status is not known. An algorithm is complete if it returns a correct formula precisely when the target query is identifiable. Problems 1–7 are special cases of problem 8 and problems 1–10 are special cases of problem 11.

Other available software for causal effect identifiability problems are only applicable to a subset of the problems presented in Table <ref>. The causaleffect R package [Tikka and Karvanen, 2017] can only be used for problems on rows 1–3 and 5–7 of the table, providing implementations of the relevant algorithms. For the problem on row 1, the generalized adjustment criterion [Perković et al., 2015] can be applied using the dagitty package [Textor and van der Zander, 2016] or the pcalg package [Kalisch et al., 2012] in R. The generalized back-door criterion [Maathuis and Colombo, 2015] is also available in the pcalg R package. The standard back-door criterion is available in the Python package DoWhy [Sharma et al., 2019]. Algorithms implemented in causaleffect run in polynomial time and can outperform dosearch in their respective restricted problem settings especially with larger graphs. For a comprehensive performance comparison between causaleffect and various adjustment criteria, see [van der Zander et al., 2019].

The paper is structured as follows. Section <ref> formulates our general identification problem and explains the scenarios in Table <ref> and previous research in detail. Section <ref> presents the search algorithm, including the rules we use, search space reduction techniques, heuristics and theoretical properties. Section <ref> shows how the search can be extended to problems that involve missing data. Section <ref> demonstrates how the search can be used in R via the dosearch package. Efficacy of the search is assessed via simulations. Section <ref> shows a number of new problems for which we can find solutions by using the search including a real-world application. These problems include combined transportability and selection bias, multiple sources of selection bias, and causal effect identification from arbitrary (experimental) distributions. This section also includes a systematic analysis of missing data problems and case-control designs. Section <ref> discusses the merits and limitations of the approach. Section <ref> offers concluding remarks.

§ THE GENERAL CAUSAL EFFECT IDENTIFICATION PROBLEM

Our presentation is based on Structural Causal Models (SCM) and the language of directed graphs. We assume the reader to be familiar with these concepts and refer them to detailed works on these topics for extended discussion and descriptions, such as [Pearl, 2009] and [Koller and Friedman, 2009].

Following the standard set-up of do-calculus [Pearl, 1995], we assume that the causal structure can be represented by a semi-Markovian causal graph $G$ over a set of vertices $\+ V$ (see Fig <ref>(fig:intro_backdoor) for example). The directed edges correspond to direct causal relations between the variables (relative to $\+ V$); directed edges do not form any cycles. Confounding of any two observed variables in $\+ V$ by some unobserved common cause is represented by a bidirected edge between the variables. This graphical representation allows us to deal with any causal structures where some variables are unmeasured. We assume a positive distribution over the variables [Huang and Valtorta, 2006] ensuring that all considered causal effects and conditional distributions are well-defined.

In a non-parametric setting, the problem of expressing a causal quantity of interest in terms of available information has been described in various ways depending on the context. When available data are affected by selection bias or missing data, a typical goal is to “recover” a joint or marginal distribution. If data are available from multiple conceptual domains, a distribution is “transported” from the source domains, from which a combination of both observational and experimental data are available, to a target domain. The aforementioned settings can be expressed in the SCM framework by equipping the graph of the model with special vertices. However, on a fundamental level these problems are simply variations of the original identifiability problem of causal effects and as such, our goal is to represent them as a single generalized identifiability problem. Formally, identifiability can be defined as follows [Pearl, 2009, Shpitser and Pearl, 2008].

Let $\+ M$ be a set of models with a description $T$ and two objects $\phi$ and $\theta$ computable from each model. Then $\phi$ is identifiable from $\theta$ in $T$ if $\phi$ is uniquely computable from $\theta$ in any model $M \in \+ M$. In other words, all models in $\+ M$ which agree on $\theta$ also agree on $\phi$.

In the simplest case, the description $T$ refers to the graph induced by causal model, $\theta$ is the joint distribution of the observed variables $P(\+ V)$ and the query $\phi$ is a causal effect $P(Y \given \doo(X))$. On the other hand, proving non-identifiability of  $\phi$ from $\theta$ can be obtained by describing two models $M^1, M^2 \in \+ M$ such that $\theta$ is the same in $M^1$ and $M^2$, but object $\phi$ in $M^1$ is different from $\phi$ in $M^2$.

The general form for a causal identifiability problem that we consider in this paper is formulated as follows.

A set of input distributions of the form $P(\+ A_i \given \doo(\+ B_i),\+ C_i)$, a query $P(\+ Y \given \doo(\+ X), \+ Z)$  and a semi-Markovian causal graph $G$ over $\+ V$.
Output a formula  for the query  $P(\+ Y \given \doo(\+ X),\+ Z)$ over the input distributions, or decide that it is not identifiable.

Here  $\+ A_i,\+ B_i, \+ C_i$ are disjoint subsets of $\+ V$ for all $i$, and $\+ X,\+ Y,\+ Z$ are disjoint subsets of $\+ V$. The causal graph $G$ may contain vertices which describe mechanisms related to transportability and selection bias. In the following sections we explain several important special cases of this problem definition, some that have been considered in the literature and some which have not been.

§.§ Previously Considered Scenarios as Special Cases

We restate the concepts of transportability and selection bias under the causal inference framework, and show that identifiability in the scenarios of rows 1–7 of Table <ref> falls under the general form on row 8. We return to problems that involve missing data on rows 9–11 later in Section <ref>.

Causal Effect Identifiability Input is restricted to a passive observational distribution $P(\+ V)$. The target is either a causal effect $P(\+ Y \given \doo ( \+ X) )$ for row 1 of Table <ref> or a conditional causal effect $P(\+ Y \given \doo ( \+ X), \+ Z )$ for row 2 of Table <ref> [Shpitser and Pearl, 2006, Shpitser and Pearl, 2006].

$z$-identifiability Similarly to ordinary causal effect identification, the input consists of the passive observational distribution $P(\+ V)$ but also of experimental distributions known as surrogate experiments intervening on a set $\+ B$ [Bareinboim and Pearl, 2012]. Two restricting assumptions, which we call nested experiments and entire distributions, apply to surrogate experiments. Experiments are called nested experiments (NE) when for each experiment intervening a set of variables $\+B$, experiments intervening on all subsets of $\+B$ are available as well. Entire distributions (ED) denote the assumption that the union of observed and intervened variables is always the set of all variables $\+V$.

$g$-identifiability Unlike $z$-identifiability, the input does not contain the passive observational distribution $P(\+ V)$ but consists instead of surrogate experiments on the sets $\{ \+ B_i \}$ [Lee et al., 2019] without the assumption of nested experiments. The assumption of entire distributions holds.

Surrogate Outcome Identifiability Surrogate outcomes generalize the notion of surrogate experiments from $z$-identifiability. For surrogate outcomes, the assumption of nested experiments still holds, but the assumption of entire distributions can be dropped. Some less strict assumptions (SO) still apply [Tikka and Karvanen, 2019]. The idea of surrogate outcomes is that data from previous experiments are available, but the target $\+ Y$ was at most only partially measured in these experiments and the experiments do not have to be disjoint from $\+ X$.

Transportability The problem of incorporating data from multiple causal domains is known as transportability [Bareinboim and Pearl, 2013]. Formally, the goal is to identify a query in a target domain $\pi^*$ using data from source domains $\pi_1, \ldots ,\pi_n$. The domains are represented in the causal graph using a special set of transportability nodes $\+ T$ which is partitioned into disjoint subsets $\+ T_1,\ldots,\+ T_n$ corresponding to each domain $\pi_i$. The causal graph contains an extra edge $T_{ij} \rightarrow V_j$ whenever a functional discrepancy in $f_{V_j}$ or in $P(u_{V_j})$ exists between the target domain $\pi^*$ and the source domain $\pi_i$. The discrepancy is active if $T_{ij} = 1$ and inactive otherwise. A distribution associated with a domain $\pi_i$ is of the form $P(\+ A \given \doo(\+ B),\+ C,\+ T_i = 1,\+ T_{-i} = 0)$, where $\+ T_{-i}$ denotes the other subsets of the partition of $\+ T$ except $\+ T_i$. In other words, only the discrepancies between the $\pi_i$ and $\pi^*$ are active. A distribution corresponding to the target domain has no active discrepancies meaning that it is of the form $P(\+ A \given \doo(\+ B),\+ C,\+ T = 0$). Any variable is conditionally independent from inactive transportability nodes since their respective edges vanish. Furthermore, since transportability nodes set to $0$ vanish, we can assume any present transportability node to have the value $1$. Thus an input distribution from a domain $\pi_i$ takes the form $P(\+ A \given \doo(\+ B),\+ C,\+ T_i)$. In the specific case of $mz$-transportability, the assumptions of entire distributions (ED) and nested experiments in different domains (NEDD) apply, which means that $P(\+ V \setminus (\+ B_i^\prime \cup \+ T_i) \given \doo(\+ B_i^\prime), \+ T_i)$ is available for every subset $\+ B_i^\prime$ of $\+ B_i$ in each domain $\pi_i$.

Selection Bias Recoverability Selection bias can be seen as a special case of missing data, where the mechanism responsible for the preferential selection is represented in the causal graph by a special sink vertex $S$ [Bareinboim and Pearl, 2012]. Typical input for the recoverability problem is $P(\+ V \given S = 1)$, the joint distribution observed under selection bias. Just as in the case of transportability nodes, selection bias nodes only appear when the mechanism has been enabled. Thus we may assume that the input is of form $P(\+ V \given S)$. More generally, we can consider input distributions of the form $P(\+ A \given \doo(\+ B),\+ C, S)$.

§.§ New Scenarios as Special Cases

The following settings are special cases of the general identifiability problem of row 8 in Table <ref> that do not fall under any of the problems of rows 1–7. They serve as interesting additions to the cases considered in the literature. Concrete examples on these new scenarios are presented in Section <ref>. Section <ref> extends the general problem of row 8 in Table <ref> to the general problem with missing data on row 11 while also showcasing the special cases of rows 9 and 10.

Multiple Data Sources with Partially Overlapping Variable Sets

The scenario where only subsets of variables are ever observed together has been extensively considered in the causal discovery literature [Danks et al., 2009, Tillman and Spirtes, 2011, Triantafillou et al., 2010], but not in the context of causal effect identification. In the basic setting the input consists of passively observed distributions $P(\+ A_i)$ such that $\+ A_i \subset \+ V$.  We may also observe experimental distributions $P(\+ A_i \given \doo (\+ B_i))$ [Hyttinen et al., 2012, Triantafillou and Tsamardinos, 2015] or even conditionals $P(\+ A_i \given \doo(\+ B_i), \+ C_i)$. Our approach sets no limitations for the number or types of input distributions.

Combining Transportability and Selection Bias To the best of our knowledge, the frameworks of transportability and selection bias have not been considered simultaneously. The combination of these scenarios fits into the general problem formulation. For example, we may have access to two observational distributions originating from different source domains, but affected by the same biasing mechanism: $P(\+ A_1 \given \+ C_1, T_1, S)$ and $P(\+ A_2 \given \+ C_2, T_2, S)$, where $T_1$ and $T_2$ are the transportability nodes corresponding to the two source domains and $S$ is the selection bias node.

Recovering from Multiple Sources of Selection Bias In recent literature on selection bias as a causal inference problem, the focus has been on settings where only a single selection bias node is present <cit.>. However, multiple sources of selection bias are typical in longitudinal studies where dropout occurs at different stages of the study. Our approach is applicable for an arbitrary number of selection bias mechanisms and input distributions affected by arbitrary combinations of these mechanisms. In other words, if $\+ S$ is the set of all selection bias nodes present in the graph, the inputs can take the form $P(\+ A \given \doo(\+ B), \+ C, \+ S^\prime)$, where $\+ S^\prime$ is an arbitrary subset of $\+ S$.

§ A SEARCH BASED APPROACH FOR CAUSAL EFFECT IDENTIFICATION

The key to identification of causal effects is that interventional expressions can be manipulated using the rules of do-calculus. We present these rules for augmented graphs
where an additional intervention variable $I_X$ such that $I_X \rightarrow X$ is added to the induced graph for each variable $X$ [Spirtes et al., 1993, Pearl, 2009, Lauritzen, 2000] (see Figure <ref>(fig:augmented)).

Now a d-separation condition (or m-separation [Richardson, 2003]) of the form $\+ Y \independent{} \+ I_{\+ Z} \given \+ X, \+ Z, \+W \doublebar \+X$ means that nodes $\+ Y$ and intervention nodes  $\+ I_{\+ Z}$ of $\+ Z$ are d-separated given $\+ X$, $\+ Z$ and $\+ W$ in a graph where edges incoming to (intervened) $\+ X$ have been removed [Hyttinen et al., 2015, Dawid, 2002].
The three rules of do-calculus [Pearl, 1995] can be expressed as follows:

\begin{equation} \label{eq:docalc}
\begin{aligned}
P(\+ Y \given \doo(\+ X), \+ Z, \+ W) &= P(\+ Y \given \doo(\+ X), \+ W), \text{ if } \+ Y \independent{} \+ Z \given \+ X, \+ W \doublebar \+ X \\
P(\+ Y \given \doo(\+ X, \+ Z), \+ W) &= P(\+ Y \given \doo(\+ X), \+ Z, \+ W), \text{ if } \+ Y \independent{} \+ I_{\+ Z} \given \+ X, \+ Z, \+W \doublebar \+X \\
P(\+ Y \given \doo(\+ X, \+ Z), \+ W) &= P(\+ Y \given \doo(\+ X), \+ W), \text{ if } \+ Y \independent{}\+  I_{\+ Z} \given \+ X, \+W \doublebar \+X
\end{aligned}
\end{equation}

The rules are often referred to as insertion/deletion of observations, exchange of actions and observations, and insertion/deletion of actions respectively. Each rule of do-calculus is only applicable if the accompanying d-separation criterion (on the right-hand side) holds in the underlying graph. In addition to these rules, most derivations require basic probability calculus.

[1]
Target $Q = P(\+ Y \given \doo(\+ X),\+ W)$, a semi-Markovian graph $G$ and a set of known input distributions $\+ P=\{P_1,\ldots,P_n\}$.
A formula for $Q$

or NA if the effect is not identifiable.
for each $P_i \in \+ P$ do

Derive new distributions
from $P_i$

such that:

• The required d-separation criteria are satisfied by $G$.
• Any possible additional input required must also be in $\+ P$.

Add the new identified distributions to $\+ P$.
If $Q$ was derived, return a formula for it.
Return NA.

An outline of a search for causal effect identification.

Do-calculus directly motivates a forwards search over its rules. The outline of this type of search is given in Algorithm <ref>. The algorithm derives new identifiable distributions based on what has been given as the input or identified in the previous steps. For each identified distribution every rule of do-calculus and standard probability manipulations of marginalization and conditioning are applied in succession, until the target distribution is found, or no new distributions can be found to be identifiable. A preliminary version of this kind of search is used by Hyttinen et al., 2015 as a part of an algorithmic solution to causal effect identifiability when the underlying graph is unavailable.

The formulas produced by Algorithm <ref> correspond to short derivations and unnecessarily complicated expressions are avoided. Also, only distributions guaranteed to be identifiable are derived and used during the search.  Formulas for intermediary queries that were identified during the search are also available as a result. Alternatively, one could also start with the target and search towards the input distributions; a search in this direction will spend time deriving a number expressions that are inevitably non-identifiable based on the input. A depth-first search would produce unnecessarily complicated expressions.

[dot = 00Xbelow left] at (0,0) ;
[dot = 00Ybelow right] at (1,0) ;
[dot = 00Zabove] at (0.5,0.66) ;

[->] (X) – (Y);
[->] (Z) – (Y);
[->] (Z) – (X);

Example graph.

[dot = 00I_Xabove] at (-0.2,0.5) ;
[dot = 00Xbelow left] at (0,0) ;
[dot = 00Ybelow right] at (1,0) ;
[dot = 00Zabove] at (0.5,0.66) ;

[->] (I_X) – (X);
[->] (X) – (Y);
[->] (Z) – (Y);
[->] (Z) – (X);

Augmented graph.

[bordered] (v1) at (0,4) $P(X,Y,Z)$;
[bordered] (v2) at (-1,2.75) $P(Z)$;
[bordered] (v3) at (1,2.75) $P(Y \given X,Z)$;
[bordered] (v4) at (-1.35,1.5) $P(Z \given \doo(X))$;
[bordered] (v5) at (1.35,1.5) $P(Y \given \doo(X), Z)$;
[bordered] (v6) at (0,0.25) $P(Y,Z \given \doo(X))$;
[bordered] (v7) at (0,-1) $P(Y \given \doo(X))$;
[->] (v1) edge node[left,xshift=-2] M (v2);
[->] (v1) edge node[right,xshift=2] C (v3);
[->] (v2) edge node[left] R3 (v4);
[->] (v3) edge node[right] R2 (v5);
[->] (v4) edge node[left,xshift=-1] P (v6);
[->] (v5) edge node[right,xshift=1] P (v6);
[->] (v6) edge node[left] M  (v7);

A derivation for $P(Y \given \doo(X))$.

The back-door criterion holds in the example graph (fig:intro_backdoor) for $Z$. The augmented graph (fig:augmented) includes the intervention node $I_X$ for $X$ explicitly. The labels M, C, P, R2 and R3 in the derivation of (fig:derivation)
refer to marginalization, conditioning, product rule and rules 2 and 3 of do-calculus respectively (see Table <ref>). The required d-separation conditions $Y \independent I_X \given Z,X$ for R2 and $Z \independent I_X$ for R3 hold in the augmented graph (fig:augmented).

The search can easily derive for example the back-door criterion in the graph of Figure <ref>(fig:intro_backdoor) as shown by the derivation in Figure <ref>(fig:derivation). The target is $Q = P(Y \given \doo(X))$ and input is $\+ P = \{P(X,Y,Z)\}$. From $P(X,Y,Z)$ the search first derives the marginal $P(Z)$ and the conditional $P(Y \given X,Z)$. Then $P(Z \given \doo(X))$ is derived by the third rule of do-calculus because $Z \independent I_X$. The second rule derives $P(Y \given \doo(X), Z)$  from $P(Y \given X,Z)$ as $Y \independent I_X \given Z,X$. The two terms can be combined via the product rule of probability calculus to get $P(Y, Z \given \doo(X))$ and finally the target is $P(Y \given \doo(X))$ is just a marginalization of this. The familiar formula $\sum_{Z} P(Y \given X,Z)P(Z)$ is thus obtained.

However, it is not straightforward to make a search over do-calculus computationally feasible. The search space in Figure <ref>(fig:derivation) shows only the parts that resulted in the identifying formula: for example all passively observed marginals and conditionals over $\+ V$ can be derived from the input $P(\+ V)$. Especially in a non-identifiable case a naive search may go through a huge space before it can return the non-identifiable verdict. The choice of rules is also not obvious: a redundant rule may make the search faster or slower; false non-identifiability may be concluded if a necessary rule is missing. Also the order in which the rules are applied can have a large impact on the performance of the search. In the following sections we will provide non-trivial solutions to these challenges.

§.§ Rules

Table <ref> lists the full set of rules used to manipulate distributions during the search, generalizing the work by Hyttinen et al., 2015.

Rule     Additional Input     Output     Description

$1+$        $P(\+ Y \given \doo(\+ X), \+ Z, \+ W)$     Insertion of observations

$1-$        $P(\+ Y \given \doo(\+ X), \+ W \setminus \+ Z)$     Deletion of observations

$2+$        $P(\+ Y \given \doo(\+ X,\+ Z), \+ W \setminus \+ Z)$     Observation to action exchange

$2-$        $P(\+ Y \given \doo(\+ X \setminus \+ Z),\+ Z, \+ W)$     Action to observation exchange

$3+$        $P(\+ Y \given \doo(\+ X,\+ Z), \+ W)$     Insertion of actions

$3-$        $P(\+ Y \given \doo(\+ X \setminus \+ Z), \+ W)$     Deletion of actions

$4$        $P(\+ Y \setminus \+ Z \given \doo(\+ X),\+ W)$     Marginalization

$5$        $P(\+ Y \setminus \+ Z \given \doo(\+ X),\+ Z,\+ W)$     Conditioning

$6+$     $P(\+ Z \given \doo(\+X),\+ W \setminus \+ Z)$     $P(\+Y, \+Z \given \doo(\+ X),\+ W \setminus \+ Z)$     Chain rule multiplication

$6-$     $P(\+ Z \given \doo(\+ X), \+ Y,\+ W)$      $P(\+Y, \+Z \given \doo(\+ X),\+W)$     Chain rule multiplication

The rules used to manipulate input distributions of the form $P(\+ Y \given \doo(\+ X), \+ W)$. The output distribution is identified if the input has been previously identified and if the corresponding d-separation criteria (<ref>) hold in the graph (for rules $1\pm,2\pm$ and $3\pm$) or if the additional input has also been identified (rules $6\pm$). The sets $\+ Y, \+ X$ and $\+ W$ are disjoint. The role of the set $\+ Z$ depends on the rule being applied (see Table <ref>).

Do-calculus Rules $1\pm, 2\pm$ and $3\pm$ correspond to the rules of do-calculus such that rules $1+, 2+, 3+$ are used to add conditional variables and interventions and rules $1-, 2-, 3-$ are used to remove them. Each rule is only valid if the corresponding d-separation criterion given in the beginning of Section <ref> holds.

Probability theory Rule $4$ performs marginalization over $\+ Z \subset \+ Y$, and produces a summation at the formula level:

\[
P(\+Y \setminus \+ Z \given \doo(\+ X),\+ W) = \sum_{\+ Z} P(\+ Y \given \doo(\+ X),\+ W).
\]

Similarly, rule $5$ conditions on a subset $\+ Z \subset \+ Y$ to obtain the following formula:

\[
P(\+ Y \setminus \+ Z \given \doo(\+ X),\+Z, \+ W) = \frac{P(\+ Y \given \doo(\+ X),\+ W)}{\sum_{\+ Y \setminus \+ Z} P(\+ Y \given \doo(\+ X),\+ W)}.
\]

Rules $6+$ and $6-$ perform multiplication using the chain rule of probability which requires two known distributions. When rule $6+$ is applied, the distribution $P(\+ Y \given \doo(\+ X),\+ W)$ is known and we check whether $P(\+ Z \given \doo(\+ X),\+ W \setminus \+ Z)$ is known as well. For rule $6-$, the roles of the distributions are reversed. In the case of rule $6+$, $\+ Z$ is a subset of $\+ W$ and we obtain

\[
P(\+ Y,\+ Z \given \doo(\+ X), \+ W \setminus \+ Z) = P(\+ Y \given \doo(\+ X),\+ W) P(\+ Z \given \doo(\+ X),\+ W \setminus \+ Z).
\]

The two version of the chain rule are needed: it may be the case that when expanding $P(\+ Y \given \doo(\+ X), \+ W)$ with rule $6+$ the additional input $P(\+ Z \given \doo(\+X),\+ W \setminus \+ Z)$ is only identified later in the search. Then, $P(\+Y, \+Z \given \doo(\+ X),\+W)$ is identified when rule $6-$ is applied to $P(\+ Y \given \doo(\+X),\+ W)$.

§.§ Improving the Efficacy of the Search

In this section, we present various techniques that improved the efficiency of the search.
These findings are implemented in a search algorithm in Section <ref>.

§.§.§ Term Expansion

Term expansion refers to the process of deriving new distributions from an input distribution using the rules of Table <ref>. By term we mean a single identified distribution. A term is considered expanded if the rules of Table <ref> have been applied to it in every possible way when the term is in the role of the input. Note that an expanded distribution may still take the role of an additional input when another term is being expanded. Consider the step of expanding the input term in Table <ref> to all possible outputs with any rule. This can be done by enumerating every non-empty subset $\+ Z$ of $\+ V$, and applying the rule with regard to it.

Rule     Validity condition     Termination condition

$1+$      $\+ Z \cap (\+ Y \cup \+ X \cup \+ W) = \emptyset$

$1-$      $\+ Z \subseteq \+ W$                                   $\+ W = \emptyset$

$2+$      $\+ Z \subseteq \+ W$                                   $\+ W = \emptyset$

$2-$      $\+ Z \subseteq \+ X$                                   $\+ X = \emptyset$

$3+$      $\+ Z \cap (\+ Y \cup \+ X \cup \+ W) = \emptyset$

$3-$      $\+ Z \subseteq \+ X$                                   $\+ X = \emptyset$
$4$       $\+ Z \subset \+ Y$                                     $|\+ Y| = 1$

$5$       $\+ Z \subset \+ Y$                                     $|\+ Y| = 1$

$6+$      $\+ Z \subseteq \+ W$                                   $\+ W = \emptyset$

$6-$      $\+ Z \cap (\+ Y \cup \+ X \cup \+ W) = \emptyset$

The conditions for the enumerated subset $\+ Z$ for applying the rules of Table <ref> to a term $P(\+ Y \given \doo(\+ X),\+ W)$. For rules $6+$ and $6-$, the conditions specify valid variables of the second required term.

Table <ref> outlines the requirements for $\+ Z$ for each rule of the search. Table <ref> tells us that when an observation $\+ Z$ is added using rule $1+$, it cannot be contained in any of the sets $\+ Y, \+ X$ or $\+ W$ since they are already present in the term. Only observations that are present can be removed, which is why $\+ Z$ has to a subset of $\+ W$ when applying rule $1-$. We may skip the application of this rule if the set of observations is empty for the current term. The exchange of observations to experiments using rule $2+$ has similar requirements for set $\+ Z$ as rule $1-$. Exchanging experiments to observations using rule $2-$ works in a similar fashion. Only experiments that are present can be exchanged which means that $\+ Z \subseteq \+ X$. This rule can be skipped if the set of experiments is empty. New experiments are added using rule $3+$ with similar requirements as rule $1+$. Well-defined subsets for using rule $3-$ are the same as for rule $2-$. For rules $4$ and $5$, the only requirement is that $\+ Z$ is a proper subset of $\+ Y$. When the chain rule is applied with rule $6+$, we require that the variables of the second product term is observed in the first term. When applied in reverse with rule $6-$, the variables of the second term must not be present in the first term.

§.§.§ Termination Conditions

Additionally, Table <ref> lists the termination condition for each rule: if it is satisfied by the current term to be expanded we know that the rule cannot be applied to it. The following simple lemma shows that when any of the termination conditions hold, no new distributions can be derived from it using the respective rule, which allows the search to directly proceed to the next rule.

Let $G$ be a semi-Markovian graph and let $\+ Y, \+ X$ and $\+ W$ be disjoint subsets of $\+ V$. Then all of the following are true:

(i)If $\+ W = \emptyset$, then rule $1-$ of Table <ref> cannot be used.

(ii)If $\+ W = \emptyset$, then rule $2+$ of Table <ref> cannot be used.

(iii)If $\+ X = \emptyset$, then rule $2-$ of Table <ref> cannot be used.

(iv)If $\+ X = \emptyset$, then rule $3-$ of Table <ref> cannot be used.

(v)If $|\+ Y| = 1$, then rule $4$ of Table <ref> cannot be used.

(vi)If $|\+ Y| = 1$, then rule $5$ of Table <ref> cannot be used.

(vii)If $\+ W = \emptyset$, then rule $6+$ of Table <ref> cannot be used.

For (i), the set $\+ W$ is empty so the application of rule $1-$ using any subset $\+ Z$ would result in $P(\+ Y \given \doo(\+ X), \+ W \setminus \+ Z) = P(\+ Y \given \doo(\+ X), \+ W)$ which is already identified. For (ii),
the set $\+ W$ is empty so no observation can be exchanged for an action using the second rule of do-calculus. For (iii), the set $\+ X$ is empty so no action can be exchanged for an observation using the second rule of do-calculus. For (iv), the  set $\+ X$ is empty so the application of rule $3-$ using any subset $\+ Z$ would result in $P(\+ Y \given \doo(\+ X \setminus \+ Z), \+ W) = P(\+ Y \given \doo(\+ X), \+ W)$ which is already identified.
For (v) and (vi), the set $\+ Y$ only has a single vertex, so it cannot have a non-empty subset. For (vii), the set $\+ W$ is empty so no subset $\+ Z \subset \+ W$ can exist for the second input.

§.§.§ Rule Necessity

The rule 1 of do-calculus can be omitted as shown by <cit.>. Instead of inserting an observation using rule 1, we can insert an intervention and then exchange it for an observation. Similarly, an observation can be removed by first exchanging it for an intervention and then deleting the intervention. It follows that rules $1+$ and $1-$ of Table <ref> are unnecessary for the search.

[dot = 00X_1left] at (0,0) ;
[dot = 00X_2right] at (2,0) ;
[dot = 00Z_1above] at (1,1) ;
[dot = 00Z_2below] at (1,-1) ;
[dot = 00Wabove right] at (1,0) ;
[dot = 00Yright] at (3,0) ;

[->] (X_1) – (Z_1);
[->] (X_1) – (Z_2);
[->] (X_2) – (Z_1);
[->] (X_2) – (Z_2);
[->] (Z_1) – (Y);
[->] (Z_2) – (Y);
[->] (Z_1) – (W);
[->] (Z_2) – (W);
[->] (X_1) – (W);
[->] (X_2) – (W);

A graph for the example where all rules of Table <ref> are required for identifying the target quantity.

The following example shows that the remaining rules of Table <ref> are all necessary. In the graph of Figure <ref>, the causal effect $P(Y, X_1 \given \doo(X_2), W)$ can be identified from the inputs $P(W \given \doo(X_2), Y, X_1)$,  $P(Y \given \doo(X_2), Z_1, Z_2, X_1)$, $P(X_1 \given \doo(X_2), W)$, $P(Z_2, X_2 \given \doo(X_1))$ and $P(Z_1 \given \doo(X_1, Y), X_2)$ when all rules are available, but not when any individual rule is omitted. This can be verified by running the search algorithm presented at the beginning of Section <ref> or the more advanced algorithm of Section <ref> with each rule switched off individually.

§.§.§ Early Detection of Non-identifiable Instances

Worst-case performance of the search can be improved by detecting non-identifiable quantities directly based on the set of inputs before launching the search. The following theorem provides a sufficient criterion for non-identifiability.

Let $G$ be a semi-Markovian graph, let $Q = P(\+ Y \given \doo(\+ X ), \+ W)$ and let

\[
\+ P = \{P(\+ A_1 \given \doo(\+ B_1), \+ C_1), \ldots, P(\+ A_n \given \doo(\+ B_n), \+ C_n) \}.
\]
Then $Q$ is not identifiable from $\+ P$ in $G$ via rules of Table <ref> if

\[
\+ Y \not\subseteq \bigcup_{i = 1}^n \+ A_i,
\]

Since $ \+ Y \not\subseteq \bigcup_{i = 1}^n \+ A_i$, there exists a variable $Y_j \in \+ Y$ such that none of the sets $\+ A_i$ contain it. No rule of Table <ref> outputs a distribution $P(\+ Y^\prime \given \doo(\+ X^\prime), \+ W^\prime)$ such that some member of $\+ Y^\prime$ would not already exist on the left-hand side of the input or additional input of the rule. Thus there is no sequence of rules that when applied to the available inputs $\+ P$ would result in a distribution of the form $P(Y_j, \cdot \given \doo(\cdot), \cdot)$. Thus there is no such sequence for $P(\+ Y \given \doo(\+ X), \+ W)$.

In other words, Theorem <ref> can be used to verify that the entire set $\+ Y$ of a target distribution $P(\+ Y \given \cdot)$ cannot be constructed from the inputs. If this is the case, the target quantity is not identifiable.

§.§.§ Heuristics

During the search, we always expand one term at a time through the rules and store the newly identified distributions. In order for the search to perform fast, we need to decide which branches are the most promising and should therefore be expanded first. We can do this by defining a proximity function relating the source terms and the target query, and by always expanding the closest term first.

Our suggestion here is motivated by the way an educated person might apply do-calculus in a manual derivation. Our chosen proximity function $h$ links the target distribution $P^{t} = P(\+ A_t \given \doo(\+ B_t),\+C_t)$ and a source distribution $P^{s} = P(\+ A_s \given \doo(\+B_s),\+ C_s)$ in the following way:

\begin{align*}
h(P^t,P^s) &= 10 |\+ A_t \cap \+ A_s| + 5 |\+ B_t \cap \+ B_s| + 3 |\+ C_t \cap \+ C_s| - 2|\+ A_t \setminus \+ A_s| - 2|\+ B_t \setminus \+ B_s| \\
&\quad - 2|\+ B_s \setminus \+ B_t| - |\+ C_t \setminus \+ C_s| - |\+ C_s \setminus \+ C_t|.
\end{align*}

Each input distribution and terms derived using the search are assigned into a priority queue, where the priority is determined by the value given by $h$. Distributions closer to the target are prioritized over other terms.

The weight 10 for the term $|\+ A_t \cap \+ A_s|$ indicates that having the correct response variables is considered as the first priority. Having the correct intervention is considered as the second priority (weight 5) and having the correct condition as the third priority (weight 3). The remaining terms in $h$ penalize variables that are in the target distribution but not in the source distribution or vice versa. Again, variables that are intervened on are considered to be more important than conditioning variables.

§.§ The Search Algorithm

We take Algorithm <ref> as our starting point and compile the results of Section <ref> into a new search algorithm called . This algorithm is capable of solving generalized identifiability problems (row 8 in Table <ref>) while streamlining the search process through a heuristic search order and elimination of redundant rules and subsets. The pseudo-code for  is shown in Algorithm <ref>.

The algorithm begins by checking whether the query can be solved trivially without performing the search. This can happen if the target $Q$ is a member of the set of inputs or if Theorem <ref> applies. Next, we note that each input distribution in the set $\+ P$ is marked as unexpanded at the beginning of the search. Distributions in $\+ P$ are expanded one at a time by applying every rule of Table <ref> in every possible way.

[1]
Target $Q = P(\+ Y \given \doo(\+ X),\+ W)$, a semi-Markovian graph $G$ and a set of known distributions $\+ P=\{P_1,\ldots,P_n\}$.
A formula $F$ for $Q$ in terms of $\+ P$ or NA
if $Q \in \+ P$, return $Q$
if target is non-identifiable by Theorem <ref>, then return NA
let $\+ U$ be the set of unexpanded distributions,
initially $\+ U := \+ P$
while $\+ U \neq \emptyset$, do

let $P^\prime$ be the unexpanded distribution closest to the target: $P^\prime = \underset{P_i \in \+ U}{\mathrm{argmax}}\;h(Q, P_i)$
let $\+ M$ be the set of rules of Table <ref>, without rules $1\pm$, such that the termination
conditions of Table <ref> do not hold with respect to $P^\prime$.

let $\+ P^*$ be the set of all distributions derived from $P^\prime$ using the rules in $\+ M$
for each new candidate distribution $P^* \in \+ P^*$, do

if $P^*$ is already in $\+ P$, then continue
if the validity conditions of Table <ref> are not satisfied by $P^*$, then continue
if an additional input is required that is not in $\+ P$, then continue

if rule $2\pm$ or $3\pm$ of Table <ref> is applied and the corresponding d-separation criterion (<ref>)
is not satisfied by $G$, then continue
if $P^* = Q$, then
Derive a formula $F$ for $Q$ by backtracking.
return $F$
Add $P^*$ to $\+ P$, add $P^*$ to $\+ U$

Mark $P^\prime$ as expanded: remove $P^\prime$ from $\+ U$

return NA

The iteration over the unexpanded distributions $\+ U$ proceeds as follows (lines 4–5). Each input distribution and terms derived from it are assigned into a priority queue, where the priority is determined by the value given by the proximity function $h$. Distributions closest to the target are expanded first. In the implementation, only the actual memory addresses of the distribution objects are placed into the queue. The set $\+ P$ is implemented as a hash table that serves as a container for all input distributions and those derived from them. Each new distribution is assigned a unique index that also serves as the hash function for this table. The distribution objects contained in the table are represented uniquely by three integers corresponding to the sets $\+ A, \+ B$ and $\+ C$ of the general form $P(\+ A \given \doo(\+ B), \+ C)$. A distribution object also contains additional auxiliary information such as which rule was used to derive it, whether it is expanded or not and from which distribution it was obtained. This information is used to construct the derivation if the target is found to be identifiable.

Multiple distributions can share the same value of the proximity function $h$. In the case that multiple candidates share the maximal value, the one that was derived the earliest takes precedence. When the unexpanded distribution currently closest to the target is determined, the rules of Table <ref> are applied sequentially for all valid subsets dictated by Table <ref>. When rules two and three of do-calculus are considered, the necessary d-separation criteria is checked from $G$ (line 12). For the chain rule, the presence of the required second input is also verified. The reverse lookup is implemented by using another hash table, where the hash function is based on the unique representation of each distribution object. The values contained in the table are the indices of the derived distributions. The same hash table is also used to ensure that we do not attempt to derive distributions again that have been previously found to be identifiable from the inputs.

We construct a set $\+ M$ of applicable rules for each unexpanded distribution $P^\prime$ using the termination conditions of Table <ref> (line 6). If all the necessary conditions have been found to hold for an applicable rule and a subset, the newly derived distribution $P^*$ is added to the set of known distributions and placed into the priority queue as an unexpanded distribution. When the applicable rules and subsets have been exhausted for the current distribution $P^\prime$, the term is marked as expanded and removed from the queue (line 17). If the target distribution is found at any point (line 13), a formula is returned for it in terms of the original inputs. Alternatively, we can also continue deriving distributions to obtain different search paths to the target that can possibly produce different formulas for it. If instead we exhaust the set of unexpanded distributions by emptying the queue, the target is deemed non-identifiable by the search (line 18).

We keep track of the rules that were used to derive each new distribution in the search. This allows us to construct a directed graph of the derivation where each root node is a member of the original input set $\+ P$ and their descendants are the distributions derived from them during the search. Each edge represents a manipulation of the parent node(s) to obtain the child node. For an identifiable target quantity, the formula $F$ is obtained by backtracking the chain of manipulations recursively until the roots are reached (line 14). The derivation of the example in the beginning of Section <ref> depicted in Figure <ref>(fig:derivation) can be efficiently found by applying this procedure.

We assess the worst case complexity of  in terms of the input graph size, which is the primary determining factor of the search time.
Checks for separation and termination and validation conditions all run in polynomial time with respect to the number of vertices, but the number of distributions grows very rapidly. In a hypothetical absolute worst case scenario, the target is not identifiable, but every other distribution is. Supposing a graph with $d$ vertices, we can determine
how many distributions would have to be derived by the search in order to exhaust the search space. For a distribution of the form $P(\+ A \given \doo(\+ B), \+ C)$, any variable in $\+ V$ can either belong to the sets $\+ A, \+ B$ or $\+ C$ or not be present in the distribution, meaning that there are $4^d$ ways to categorize every variable. However, we must account for the fact that there must always be at least one variable in the set $\+ A$. There are $3^d$ ways to assign all the variables in such a way that no variable is a member of set $\+ A$.

Thus the total number of distributions considered by  grows as $O(4^d)$. See the simulations of Section <ref> for running time performance in more realistic settings.

§.§ Soundness and Completeness Properties

We are ready to establish some key theoretical properties of . The first theorem considers the correctness of the search.

always terminates: if it returns an expression for the target $Q$, it is correct, if it returns NA then $Q$ is not identifiable with respect to the rules of do-calculus and standard probability manipulations (in Table <ref>).

Each new distribution is derived by using only well-defined manipulations as outlined by Table <ref> and by ensuring that the required separation criteria hold in $G$ when rules of do-calculus are concerned. It follows that if the search terminates and returns a formula for the target distribution, it was reached from the set input distributions through a sequence of valid manipulations. If  terminates as a result of Theorem <ref>, we are done.
Suppose now that Theorem <ref> does not apply. By definition,  enumerates every rule of Table <ref> for every well-defined subset of Table <ref>. By Lemma <ref>, no distributions are left out by applying the termination criteria of Table <ref>. We know that if rules $1\pm$ of Table <ref> are omitted, the distributions generated by these rules can be obtained by a combination of rules $2\pm$ and $3\pm$. Furthermore, the order in which the distributions are expanded does not matter, as every possible manipulation is carried out nonetheless. The search will eventually terminate, since distributions that have already been derived are not added again to the set of unexpanded distributions and there are only finitely many ways to apply the rules of Table <ref>.

The following theorem provides a completeness result in connection to existing identifiability results. Since do-calculus has been shown to be complete with respect to (conditional) causal effect identifiability, $z$-identifiability, $g$-identifiability and transportability, it follows that  is complete for these problems as well.
If  returns NA in the settings in rows 1–4 and 6 in Table <ref>, then the query is non-identifiable.

Do-calculus has been shown to be complete in these settings. The rules of probability calculus encode what is used in the algorithms as can be seen for example from the proofs of Theorem 7 and Lemmas 4–8 of Shpitser and Pearl, 2006.

It is not known whether the rules implemented in  are sufficient for other more general identifiability problems since it is conceivable that some additional rules might exist that would be required to achieve completeness. One such generalization is the inclusion of missing data in the causal model, which we present in Section <ref>. However, if one were to show that do-calculus (or any other set of rules included in ) is complete for some special case of the generalized identifiability problem, then  would be complete for this problem as well. In the following sections we will use the term “identifiable by ” to refer to causal queries that can be identified by .

§ EXTENSION TO MISSING DATA PROBLEMS

The SCM framework can be extended to describe missing data mechanisms. For each variable $V_i$, two special vertices are added to the causal graph. The vertex $V_i^*$ is the observed proxy variable which is linked to the true variable $V_i$ via the missingness mechanism [Little and Rubin, 1986, Mohan et al., 2013]:

\begin{equation} \label{eq:missingness}
V_i^* = \begin{cases}
V_i, & \mathrm{if}\; R_{V_i} = 1, \\
\NA,  & \mathrm{if}\; R_{V_i} = 0,
\end{cases}
\end{equation}

where $\NA$ denotes a missing value and $R_{V_i}$ is called the response indicator (of $V_i$). In other words, the variable $V_i^*$ that is actually observed matches the true value $V_i$ if it is not missing ($R_{V_i} = 1$). We note that in this formulation, each true variable has its own response indicator, meaning that we do not consider shared indicators between variables or multiple indicators for a single variable. Figure <ref> in Section <ref> depicts some examples of graphs containing missing data mechanisms. Furthermore, if there is no missingness associated with a given variable $V_i$ meaning that it is fully observed, the corresponding response indicator $R_{V_i}$ always has the value $1$. The omission of a proxy variable and a response indicators of a specific variable from a graph encodes the assumption that the variable in question if fully observed. Note that intervention nodes are added for true variables and response indicators but not for proxy variables. On a symbolic level one could intervene on proxy variables, however we are only interested in interventions that keep equation (<ref>) intact.

The observed vertices of the causal diagram can be partitioned into three categories

\[
\+ V = \+ V^\true \cup \+ V^* \cup \+ V^r,
\]

where $\+ V^\true$ is the set of true variables, $\+ V^*$ is the set of proxy variables and $\+ V^r$ is the set of response indicators. For any subset $\+ Z \subset \+ V$ we define the same partition via $\+ Z^\true = \+ Z \cap \+ V^\true$, $\+ Z^* = \+ Z \cap \+ V^*$ and $\+ Z^r = \+ Z \cap \+ V^r$.

The definition of the response indicator connects a proxy variable and a true variable. Typically this connection is only of interest for those variables $V_i$ that have missing data, meaning that $P(R_{V_i} = 1) < 1$. Furthermore, we often utilize a proxy variable corresponding to a specific true variable and conversely, a true variable corresponding to a specific proxy variable. We define this correspondence explicitly in the following way

\begin{align*}
\+ Z^\parprox &= \{V_i^* \in \+ V^* \given V_i \in \+ Z^\true, P(R_{V_i} = 1) < 1\}, \\
\+ Z^\proxpar &= \{V_i \in \+ V^\true \given V_i^* \in \+ Z^*, P(R_{V_i} = 1) < 1\}.
\end{align*}

Similarly, given a set $\+ Z$ we often require the set of the response indicators that define the missingness mechanism for the true variables that are member of $\+ Z$. This set is defined as follows

\[
\+ R_{\+Z} = \{R_{V_i} \in \+ V^r \given V_i \in \+ Z^\true, P(R_{V_i} = 1) < 1 \}.
\]

It is important to note the difference between the sets $\+ Z^r$ and $\+ R_{\+Z}$; the first set denotes the set of response indicators that are members of $\+ Z$ while the second gives the corresponding response indicators for the true variables that are member of $\+ Z$.

Our method is also capable of processing queries when the causal graph contains missing data mechanisms where the sets $\+ A_i, \+ B_i$ and $\+ C_i$ of some of the input distributions may be restricted to contain observed variables in $\+ V^* \cup \+ V^r$. An active response indicator $R_{V_i} = 1$ is denoted by $R_{V_i}^1$. Similarly, for sets of response indicators $\+ R_{\+Z}^1$ denotes that all indicators in the set are active. Proxy variables are not explicitly shown in graphs for clarity.

Determining identifiability is challenging under missing data. As evidence of this, even some non-interventional queries require the application of do-calculus [Mohan and Pearl, 2018]. Furthermore, the rules used in the search of Table <ref> are no longer sufficient and deriving the desired quantity necessitates the use of additional rules that stem from the definition of the proxy variables and the response indicator. Each new true variable also has a higher impact on computational complexity, since the corresponding response indicator and proxy variable are always added to the graph as well.

Rule     Additional Input     Output     Description

$1+$        $P(\+ Y \given \doo(\+ X), \+ Z, \+ W)$     Insertion of observations

$1-$        $P(\+ Y \given \doo(\+ X), \+ W \setminus \+ Z)$     Deletion of observations

$2+$        $P(\+ Y \given \doo(\+ X,\+ Z), \+ W \setminus \+ Z)$     Obs. to action exchange

$2-$        $P(\+ Y \given \doo(\+ X \setminus \+ Z),\+ Z, \+ W)$     Action to obs. exchange

$3+$        $P(\+ Y \given \doo(\+ X,\+ Z), \+ W)$     Insertion of actions

$3-$        $P(\+ Y \given \doo(\+ X \setminus \+ Z), \+ W)$     Deletion of actions

$ 4$        $P(\+ Y \setminus \+ Z \given \doo(\+ X),\+ W)$     Marginalization

$ 5$        $P(\+ Y \setminus \+ Z \given \doo(\+ X),\+ Z,\+ W)$     Conditioning

$6+$     $P(\+ Z \given \doo(\+X),\+ W \setminus \+ Z)$     $P(\+Y, \+Z \given \doo(\+ X), \+ W \setminus \+ Z)$     Chain rule multiplication

$6-$     $P(\+ Z \given \doo(\+ X), \+ Y,\+ W)$      $P(\+Y, \+Z \given \doo(\+ X),\+W)$     Chain rule multiplication

$7+$      $P(\+Z \given \doo(\+X),\+W)$                          $P(\+ Y \setminus \+ Z \given \doo(\+X),\+Z,\+W)$     Chain rule conditioning (numerator)

$7-$      $P(\+Z \given \doo(\+X),\+W, \+ Y \setminus \+ Z)$     $P(\+ Y \setminus \+ Z \given \doo(\+X), \+W)$     Chain rule conditioning (numerator)

$8+$      $P(\+Y, \+ Z \given \doo(\+X),\+W)$                    $P(\+ Z \given \doo(\+X),\+W, \+ Y)$     Chain rule conditioning (denominator)

$8-$      $P(\+Y, \+ Z \given \doo(\+X),\+W \setminus \+ Z)$     $P(\+ Z \given \doo(\+X),\+W)$     Chain rule conditioning (denominator)

$9+$         $P(\+ Y \given \doo(\+X),\+W \setminus \+ R_{\+ Z}, \+ R_{\+ Z}^1)$     Enable response indicators

$9-$         $P(\+ Y \setminus \+ R_{\+ Z}, \+ R_{\+ Z}^1 \given \doo(\+ X),\+ W)$     Enable response indicators

$10+$        $P(\+ Y \given \doo(\+X),\+ W \setminus \+ Z^*,\+ Z^\proxpar)$     Proxy variable exchange

$10-$        $P(\+ Y \setminus \+ Z^*,\+Z^\proxpar \given \doo(\+X),\+W)$     Proxy variable exchange

Extended set of rules for missing data problems used to manipulate input distributions of the form $P(\+ Y \given \doo(\+ X), \+ W)$.  Rules $1\pm, 2\pm, 3\pm, 4, 5$ and $6\pm$ are the same as in Table <ref>. For rules $7\pm$ and $8\pm$, the additional input has also been identified. The sets $\+ Y, \+ X$ and $\+ W$ are disjoint. Sets $\+ Y$ and $\+ W$ may contain true variables, proxy variables and response indicators. Set $\+ X$ may only contain true variables and response indicators. The roles of the sets $\+ Z$ and $\+ R_{\+ Z}$ depend on the rule being applied (see Table <ref>).

Table <ref> extends the set of rules of Table <ref> to missing data problems by providing manipulations related to the missingness mechanism. Rules $7\pm$ and $8\pm$ perform conditioning using the chain rule. These rules are necessary in the case that set $\+ Y$ contains missing data mechanisms that have been enabled and thus cannot be marginalized over by using rule $5$.

Rules $9\pm$ are used to enable response indicators, which then facilitates the use of rules $10\pm$. These last two rules exchange proxy variables to their true counterparts when the corresponding response indicators are enabled. For example, under the conditions specified in Table <ref>, rule $9+$ can be applied on $P(Y,X^*\given R_X)$ to first obtain $P(Y,X^* \given R^1_X)$ by enabling $R_X$. Then, rule $10+$ can applied to this distribution to obtain $P(Y,X \given R^1_X)$ by exchanging $X^*$ for $X$.

Similarly to Table <ref>, Table <ref> outlines the valid subsets $\+ Z$ for applying the extended rules of Table <ref>. A major difference to the original validity and termination conditions is the addition of the missing data condition that outlines the additional requirements that must be satisfied when missingness mechanisms are present. For the rules that are shared by Tables <ref> and <ref>, the missing data condition ensures that a true variable and its proxy counterpart never appear in the same term at the same time. For example, we cannot add an intervention on $X$ to $P(X^*)$. It also ensures that we do not carry out summation over enabled response indicators in the case of rules 4 and 5. When applying rules $9\pm$, the condition also ensures that we do not attempt to enable a response indicator that is already enabled. For rules $10\pm$, the conditions guarantee that a proxy can only be exchanged to its true counterpart if its corresponding response indicator is enabled and present in the input term.

Rule     Validity cond.     Missing data condition     Term. cond.

$1+$       $\+ Z \cap \+ T = \emptyset$                                  $\+ Z \cap \left(\+ T^\proxpar \cup \+ T^\parprox \cup \+ Z^\proxpar \cup \+ Z^\parprox\right) = \emptyset$

$1-$       $\+ Z \subseteq \+ W$                                                                           $\+ W = \emptyset$

$2+$       $\+ Z \subseteq \+ W$                                   $\+ Z \cap \+ W^* = \emptyset$          $\+ W = \emptyset$

$2-$       $\+ Z \subseteq \+ X$                                                                           $\+ X = \emptyset$

$3+$       $\+ Z \cap \+ T = \emptyset$       $\+ Z \cap \left(\+ T^\proxpar \cup \+ T^\parprox \cup \+ Z^\proxpar \cup \+ Z^\parprox\right) = \emptyset$

$3-$       $\+ Z \subseteq \+ X$                                                                           $\+ X = \emptyset$
$4$        $\+ Z \subset \+ Y$                                                                 $\+ Z \cap (\+ R^a \cap \+ Y) = \emptyset$            $|\+ Y| = 1$

$5$        $\+ Z \subset \+ Y$                                                                 $ (\+ Y \setminus \+ Z) \cap (\+ R^a \cap \+ Y) = \emptyset$     $|\+ Y| = 1$

$6+$       $\+ Z \subseteq \+ W$                                               $\+ W = \emptyset$

$6-$       $\+ Z \cap \+ T= \emptyset$      $\+ Z \cap \left(\+ T^\proxpar \cup \+ T^\parprox \cup \+ Z^\proxpar \cup \+ Z^\parprox\right) = \emptyset$
$7+$         $\+ Z \subset \+ Y$                                              $|\+ Y| = 1$

$7-$         $\+ Z \subset \+ Y$                                              $|\+ Y| = 1$

$8+$         $\+ Z \cap \+ T = \emptyset$

$8-$         $\+ Z \subseteq \+ W$                                            $\+ W = \emptyset$

$9+$         $\+ R_{\+ Z} \subseteq \+ W^r, \+ R_{\+ Z} \cap \+ R^a = \emptyset$                                $\+ W^r = \emptyset$

$9-$         $\+ R_{\+ Z} \subseteq \+ Y^r,\; \+ R_{\+ Z} \cap \+ R^a = \emptyset$                                $\+ Y^r = \emptyset$

$10+$        $\+ Z^* \subseteq \+ W^*, \+ R_{\+ Z^\proxpar} \subseteq \+ R^a, \+ R_{\+ Z^\proxpar} \subseteq \+ W^r$                     $\+ R^a = \emptyset$

$10-$        $\+ Z^* \subseteq \+ Y^*,\; \+ R_{\+ Z^\proxpar} \subseteq \+ R^a, (\+ R_{\+ Z^\proxpar} \subseteq \+ W^r \textrm{ or } \+ R_{\+ Z^\proxpar} \subseteq\+ Y^r) $                     $\+ R^a = \emptyset$

The conditions for the enumerated subset $\+ Z$ for applying the rules of Table <ref> to a term in the input column. Here $\+ T = \+ Y \cup \+ X \cup \+ W$ and the sets $\+ Y, \+ X$ and $\+ W$ are those present in the input term $P(\+ Y \given \doo(\+ X), \+ W)$. Active response indicators of the input are denoted by $\+ R^a$. For rules $6\pm, 7\pm$ and $8\pm$, the conditions specify valid variables of the second required term. Validity conditions for rules $1\pm,2\pm,3\pm,4,5$ and $6\pm$ are the same as in Table <ref>.

Additional termination conditions also apply to the new rules and their correctness is easily verified.
Let $G$ be a semi-Markovian graph and let $\+ Y, \+ X$ and $\+ W$ be disjoint subsets of $\+ V$. Then all of the following are true:

(i)If $\+ |\+ Y| = 1$, then rules $7\pm$ of Table <ref> cannot be used.

(ii)If $\+ W = \emptyset$, then rule $8-$ of Table <ref> cannot be used.

(iii)If $\+ W^r = \emptyset$ then rule $9+$ of Table <ref> cannot be used.

(iv)If $\+ Y^r = \emptyset$ then rule $9-$ of Table <ref> cannot be used.

(v)If $\+ R^a = \emptyset$, then rules $10\pm$ of Table <ref> cannot be used.

For (i), the set $\+ Y$ only has a single vertex, so it cannot have a non-empty subset. For (ii), the set $\+ W$ is empty so no subset $\+ Z \subset \+ W$ can exist for the second input. For (iii), and the set $\+ W^r$ is empty so no assignment to value $1$ can be performed. Similarly for (iv), the set $\+ Y^r$ is empty so no assignment to value $1$ can be performed. For (v) the set of active response indicators $\+ R^a$ is empty, so no transformation from proxy variables to true variables via the missingness mechanism in equation (<ref>) can take place.

The task of selecting a suitable heuristic becomes more difficult when missing data are involved with the identifiability problem. The heuristic approach of Section <ref> is no longer directly applicable due to the relation between proxy variables, response indicators and true variables. The proximity function considers $X$ and $X^*$ as entirely different variables despite their connection and does not prefer the inclusion of response indicators. If the heuristic is applied as such, the search path will often involve a large number of manipulations which in turn leads to complicated expressions. For these reasons we do not apply a heuristic to missing data problems, but expand terms in the order in which they were identified. The improvements described in Section <ref> still apply.

It is straightforward to adapt  to the new extended set of rules. In the pseudo-code shown in Algorithm <ref>, we simply replace all references to Tables <ref> and <ref> by references to Tables <ref> and Tables <ref>, respectively. When the validity condition is checked, we also verify that the missing data condition holds. Lemma <ref> guarantees the correctness of the new termination criteria. Theorem <ref> is also valid when the sets $\+ A_i$ are replaced by $\+ A_i \cup \+ A_i^\proxpar$, since it may be possible to exchange some proxy variable to a true variable that is present in the set $\+ Y$ of the target $P(\+ Y \given \doo(\+ X), \+ W)$.

§ THE DOSEARCH PACKAGE

We implemented  (Algorithm <ref>) in C++ and constructed an R interface using the Rcpp package [Eddelbuettel and Francois, 2011]. This interface is provided by the R package dosearch. Calling the search from R is straightforward via the primary function that carries the name of package.

dosearch(data, query, graph,
transportability, selection_bias, missing_data,

The required inputs of the function are data, query and graph. Parameter data is used to encode the set $\+ P$ of known input distributions of Algorithm <ref> as a character string, where each distribution is separated by a new line. For example, if we have access to a set of  distributions $\+ P = \{P(W), P(Y \given X), P(Z \given \doo(X), W)\}$, we would write

R> data <- "
+    P(W)
+    P(Y|X)
+    P(Z|do(X),W)
+  "

The $\doo(\cdot)$-operator can either precede or succeed conditioning variables, but it must appear only once in a given term, meaning that expressions such as P(Y|do(A),B,do(C)) are not allowed, but should instead be given as P(Y|B,do(A,C)) or P(Y|do(A,C),B). If variable sets are desired, each member of the set has to be included explicitly.

Parameter query is used to describe the target $Q$ of Algorithm <ref> as a character string, similarly as the data. If we are interested in identifying $P(Y \given \doo(X), W)$ we would write

R> query <- "P(Y|do(X),W)"

Instead of describing distributions via text, it is also possible to use the following structure that encodes the role of each variable via a numeric vector:

R> query <- c(Y = 0, X = 1, W = 2)

Given a distribution of the form $P(\+ A \given \doo(\+B),\+ C)$ and a variable $V$, a value 0 means that $V \in \+ A$, value 1 means that $V \in \+ B$ and value 2 means that $V \in \+ C$. This format can also be used to input data as a list of numeric vectors:

R> data <- list(
+    c(W = 0),
+    c(Y = 0, X = 2),
+    c(Z = 0, X = 1, W = 2)
+  )

Finally, graph encodes the semi-Markovian graph $G$ of the causal model as a character string with each edge on its own line. A directed edge from $X$ to $Y$ is given as X -> Y and a bidirected edge between $X$ and $Y$ is given as X <-> Y. Intervention nodes should not be given explicitly, since they are added automatically after calling dosearch. Furthermore, only vertices with incoming or outgoing edges should be included in graph. A variable with no connected edges can still appear in the input distributions and is automatically added to the graph as well. As an example, we can encode the graph of Figure <ref>(fig:intro_backdoor) with an added bidirected edge between $X$ and $Y$ as follows

R> graph <- "
+    X -> Y
+    Z -> X
+    Z -> Y
+    X <-> Y
+  "

Alternatively, one may use igraph graphs [Csardi and Nepusz, 2006] in the syntax of the causaleffect package or DAGs created using the dagitty package.

R> library("igraph")
R> graph <- graph.formula(X -+ Y, Z -+ X, Z -+ Y, X -+ Y, Y -+ X)
R> graph <- set.edge.attribute(graph, "description", 4:5, "U")

R> library("dagitty")
R> graph <- dagitty("dagX -> Y; Z -> X; Z -> Y; X <-> Y")

The next two optional parameters of dosearch, transportability and selection_bias, are used to denote those vertices of $G$ that should be understood as either transportability nodes or selection bias nodes, respectively. Providing these parameters may increase search performance in relevant problems. Both of these parameters should be given as character strings, where individual variables are separated by a comma, for example transportability = "S,T". Parameter missing_data, as the name suggests, is used to define missingness mechanisms (<ref>) as a character string, where individual mechanisms are separated by a comma. In order to describe that $R_X$ is the response indicator of $X$ we would write R_X : X, which also implicitly defines that X* is the proxy variable of X. Proxy variables do not need to be manually described in graph as they are automatically constructed based on the missing_data argument.

The list control can be used to set various additional parameters that are not directly related to the identifiability problem itself, but more so to the output of the search and other auxiliary details, such as benchmarking and obtaining derivations such as Figure <ref>(fig:derivation). One such control parameter determines whether to use the search heuristic or not (heuristic = TRUE by default for problems without missing data). Documentation of the dosearch package contains detailed information on the full list of control parameters.

The return object of dosearch is a list with three components by default. The first component, identifiability, is a logical value that takes the value TRUE when the target distribution described by query is identifiable from the inputs of data. The second component, formula, is a character string describing the target distribution in terms of the inputs in LaTeX syntax if the target is identifiable. Otherwise this component is just an empty character string. The third component call contains the arguments of the original function call.

§.§ Simulations

Here we report the results of a simulation study to assess the running time performance of  and the impact of the search space reduction techniques as well as the search heuristic outlined in Section <ref>.

Our synthetic simulation scenario consisted of $1071$ semi-Markovian causal graphs of $10$ vertices that were generated at random by first generating a random topological order of the vertices followed by a random lower triangular adjacency matrices for both directed and bidirected edges[The computation took place in a cluster with $9$ parallel jobs, each with $120$ tasks where the supervisor task only collated the results from the 119 workers resulting in $9 \cdot (120 - 1) = 1071$ graphs in total.]. Graphs without a directed path from $X$ to $Y$ were discarded. We sampled sequentially input distributions of the form $P(\+ A \given \doo(\+ B), \+ C)$ at random by generating disjoint subsets such that $\+ A$ is always non-empty. This was continued until the target quantity $P(Y \given \doo(X))$ was found to be identifiable by the search. Then for each graph, we recorded the search times for the specific set of inputs that first resulted in the query to be identified and for the last set such that the target was non-identifiable. In other words, each graph generates two simulation instances, one for an identifiable query and one for a non-identifiable query. This setting directly corresponds to the setting of partially overlapping experimental data sets discussed in Section <ref> for which no other algorithmic solutions exist.

To understand the impact of the search heuristic and the various improvements, we compare four different search configurations: the basic  without the search heuristic or improvements[In this configuration, terms are expanded in the order they were identified; the conditions in Table <ref> are not checked.], one that only uses the search heuristic, one that only uses the improvements of Section <ref> and one that uses them both.

Scatter plots of the search times from identifiable instances under different search configurations compared to the basic  without a heuristic or improvements.

Figure <ref> shows the search times of the configurations compared to the basic configuration for identifiable instances. Most importantly, a vast majority of instances (93 %) are solved faster than the basic configuration when both heuristics and improvements are used. The average search time with both heuristics and improvements enabled was 32.7 seconds and 75.2 seconds for the basic configuration. The search heuristic provides the greatest benefit for these instances as can be seen from Figure <ref>(fig:scatter_id_h). Using a heuristic can sometimes hinder performance by leading the search astray and by causing additional computational steps through the evaluation of the proximity function. For example, there is a small number of instances where the search is over ten times slower than the basic configuration when using a heuristic. Fortunately, there are several instances in the opposite direction, where the heuristic provides over one hundred fold reduction in search time. Curiously, even using the improvements sometimes results in slower search times. This is most likely due to the elimination of rule 1 of do-calculus, since it may be the case that the basic search is able to use this rule to reach the target distribution faster. More importantly, Figure <ref>(fig:scatter_id_i) shows that the improvements clearly benefit the search. Furthermore, the benefit tends to increase as the instances get harder.

Figure <ref> shows the search times of the configurations for non-identifiable instances. Relying only on a search heuristic provides no benefit here, as expected. The improvements to the search are most valuable for these instances, and in this scenario every non-identifiable instance was solved faster than baseline using the improvements, and when applied with the heuristic only three instance were slower than baseline. The average search time with both heuristic and improvements enabled was 105.2 seconds and 139.7 seconds for the basic configuration. The almost zero second instances are a result of Theorem <ref> when no search has to be performed in order to determine the instance to be non-identifiable. The benefit of the improvements tends to increase as the instances get harder also for these instances.

Scatter plots of the search times from non-identifiable instances under different search configurations compared to the baseline configuration.

Boxplots of search times for both identifiable and non-identifiable instances in graphs of $n = 4,\ldots,10$ vertices. The vertical axis uses a logarithmic scaling. Instances where the search time was less than $10^{-5}$ seconds were omitted for clarity.

Finally we examined the average run time performance of , with all improvements and heuristics enabled. We replicated the previously described simulation scenario with the same number of instances (1071) for graphs up to 10 vertices. Figure <ref> shows the boxplots of search times on a log-scale for graphs of different size, including both identifiable and non-identifiable instances. Note that for every graph size there are a number of easily solvable instances that show up as outliers in this plot. Instances with graphs of 10 vertices are solved routinely under 100 seconds. In this plot, the running times increase exponentially with increasing graph size (or number of variables).

§ NEW CAUSAL EFFECT IDENTIFICATION RESULTS

We present a number of results for various identifiability problems to showcase the versatility of  with the accompanying R code for some specific examples.

§.§ Multiple Data Sources with Partially Overlapping Variable Sets

Earlier generalizations of the identifiability problem assume nested experiments or entire distributions with the exception of surrogate outcome identifiability [Tikka and Karvanen, 2019] which also has its own intricate set of assumptions regarding the available distributions. None of these assumptions are needed in  and it can be used to solve identifiability problems from completely arbitrary collections of input distributions.

We showcase identifiability from multiple data sources by three examples. The first example is the human resource management problem presented in the introduction and shown in Figure <ref>. The question of interest was the identifiability of $P(Y \given \doo(X))$ from the data sources $P(Y,B,E,X)$ and $P(A,B,X)$. The answer can be provided with the following lines of R code:

R> library("dosearch")
R> data <- "
+    p(y,b,e,x)
+    p(a,b,x)
+  "
R> query <- "p(y|do(x))"
R> graph <-  "
+    e -> x
+    e -> y
+    a -> b
+    a -> x
+    x -> b
+    x -> y
+    b -> y
+  "
R> dosearch(data, query, graph)

The result means that the causal effect is identifiable and the returned formula is

\[
\sum_{B,A}P(A)P(B \given X,A)\sum_{E} P(E)P(Y|X,B,E).
\]

By running an additional line of code

R> dosearch(data, "p(y,b,e,x,a)", graph)

The query p(y,b,e,x,a) is non-identifiable.

we learn that the joint distribution $P(Y,B,E,X,A)$ is not identifiable. This rules out the possibility to solve the problem with causaleffect which requires the joint distribution as an input.

In the second example we consider identifiability of $P(Y_1,Y_2 \given \doo(X_1,X_2))$ in the graph of Figure <ref>(fig:obsandexp_1) from $P(\+ V)$, $P(Y_1,Y_2 \given \doo(X_1),Z,W,X_2)$, $P(W \given \doo(X_1, X_2))$ and $P(Z \given \doo(X_2))$. The target quantity is identifiable and  produces the following formula for it

\[
\sum_{Z,W} P(Y_1,Y_2 \given \doo(X_1), Z, W, X_2)P(Z \given \doo(X_2))P(W \given \doo(X_2, X_1))
\]

In the third example we consider identifiability of $P(Y_1,Y_2 \given \doo(X_1,X_2))$ in the graph of Figure <ref>(fig:obsandexp_2) from $P(\+ V)$, $P(Y_1 \given \doo(X_1),Y_2,W,Z,X_2)$, $P(X_2,W \given \doo(X_1))$, $P(X_2 \given \doo(X_1,W))$, $P(Y_2 \given \doo(X_1),Z,W,X_2)$, $P(Y_2 \given \doo(Z),X_1,W,X_2)$, and $P(Y_1,Y_2 \given \doo(Z),W,X_1,X_2)$. Again, the target quantity is identifiable and  outputs the following formula

\begin{align*}
\sum_{W}&\left(P(W \given \doo(X_1),X_2) \sum_{X_2} P(X_2 \given \doo(X_1,W)) \times \right. \\
&\left.\frac{\sum_{Z}\! P(X_2,W,Z \given X_1)P(Y_1,Y_2 \given \doo(X_1),X_2,W,Z)}{\sum_{Y^\prime_1,Y^\prime_2,Z}\! P(X_2,W,Z \given X_1)P(Y^\prime_1,Y^\prime_2 \given \doo(X_1),X_2,W,Z)} \right).
\end{align*}

This example shows that a heuristic approach can also help us to find shorter formulas. If we run  again without the heuristic in this instance, the output formula is instead

\begin{align*}
&\sum_{W,Z} \left(P(Z)P(W \given X_2,X_1,Z) \sum_{X_2} P(X_2 \given X_1,Z) \sum_{Y_2} P(Y_2 \given \doo(X_1),X_2,W,Z)\right. \, \times \\
&\quad\left. P(Y_1 \given \doo(X_1),X_2,Y_2,W,Z) \frac{P(Y_2 \given \doo(X_1),X_2,W,Z)P(Y_1 \given \doo(X_1),X_2,Y_2,W,Z)}{\sum_{Y^\prime_2} P(Y^\prime_2 \given \doo(X_1),X_2,W,Z)P(Y_1 \given \doo(X_1),X_2,Y^\prime_2,W,Z)} \right).
\end{align*}

We can run these examples in R by writing

R> data <- "
+    p(x_1,y_1,x_2,y_2,z,w)
+    p(y_1,y_2|z,w,x_2,do(x_1))
+    p(y_2|y_1,z,w,x_2,do(x_1))
+    p(w|do(x_1,x_2))
+    p(z|do(x_2))
+  "
R> query <- "p(y_1,y_2|do(x_1,x_2))"
R> graph <- "
+    z -> y_1
+    w -> y_1
+    y_1 -> y_2
+    x_2 -> z
+    x_1 -> w
+    y_1 <-> x_1
+    y_1 <-> y_2
+    y_2 <-> z
+    y_1 <-> w
+    y_2 <-> w
+  "
R> dosearch(data, query, graph)

R> data <- "
+    p(x_1,y_1,x_2,y_2,z,w)
+    p(y_1,y_2|w,x_1,x_2,do(z))
+    p(y_1|y_2,w,z,x_2,do(x_1))
+    p(y_2|x_1,w,x_2,do(z))
+    p(x_2,w|do(x_1))
+    p(x_2|do(x_1,w))
+    p(y_2|z,w,x_2,do(x_1))
+  "
R> query <- "p(y_1,y_2|do(x_1,x_2))"
R> graph <- "
+    y_2 -> y_1
+    w -> y_1
+    x_1 -> x_2
+    x_1 -> y_2
+    z -> y_2
+    w -> y_2
+    x_2 -> w
+    x_1 <-> y_1
+    x_1 <-> y_2
+    y_1 <-> z
+    x_2 <-> z
+  "
R> dosearch(data, query, graph)

∑_z(p(x_2,w,z|x_1)p(y_1,y_2|do(x_1),x_2,w,z))/∑_y_1,y_2 ∑_z(p(x_2,w,z|x_1)p(y_1,y_2|do(x_1),x_2,w,z)

[dot = 00X_1above left] at (0,0) ;
[dot = 00X_2above left] at (1,0) ;
[dot = -0.1-0.1Y_1above right] at (2,-1) ;
[dot = 00Y_2above right] at (3,-1.5) ;
[dot = 00Zabove] at (2,-0.5) ;
[dot = 0.10.1Wbelow left] at (1,-0.5) ;

[->] (X_1) – (W);
[->] (X_2) – (Z);
[->] (W) – (Y_1);
[->] (Z) – (Y_1);
[->] (Y_1) – (Y_2);

[<->,dashed] (X_1) to [bend right=50] (Y_1);
[<->,dashed] (W) to [bend right=50] (Y_2);
[<->,dashed] (W) to [bend left=30] (Y_1);
[<->,dashed] (Y_1) to [bend right=30] (Y_2);
[<->,dashed] (Y_2) to [bend right=30] (Z);

[dot = 00X_1below left] at (0,0) ;
[dot = 00X_2above left] at (1,1.5) ;
[dot = 00Y_1above left] at (2,1) ;
[dot = 00Y_2below right] at (2,0) ;
[dot = 00Zright] at (2.5,0.5) ;
[dot = 00Wleft] at (1.5,0.5) ;

[->] (X_2) – (W);
[->] (X_1) – (Y_2);
[->] (X_1) – (X_2);
[->] (W) – (Y_1);
[->] (W) – (Y_2);
[->] (Y_2) – (Y_1);
[->] (Z) – (Y_2);

[<->,dashed] (Y_2) to [bend right=30] (X_1);
[<->,dashed] (X_2) to [bend left=50] (Z);
[<->,dashed] (Y_1) to (Z);
[<->,dashed] (Y_1) to [bend right=20] (X_1);

Graphs for the examples on identifiability problems combining both observational and experimental distributions.

§.§ Combining Transportability and Selection Bias

Input distributions that originate from multiple sources while being simultaneously affected by selection bias can be considered with . This kind of problem cannot be solved with algorithms RC or TR$^{\textrm{mz}}$ of Table <ref>. As an example we consider one source domain and a target domain with two input data sets: a biased distribution $P(X,Y,Z \given S)$ from the target domain and an unbiased experimental distribution $P(Y,Z \given \doo(X), T)$ from the source domain. We evaluate the query $P(Y \given \doo(X))$ in the graph of Figure <ref> using these inputs. In the figure transportability node $T$ is depicted as a gray square and selection bias node $S$ is depicted as an open double circle.

[dot = 00Xbelow left] at (0,0) ;
[dot = 00Zabove] at (1,0) ;
[dot = 00Ybelow right] at (2,0) ;
[sb = 00Sbelow] at (0,-0.5) ;
[black, inner sep = 0.5pt] (0,-0.5) circle (1.15pt);
[tr = 00Tbelow] at (1,-0.5) ;
[->] (X) – (Z);
[->] (Z) – (Y);
[->] (X) – (S);
[->] (T) – (Z);
[<->,dashed] (X) to [bend left=50]  (Y);

Graph that contains both selection bias and transportability nodes.

The query is identifiable and  outputs the following formula for it

\[
P(Y \given \doo(X)) = \sum_{Z} P(Y \given \doo(X),Z,T) \sum_{Y^\prime} P(Z,Y^\prime \given X,S).
\]

In R we may write

R> data <- "
+    p(x,z,y|s)
+    p(y,z|t,do(x))
+  "
R> query <- "p(y|do(x))"
R> graph <- "
+    x -> z
+    z -> y
+    x -> s
+    t -> z
+    x <-> y
+  "
R> dosearch(data, query, graph,
+    transportability = "t", selection_bias = "s")

§.§ Recovering from Multiple Sources of Selection Bias

We present an example where bias originates from two sources with two input data sets: a distribution affected by both biasing mechanisms $P(X,Y,Z,W_1,W_2 \given S_1,S_2)$ and a distribution affected only by a single bias source $P(Z \given S_1)$. We evaluate the query $P(Y \given \doo(X))$ in the graph of Figure <ref> using the inputs.

[dot = 00Xbelow] at (0,0) ;
[dot = 00Zabove] at (2,2) ;
[dot = 00Ybelow right] at (3,0) ;
[dot = 00W_2below right] at (1,1) ;
[dot = 00W_1above left] at (-0.5,1.5) ;
[sb = 00S_1below] at (-1,0) ;
[black, inner sep = 0.5pt] (-1,-0) circle (1.40pt);
[sb = 00S_2below] at (2,1) ;
[black, inner sep = 0.5pt] (2,1) circle (1.40pt);
[->] (X) – (Y);
[->] (W_2) – (X);
[->] (Z) – (W_2);
[->] (Z) – (Y);
[->] (Z) – (S_2);
[->] (W_1) – (X);
[->] (W_1) – (W_2);
[->] (W_1) – (S_1);

Graph where two selection bias nodes are present.

The query is identifiable and the following formula is obtained by

\[
\sum_{Z} P(Z \given S_1)P(Y \given X,Z,W_1,W_2,S_1,S_2).
\]

The result can be obtained in R as follows. In this case, slightly nicer expression is produced when not using the heuristic (control = list(heuristic = FALSE)) than when using the heuristic.

R> data <- "
+    p(x,y,z,w_1,w_2|s_1,s_2)
+    p(z|s_1)
+  "
R> query <- "p(y|do(x))"
R> graph <- "
+    w_1 -> w_2
+    z -> w_2
+    x -> y
+    z -> y
+    z -> s_2
+    w_1 -> x
+    w_2 -> x
+    w_1 -> s_1
+  "
R> dosearch(data, query, graph, selection_bias = "s_1, s_2",
+    control = list(heuristic = FALSE))

§.§ Systematic Analysis of Bivariate Missing Data Problems

We apply  using the extended rule set of Table <ref> for all identifiability problems in bivariate missingness graphs. By bivariate missingness graphs we mean semi-Markovian graphs for two variables, $X$ and $Y$, and their missingness indicators, $R_X$ and $R_Y$. Noting that edges from $\{R_X, R_Y \}$ to $\{X, Y \}$ are not allowed, there are 9216 such graphs. We consider only 6144 graphs of which 3072 have the edge $X \rightarrow Y$ and 3072 do not have an edge between $X$ and $Y$. Graphs with the edge  $Y \rightarrow X$ are obtained from the studied graphs by swapping the roles of $X$ and $Y$. The maximum number of edges in a bivariate missingness graph is 12 (when a bidirected edge is counted as a single edge).

The available theoretical results for missing data problems include a theorem by Mohan et al., 2013 that gives a sufficient and necessary condition for the identifiability of the joint distribution $P(\+ V)$ but is restricted to graphs that do not have edges between the missingness indicators (row 9 of Table <ref>). In our example, 5120 graphs out of 6144 have such edges. The algorithm by Shpitser et al., 2015 does not have this restriction but it is not complete as shown by Bhattacharya et al., 2019 (row 10 of Table <ref>). It follows from the results of Bhattacharya et al., 2019 that the rules of Table <ref> are not complete for missing data problems. Differently from the theorem by Mohan et al., 2013 and the algorithms by Shpitser et al., 2015 and Bhattacharya et al., 2019,  can however address missing data problems where we consider identification of a marginal or a conditional distribution. In addition,  can address missing data problems with multiple input distributions.

The queries $P(X,Y)$, $P(X)$, $P(Y)$, $P(Y \given X)$ and $P(Y \given \doo(X))$ were evaluated using  in these 6144 graphs with the input distribution $P(X^*, Y^*, R_X, R_Y)$. The results are summarized by Venn diagrams in Figure <ref>. The results are also available as a data set  in the R package dosearch. Using this data set we are able to showcase examples on non-identifiability and find interesting special cases by direct evaluation of all possible bivariate missingness graphs. The first example relates non-identifiability to the of the number of edges present in the graph.

Let $K$ denote the number of edges in a bivariate missingness graph that does not have edge $Y \rightarrow X$. The joint distribution $P(X,Y)$ is not identifiable by  if $K>5$, marginal distribution $P(X)$ is not identifiable by  if , marginal distribution $P(Y)$ and conditional distribution $P(Y \given X)$ are not identifiable by  if .

The next example specifies the graph with the largest number of edges where both the joint distribution of $X$ and $Y$ and the causal effect of $X$ on $Y$ can be identified.

The graph in Figure <ref>(fig:jointanddo5edges) is the
only bivariate missingness graph that (i) has edge $X \rightarrow Y$, (ii) has five edges, and (iii) allows for the identification of $P(X,Y)$ and $P(Y \given \doo(X))$ by .

The third example specifies the graph with the largest number of edges where the marginal distributions are identifiable while the joint distribution and the causal effect of $X$ on $Y$ are non-identifiable.

The graph in Figure <ref>(fig:marginalsonly5edges) is the
only bivariate missingness graph that (i) has five edges, and (ii) allows for the identification of $P(X)$ and $P(Y)$, and (iii) does not allow for the identification of $P(X,Y)$ or $P(Y \given \doo(X))$ by . No bivariate missingness graph that has more than five edges fulfills the conditions (ii) and (iii).

Graphs with arrow $X \rightarrow Y$.

Graphs where there is neither $X \rightarrow Y$ nor $Y \rightarrow X$.

Venn diagrams indicating the number of graphs were different distributions can be identified . The intersection of $P(X)$ and $P(Y \given X)$ shows the number of graphs were $P(X,Y)$ can be identified. The total number of possible graphs is 3072 in both cases.

Some interesting examples are shown in Figure <ref>. Graphs (a) and (b) are the unique graphs that fulfill the conditions specified in Examples <ref> and <ref>, respectively. Graph (c) is the graph with the smallest number of edges where marginals $P(X)$ and $P(Y)$ can be identified but the joint distribution $P(X,Y)$ or causal effect $P(Y \given \doo(X))$ cannot be identified by . In graph (d), $P(X)$, $P(Y)$, $P(X,Y)$ and $P(Y \given \doo(X))$ are not identifiable by  but the conditional distribution $P(Y \given X)$ can be identified as follows

\begin{equation} \label{eq:cross1and2}
P(Y \given X) = \frac{ P(Y \given R_Y=1) P(X \given Y, R_X=1, R_Y=1) } {\sum_{Y^\prime} P(Y^\prime \given R_Y=1) P(X \given Y^\prime, R_X=1, R_Y=1)}.
\end{equation}

In equation (<ref>), the numerator resembles the joint distribution $P(X,Y \given R_X = 1,R_Y = 1)$ but is different because $Y$ and $R_X$ are not independent. The denominator is the marginal of this pseudo joint distribution.
In graph (e), $P(X)$, $P(Y)$ and $P(X,Y)$ are not identifiable by  but $P(Y \given X)$ and $P(Y \given \doo(X))$ are identifiable and can be both estimated with equation (<ref>). In graph (f), $P(X,Y)$, $P(X)$ and $P(Y \given \doo(X))$ are not identifiable by  but $P(Y)$ and $P(Y \given X)$ can be identified as follows

\begin{align}
P(Y) &= \sum_{R_X,X^*} P(Y\given X^*, R_X, R_Y=1) P(R_X,X^*), \label{eq:openboxY}\\
P(Y \given X) &= P(Y \given X, R_X=1, R_Y=1) \nonumber
\end{align}

In equation (<ref>), the summation also goes over the cases where $X^* = \textrm{NA}$ and the distribution of $Y$ must be estimated also on the condition that $X$ is not observed.

[dot = 00Xbelow left] at (0,0) ;
[dot = 00Ybelow right] at (2,0) ;
[dot = 00R_Xabove left] at (0,2) ;
[dot = 00R_Yabove right] at (2,2) ;
[->] (R_X) – (R_Y) ;
[->] (X) – (R_Y) ;
[->] (X) – (Y) ;
[<->,dashed] (X) to [bend right=30] (R_Y) ;
[<->,dashed] (R_X) to [bend left=30] (R_Y) ;

[dot = 00Xbelow left] at (0,0) ;
[dot = 00Ybelow right] at (2,0) ;
[dot = 00R_Xabove left] at (0,2) ;
[dot = 00R_Yabove right] at (2,2) ;
[->] (R_Y) – (R_X) ;
[->] (X) – (Y) ;
[<->,dashed] (X) to [bend right=30] (Y) ;
[<->,dashed] (Y) to [bend left=30] (R_X) ;
[<->,dashed] (R_X) to [bend left=30] (R_Y) ;

[dot = 00Xbelow left] at (0,0) ;
[dot = 00Ybelow right] at (2,0) ;
[dot = 00R_Xabove left] at (0,2) ;
[dot = 00R_Yabove right] at (2,2) ;
[->] (X) – (Y) ;
[<->,dashed] (Y) to [bend left=30] (R_X) ;

[dot = 00Xbelow left] at (0,0) ;
[dot = 00Ybelow right] at (2,0) ;
[dot = 00R_Xabove left] at (0,2) ;
[dot = 00R_Yabove right] at (2,2) ;
[->] (R_Y) – (R_X) ;
[->] (X) – (R_Y) ;
[->] (Y) – (R_X) ;
[->] (X) – (Y) ;
[<->,dashed] (X) to [bend right=30] (Y) ;

[dot = 00Xbelow left] at (0,0) ;
[dot = 00Ybelow right] at (2,0) ;
[dot = 00R_Xabove left] at (0,2) ;
[dot = 00R_Yabove right] at (2,2) ;
[->] (R_Y) – (R_X) ;
[->] (X) – (R_Y) ;
[->] (Y) – (R_X) ;
[->] (X) – (Y) ;
[<->,dashed] (X) to [bend right=30] (R_Y) ;

[dot = 00Xbelow left] at (0,0) ;
[dot = 00Ybelow right] at (2,0) ;
[dot = 00R_Xabove left] at (0,2) ;
[dot = 00R_Yabove right] at (2,2) ;
[->] (R_X) – (R_Y) ;
[->] (X) – (R_X) ;
[->] (X) – (Y) ;
[<->,dashed] (X) to [bend right=30] (Y) ;

Missingness graphs used as example cases. Proxy variables are omitted for clarity.

§.§ Causal Inference under Case-control Design

Case-control design [Breslow, 1996] is commonly used in epidemiology to study risk factors of rare diseases. In the basic setup, a fixed number of disease cases and a fixed number of controls are selected for the risk factor measurements. When the disease is rare, this design leads to substantial savings in the sample size compared to simple random sampling. Figure <ref>(fig:simplecc) shows the missingness graph for a situation where the inclusion to the study (indicator $R_Y$) depends on the disease endpoint $Y$. The risk factors $X$ are measured for the subset $R_Y=1$ but occasionally the values are missing (indicator $R_X$). It is immediately seen that neither the causal effect $P(Y \given \doo(X))$ nor conditional distribution $P(Y \given X)$ can be identified because of the edge $Y \rightarrow R_Y$. However, if the prevalence of the disease in the population, i.e., the marginal distribution $P(Y)$, is known, the causal effect $P(Y \given \doo(X))$ can be identified. The result is provided by

\begin{equation} \label{eq:simplecc}
P(Y \given \doo(X)) = \frac{P(Y)P(X \given Y, R_Y=1, R_X=1)}{\sum_{Y^\prime} P(Y^\prime)P(X \given Y^\prime, R_Y=1, R_X = 1)}.
\end{equation}

In typical applications response $Y$ is binary but in the non-parametric formula of equation (<ref>) response can be discrete or continuous.

[dot = 00Xbelow] at (0,0) ;
[dot = 00Ybelow] at (1.6,0) ;
[dot = 00R_Yleft] at (1.8,1) ;
[dot = 00R_Xleft] at (2.0,2) ;
[->] (X) – (Y) ;
[->] (Y) – (R_Y) ;
[->] (R_Y) – (R_X) ;

Basic case-control design.

[dot = 00Xbelow] at (0,0) ;
[dot = 00Zbelow] at (1.6,0) ;
[dot = 00Ybelow] at (3.2,0) ;
[dot = -0.10R_Yleft] at (3.4,1) ;
[dot = -0.10R_Xleft] at (3.6,2) ;
[dot = 00R_Zright] at (4.6,2) ;
[->] (X) – (Z) ;
[->] (Z) – (Y) ;
[->] (Y) – (R_Y) ;
[->] (R_Y) – (R_X) ;
[->] (R_Y) – (R_Z) ;
[<->,dashed] (X) to [bend left=30] (Y) ;
[<->,dashed] (R_X) to [bend right=30] (R_Y) ;
[<->,dashed] (R_X) to [bend left=30] (R_Z) ;
[<->,dashed] (R_Z) to [bend left=30] (R_Y) ;

Case-control design for the front-door situation.

Missingness graph for the case-control examples.

A more complicated example is shown in Figure <ref>(fig:frontdoorcc) where the causal effect of risk factor $X$ on disease endpoint $Y$ fulfills the front-door criterion [Pearl, 1995] with respect to mediator $Z$ and the data are collected from a case-control design where the selection depends $Y$ and there is occasional item non-response in $X$ and $Z$. We observe data $P(Y^*,X^*,Z^*,R_Y,R_X,R_Z)$ and know the marginal distribution $P(Y)$ from other sources. Applying  we obtain the result

\begin{equation}
\begin{aligned}
&  P(Y \given \doo(X)) = \\
&\quad \sum_Z \left[ \frac{\sum_{Y^\prime} P(Y^\prime)P(X,Z \given Y^\prime,R_X = 1, R_Y = 1, R_Z = 1)}{\sum_{Z^\prime,Y^\prime} P(Y^\prime)P(X,Z^\prime \given Y^\prime,R_X = 1, R_Y = 1, R_Z = 1)}\, \right. \times \\
&\qquad \sum_{X^\prime} \left( \sum_{Y^\prime,Z^\prime} P(Y^\prime)P(X^\prime,Z^\prime \given Y^\prime,R_X = 1, R_Y = 1, R_Z = 1)\, \right. \times \\
&\quad\qquad \left. \left. \frac{P(Y)P(X^\prime,Z \given Y, R_X = 1, R_Y = 1,R_Z = 1)}{\sum_{Y^\prime} P(Y^\prime)P(X^\prime,Z \given Y^\prime,R_X = 1,R_Y = 1,R_Z = 1)} \vphantom{\sum_{Z^\prime}}\right) \right].
\end{aligned}
\label{eq:frontdoorcc}
\end{equation}

Expression (<ref>) follows the general structure of the front-door adjustment

\[
P(Y \given \doo(X)) = \sum_Z P(Z \given X) \sum_{X^\prime} P(X^\prime) P(Y \given X^\prime,Z),
\]

\begin{align*}
P(Z \given X) &= \frac{\sum_{Y^\prime} P(Y^\prime)P(X,Z \given Y^\prime,R_X = 1, R_Y = 1, R_Z = 1)}{\sum_{Z^\prime,Y^\prime} P(Y^\prime)P(X,Z^\prime \given Y^\prime,R_X = 1, R_Y = 1, R_Z = 1)}, \\
P(X) &= \sum_{Y^\prime,Z^\prime} P(Y^\prime)P(X,Z^\prime \given Y^\prime,R_X = 1, R_Y = 1, R_Z = 1), \\
P(Y \given X,Z) &= \frac{P(Y)P(X,Z \given Y,R_X = 1,R_Y = 1,R_Z = 1)}{\sum_{Y^\prime} P(Y^\prime)P(X,Z \given Y^\prime,R_X = 1,R_Y = 1,R_Z = 1)}.
\end{align*}

Note that $P(X,Y,Z) = P(Y)P(X,Z \given Y,R_X = 1,R_Y = 1,R_Z = 1)$. In [Karvanen, 2015], a similar example was studied assuming that $X$, $Z$ and $Y$ are binary but in expression (<ref>) there are no such restrictions. This factorization can be obtained in R as follows

R> data <- "
+    p(x*,y*,z*,r_x,r_y,r_z)
+    p(y)
+  "
R> graph <- "
+    x -> z
+    z -> y
+    y -> r_y
+    x <-> y
+    r_y -> r_x
+    r_y -> r_z
+    r_y <-> r_x
+    r_y <-> r_z
+    r_z <-> r_x
+  "
R> md <- "r_x : x, r_y : y, r_z : z"
R> query1 <- "p(z|x)"
R> query2 <- "p(x)"
R> query3 <- "p(y|x,z)"
R> dosearch(data, query1, graph, missing_data = md)
R> dosearch(data, query2, graph, missing_data = md)
R> dosearch(data, query3, graph, missing_data = md)

∑_y(p(y)p(x,z|r_x = 1,y,r_y = 1,r_z = 1))/∑_z ∑_y(p(y)p(x,z|r_x = 1,y,r_y = 1,r_z = 1))

∑_y,z(p(y)p(x,z|r_x = 1,y,r_y = 1,r_z = 1))

(p(y)p(x,z|r_x = 1,y,r_y = 1,r_z = 1))/∑_y (p(y)p(x,z|r_x = 1,y,r_y = 1,r_z = 1))

§ DISCUSSION

The presented algorithm, , removes the need for manual application of do-calculus, which is time-consuming and prone to errors. Systematic analyses such as the one in Section <ref> are practically unreachable with manual application of do-calculus. Superiority of  over a simple forwards breadth-first search was attained through a combination of a search heuristic and a reduction of the search space. Some further approaches were attempted but later discarded as non-beneficial. These include caching separation criteria that hold in the graph after they are first evaluated, pre-computing valid subsets for each subset size and enumerating subsets in an order of increasing cardinality.

As the simulations showed, our intuitive heuristic yielded significant improvements in search performance. The proximity function defined in Section <ref> uses only the information contained in the distributions themselves. One approach could be to also take the structure of the graph into account in the proximity function. Further study is needed for finding a heuristic that performs well when missing data mechanisms are present in the graph.

The scalability of  is limited due to vast search space of possibly identified causal effects. Currently, algorithms with polynomial complexity currently exist only for the simpler problems (see Table <ref>). However, based on the simulation results,  solves identifiability problems in graphs of ten vertices in under two minutes on average. By our observation, graphs typically analyzed in literature related to identifiability problems have fewer vertices. The theoretical computational complexity of the general form of the causal identifiability problem defined in Section <ref> remains an important and interesting question.

The search could also be used to obtain formulas that are in some sense simpler than those produced by existing identifiability algorithms. A simplification algorithm by Tikka and Karvanen, 2017 functions as a post-processing step after the identifying formula has already been obtained by the ID algorithm. Given a measure of simplicity, the search heuristic could be adjusted to find simple formulas directly without resorting to separate simplification procedures. In some specific scenarios, such as the standard causal effect identifiability problem, an approach known as pruning [Tikka and Karvanen, 2018] could be incorporated into the search. Pruning refers to the removal of vertices from the graph, that are not required for determining identifiability.

Finally we note that identifiability has also been studied under the assumption that the functional relationships depicted by the causal model are linear [Angrist et al., 1996, van der Zander and Liskiewicz, 2016, Chen et al., 2017] or non-parametric with additive error terms [Peters et al., 2014, Peña and Bendtsen, 2017] and when the causal graph is not completely known [Maathuis et al., 2009, Entner et al., 2013, Hyttinen et al., 2015, Perković et al., 2015, Malinsky and Spirtes, 2017, Jaber et al., 2018]. Extending the search in these directions is an interesting line of future research.

§ CONCLUSION

We presented : a do-calculus based search capable of solving identifiability problems for which no known solutions exist. This contribution is especially useful for researchers working in the field of causal inference to confirm theoretical results or to find counterexamples to identifiability claims. In practical terms, the search can also provide solutions to complicated problems such as combining transportability and selection bias, recovering from multiple bias sources or identifying causal quantities in the presence of missing data that cannot be solved by any other existing method. The R package dosearch providing an implementation of  is available on CRAN.

§ ACKNOWLEDGMENTS

This work belongs to the thematic research area “Decision analytics utilizing causal models and multiobjective optimization” (DEMO) supported by Academy of Finland (grant number 311877). AH was supported by Academy of Finland through grant 295673. The authors wish to acknowledge CSC – IT Center for Science, Finland, for computational resources.

[Angrist et al., 1996]
Angrist JD, Imbens GW, Rubin DB (1996).
Identification of Causal Effects Using Instrumental
Journal of the American Statistical Association,
91(434), 444–455.

[Bareinboim and Pearl, 2012]
Bareinboim E, Pearl J (2012a).
Causal Inference by Surrogate Experiments:
In N de Freitas, K Murphy (eds.), Proceedings of the 28th
Conference on Uncertainty in Artificial Intelligence, pp. 113–120. AUAI

[Bareinboim and Pearl, 2012]
Bareinboim E, Pearl J (2012b).
Controlling Selection Bias in Causal Inference.
In ND Lawrence, M Girolami (eds.), Proceedings of the 15th
International Conference on Artificial Intelligence and Statistics,
volume 22, pp. 100–108.

[Bareinboim and Pearl, 2013]
Bareinboim E, Pearl J (2013).
A General Algorithm for Deciding Transportability of
Experimental Results.
Journal of Causal Inference, 1, 107–134.

[Bareinboim and Pearl, 2014]
Bareinboim E, Pearl J (2014).
Transportability from Multiple Environments with Limited
Experiments: Completeness Results.
In Proceedings of the 27th Annual Conference on Neural
Information Processing Systems, pp. 280–288.

[Bareinboim and Tian, 2015]
Bareinboim E, Tian J (2015).
Recovering Causal Effects from Selection Bias.
In Proceedings of the 29th AAAI Conference on Artificial
Intelligence, pp. 3475–3481.

[Bareinboim et al., 2014]
Bareinboim E, Tian J, Pearl J (2014).
Recovering from Selection Bias in Causal and Statistical
In Proceedings of the 28th AAAI Conference on Neural
Information Processing Systems.

[Bhattacharya et al., 2019]
Bhattacharya R, Nabi R, Shpitser I, Robins JM (2019).
Identification in Missing Data Models Represented by
Directed Acyclic Graphs.
In Proceedings of the 35th Conference on Uncertainty in
Artificial Intelligence.

[Breslow, 1996]
Breslow NE (1996).
Statistics in Epidemiology: The Case-Control Study.
Journal of the American Statistical Association,
91(433), 14–28.

[Chen et al., 2017]
Chen B, Kumor D, Bareinboim E (2017).
Identification and Model Testing in Linear Structural
Equation Models using Auxiliary Variables.
In D Precup, YW Teh (eds.), Proceedings of the 34th
International Conference on Machine Learning, volume 70, pp. 757–766.

[Correa and Bareinboim, 2017]
Correa J, Bareinboim E (2017).
Causal Effect Identification by Adjustment under Confounding
and Selection Biases.
In Proceedings of the 31st AAAI Conference on Artificial

[Correa et al., 2018]
Correa J, Tian J, Bareinboim E (2018).
Generalized Adjustment under Confounding and Selection
In Proceedings of the 32nd AAAI Conference on Artificial

[Csardi and Nepusz, 2006]
Csardi G, Nepusz T (2006).
The igraph Software Package for Complex Network
InterJournal, Complex Systems, 1695.

[Danks et al., 2009]
Danks D, Glymour C, Tillman RE (2009).
Integrating Locally Learned Causal Structures with
Overlapping Variables.
In Advances in Neural Information Processing Systems, pp.

[Dawid, 2002]
Dawid AP (2002).
Influence Diagrams for Causal Modelling and Inference.
International Statistical Review, 70(2), 161–189.

[Eddelbuettel and Francois, 2011]
Eddelbuettel D, Francois R (2011).
Rcpp: Seamless R and C++
Journal of Statistical Software, 40(8), 1–18.

[Entner et al., 2013]
Entner D, Hoyer P, Spirtes P (2013).
Data-driven Covariate Selection for Nonparametric Estimation
of Causal Effects.
In Proceedings of the 16th International Conference on
Artificial Intelligence and Statistics, volume 31, pp. 256–264. PMLR.

[Greenland et al., 1999]
Greenland S, Robins JM, Pearl J (1999).
Confounding and Collapsibility in Causal Inference.
Statistical Science, 14(1), 29–46.

[Huang and Valtorta, 2006]
Huang Y, Valtorta M (2006a).
Identifiability in Causal Bayesian Networks: A Sound and
Complete Algorithm.
In Proceedings of the 21st National Conference on Artificial
intelligence – Volume 2, pp. 1149–1154. AAAI Press.

[Huang and Valtorta, 2006]
Huang Y, Valtorta M (2006b).
Pearl's Calculus of Intervention Is Complete.
In Proceedings of the 22nd Conference on Uncertainty in
Artificial Intelligence, pp. 217–224. AUAI Press.

[Hyttinen et al., 2012]
Hyttinen A, Eberhardt F, Hoyer PO (2012).
Causal Discovery of Linear Cyclic Models from Multiple
Experimental Data Sets with Overlapping Variables.
In Proceedings of the 28th Conference on Uncertainty in
Artificial Intelligence, pp. 387–396.

[Hyttinen et al., 2015]
Hyttinen A, Eberhardt F, Järvisalo M (2015).
Do-calculus when the True Graph is Unknown.
In Proceedings of the 31st Conference on Uncertainty in
Artificial Intelligence, pp. 395–404. AUAI Press.

[Jaber et al., 2018]
Jaber A, Zhang J, Bareinboim E (2018).
Causal Identification under Markov Equivalence.
In Proceedings of the 34th Conference on Uncertainty in
Artificial Intelligence, pp. 978–987. AUAI Press.

[Kalisch et al., 2012]
Kalisch M, Mächler M, Colombo D, Maathuis MH, Bühlmann P (2012).
Causal Inference Using Graphical Models with the
R Package pcalg.
Journal of Statistical Software, 47(11), 1–26.

[Karvanen, 2015]
Karvanen J (2015).
Study Design in Causal Models.
Scandinavian Journal of Statistics, 42(2), 361–377.

[Koller and Friedman, 2009]
Koller D, Friedman N (2009).
Probabilistic Graphical Models: Principles and Techniques.
MIT Press.

[Lauritzen, 2000]
Lauritzen SL (2000).
Causal Inference from Graphical Models.
In OE Barndorff-Nielsen, DR Cox, C Klüppelberg (eds.),
Complex Stochastic Systems, pp. 67–107. CRC Press.

[Lee et al., 2019]
Lee S, Correa J, Bareinboim E (2019).
General Identifiability with Arbitrary Surrogate
In Proceedings of the 35th Conference on Uncertainty in
Artificial Intelligence.

[Little and Rubin, 1986]
Little RJA, Rubin DB (1986).
Statistical Analysis with Missing Data.
John Wiley & Sons, Inc., New York, NY, USA.

[Maathuis and Colombo, 2015]
Maathuis MH, Colombo D (2015).
A Generalized Backdoor Criterion.
The Annals of Statistics, pp. 1060–1088.

[Maathuis et al., 2009]
Maathuis MH, Kalisch M, Bühlmann P (2009).
Estimating High-dimensional Intervention Effects from
Observational Data.
The Annals of Statistics, 37(6A), 3133–3164.

[Malinsky and Spirtes, 2017]
Malinsky D, Spirtes P (2017).
Estimating Bounds on Causal Effects in High-dimensional and
Possibly Confounded Systems.
International Journal of Approximate Reasoning, 88,

[Mohan and Pearl, 2018]
Mohan K, Pearl J (2018).
Graphical Models for Processing Missing Data.
Forthcoming, https://arxiv.org/abs/1801.03583.

[Mohan et al., 2013]
Mohan K, Pearl J, Tian J (2013).
Graphical Models for Inference with Missing Data.
In Proceedings of the 26th International Conference on Neural
Information Processing Systems, pp. 1277–1285.

[Peña and Bendtsen, 2017]
Peña JM, Bendtsen M (2017).
Causal Effect Identification in Acyclic Directed Mixed
Graphs and Gated Models.
International Journal of Approximate Reasoning, 90,

[Pearl, 1995]
Pearl J (1995).
Causal Diagrams for Empirical Research.
Biometrika, 82(4), 669–688.

[Pearl, 2009]
Pearl J (2009).
Causality: Models, Reasoning, and Inference.
Second edition. Cambridge University Press.

[Perković et al., 2015]
Perković E, Textor J, Kalisch M, Maathuis M (2015).
A Complete Generalized Adjustment Criterion.
In Proceedings of the 31st Conference on Uncertainty in
Artificial Intelligence, pp. 682–691. AUAI Press.

[Peters et al., 2014]
Peters J, Mooij JM, Janzing D, Schölkopf B (2014).
Causal Discovery with Continuous Additive Noise Models.
Journal of Machine Learning Research, 15, 2009–2053.

[R Core Team, 2019]
R Core Team (2019).
R: A Language and Environment for Statistical
R Foundation for Statistical Computing, Vienna, Austria.

[Richardson, 2003]
Richardson T (2003).
Markov Properties for Acyclic Directed Mixed Graphs.
Scandinavian Journal of Statistics, 30(1), 145–157.

[Sharma et al., 2019]
Sharma A, Kiciman E, et al. (2019).
DoWhy: A Python Package for Causal

[Shpitser et al., 2015]
Shpitser I, Mohan K, Pearl J (2015).
Missing Data as a Causal and Probabilistic Problem.
In M Meila, T Heskes (eds.), Proceedings of the 31st Conference
on Uncertainty in Artificial Intelligence, pp. 802–811. AUAI Press.

[Shpitser and Pearl, 2006]
Shpitser I, Pearl J (2006a).
Identification of Conditional Interventional Distributions.
In Proceedings of the 22nd Conference on Uncertainty in
Artificial Intelligence, pp. 437–444. AUAI Press.

[Shpitser and Pearl, 2006]
Shpitser I, Pearl J (2006b).
Identification of Joint Interventional Distributions in
Recursive Semi-Markovian Causal Models.
In Proceedings of the 21st National Conference on Artificial
Intelligence – Volume 2, pp. 1219–1226. AAAI Press.

[Shpitser and Pearl, 2008]
Shpitser I, Pearl J (2008).
Complete Identification Methods for the Causal Hierarchy.
Journal of Machine Learning Research, 9, 1941–1979.

[Spirtes et al., 1993]
Spirtes P, Glymour C, Scheines R (1993).
Causation, Prediction, and Search.
(2nd ed. MIT Press 2000).

[Textor and van der Zander, 2016]
Textor J, van der Zander B (2016).
dagitty: Graphical Analysis of Structural Causal Models.
R package version 0.2-2,

[Tikka et al., 2020]
Tikka S, Hyttinen A, Karvanen J (2020).
dosearch: Causal Effect Identification from Multiple
Incomplete Data Sources.
R package version 1.0.5,

[Tikka and Karvanen, 2017]
Tikka S, Karvanen J (2017a).
Identifying Causal Effects with the R Package
Journal of Statistical Software, 76(12), 1–30.

[Tikka and Karvanen, 2017]
Tikka S, Karvanen J (2017b).
Simplifying Probabilistic Expressions in Causal Inference.
Journal of Machine Learning Research, 18(36), 1–30.

[Tikka and Karvanen, 2018]
Tikka S, Karvanen J (2018).
Enhancing Identification of Causal Effects by Pruning.
Journal of Machine Learning Research, 18(194), 1–23.

[Tikka and Karvanen, 2019]
Tikka S, Karvanen J (2019).
Surrogate Outcomes and Transportability.
International Journal of Approximate Reasoning, 108,

[Tillman and Spirtes, 2011]
Tillman R, Spirtes P (2011).
Learning Equivalence Classes of Acyclic Models with Latent
and Selection Variables from Multiple Datasets with Overlapping Variables.
In Proceedings of the 14th International Conference on
Artificial Intelligence and Statistics, pp. 3–15.

[Triantafillou and Tsamardinos, 2015]
Triantafillou S, Tsamardinos I (2015).
Constraint-based Causal Discovery from Multiple
Interventions over Overlapping Variable Sets.
Journal of Machine Learning Research, 16, 2147–2205.

[Triantafillou et al., 2010]
Triantafillou S, Tsamardinos I, Tollis I (2010).
Learning Causal Structure from Overlapping Variable Sets.
In Proceedings of the 13th International Conference on
Artificial Intelligence and Statistics, pp. 860–867.

[van der Zander and Liskiewicz, 2016]
van der Zander B, Liskiewicz M (2016).
On Searching for Generalized Instrumental Variables.
In Proceedings of the 19th International Conference on
Artificial Intelligence and Statistics.

[van der Zander et al., 2019]
van der Zander B, Liśkiewicz M, Textor J (2019).
Separators and Adjustment Sets in Causal Graphs: Complete
Criteria and An Algorithmic Framework.
Artificial Intelligence, 270, 1–40.