# Causal Markov Boundaries

Sofia Triantafillou  Fattaneh Jabbari  Greg Cooper

###### Abstract

Feature selection is an important problem in machine learning, which aims to
select variables that lead to an optimal predictive model. In this paper, we
focus on feature selection for post-intervention outcome prediction from pre-
intervention variables. We are motivated by healthcare settings, where the
goal is often to select the treatment that will maximize a specific patient’s
outcome; however, we often do not have sufficient randomized control trial
data to identify well the conditional treatment effect. We show how we can use
observational data to improve feature selection and effect estimation in two
cases: (a) using observational data when we know the causal graph, and (b)
when we do not know the causal graph but have observational and limited
experimental data. Our paper extends the notion of Markov boundary to
treatment-outcome pairs. We provide theoretical guarantees for the methods we
introduce. In simulated data, we show that combining observational and
experimental data improves feature selection and effect estimation.

##  1 Introduction

Feature selection is a fundamental problem in machine learning that aims to
select the minimal set of features that lead to the optimal prediction of a
target variable $Y$. For observational distributions, this set is the Markov
boundary of $Y$, $\textnormal{MB}(Y)$. In causal graphical models, this set
can be identified from the causal graph $\mathcal{G}$ [Pearl, 2000]. This set
exhausts the predictive information for the state of a variable $Y$, and can
be used to obtain the best (and minimal) predictive model
$P(Y|\textnormal{MB}(Y))$ for $Y$.

Observational Markov Boundary (OMB) of $Y$: MB(Y) |  The Markov boundary of $Y$. Leads to optimal prediction of $Y$ from observational data.
---|---
Interventional Markov Boundary (IMB) of $Y$ relative to $X$: $\textnormal{MB}_{\overline{X}}(Y)$ |  The Markov boundary of $Y$ in the post-intervention distribution $P_{X}$. Leads to optimal prediction of $Y_{x}$ from experimental data.
Causal Markov Boundaries (CMB) of $Y$ relative to $X$: $\textnormal{{CMB}}_{\overline{X}}(Y)$ |  Sets of measured variables that satisfy Definition 3.2. Possibly not unique, and possibly empty. If not empty, one of the CMBs leads to the optimal prediction of $Y_{x}$ from observational data.
Table 1: Different Markov boundaries discussed in this paper.

In decision making, we are often interested in finding the optimal predictive
model for the post-intervention distribution of an outcome $Y$ after we
intervene on a treatment $X$, when we only have observational data. Ideally,
we would like to include in our model the Markov boundary $\mathbf{Z}$ of $Y$
in the post-intervention causal graph that is parameterized with the post-
intervention distribution. However, under causal insufficiency in which latent
confounding may exist, the conditional post-interventional distribution
$P(Y|do(X),\mathbf{Z})$ may not be identifiable. For example, in Fig. 1,
$P(Y|do(X),A,B)$ is not identifiable from the observational distribution
alone. In this case, we are interested in identifying the optimal set
$\mathbf{Z}$ for which the post-intervention distribution
$P(Y|do(X),\mathbf{Z})$ is identifiable from observational data, which we call
the causal Markov boundary.

Moreover, even when experimental data are available, they typically have much
smaller sample sizes and are not powered to identify conditional
distributions. In that case, we would like to combine large observational data
with limited experimental data to improve interventional feature selection and
effect estimation.

Our methods are heavily motivated by embedded clinical trials [Angus, 2015,
Angus et al., 2020], which take place within usual clinical care. In these
trials, patients who agree to participate are randomized to receive a
treatment from among those considered effective for that patient. The
electronic health records (EHRs) of the health system in which the trial is
being conducted contains both experimental data from the trial, and
observational data obtained outside (e.g., before/after) the trial, all
measuring the same variables. Combining observational and experimental data
has the potential to better predict the most effective treatments for
individual patients, than either type of data alone.

Our contributions are the following:

  * •

We define the interventional Markov boundary
${\textnormal{MB}_{\overline{X}}(Y)}\;$ and the causal Markov boundaries
$\textnormal{{CMB}}_{\overline{X}}(Y)$ for an outcome $Y$ and a treatment $X$.
These sets correspond to the minimal set of covariates $\mathbf{Z}$ that are
maximally informative for $Y|do(X)$, from experimental and observational data,
respectively (Sec 3). Table 1 summarizes the types of Markov boundaries
discussed in this paper.

  * •

We present a Bayesian method that combines observational and experimental data
to learn interventional Markov boundaries. The method provides estimates of
the post-interventional distribution that are based on both observational and
experimental data, when possible (Sec. 4), in which case the IMB is a CMB. In
simulated data, we show that our method improves causal effect estimation
(Sec. 6).

$X$A$B$$Y$ Figure 1:  An example SMCM. $\\{X,A,B\\}$ is the Markov boundary of
$Y$ in $\mathcal{G}_{\overline{X}}$ . $P(Y|do(X),A,B)$ is not identifiable
from the observational distribution $P$, but $P(Y|do(X),A)$ and $P(Y|do(X),B)$
are. $\\{X,B\\}$ is the causal Markov boundary for $Y_{X}$.

##  2 Preliminaries

We use the framework of semi-Markovian causal models [SMCMs, Tian and
Shpitser, 2003], and assume the reader is familiar with related terminology.
Variables are denoted in uppercase, their values in lowercase, and variable
sets in bold. We use $\mathcal{G}$ to denote a causal graph, and say
$\mathcal{G}$ induces a probability distribution $P$ if $P$ factorizes
according to $\mathcal{G}$ and the causal Markov condition.

We use $Y|do(X)$ or $Y_{X}$ to denote a variable $Y$ after the hard
intervention on variable $X$. If we know the causal SMCM $\mathcal{G}$, a hard
intervention of where a treatment $X$ is set to $x$ can be represented with
the do-operator, $do(X\texttt{=}x)$. We use $P_{x}$ or to denote the
interventional distribution over the same variables for $do(X\texttt{=}x)$. In
the corresponding graph, this is equivalent to removing all incoming edges
into $X$, while keeping all other mechanisms intact. We use
${\mathcal{G}_{\overline{X}}}\;$ to denote the graph stemming from
$\mathcal{G}$ after removing edges into $X$. We use
$\mathcal{G}_{\underline{X}}$ to denote the graph stemming from $\mathcal{G}$
after removing edges out of $X$. We use the terms
$Pa_{\mathcal{G}}(Z),Ch_{\mathcal{G}}(Z)$ to denote the set of parents and
children of $Z$ in $\mathcal{G}$, respectively. The set of variables that are
connected with a variable $Y$ through a bidirected path (i.e., a path that
only has bidirected edges) is called the district of $Y$, and denoted
$Dis_{\mathcal{G}}(Y)$.

##  3 Markov Boundaries

A Markov blanket of a variable $Y$ in a set of variables $\mathbf{V}$ is a
subset $\mathbf{Z}$ of $\mathbf{V}$ conditioned on which other variables are
independent of $Y$:
$Y{\,\perp\mkern-12.0mu\perp\,}\mathbf{V}\setminus\mathbf{Z}|\mathbf{Z}$. The
Markov boundary of $Y$ is the Markov blanket that is also minimal (i.e., no
subset of the Markov boundary is a Markov blanket) [Pearl, 2000]. In
distributions that satisfy the intersection property (including faithful
distributions), the Markov boundary of a variable $Y$ is unique [Pearl, 1988].
To distinguish from other types of Markov boundaries defined in this work, we
often use the terminology Observational Markov Boundary (OMB) to denote the
Markov boundary of a variable.

For a DAG $\mathcal{G}$, the OMB of a variable $Y$ in any distribution
faithful to $\mathcal{G}$ is the set parents, children, and spouses of $Y$:
$\textnormal{MB}(Y)=Pa_{\mathcal{G}}(Y)\cup Ch_{\mathcal{G}}(Y)\cup
Pa_{\mathcal{G}}(Ch_{\mathcal{G}}(Y))$. For SMCMs, it has been shown that the
OMB of a variable $Y$ is the set of parents, children, children’s parents
(spouses) of $Y$, district of $Y$ and districts of the children of $Y$, and
the parents of each node of these districts [Richardson, 2003, Pellet and
Elisseeff, 2008]111Pellet and Elisseeff [2008] prove this for maximal
ancestral graphs, but the proof can be readily adapted to SMCMs..

The OMB has been shown to be the minimal set of variables with optimal
predictive performance for a given distribution and response variable, given
some assumptions on the learner and the loss function [Tsamardinos and
Aliferis, 2003]. In this work, we are interested in the model that gives the
optimal prediction of the post-intervention distribution, with the goal of
designing optimal policies. For this reason, we are not interested in
including post-intervention covariates in this model, because these variables
are not known prior to treatment assignment, and thus, cannot affect the
assignment. In the rest of this document, we make the following assumption:

######  Assumption 3.1.

Covariates $\mathbf{V}$ are pre-treatment.

This simplifies the expressions for the OMBs, because we no longer need to
consider children of $Y$ and their districts. Knowing the OMB allows a more
efficient representation of the conditional distribution of $Y$ given
$\mathbf{V}$, since the following equation holds:

$P(Y|\mathbf{V})=P(Y|MB(Y)).$ (1)

###  3.1 Interventional Markov Boundary

Our goal is to identify the set of variables that lead to the optimal model
for the post-intervention distribution of a target $Y$ relative to a specific
treatment $X$. We call this set the interventional Markov boundary (IMB) of
$Y$ relative to $X$, and denote it $\textnormal{MB}_{\overline{X}}(Y)$ .
Obviously, ${\textnormal{MB}_{\overline{X}}(Y)}\;\subseteq\textnormal{MB(Y)}$.
When we have data from the post-intervention distribution, we can apply
statistical methods for OMB identification to obtain the IMB of $Y$ relative
to $X$. However, experimental data are often limited in sample sizes, while
OMB identification methods may require large sample sizes.

If we know the causal graph $\mathcal{G}$, the post-intervention distribution
with respect to $X$ is induced by the manipulated graph
${\mathcal{G}_{\overline{X}}}\;$. The IMB of $Y$ is then the OMB of $Y$ in
${\mathcal{G}_{\overline{X}}}\;$, and can be identified using the definition
of the Markov Boundary above. However, the post-intervention distribution
$P(Y|do(X),{\textnormal{MB}_{\overline{X}}(Y)}\;\setminus X)$, may not be
identifiable from the observational distribution. For example, in Fig. 1,
${\textnormal{MB}_{\overline{X}}(Y)}\;=\\{X,A,B\\}$, but $P(Y|do(X),A,B)$ is
not identifiable from observational data. We then want to answer the following
question: _What is the best model for predicting $Y_{X}$ from the
observational distribution, when the causal graph is known?_

###  3.2 Causal Markov Boundaries

To answer this question, we define the causal Markov boundaries of an outcome
$Y$ relative to a treatment $X$ as follows:

######  Definition 3.2.

Let $\mathbf{Z}\subseteq(\mathbf{V}\cup X)$, and
$\mathbf{W}=\mathbf{Z}\setminus X$. Then $\mathbf{Z}$ is a causal Markov
boundary (CMB) for $Y$ relative to $X$ if it satisfies the following
properties:

  1. 1.

$P(Y|do(X),\mathbf{W})$ is identifiable from $P(X,Y,\mathbf{V})$.

  2. 2.

For every subset $\mathbf{W}^{\prime}$ of $\mathbf{V}\setminus\mathbf{W}$
either $P(Y|do(X),\mathbf{W},\mathbf{W}^{\prime})$ $=P(Y|do(X),\mathbf{W})$ or
$P(Y|do(X),\mathbf{W},\mathbf{W}^{\prime})$ is not identifiable from
$P(X,Y,\mathbf{V})$.

  3. 3.

$\nexists\mathbf{W}^{\prime}\subset\mathbf{W}$ s.t.
$P(Y|do(X),\mathbf{W}^{\prime})=P(Y|do(X),\mathbf{W})$.

Condition (1) ensures that the post-intervention conditional probability of
$Y_{X}$ given a CMB is identifiable. Condition (2) states that the covariates
that are not in that CMB are either redundant for the prediction of $Y_{X}$
given the CMB, or they make the post-intervention distribution non-
identifiable. Condition (3) ensures that $\mathbf{Z}$ is additionally
maximally informative for $Y_{X}$ in the sense that you cannot remove any
variable from $\mathbf{Z}$ without losing some information for $Y_{X}$. This
condition rules out sets like $\\{X,A\\}$ in Fig. 1, where, while
$P(Y|do(X),A)$ is identifiable from $P$, it is equal to $P(Y|do(X))$. Thus,
conditioning on $A$ does not improve the prediction of $Y_{X}$ compared to its
subset $\emptyset$.

Notice that this definition does not capture the spirit of Markov boundaries
precisely: Markov boundaries make all remaining variables redundant for
predicting $Y$; however, this does not necessarily hold with CMBs. For
example, in Fig. 1, $\\{X,B\\}$ is a CMB according to the definition above,
but $\\{A\\}$ remains relevant for predicting $Y_{X}$; however, including it
with $B$ in the CMB leads to non-identifiability.

CMB is not necessarily unique; it is possible that multiple sets satisfy
Definition 3.2. For example, assume the distribution $P$ is induced by the
SMCM shown in Fig. 2. Both $\\{X,B,C\\}$ and $\\{X,A,D\\}$, satisfy Definition
3.2. The best predictive CMB for predicting $Y_{X}$ will depend on the
parameters in $P$. We use the notation $\textnormal{{CMB}}_{\overline{X}}(Y)$
to denote the set of causal Markov boundaries of $Y$ relative to $X$. Thus, we
will generally need to find all CMBs and then determine which of them leads to
the best prediction of $Y_{X}$. Also, notice that the
$\textnormal{{CMB}}_{\overline{X}}(Y)$ can be empty; thus, no subset of
$\mathbf{V}$ satisfies the Definition 3.2. This can happen for example if
$X\rightarrow Y$ and $X\leftrightarrow Y$ in $\mathcal{G}$.

The CMB is useful in determining a minimal set of maximally predictive
variables for which we can use observational data to predict post-
interventional distributions. In the next section we show that, for pre-
treatment covariates, CMBs satisfy the back-door criterion and are subsets of
the observational Markov boundary. These results enable more efficient
algorithms for finding CMBs, limiting the types of estimators and the number
of variable sets we need to consider.

$X$$A$$B$$C$$D$$Y$ Figure 2:  Causal Markov boundaries are not necessarily
unique. Both $\\{X,B,C\\}$ and $\\{X,A,D\\}$ are causal Markov boundaries for
$Y$ relative to $X$.

####  3.2.1 Characterization

Given a graph $\mathcal{G}$, Shpitser and Pearl [2006a] provide a sound and
complete algorithm (IDC) for estimating conditional post-intervention
distributions from observational distributions induced by $\mathcal{G}$. The
output of this algorithm is an expression for $P(Y|do(X),\mathbf{W})$ if the
distribution is identifiable from distribution $P$ and $\mathcal{G}$, or N/A
otherwise. Thus, we can identify CMBs in a brute-force way by running IDC for
every possible subset of $\mathbf{V}$, and then check for sets that satisfy
the conditions in Def. 3.2. This process is computationally expensive and
would not be possible for graphs with more than a few variables.

In this section, we provide theoretical results that lead to a much easier
process when all candidate conditioning variables are pre-treatment (all
proofs can be found in the supplementary). For pre-treatment covariates, one
obvious family of sets for which the conditional post-intervention
distributions are identifiable are sets that m-separate $X$ and $Y$ in
$\mathcal{G}_{\underline{X}}$. These sets satisfy Rule 2 of do-calculus
[Pearl, 2000], so the conditional interventional distribution
$P(Y|do(X),\mathbf{W})$ is equal to the observational distribution
$P(Y|X,\mathbf{W})$. Sets of pre-treatment covariates that m-separate $X$ and
$Y$ in $\mathcal{G}_{\underline{X}}$ are also known to satisfy the backdoor
criterion [Van der Zander et al., 2014] and the adjustment criterion [Shpitser
et al., 2012]. However, these definitions are more general to include possible
post-treatment covariates, and are intended for estimating marginal post-
intervention distributions (or average effects). For brevity, we will call
sets that m-separate $X$ and $Y$ in $\mathcal{G}_{\underline{X}}$ backdoor
sets, since they block all back-door paths between $X$ and $Y$.

One question that arises is if there are sets that are not backdoor sets that
may satisfy the conditions in Definition 3.2. In that case, identifiability
could stem from some sequential application of do-calculus rules. As we show
next, this is not possible for pre-treatment covariates. This ensures that we
only need to check CMB Conditions (2) and (3) for sets for which
$P(Y|do(X),\mathbf{W})=P(Y|X,\mathbf{W})$. This makes the identification of
$P(Y|do(X),\mathbf{W})$ more straightforward than having to compute more
complex probability expressions.

######  Theorem 3.3.

We assume that $P_{x}$ and $\mathcal{G}_{\overline{X}}$ are faithful to each
other. If $\mathbf{Z}$ is a CMB for $Y$ relative to $X$, then
$\mathbf{Z}\setminus X$ is a backdoor set for $X$ relative to $Y$.

The second theoretical result is that any CMB of $Y$ relative to $X$ is a
subset of the OMB of $Y$. While this sounds intuitive, it is not completely
straightforward. It could be the case that conditioning on every subset of the
OMB opens some m-connecting path between $X$ and $Y$, that can only be blocked
by a variable that is not a member of the Markov boundary. The following
theorem proves that this is not possible, allowing for more efficient search
algorithms:

######  Theorem 3.4.

We assume that $P_{x}$ and $\mathcal{G}_{\overline{X}}$ are faithful to each
other. Every CMB $\mathbf{Z}$ of an outcome variable $Y$ w.r.t a treatment
variable $X$ is a subset of the OMB $\textnormal{MB}(Y)$.

Based on Theorem 3.4, we only need to look for CMBs within subsets of the OMB
of $Y$. So far, we have shown that both the IMB and any CMB are subsets of the
OMB. We can also show that when the IMB is a CMB, then it also coincides with
the OMB:

######  Theorem 3.5.

If $\textnormal{MB}_{\overline{X}}(Y)$ is a causal Markov boundary, then
${\textnormal{MB}_{\overline{X}}(Y)}\;=\textnormal{MB}(Y)$.

##  4 Combining observational and experimental data

When the causal graph is known, we can obtain CMBs by looking for subsets of
$\textnormal{MB}(Y)$ that satisfy Def. 3.2. Unfortunately, in most real-world
applications, the true graph is unknown, and selecting the
causal/interventional Markov boundary is not possible from observational data
alone. Experimental data may exist, but are typically much fewer than
observational data, due to expense or ethical concerns. This scenario is
common in embedded trials, where non-randomized patients are much more common
than trial participants. In such cases, the experimental data may be under-
powered to accurately estimate conditional effects. As a result, the
conditional effects that can be derived from the experimental data have high
variance and may not be reliable. In this case, combining all data
(observational and experimental) in a Bayesian manner may help improve the
prediction of $Y_{x}$.

We assume that we have observational data $D_{o}$ and experimental data
$D_{e}$ measuring treatment $X$, outcome $Y$, and pre-treatment covariates
$\mathbf{V}$. We use $N_{o},N_{e}$ to denote the number of samples in
$D_{o},D_{e}$, respectively.

We present a Bayesian method, called FindIMB, that uses both $D_{e}$ and
$D_{o}$ to estimate the probability of a set being the
${\textnormal{MB}_{\overline{X}}(Y)}\;\\!$, and estimate
$P(Y|do(X),\mathbf{V})=P(Y|do(X),{\textnormal{MB}_{\overline{X}}(Y)}\;\setminus
X)$. The method is presented in Alg. 1. The method first estimates the OMB of
$Y$ in observational data $\textnormal{MB}(Y)$ (Line 1), and then looks among
subsets of $\textnormal{MB}(Y)$ for sets that are IMBs (Line 2). It uses
$D_{e}$ and $D_{o}$ to evaluate the probability that a set is an IMB (Line 3),
and then returns a weighted average for $P(Y|do(X),\mathbf{V})$ based on these
probabilities (Line 5).

The enabling idea of the method is that, when the IMB is a CMB, we can use
both the $D_{o}$ and $D_{e}$ to estimate the conditional post-intervention
distribution. Otherwise, we use only $D_{e}$ to derive the estimate. We use
the following notation to express these hypotheses:

  * •

${H_{\mathbf{\mathbf{Z}}}^{c}}\;$ is a binary variable denoting the hypothesis
that $\mathbf{Z}$ is the IMB ${\textnormal{MB}_{\overline{X}}(Y)}\;$, and it
is also a CMB:
$\mathbf{Z}={\textnormal{MB}_{\overline{X}}(Y)}\;\wedge\mathbf{Z}\in{\textnormal{{CMB}}_{\overline{X}}(Y)}\;$.

  * •

${H_{\mathbf{\mathbf{Z}}}^{\overline{c}}}\;$ is a binary variable denoting the
hypothesis that $\mathbf{Z}$ is the IMB
${\textnormal{MB}_{\overline{X}}(Y)}\;$, but it is not a CMB:
$\mathbf{Z}={\textnormal{MB}_{\overline{X}}(Y)}\;\wedge\mathbf{Z}\not\in{\textnormal{{CMB}}_{\overline{X}}(Y)}\;$.

For a set $\mathbf{Z}^{\star}$, if either
${H_{\mathbf{\mathbf{Z}^{\star}}}^{c}}\;$ or
${H_{\mathbf{\mathbf{Z}^{\star}}}^{\overline{c}}}\;$ is true,
$\mathbf{Z}^{\star}$ is an IMB and therefore
$P(Y|do(X),\mathbf{V})=P(Y|do(X),\mathbf{Z}^{\star}\setminus X)$. Under
${H_{\mathbf{\mathbf{Z}^{\star}}}^{c}}\;$ though, $\mathbf{Z}^{\star}$ is also
a CMB and therefore the pre- and post- intervention distributions are the
same, i,e,

$P(Y|do(X),\mathbf{Z}^{\star}\setminus
X,{H_{\mathbf{Z^{\star}}}^{c}}\;)=P(Y|X,\mathbf{Z}^{\star}\setminus X)$ (2)

In contrast, under $H_{\mathbf{Z^{\star}}}^{\overline{c}}$ , Eq. 2 does not
hold: $\mathbf{Z}^{\star}$ is the IMB, but not a CMB. This means that
$P(Y|do(X),\mathbf{Z}^{\star}\setminus X)$ is not identifiable from
observational data, and we cannot use $D_{o}$ to estimate
$P(Y|do(X),\mathbf{Z}^{\star}\setminus X)$222Notice however that $D_{o}$ may
still place some constraints on $P(Y|do(X),\mathbf{Z}^{\star}\setminus X)$,
like for example provide bounds.. _In summary, if either of
${H_{\mathbf{Z^{\star}}}^{c}}\;or{H_{\mathbf{Z^{\star}}}^{\overline{c}}}\;$
holds, $\mathbf{Z}^{\star}$ is the IMB and
$P(Y|do(X),\mathbf{V})=P(Y|do(X),\mathbf{Z}^{\star}\setminus X)$. If
${H_{\mathbf{Z^{\star}}}^{c}}\;$ holds, we can use both $D_{o}$ and $D_{e}$ in
our estimation of $P(Y|do(X),\mathbf{Z}^{\star}\setminus X)$, while if
${H_{\mathbf{Z^{\star}}}^{\overline{c}}}\;$ holds we can only use $D_{e}$._

Based on this observation, we want to compute
$P({H_{\mathbf{Z}}^{c}}\;|D_{e},D_{o})$ and
$P({H_{\mathbf{Z}}^{\overline{c}}}\;|D_{e},D_{o})$ for possible IMBs
$\mathbf{Z}$. These probabilities tell us both how likely it is that
$\mathbf{Z}$ is an IMB (their sum), and if we can include observational data
in the estimation of $P(Y|\do(X),\mathbf{V})$. Using Bayes rule, we obtain:

$\begin{split}&P({H_{\mathbf{Z^{\prime}}}^{c}}\;|D_{e},D_{o})=\\\
&\frac{P(D_{e}|D_{o},{H_{\mathbf{Z^{\prime}}}^{c}}\;)P(D_{o}|{H_{\mathbf{Z^{\prime}}}^{c}}\;)P({H_{\mathbf{\mathbf{Z}{{}^{\prime}}}}^{c}}\;)}{\displaystyle\sum_{\mathbf{Z}}{\sum_{C=c,\overline{c}}P(D_{e}|D_{o},{H_{\mathbf{\mathbf{Z}}}^{C}}\;)P(D_{o}|{H_{\mathbf{\mathbf{Z}}}^{C}}\;)P({H_{\mathbf{\mathbf{Z}}}^{C}}\;)}}.\end{split}$
(3)

We can similarly derive
$P({H_{\mathbf{\mathbf{Z}}}^{\overline{c}}}\;|D_{e},D_{o})$ by replacing each
appearance of $c$ with $\overline{c}$ in the numerator. The denominator is the
same for all sets. $P({H_{\mathbf{\mathbf{Z}{{}^{\prime}}}}^{c}}\;)$ and
$P({H_{\mathbf{\mathbf{Z}{{}^{\prime}}}}^{\overline{c}}}\;)$ is our prior that
${H_{\mathbf{\mathbf{Z}{{}^{\prime}}}}^{c}}\;$ or
${H_{\mathbf{\mathbf{Z}{{}^{\prime}}}}^{\overline{c}}}\;$ holds. We set this
to be uniform over both values of $C$ and all $\mathbf{Z}$.

As Eq. 3 shows, using Bayes rule, we can estimate the posterior probabilities
for the set of hypotheses ${H_{\mathbf{Z}}^{c}}\;$ and
${H_{\mathbf{Z}}^{\overline{c}}}\;$ using marginal likelihoods of the
experimental and observational data. In the next sections, we show how we can
compute each term in Eq. 3.

We present our results for multinomial distributions, but we believe the
method can be readily extended to any type of distribution with closed-form
marginals. Due to space constraints, the closed-form solution for each
equation appearing in remainder of the paper is presented in Supplementary
Table 1.

Estimating
$P(D_{e}|D_{o},{H_{\mathbf{Z^{\prime}}}^{c}}\;),P(D_{e}|D_{o},{H_{\mathbf{Z^{\prime}}}^{\overline{c}}}\;)$:
Let $\mathbf{W}=\mathbf{Z}\setminus X$, and let $\theta_{Y_{x}|\mathbf{W}}$ be
a set of parameters expressing the conditional probabilities for
$P(Y|do(X),\mathbf{W})$. Also, let $\theta_{Y|X,\mathbf{W}}$ denote the
observational parameters for $P(Y|X,\mathbf{W})$. By integrating over all
$\theta_{Y_{x}|\mathbf{W}}$, we obtain

$\begin{split}&P(D_{e}|D_{o},{H_{\mathbf{Z}}^{c}}\;)=\\\
&\int_{Y_{x}|\mathbf{W}}P(D_{e}|\theta_{Y_{x}|\mathbf{W}})f(\theta_{Y_{x}|\mathbf{W}}|D_{o},{H_{\mathbf{Z}}^{c}}\;)d\theta_{Y_{x}|\mathbf{W}}\end{split}$
(4)

$f(\theta_{Y_{x}|\mathbf{W}}|D_{o},{H_{\mathbf{Z}}^{c}}\;)$ is the posterior
for $\theta_{Y_{x}|\mathbf{W}}$ given the observational data, when
$\mathbf{Z}$ is the IMB and the CMB. In this case,
$P(Y|do(X),\mathbf{W})=P(Y|X,\mathbf{W})$, and therefore
$f(\theta_{Y_{x}|\mathbf{W}}|D_{o},{H_{\mathbf{Z}}^{c}}\;)=f(\theta_{Y|X,\mathbf{W}}|D_{o})$.
Eq. 4 can then be rewritten in terms of the observational parameters as

$\begin{split}&P(D_{e}|D_{o},{H_{\mathbf{Z}}^{c}}\;)=\\\
&\int_{\theta_{Y|X,\mathbf{W}}}P(D_{e}|\theta_{Y|X,\mathbf{W}})f(\theta_{Y|X,\mathbf{W}}|D_{o})d\theta_{Y|X,\mathbf{W}}\end{split}$
(5)

Eq. 5 is the marginal likelihood of $Y$ in experimental data, with parameter
density $f(\theta_{Y|X,\mathbf{W}}|D_{o})$ being equal to the parameter
posterior given $D_{o}$. In other words, under ${H_{\mathbf{Z}}^{c}}\;$, the
observational and experimental parameters coincide. Therefore, $D_{o}$ gives
us a strong ”prior” for $D_{e}$. Eq. 5 can be computed in closed-form for
distributions with conjugate priors.

Under $H_{\mathbf{Z}}^{\overline{c}}$ , the equality of the observational and
experimental parameters does not hold, and we cannot use the
$\theta_{Y|X,\mathbf{W}}$ to inform $\theta_{Y_{X}|\mathbf{W}}$, at least not
in a straightforward way. Instead, we impose that
$f(\theta_{Y_{X}|\mathbf{W}}|D_{o})=f(\theta_{Y_{X}|\mathbf{W}})$. Then
$P(D_{e}|D_{o},{H_{\mathbf{Z}}^{\overline{c}}}\;)$ corresponds to the marginal
likelihood of $Y$ in the experimental data, using a prior that we model as
being non-informative.

Estimating
$P(D_{o}|{H_{\mathbf{Z}}^{c}}\;),P(D_{o}|{H_{\mathbf{Z}}^{\overline{c}}}\;)$:
These probabilities score how well the observational data fit with the
hypotheses ${H_{\mathbf{Z}}^{c}}\;,{H_{\mathbf{Z}}^{\overline{c}}}\;$. We can
derive these terms on the basis of the OMB and its connection to the IMB and
the CMBs. We first need to express the hypothesis that a set $\mathbf{U}$ is
the OMB of $Y$: Let ${H_{\mathbf{U}}^{o}}\;$ denote this hypothesis; thus, for
any $\mathbf{U}\subseteq\mathbf{V}\cup X$, ${H_{\mathbf{\mathbf{U}}}^{o}}\;$
is true iff $\mathbf{U}$ is the OMB for $Y$. Then we can write

$P(D_{o}|{H_{\mathbf{Z}}^{C}}\;)=\sum_{\mathbf{U}\subseteq\mathbf{V}\cup
X}P(D_{o}|{H_{\mathbf{U}}^{o}}\;)P({H_{\mathbf{U}}^{o}}\;|{H_{\mathbf{Z}}^{C}}\;),$
(6)

for $C=c,\overline{c}$. Under ${H_{\mathbf{Z}}^{c}}\;\\!$, Theorem 3.5 implies
that $P({H_{\mathbf{U}}^{o}}\;|{H_{\mathbf{Z}}^{c}}\;)=1$ if
$\mathbf{U}=\mathbf{Z}$, and zero otherwise. Under
${H_{\mathbf{Z}}^{\overline{c}}}\;$, the IMB is not a CMB. Instead, the IMB
has to be a subset of $\mathbf{U}$, therefore
$P({H_{\mathbf{U}}^{o}}\;|{H_{\mathbf{Z}}^{C}}\;)=0$ for any
$\mathbf{Z}\supset\mathbf{U}$.

$P(D_{o}|{H_{\mathbf{U}}^{o}}\;)$ is the marginal likelihood of $Y$ in
$D_{o}$, under the hypothesis that $\mathbf{U}$ is the data-generating OMB for
$Y$ in the observational data. We can obtain this likelihood using a Bayesian
scoring algorithm like FGES, by scoring a DAG where $Y$ is a child of
variables $\mathbf{U}$ (and no other edges are in the graph). We call this
algorithm FGESMB. For discrete variables, in the large sample limit this
probability will be maximum only for the true OMB:

######  Theorem 4.1.

Given dataset $D_{o}$ that contains samples from a strictly positive
distribution $P$, which is a perfect map for a SMCM $\mathcal{G}$, the BD
score [Heckerman et al., 1995] will assign the highest score to the OMB of $Y$
in the large sample limit.

Eq. 6 needs to be computed for all subsets of the covariate sets. In practice,
however, for large $N_{o}$, these probabilities are dominated by the true OMB.
Assuming our sample is large enough, we can use an algorithm with asymptotic
guarantees for identifying the true OMB. In fact, once we commit to the an
observational Markov boundary $\mathbf{U}^{\star}$, Eq. 6 leads to the
following equations

$\begin{split}&P(D_{o}|{H_{\mathbf{Z}}^{c}}\;)\texttt{=}P(D_{o}|{H_{\mathbf{U^{\star}}}^{o}}\;)\textnormal{
for }\mathbf{Z}\texttt{=}\mathbf{U}^{\star}\\\
&P(D_{o}|{H_{\mathbf{Z}}^{\overline{c}}}\;)=P(D_{o}|{H_{\mathbf{U^{\star}}}^{o}}\;)\textnormal{
for all }\mathbf{Z}\subseteq\mathbf{U}^{\star}\\\ \end{split}$ (7)

For the remaining cases (${H_{\mathbf{Z}}^{c}}\;$ and
$\mathbf{Z}\texttt{=}\mathbf{U}^{\star}$, or ${H_{\mathbf{Z}}^{c}}\;$ and
$\mathbf{Z}\supset\mathbf{U}^{\star}$), the corresponding probabilities are
zero.

Bayesian estimation of $P(Y|do(X),\mathbf{V},D_{e},D_{o})$:
We now compute the $P(Y|do(X),\mathbf{V})$ using Bayesian model averaging over
the hypotheses ${H_{\mathbf{Z}}^{C}}\;$. Let $x,y,\mathbf{V}=\mathbf{v}$,
denote given instances of $X,Y$ and $\mathbf{V}$, respectively. When
$\mathbf{V}=\mathbf{v}$, we use $\mathbf{W}=\mathbf{w}_{v}$ to denote the
corresponding values of a set $\mathbf{W}\subset\mathbf{V}$. Recall that under
${H_{\mathbf{Z}}^{c}}\;$ and under , ${H_{\mathbf{Z}}^{\overline{c}}}\;$,
$P(Y|do(X),\mathbf{V})=P(Y|X,\mathbf{W})$, where
$\mathbf{W}=\mathbf{Z}\setminus X$.Then for a given instance of
$\mathbf{V}=\mathbf{v}$, we have

$\begin{split}&P(y|do(x),\mathbf{v},D_{e},D_{o})=\\\
&\sum_{\mathbf{Z}\subset\mathbf{V}}\sum_{C=c,\overline{c}}P(y|do(x),\mathbf{w}_{v},D_{e},D_{o},{H_{\mathbf{Z}}^{C}}\;)P({H_{\mathbf{Z}}^{C}}\;|D_{o},D_{e})\end{split}$
(8)

This equation computes the _expectation_ of the conditional probability
parameter. The individual probabilities
$P(Y|do(X),\mathbf{W},D_{e},D_{o},{H_{\mathbf{Z}}^{C}}\;)$ can be estimated as
posterior expectations of $P(Y|do(X),\mathbf{W})$ from the data. Specifically,
under given $H_{\mathbf{Z}}^{c}$ , $P(Y|do(X),\mathbf{W})=P(Y|X,\mathbf{W})$,
and therefore we can use both $D_{e}$ and $D_{o}$ for the posterior
expectation. In contrast, under $H_{\mathbf{Z}}^{\overline{c}}$ , we only use
$D_{e}$. Analytical equations for these probabilities for multinomial
distributions can be found in Supplementary Table 1.

input : $D_{o},D_{e}$, treatment $X$, outcome $Y$, pre-treatment covariates
$\mathbf{V}$

output : Post-intervention distribution $P(Y|do(X),\mathbf{V})$

1 $\textnormal{MB}(Y)\leftarrow\textnormal{{MarkovBoundary}}(Y,D_{o})$;

2 foreach  _subset $\mathbf{Z}$ of $\textnormal{MB}(Y)$ and
$C=c,\overline{c}$_ do

3 Compute $P({H_{\mathbf{\mathbf{Z}}}^{C}}\;|D_{e},D_{o})$ using Eq. 3;

4 Compute $P(Y|do(X),\mathbf{V},D_{e},D_{o},{H_{\mathbf{\mathbf{Z}}}^{C}}\;)$
using Eq. 8;

5$P(Y|do(X),\mathbf{V})\leftarrow\sum_{\mathbf{Z}}\sum_{C=c,\overline{c}}P(Y|do(X),\mathbf{V},D_{e},D_{o},{H_{\mathbf{\mathbf{Z}}}^{C}}\;)P({H_{\mathbf{\mathbf{Z}}}^{C}}\;|D_{e},D_{o})$;

Algorithm 1 FindIMB |  |
---|---|---
Figure 3:  Boxplots of absolute bias in the estimation of
$P(Y|do(X),\mathbf{V})$ using (a) FindIMB (b) IMB (c) OMB and (d) FCIt-IMB.
Data were simulated from random DAGs with 10 observed and 5 latent variables.
$D_{o}$ included 10,000 samples, and $D_{e}$ included 100 (left), 200 (middle)
or 1000 (right) samples. FindIMB improves the estimation of
$P(Y|do(X),\mathbf{V})$ particularly for smaller experimental sample sizes.
Black asterisks denote statistical significance, assessed with the Wilcoxon
signed-rank test. Three stars correspond to $p<0.001$.

##  5 Related work

We are not aware of other methods that try to identify causal and
interventional Markov boundaries. Our work has connections and builds on work
from many different areas. Due to space constraints, we only focus on methods
that do not require causal sufficiency. Markov boundaries: Several algorithms
learn OMBs from data under causal insufficiency [Yu et al., 2018, 2020]. In
addition, FGESMB presented in Sec. 4 is also a sound and complete method for
learning OMBs from data. These methods can be used to identify the IMBs from
the experimental data, but they do not combine observational and experimental
data to learn IMBs. Identifiability: Shpitser and Pearl [2006a, b] and Tian
and Shpitser [2003] provide sound and complete identifiability results for
post-intervention distributions from observational data when the causal graph
is known. These methods can answer queries for a specific marginal or
conditional probability of interest. Hyttinen et al. [2015] and Jaber et al.
[2019] provide similar identifiability results when the graph is unknown,
using the Markov equivalence class of graphs that are consistent with the
observational data. Hyttinen et al. [2015] can provide identifiability results
for graphs that are consistent with conditional independencies in both $D_{e}$
and $D_{o}$. However, the method is not proven to be complete for these
settings. These methods are not directly comparable with our method because
they do not select features for optimal prediction. Moreover, they provide
expressions for the post-intervention distributions that are based on
observational data alone, not by combining $D_{o}$ and $D_{e}$ like FindIMB.
Combining observational and experimental data to learn causal graphs: Several
causal discovery methods combine observational and experimental data to learn
causal structure [Triantafillou and Tsamardinos, 2015, Hyttinen et al., 2014,
Mooij et al., 2020, Andrews et al., 2020]. These methods return a summarized
version of all the causal graphs that are consistent with all the independence
constraints in all the data sets, observational and experimental. While these
methods can be used to improve the estimation of IMBs, it is not clear that
they can always provide a unique solution in this setting. Two additional
drawbacks they have for the purpose of optimized target prediction are that
(a) they rely on conditional independence tests that are unreliable when
$N_{e}$ is low, and (b) they learn the entire graph and do not focus on
finding the neighborhood of the target variable. This can result in unreliable
orientations due to error propagation. The method that has the closest setting
to ours is by Andrews et al. [2020] (Mooij et al. [2020] is also related, but
more general, and the two are equivalent for our setting). The authors propose
a method called FCItiers that can learn a family of SMCMs from $D_{e}$ and
$D_{o}$ when (a) the target of the intervention is known and (b) we specify
”tiered knowledge” on the variables (e.g., we know which variables are pre-
treatment). The method is complete in these settings. In the experimental
section, we develop a baseline comparison method based on FCItiers. Selecting
optimal adjustment sets. Some methods seek to select optimal adjustment sets
for efficient average treatment effect estimation [Perkovic et al., 2017,
Rotnitzky and Smucler, 2019, 2020]. While these methods have a different
purpose than ours, they have some connections with our work since, for pre-
treatment variables, any CMB is also an adjustment set. We point out that
while optimal adjustment sets and CMBs may often coincide (for example in
DAGs), they are not always the same (see example Fig S1 in the Supplementary).
Moreover, these methods are not directly comparable to ours since they focus
on identifying average treatment effects. Potential Outcomes Approaches:
Kallus et al. [2018] present a method for estimating conditional average
treatment effects (CATEs) by combining $D_{o}$ and $D_{e}$. The method assumes
a binary treatment and uses the experimental data to model the effect of
possibly unmeasured confounders as a function of the measured covariates. The
CATE is obtained from the $D_{o}$ by adding the modeled correction. The main
assumption of the method is that the hidden confounding has an identifiable
parametric structure. The method is implemented for continuous covariates and
outcome and a linear correction function, obtained by solving a least squares
optimization problem. It is not directly applicable to our settings of
categorical covariates, and extending the optimization problem in these
settings is not straightforward. Transportability: Finally, our work has some
connections with the field of transportability Bareinboim and Pearl [2013],
where knowledge of the causal graph is used to determine if the results of an
experimental trial apply to a different population. However, the methods
require knowing the causal graph, and focus on transferring estimators across
distributions rather than combining data to improve estimators.

|  |
---|---|---
Figure 4:  Boxplots of areas under the ROC curve for predicting $Y_{x}$ using
(a) FindIMB (b) IMB and (c) OMB. Data were simulated from random DAGs with 40
observed and 20 latent variables. $D_{o}$ included 10,000 samples, and $D_{e}$
included 100 (left), 200 (middle) or 1000 (right) samples. Black asterisks
denote statistical significance, assessed with the Wilcoxon signed-rank test.
Three stars correspond to $p<0.001$, no stars denote non-significance. Figure
5:  Performance of FindIMB, IMB and OMB for estimating
$P(Y|do(X\texttt{=}1),M)$ with increasing m-bias.

##  6 Experiments

In this section, we show the performance of FindIMB using simulated data. We
simulated random DAGs with a varying number of discrete variables, with mean
in-degree 2. Each DAG includes a binary treatment $X$ and outcome $Y$, where
$X\rightarrow Y$. The remaining covariates $\mathbf{V}$ are pre-treatment and
are binary or ternary, and $1/3$ of the variables are set to be latent. The
observational data $D_{o}$ consist of 10,000 simulated samples from the ground
truth DAGs, and do not include values for the latent variables. Comparison to
other approaches. We compared FindIMB to the following approaches: (a) IMB:
using only experimental data. We used $D_{e}$ to identify the
$\textnormal{MB}_{\overline{X}}(Y)$ using FGESMB. After identifying
${\textnormal{MB}_{\overline{X}}(Y)}\;$, and we used the posterior expectation
$P(Y|do(X),{\textnormal{MB}_{\overline{X}}(Y)}\;\setminus X,D_{e})$ as the
estimator for $P(Y|do(X),\mathbf{V})$. (b) OMB: using only observational data.
We used $\textnormal{{FGESMB}}(D_{o})$ to identify the OMB of $Y$, $MB(Y)$,
and used the posterior expectation $P(Y|do(X),MB(Y)\setminus X,D_{o})$
estimated on $D_{o}$ as the estimator of $P(Y|do(X),\mathbf{V})$. This
estimator is unbiased when conditional ignorability holds for the OMB of $Y$.
(c) FCIt-IMB: Using both observational and experimental data based on
FCItiers: We use FCItiers using as input a data set $D$ constructed by
concatenating $D_{e}$ and $D_{o}$, and adding a binary variable
$I_{e}\rightarrow X$ that corresponds to the presence or absence of
manipulation of $X$. So, $I_{e}=1$ for samples in $D_{e}$ and $I_{e}=0$
otherwise. FCIt-IMB outputs a PAG $\mathcal{P}$ representing all possible
underlying SMCMs. Let $\mathcal{P}_{\overline{X}}$ denote the corresponding
manipulated PAG. We then take the Markov boundary of $Y$ in $P_{\overline{X}}$
to be the $IMB_{X}(Y)$. After identifying $IMB_{X}(Y)$, we test if it is a
backdoor set in $\mathcal{P}_{\overline{X}}$. If so, we used both $D_{e}$ and
$D_{o}$ pooled together to estimate $P(Y|do(X),IMB_{X}(Y)\setminus X)$.
Otherwise, we only used $D_{e}$.

First, we tested if our FindIMB method improves estimation of the probability
$P(Y|do(X),\mathbf{V})$. We simulated DAGs with 10 observed and 5 latent
variables, and applied the methods described above. Each method outputs a set
of variables $\mathbf{Z}$, that is used as an estimate for
$P(Y|do(X),\mathbf{Z},\mathbf{V}\setminus Z).$ Notice that even for 10
variables, the number of possible configurations of $\mathbf{Z}$ can be very
large, and some of these configurations may be very rare. To avoid computing
these parameters for all possible configurations, we tested the methods in a
test dataset $D_{e}^{test}$, that includes 1000 treatment and 1000 control
cases simulated from the manipulated ground truth graph. For each sample in
$D_{e}^{test}$, we obtained an estimate $\hat{P}(Y|do(X),\mathbf{V})$ with the
four methods. The ground truth probability was estimated from the original
manipulated Bayesian network with the junction tree algorithm. We then
computed the average absolute bias
$|\hat{P}(Y|do(X),\mathbf{V})-P(Y|do(X),\mathbf{V})|$ over all test samples.
Fig. 3 shows that FindIMB produced the most accurate probabilities, compared
to using only observational or only experimental data. Moreover, FindIMB
outperforms FCIt-IMB ($p<10^{-3}$ in all cases for a left-tailed t-test). One
reason is that FCIt-IMB selects much larger IMBs than FindIMB, possibly due to
error propagation that results in many bidirected edges. Thus, the resulting
parameters are estimated based on much fewer samples.

We also tested the scalability of FindIMB, using DAGs with 40 observed and 20
latent variables. For this experiment, we could not test against FCIt-IMB
because the method results in very large IMBs. For the same reason, we could
not estimate the true parameters $P(Y|do(X),\mathbf{V})$ for computational
reasons, since the ground truth IMBs can also be very large and include rare
configurations. For this reason, instead of measuring the average absolute
bias, we measured the performance of the methods to classify the test samples
correctly. Fig. 4 shows the area under the ROC curve (AUC) of the models based
FindIMB, IMB, and OMB. FindIMB performs on par or better than the two
alternatives. Average running time for FindIMB (learning the model) was
$1.71\pm 2.46$ seconds per iteration.

One interesting finding is that in all experiments, using the observational
data only, performs better than using experimental data and often is close to
the performance when combining $D_{e}$ and $D_{o}$ using FindIMB. This happens
because in random graphs, the effect of variables inducing bias is often
negligible [Greenland, 2003], and proxies of the unmeasured confounders are
often included in the observed covariates. However, there are cases where the
conditioning on an observed variable in observational data can produce heavily
biased post-intervention probability estimates. A very simple example is the
graph in supplementary Fig. S1), which we call the ”m-bias” graph. To
illustrate how m-bias can affect the prediction of $Y_{x}$ from observational
data, we simulated data from the m-bias graph with binary variables. We set
$Y,X$ and $M$ to be noisy-AND functions of their parents with a parameter
$\alpha$. $\alpha$ has a monotonic relationship with the bias in estimating
$P(Y|do(X=1),M)$ using observational data: Larger $\alpha$ leads to a larger
bias. We then varied alpha from $0.5$ to $1$, and we simulated $D_{o}$ and
$D_{e}$ with 10,000 and 1000 samples, respectively. We used FindIMB, IMB and
OMB to estimate $P(Y|do(X=1),M)$ in a test data set. Fig. 5 shows the bias in
the estimated parameter. We can see that while using $D_{o}$ to estimate
$P(Y|do(X=1),M)$ leads to increasing bias, combining $D_{e}$ and $D_{o}$ can
identify the situations where the parameter is not identifiable from
observational data. We believe that noisy-AND types of distributions are not
rare in biomedical data.

##  7 Discussion

Our paper extends the concepts of Markov boundaries for predicting post-
intervention distributions, and presents a method for learning such Markov
boundaries from mixtures of observational and experimental data. The method
could be useful in settings like embedded trials, where we have abundant
observational and limited experimental data. Future work includes extensions
of the method to mixed data, overlapping covariate sets in observational and
experimental data, and to non-singleton treatments and outcomes.

## References

  * Andrews et al. [2020] Bryan Andrews, Peter Spirtes, and Gregory F Cooper.  On the completeness of causal discovery in the presence of latent confounding with tiered background knowledge.  In _International Conference on Artificial Intelligence and Statistics (AISTATS)_ , pages 4002–4011. PMLR, 2020.
  * Angus [2015] Derek C Angus.  Fusing randomized trials with big data: The key to self-learning health care systems?  _Journal of American Medical Association (JAMA)_ , 314(8):767–768, 2015.
  * Angus et al. [2020] Derek C Angus, Scott Berry, Roger J Lewis, Farah Al-Beidh, Yaseen Arabi, Wilma van Bentum-Puijk, Zahra Bhimani, Marc Bonten, Kristine Broglio, Frank Brunkhorst, et al.  The REMAP-CAP (randomized embedded multifactorial adaptive platform for community-acquired pneumonia) study rationale and design.  _Annals of the American Thoracic Society_ , 17(7):879–891, 2020.
  * Bareinboim and Pearl [2013] Elias Bareinboim and Judea Pearl.  A general algorithm for deciding transportability of experimental results.  _Journal of Causal Inference_ , 1(1):107–134, 2013.
  * Cover [1999] Thomas M Cover.  _Elements of Information Theory_.  John Wiley & Sons, 1999.
  * Greenland [2003] Sander Greenland.  Quantifying biases in causal models: Classical confounding vs collider-stratification bias.  _Epidemiology_ , 14(3):300–306, 2003.
  * Heckerman et al. [1995] David Heckerman, Dan Geiger, and David M Chickering.  Learning Bayesian networks: The combination of knowledge and statistical data.  _Machine Learning_ , 20(3):197–243, 1995.
  * Henckel et al. [2019] Leonard Henckel, Emilija Perković, and Marloes H. Maathuis.  Graphical criteria for efficient total effect estimation via adjustment in causal linear models.  _arXiv preprint arXiv:1907.02435_ , 2019.
  * Hyttinen et al. [2014] Antti Hyttinen, Frederick Eberhardt, and Matti Järvisalo.  Constraint-based causal discovery: Conflict resolution with answer set programming.  In _Proceedings of the 30th Conference on Uncertainty in Artificial Intelligence (UAI)_ , pages 340–349, 2014.
  * Hyttinen et al. [2015] Antti Hyttinen, Frederick Eberhardt, and Matti Järvisalo.  Do-calculus when the true graph is unknown.  In _Proceedings of the 31st Conference on Uncertainty in Artificial Intelligence (UAI)_ , pages 395–404, 2015.
  * Jaber et al. [2019] Amin Jaber, Jiji Zhang, and Elias Bareinboim.  Causal identification under Markov equivalence: Completeness results.  In _Proceedings of the 36th International Conference on Machine Learning (ICML)_ , pages 2981–2989, 2019.
  * Kallus et al. [2018] Nathan Kallus, Aahlad Manas Puli, and Uri Shalit.  Removing hidden confounding by experimental grounding.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , pages 10888–10897, 2018.
  * Mooij et al. [2020] JM Mooij, S Magliacane, and T Claassen.  Joint causal inference from multiple contexts.  _Journal of Machine Learning Research_ , 21(99):1–108, 2020.
  * Pearl [2000] J Pearl.  _Causality: Models, Reasoning and Inference_ , volume 113 of _Hardcover_.  Cambridge University Press, 2000.
  * Pearl [1988] Judea Pearl.  _Probabilistic Reasoning in Intelligent Systems: Networks of Plausible Inference_.  Morgan Kaufmann Publishers Inc., 1988.
  * Pellet and Elisseeff [2008] Jean-Philippe Pellet and André Elisseeff.  Finding latent causes in causal networks: An efficient approach based on Markov blankets.  In _Advances in Neural Information Processing Systems (NeurIPS)_ , pages 1249–1256, 2008.
  * Perkovic et al. [2017] Emilija Perkovic, Johannes Textor, Markus Kalisch, and Marloes H Maathuis.  Complete graphical characterization and construction of adjustment sets in Markov equivalence classes of ancestral graphs.  _Journal of Machine Learning Research_ , 18(1):8132–8193, 2017.
  * Richardson [2003] Thomas Richardson.  Markov properties for acyclic directed mixed graphs.  _Scandinavian Journal of Statistics_ , 30(1):145–157, 2003.  ISSN 03036898.
  * Richardson et al. [2002] Thomas Richardson, Peter Spirtes, et al.  Ancestral graph Markov models.  _The Annals of Statistics_ , 30(4):962–1030, 2002\.
  * Rotnitzky and Smucler [2019] Andrea Rotnitzky and Ezequiel Smucler.  Efficient adjustment sets for population average treatment effect estimation in non-parametric causal graphical models, 2019.
  * Rotnitzky and Smucler [2020] Andrea Rotnitzky and Ezequiel Smucler.  Efficient adjustment sets for population average causal treatment effect estimation in graphical models.  _Journal of Machine Learning Research_ , 21(188):1–86, 2020.  URL http://jmlr.org/papers/v21/19-1026.html.
  * Shpitser and Pearl [2006a] Ilya Shpitser and Judea Pearl.  Identification of joint interventional distributions in recursive semi-Markovian causal models.  In _In proceedings of the 21st National Conference on Artificial Intelligence_ , pages 1219–1226, 2006a.
  * Shpitser and Pearl [2006b] Ilya Shpitser and Judea Pearl.  Identification of conditional interventional distributions.  In _Proceedings of the 22nd Conference on Uncertainty in Artificial Intelligence (UAI)_ , pages 437–444, 2006b.
  * Shpitser et al. [2012] Ilya Shpitser, Tyler VanderWeele, and James M Robins.  On the validity of covariate adjustment for estimating causal effects.  _arXiv preprint arXiv:1203.3515_ , 2012.
  * Tian and Shpitser [2003] Jin Tian and Ilya Shpitser.  On the identification of causal effects.  Technical report, Cognitive Systems Laboratory, University of California at Los Angeles, 2003.
  * Triantafillou and Tsamardinos [2015] Sofia Triantafillou and Ioannis Tsamardinos.  Constraint-based causal discovery from multiple interventions over overlapping variable sets.  _Journal of Machine Learning Research_ , 16(66):2147–2205, 2015.
  * Tsamardinos and Aliferis [2003] Ioannis Tsamardinos and Constantin F Aliferis.  Towards principled feature selection: Relevancy, filters and wrappers.  In _International Conference on Artificial Intelligence and Statistics (AISTATS)_. Citeseer, 2003.
  * Van der Zander et al. [2014] Benito Van der Zander, Maciej Liskiewicz, and Johannes Textor.  Constructing separators and adjustment sets in ancestral graphs.  In _Proceedings of the 30th Conference on Uncertainty in Artificial Intelligence (UAI)_ , pages 11–24, 2014.
  * Yu et al. [2018] K. Yu, L. Liu, J. Li, and H. Chen.  Mining Markov blankets without causal sufficiency.  _IEEE Transactions on Neural Networks and Learning Systems_ , 29(12):6333–6347, 2018.  doi: 10.1109/TNNLS.2018.2828982.
  * Yu et al. [2020] K. Yu, L. Liu, and J. Li.  Learning Markov blankets from multiple interventional data sets.  _IEEE Transactions on Neural Networks and Learning Systems_ , 31(6):2005–2019, 2020.  doi: 10.1109/TNNLS.2019.2927636.

Supplementary Materials for: Causal Markov Boundaries

Eq. number | Analytical Expression
---|---
Eq. 5 | $\displaystyle P(D_{e}|D_{o},{H_{\mathbf{Z}}^{c}}\;)=\prod_{j=1}^{q}\frac{\Gamma(\alpha_{j}+N_{j}^{o})}{\Gamma(\alpha_{j}+N_{j}^{o}+N_{j}^{e})}\prod_{k=1}^{r}\frac{\Gamma(\alpha_{jk}+N_{jk}^{o}+N_{jk}^{e})}{\Gamma(\alpha_{jk}+N_{jk}^{o})}$
- | $\displaystyle P(D_{e}|D_{o},{H_{\mathbf{Z}}^{\overline{c}}}\;)=\prod_{j=1}^{q}\frac{\Gamma(\alpha_{j})}{\Gamma(\alpha_{j}+N_{j}^{e})}\prod_{k=1}^{r}\frac{\Gamma(\alpha_{jk}+N_{jk}^{e})}{\Gamma(\alpha_{jk})}$
Eq. 7 | $\displaystyle P(D_{o}|{H_{\mathbf{Z}}^{c}}\;)=\prod_{j=1}^{\tilde{q}}\frac{\Gamma(\tilde{\alpha}_{j})}{\Gamma(\tilde{\alpha}_{j}+\tilde{N}_{j}^{o})}\prod_{k=1}^{r}\frac{\Gamma(\tilde{\alpha}_{jk}+\tilde{N}_{jk}^{o})}{\Gamma(\tilde{\alpha}_{jk})}$
Eq. 7 | $\displaystyle P(D_{o}|{H_{\mathbf{Z}}^{\overline{c}}}\;)=\prod_{j=1}^{\tilde{q}}\frac{\Gamma(\tilde{\alpha}_{j})}{\Gamma(\tilde{\alpha}_{j}+\tilde{N}_{j}^{o})}\prod_{k=1}^{r}\frac{\Gamma(\tilde{\alpha}_{jk}+\tilde{N}_{jk}^{o})}{\Gamma(\tilde{\alpha}_{jk})}$
Terms in Eq. 8 | $\displaystyle P(Y=k|x,Z=j,D_{e},D_{o},{H_{\mathbf{Z}}^{c}}\;)=\frac{N_{jk}^{o}+N_{jk}^{e}+\alpha_{jk}}{N_{j}^{o}+N_{j}^{e}+\alpha_{j}}$
Terms in Eq. 8 | $\displaystyle P(Y=k|x,Z=j,D_{e},D_{o},{H_{\mathbf{Z}}^{\overline{c}}}\;)=\frac{N_{jk}^{e}+\alpha_{jk}}{N_{j}^{e}+\alpha_{j}}$

Table 1: Closed-form solutions for Eq. 4 -8

in the main paper, for multinomial distributions with Dirichlet priors.
Subscript $jk$ refers variable $Y$ taking its $k$-th configuration, and
variable set $\mathbf{Z}$ taking its $j$-th configuration. $\alpha_{jk}$ is
the prior for the Dirichlet distribution. We set $\alpha_{jk}=1$ in all
experiments. $N^{o}_{jk},N^{e}_{jk}$ corresponds to counts in the data where
$Y=k$ and $\mathbf{Z}=j$ in $D_{o}$ and $D_{e}$, respectively.
$N^{o}_{j},N^{e}_{j}$ corresponds to counts in the data where $Z=j$. Tilde
notation corresponds to the OMB $\mathbf{U}$.

$X$$A$$B$$Y$ Figure S1: An example where a CMB does not necessarily correspond
to an optimal adjustment set [Henckel et al., 2019].
${\textnormal{{CMB}}_{\overline{X}}(Y)}\;=\\{\\{X,A,B\\}\\}$, but the optimal
adjustment set depends on the parameters. $X$$A$$B$$M$$Y$ Figure S1: The
m-bias graph used to simulate data for Fig 5. $A$ and $B$ are unobserved. All
variables were binary. Parameters were as follows:
$P(A=1)=.8,P(B=1)=.8,P(M=1|A=1,B=1)=\alpha,P(X=1|A=1)=\alpha,P(Y=1|X=1,B=1)=\alpha$.
All other parameters $P(Y=1|\dots),P(M=1|\dots),P(X=1|A=0)$ were set to zero.

##  8 Proofs

In this section, we provide a proof that every causal Markov boundary is
backdoor set, which is defined below (Definition 8.1). We make the following
assumptions throughout the entire document:

  * •

$X$ causes $Y$

  * •

all variables $\mathbf{V}$ are pre-treatment.

######  Definition 8.1 (Backdoor Set).

$\mathbf{Z}$ is a backdoor set for $X$, $Y$ if and only if $\mathbf{Z}$
m-separates $X$ and $Y$ in $\mathcal{G}_{\underline{X}}$.

We use the following definitions from [Shpitser and Pearl, 2006a]:

######  Definition 8.2 (C-component).

A C-component is as set of nodes $S$ in $\mathcal{G}$ where every two nodes
are connected by a bidirected path.

######  Definition 8.3 (C-forest).

A graph $\mathcal{G}$ where the set of all of its nodes is a C-component, and
each node has at most one child is a C-forest. The set of nodes $\mathbf{R}$
without children in the C-forest is called the root, and we say that
$\mathcal{G}$ is an $\mathbf{R}$-rooted C-forest.

C-forests are useful for defining hedges:

######  Definition 8.4 (hedge).

Let $\mathbf{X}$,$\mathbf{Y}$ be sets of variables in $\mathcal{G}$. Let
$F,F^{\prime}$ be $\mathbf{R}$-rooted C-forests in $\mathcal{G}$ such that
$F^{\prime}$ is a subgraph of $F$, $\mathbf{X}$ only occurs in $F$, and
$\mathbf{R}\in An(\mathbf{Y})_{\mathcal{G}_{\overline{X}}}$. Then
$F,F^{\prime}$ form a hedge for $P(\mathbf{Y}|do(\mathbf{X}))$.

The existence of a hedge for $P(\mathbf{Y}|do(\mathbf{X}))$ in $\mathcal{G}$
is equivalent to the non-identifiability of $P(\mathbf{Y}|do(\mathbf{X}))$
(see Theorem 4 in [Shpitser and Pearl, 2006a]).

######  Lemma 8.5.

Let $\mathbf{Z}$ be a set that is not a subset of any backdoor set (i.e.,
there exists no set $\mathbf{Q}\subseteq(\mathbf{V}\setminus\mathbf{Z})$ such
that $\mathbf{Q}\cup\mathbf{Z}$ m-separate $X$ and $Y$ in
$\mathcal{G}_{\underline{X}}$). Then there exists in $\mathcal{G}$ a bi-
directed path from $X$ to $Y$ where every collider has a descendant in
$\mathbf{Z}\cup Y$.

###### Proof.

The proof is a special case of Theorem 4.2 $(iv)\Rightarrow(ii)$ in
[Richardson et al., 2002] with
$\mathbf{S}\leftarrow\mathbf{Z},\mathbf{L}\leftarrow\emptyset,\mathcal{G}\leftarrow\mathcal{G}_{\underline{X}}$.
The proof is for ancestral graphs, but it is straightforward to show that it
holds for SMCMs, given that every SMCM $\mathcal{G}$ can be transformed to a
maximal ancestral graph $\mathcal{M}$ over the same nodes (by adding some
edges) such that (a) $\mathcal{G}$ and $\mathcal{M}$ entail the exact same
m-separations and m-connections and (b) the exact same ancestral relationships
hold in both graphs. The theorem proves that if
$\forall\mathbf{Q}\subseteq\mathbf{(}V\setminus\mathbf{Z}),\mathbf{Z}\cup\mathbf{Q}$
do not m-separate $X$ and $Y$ in $\mathcal{G}_{\underline{X}}$, then there
exists a bidirected path between $X$ and $Y$ in $\mathcal{G}_{\underline{X}}$
where every variable is an ancestor of some variables in
$\mathbf{Z}\cup\\{X,Y\\}$, which means that there exists a path in
$\mathcal{G}$ a bi-directed path from $X$ to $Y$ where every collider has a
descendant in $\mathbf{Z}\cup Y$ (since $X\rightarrow Y$ by assumption). ∎

######  Lemma 8.6.

Let $\mathbf{Z}$ be a set for which $P(Y|do(X),\mathbf{Z})$ is identifiable
from $P(Y|X,\mathbf{Z})$, then $\mathbf{Z}$ is a subset of a backdoor set.

###### Proof.

First, notice that
$P(Y|do(X),\mathbf{Z})=\frac{P(Y,\mathbf{Z}|do(X))}{P(Z|do(X))}=\frac{P(Y,\mathbf{Z}|do(X))}{P(Z)}$.
Therefore $P(Y|do(X),\mathbf{Z})$ is only identifiable if
$P(Y,\mathbf{Z}|do(X))$ is identifiable. If $\mathbf{Z}$ is not a subset of a
backdoor set, then there exists a bidirected path where every variable has a
descendant in $\mathbf{Z}\cup Y$ in $\mathcal{G}$ by Lemma 8.5. Let
$\mathcal{F}$ be the graph consisting of the bidirected path, and
$\mathcal{F}$’ be the same graph without $X$. Then $\mathcal{F}$,
$\mathcal{F}$’ are $\\{Y,\mathbf{Z}\\}$ rooted C-forests, and
$\\{Y,\mathbf{Z}\\}\in An(\\{Y,\mathbf{Z}\\})$, so $\mathcal{F}$,
$\mathcal{F}$’ form a hedge for $\\{Y,\mathbf{Z}\\}$. Therefore,
$P(Y,\mathbf{Z}|do(X))$ is not identifiable, and $P(Y|do(X),\mathbf{Z})$ is
not identifiable. ∎

######  Theorem 3.3.

We assume that $P_{x}$ and $\mathcal{G}_{\overline{X}}$ are faithful to each
other. If $\mathbf{Z}$ is a causal Markov boundary for $Y$ relative to $X$,
then $\mathbf{W}=\mathbf{Z}\setminus X$ is a backdoor set.

###### Proof.

Assume $\mathbf{Z}$ is a causal Markov boundary, but $\mathbf{W}$ is not a
backdoor set. Since $P(Y|do(X),\mathbf{W})$ is identifiable, by Lemma 8.6
$\mathbf{W}$ is a subset of a backdoor set $\mathbf{W}\cup\mathbf{Q}$, where
$\mathbf{Q}\subseteq(\mathbf{V}\setminus\mathbf{W})$. Since by assumption
$\mathbf{W}$ is not a backdoor set, $\mathbf{Q}$ is not the empty set (i.e.,
$\mathbf{W}$ is a proper subset of a backdoor set). We will show that
$P(Y|do(X),\mathbf{W},\mathbf{Q})\neq P(Y|do(X),\mathbf{W})$. To show that, we
only need to show that $\mathbf{Q}$ is not independent from $\mathbf{W}$ in
$\mathcal{G}_{\overline{X}}$ . Since $\mathbf{W}$ is not a backdoor set, there
exists a backdoor path from $X$ to $Y$ that is m$-$connecting given
$\mathbf{W}$, but blocked given $\mathbf{W}\cup\mathbf{Q}$. Thus, some
$Q\in\mathbf{Q}$ is a non-collider on that path, therefore $\mathbf{Q}$ are
not independent with $Y$ given $\mathbf{W}$. Hence,
$P(Y|do(X),\mathbf{W},\mathbf{Q})\neq P(Y|do(X),\mathbf{W})$ and therefore
$\mathbf{Z}$ does not satisfy Condition (2), and $\mathbf{Z}$ is not a causal
Markov boundary (Contradiction). ∎

######  Lemma 8.7.

Let $\mathbf{Z}\subseteq\mathbf{V}$ be a backdoor set for $X,Y$, and let
$Q\in(\mathbf{Z}\setminus\textnormal{MB}(Y))$ that has an m-connecting path
$Q\pi_{QY}Y$ with $Y$ given $\mathbf{Z}\setminus Q$. Then there exists a
variable $W\in(\textnormal{MB}(Y)\setminus\mathbf{Z})$ such that:
$W\cup\mathbf{Z}$ is a backdoor set and
$W{\,\not\mkern-1.0mu\perp\mkern-12.0mu\perp\,}Y|\mathbf{Z}$ in
$\mathcal{G}_{\overline{X}}$.

###### Proof.

Let $Q$ be a variable as described above. Then there exists a variable
$W\in\textnormal{MB}(Y)$ between $Q$ and $Y$ that is a non-collider on $\pi$,
otherwise $Q\in Pa(Dis(Y))$, and therefore $Q\in\textnormal{MB}(Y)$. In
addition, $W\not\in\mathbf{Z}$, otherwise $Q\pi_{QY}Y$ would be blocked given
$\mathbf{Z}\setminus Q$. We will now show, by contradiction, that adding $W$
to the conditioning set $\mathbf{Z}$ does not open any backdoor paths from $X$
to $Y$; hence, $\mathbf{Z}\cup W$ is a backdoor set.

Assume that conditioning on $W$ opens a path $\pi_{XY}$ between $X$ and $Y$
that is blocked given just $\mathbf{Z}$. Then $W$ must be a descendant of one
or more colliders on that path. Let $C$ be the collider closest to $X$ on
$\pi_{XY}$ such that $C$ is blocked on $\pi_{XY}$ given $\mathbf{Z}$, but open
given $\mathbf{Z}\cup W$. Then $X\pi_{XC}C$ is open given $\mathbf{Z}$, and
$W$ is a descendant of $C$. Let $C\pi_{CW}W$ be the (possibly empty) directed
path from $C$ to $W$, and let $W\pi_{WY}Y$ be the subpath of $\pi_{CY}$ from
$W$ to $Y$. Since $C$ is blocked on $\pi_{XY}$ given $\mathbf{Z}$, no variable
on $\pi_{CW}$ can be in $\mathbf{Z}$. But then $X\pi_{XC}C\pi_{CW}W\pi_{WY}\
Y$ is an open path from $X$ and $Y$ given $\mathbf{Z}$ in
$\mathcal{G}_{\overline{X}}$. Contradiction, since $\mathbf{Z}$ is a backdoor
set. Thus, $W$ does not open any backdoor paths, and $\mathbf{Z}\cup W$ is
also a backdoor set.

Finally, $W$ is not independent of $Y$ given $\mathbf{Z}$ in
$\mathcal{G}_{\overline{X}}$, since $W\pi_{WY}Y$ is open given $\mathbf{Z}$. ∎

######  Theorem 3.4.

We assume that $P_{x}$ and $\mathcal{G}_{\overline{X}}$ are faithful to each
other. Every causal Markov boundary $\mathbf{Z}$ of an outcome variable $Y$
w.r.t a treatment variable $X$ is a subset of the Markov boundary
$\textnormal{MB}(Y)$.

###### Proof.

We will show this by contradiction. Specifically, we will show that any set
$\mathbf{Z}$ that includes variables $\mathbf{Q}$ not in the Markov boundary
of $Y$ cannot satisfy one of the Conditions (2) or (3) of the causal Markov
boundary.

Assume that $\mathbf{Z}$ is a causal Markov boundary for $Y$ with respect to
$X$. and let $\mathbf{W}=\mathbf{Z}\setminus X$. Let
$\mathbf{Q}=\mathbf{W}\setminus\textnormal{MB}(Y)$ be the non-empty subset of
$\mathbf{W}$ that is not a part of the Markov boundary of $Y$.

If there exists no $Q\in\mathbf{Q}$ that has an m-connecting path $Q\pi_{QY}Y$
to $Y$ given $\mathbf{W}\setminus Q$, then
$\mathbf{Q}{\,\perp\mkern-12.0mu\perp\,}Y|(\mathbf{W}\setminus\mathbf{Q})$ in
$\mathcal{G}_{\overline{X}}$. Conditioning on $X$ cannot open any paths from
$X$ to $Y$; therefore,
$\mathbf{Q}{\,\perp\mkern-12.0mu\perp\,}Y|X,(\mathbf{W}\setminus\mathbf{Q})$
in $\mathcal{G}_{\overline{X}}$. Then by Rule 1 of the do-calculus [Pearl,
2000], $P(Y|do(X),\mathbf{W})=P(Y|do(X),\mathbf{W}\setminus Q)$, and
$\mathbf{Z}$ does not satisfy Condition (3) of the causal Markov boundary
definition (Contradiction).

If there exists a $Q\in(\mathbf{W}\setminus\textnormal{MB}(Y))$ that has an
m-connecting path $Q\pi_{QY}Y$ with $Y$ given $\mathbf{Z}\setminus Q$, then by
Lemma 8.7, there exists a variable $W$ in
$\textnormal{MB}(Y)\setminus\mathbf{Z}$ such that $\mathbf{Z}\cup W$ is also a
backdoor set, and
$W{\,\not\mkern-1.0mu\perp\mkern-12.0mu\perp\,}Y|X,\mathbf{Z}$ in
$\mathcal{G}_{\overline{X}}$. Then $P(Y|do(X),\mathbf{Z},W)\neq
P(Y|do(X),\mathbf{Z})$. Thus, $\mathbf{Z}$ does not satisfy Condition (2) of
the Causal Markov boundary definition (Contradiction).

Thus, $\mathbf{Z}$ cannot include any variables that are not in the Markov
boundary of $Y$. ∎

######  Theorem 3.5.

Let $\mathcal{G}$ be a SMCM over $X$, $Y$,$\mathbf{V}$ with $\mathbf{V}$
occurring before $X$ and $Y$. Let $\mathbf{Z}\subseteq\mathbf{V}\cup X$ be the
IMB of $Y$ relative to $X$. If $\mathbf{Z}$ is a causal Markov boundary, then
$\textnormal{MB}(Y)=\mathbf{Z}$.

###### Proof.

${\textnormal{MB}_{\overline{X}}(Y)}\;\subseteq{\textnormal{MB}(Y)}\;$, so we
need to show that
${\textnormal{MB}(Y)}\;\subseteq{\textnormal{MB}_{\overline{X}}(Y)}\;$ when
${\textnormal{MB}_{\overline{X}}(Y)}\;\in{\textnormal{{CMB}}_{\overline{X}}(Y)}\;$.
Assume that $\mathbf{Z}$ is both the $\textnormal{MB}_{\overline{X}}(Y)$ and a
causal Markov boundary, but there exists a variable $Q$ in $\mathbf{Z}$ that
is not in $\textnormal{MB}(Y)$ . Then $Q$ is reachable from $Y$ through a
bidirected path in $\mathcal{G}$ but not in $\mathcal{G}_{\overline{X}}$.
Since $\mathcal{G}$ and $\mathcal{G}_{\overline{X}}$ only differ in edges that
are into $X$, this path must be going through an edge that is incoming into X.
Thus, $\mathcal{G}$ includes a bidirected path
$Y\leftrightarrow\dots\leftrightarrow X$, and every variable on this path is
in $\textnormal{MB}_{\overline{X}}(Y)$ =$\mathbf{Z}$. But then
$\mathbf{Z}\setminus X$ cannot be a backdoor set, and by Theorem 3.3
$\mathbf{Z}$ cannot be a causal Markov boundary. Contradiction. Thus, the
Markov boundary of $Y$ cannot include any more variables than $\mathbf{Z}$. ∎

##  9 Convergence Proof for Observational Markov Boundary (OMB)

######  Definition 9.1 (Conditional Entropy).

Let $P$ be the full joint probability distribution over a set of variables
$\mathbf{V}$, let $Y\in\mathbf{V}$ be a variable, and let
$\mathbf{Z}\subseteq\mathbf{V}\setminus\\{Y\\}$ be a set of variables. Then,
the conditional entropy of $Y$ given $\mathbf{Z}$ is defined as follows
[Cover, 1999]:

$H(Y|\mathbf{Z})=-\sum_{y}{\sum_{z}{P(y,z)\cdot\log P(y|z)}}$ (S1)

where $y$ and $z$ denote the values of $Y$ and $\mathbf{Z}$, respectively.

######  Lemma 9.2.

Let $X,Y\in\mathbf{V}$ be two variables and
$\mathbf{Z}\subseteq\mathbf{V}\setminus\\{X,Y\\}$ be a set of variables. Then,
$H(Y|\mathbf{Z})\geq H(Y|X,\mathbf{Z})$, where the entropies are defined by
Definition 9.1, and the equality holds if and only if
$Y{\,\perp\mkern-12.0mu\perp\,}X|\mathbf{Z}$.

###### Proof.

Applying the chain rule of entropy, the conditional mutual information can be
computed as follows [Cover, 1999]:

$I(X;Y|\mathbf{Z})=H(Y|\mathbf{Z})-H(Y|X,\mathbf{Z})\,.$ (S2)

Given that the mutual information is nonnegative (i.e., $I(X;Y|\mathbf{Z})\geq
0$) and $I(X;Y|\mathbf{Z})=0$ if and only if
$Y{\,\perp\mkern-12.0mu\perp\,}X|\mathbf{Z}$ (see [Cover, 1999], page 29), it
follows that:

$\begin{split}H(Y|\mathbf{Z})-H(Y|X,\mathbf{Z})\geq 0\\\ H(Y|\mathbf{Z})\geq
H(Y|X,\mathbf{Z})\,,\end{split}$ (S3)

where the equality holds if and only if
$Y{\,\perp\mkern-12.0mu\perp\,}X|\mathbf{Z}$. ∎

For brevity, let $\mathbf{V}=\\{\mathbf{V}\cup X\\}$, where $X$ is a treatment
variable, and let $Y$ be an outcome variable in the remainder of this section.

######  Lemma 9.3.

All Markov blankets of $Y$ have the same entropy.

###### Proof.

By definition, $\mathbf{Z}^{\prime}$ is the Markov blanket of $Y$ if and only
if $P(Y|\mathbf{Z}^{\prime},\mathbf{W})=P(Y|\mathbf{Z}^{\prime})$ for any
$\mathbf{W}\subseteq\mathbf{V}\setminus\mathbf{Z}^{\prime}$, which indicates
that $Y{\,\perp\mkern-12.0mu\perp\,}\mathbf{W}|\mathbf{Z}^{\prime}$. Also,
according to Lemma 9.2,
$H(Y|\mathbf{Z}^{\prime})=H(Y|\mathbf{Z}^{\prime},\mathbf{W})$ for any
$\mathbf{W}\subseteq\mathbf{V}\setminus\mathbf{Z}^{\prime}$. Let $\mathbf{Z}$
also be a Markov blanket of $Y$. By multiple applications of Lemma 9.2, we
obtain:

$H(Y|\mathbf{Z}^{\prime})=H(Y|\mathbf{Z}^{\prime},\mathbf{V}\setminus\mathbf{Z}^{\prime})=H(Y|\mathbf{V})=H(Y|\mathbf{Z},\mathbf{V}\setminus\mathbf{Z}^{\prime})=H(Y|\mathbf{Z})$
(S4)

∎

######  Lemma 9.4.

Let $\mathbf{Z}^{\prime}$ be a Markov blanket of $Y$ and let $\mathbf{Z}$ be a
set of variables that is not a Markov blanket of $Y$. Then,
$H(Y|\mathbf{Z}^{\prime})<H(Y|\mathbf{Z})$, where the entropies are defined by
Definition 9.1.

###### Proof.

Assume there is exists a set
$\mathbf{W}\subseteq\mathbf{V}\setminus\mathbf{Z}$ such that
$P(Y|\mathbf{Z},\mathbf{W})\neq P(Y|\mathbf{Z})$. According to Lemma 9.2 we
have:

$H(Y|\mathbf{Z},W)<H(Y|\mathbf{Z}).$ (S5)

Also, given that $\mathbf{V}$ is a superset of $(\mathbf{Z}\cup\mathbf{W})$,
we have:

$H(Y|\mathbf{V})\leq H(Y|\mathbf{Z},\mathbf{W}).$ (S6)

Therefore,

$H(Y|\mathbf{V})<H(Y|\mathbf{Z}).$ (S7)

Also, since $\mathbf{Z}^{\prime}$ is a Markov blanket of $Y$, by Lemma 9.3 we
have:

$H(Y|\mathbf{Z}^{\prime})=H(Y|\mathbf{V}).$ (S8)

Combining Equations (S7) and (S8), we obtain:

$H(Y|\mathbf{Z}^{\prime})<H(Y|\mathbf{Z}).$ (S9)

∎

######  Lemma 9.5.

Given dataset $D_{o}$ that contains samples from a strictly positive
distribution $P$, which is a perfect map for a SMCM $\mathcal{G}$, the BD
score [Heckerman et al., 1995] for $\log P(D_{o}|\mathbf{Z})$ is defined as
follows in the large sample limit:

$\lim_{N\rightarrow\infty}\log
P(D_{o}|\mathbf{Z})=\lim_{N\rightarrow\infty}-N\cdot
H(Y|\mathbf{Z})-\frac{q\cdot(r-1)}{2}\log N+const.,$ (S10)

###### Proof.

The BD score for $P(D_{o}|\mathbf{Z})$ is calculated as follows [Heckerman et
al., 1995]:

$P(D_{o}|\mathbf{Z})=\prod_{j=1}^{q}{\frac{\mathrm{\Gamma}(\alpha_{j})}{\mathrm{\Gamma}(\alpha_{j}+N_{j})}\cdot\displaystyle\prod_{k=1}^{r}\frac{\mathrm{\Gamma}(\alpha_{jk}+N_{jk})}{\mathrm{\Gamma}(\alpha_{jk})}}\,,$
(S11)

where $q$ denotes instantiations of variables in $\mathbf{Z}$ and $r$ denotes
values of variable $Y$. The term $N_{jk}$ is the number of cases in data in
which variable $Y=k$ and its parent $\mathbf{Z}=j$; also,
$N_{j}=\sum_{k=1}^{r}{N_{jk}}$. The term $\alpha_{jk}$ is a finite positive
real number that is called Dirichlet prior parameter and may be interpreted as
representing “pseudo-counts”, where $\alpha_{j}=\sum_{k=1}^{r}{\alpha_{jk}}$.
BD can be re-written in $log$ form as follows:

$\log
P(D_{o}|\mathbf{Z})=\sum_{j=1}^{q}{\left[{\log\Gamma(\alpha_{j})}-{\log\Gamma(\alpha_{j}+N_{j})}+\sum_{k=1}^{r}{\left[\log\Gamma(\alpha_{jk}+N_{jk})-\log\Gamma(\alpha_{jk})\right]}\right]}.$
(S12)

We can re-arrange the terms in Eq. (S12) to gather the constant terms as
follows:

$\begin{split}\log
P(D_{o}|\mathbf{Z})&=\sum_{j=1}^{q}{\left[-{\log\Gamma(\alpha_{j}+N_{j})}+\sum_{k=1}^{r}{\log\mathrm{\Gamma}(\alpha_{jk}+N_{jk})}\right]}+\sum_{j=1}^{q}{\left[\log\mathrm{\Gamma}(\alpha_{j})-\sum_{k=1}^{r}{\log\mathrm{\Gamma}(\alpha_{jk})}\right]}\\\
&=\sum_{j=1}^{q}{\left[-{\log\Gamma(\alpha_{j}+N_{j})}+\sum_{k=1}^{r}{\log\mathrm{\Gamma}(\alpha_{jk}+N_{jk})}\right]}+const.\end{split}$
(S13)

Using the Stirling’s approximation of
$\lim_{n\rightarrow\infty}\log\mathrm{\Gamma}(n)=(n-\frac{1}{2})\log(n)-n+const.$,
we can re-write Eq. (S13) as follows:

$\begin{split}&\lim_{N\rightarrow\infty}\log P(D_{o}|\mathbf{Z})\\\
&=\lim_{N\rightarrow\infty}\sum_{j=1}^{q}\Biggl{[}-(\alpha_{j}+N_{j}-\frac{1}{2})\log(\alpha_{j}+N_{j})+(\alpha_{j}+N_{j})+\sum_{k=1}^{r}{\left((\alpha_{jk}+N_{jk}-\frac{1}{2})\log(\alpha_{jk}+N_{jk})-(\alpha_{jk}+N_{jk})\right)}\Biggl{]}+const.\\\
&=\lim_{N\rightarrow\infty}\sum_{j=1}^{q}\Biggl{[}-\alpha_{j}\log(\alpha_{j}+N_{j})-N_{j}\log(\alpha_{j}+N_{j})+\frac{1}{2}\log(\alpha_{j}+N_{j})+\alpha_{j}+N_{j}\\\
&+\sum_{k=1}^{r}{\left(\alpha_{jk}\log(\alpha_{jk}+N_{jk})+N_{jk}\log(\alpha_{jk}+N_{jk})-\frac{1}{2}\log(\alpha_{jk}+N_{jk})-\alpha_{jk}-N_{jk}\right)}\Biggl{]}+const.\\\
&=\lim_{N\rightarrow\infty}\sum_{j=1}^{q}\Biggl{[}-N_{j}\log(\alpha_{j}+N_{j})+\sum_{k=1}^{r}N_{jk}\log(\alpha_{jk}+N_{jk})\Biggl{]}+\sum_{j=1}^{q}\Biggl{[}-\alpha_{j}\log(\alpha_{j}+N_{j})+\sum_{k=1}^{r}\alpha_{jk}\log(\alpha_{jk}+N_{jk})\Biggl{]}\\\
&+\frac{1}{2}\sum_{j=1}^{q}\Biggl{[}\log(\alpha_{j}+N_{j})-\sum_{k=1}^{r}\log(\alpha_{jk}+N_{jk})+\alpha_{j}+N_{j}-\sum_{k=1}^{r}\left(\alpha_{jk}+N_{jk}\right)\Biggl{]}+const.\\\
&=\lim_{N\rightarrow\infty}\sum_{j=1}^{q}\Biggl{[}-N_{j}\log(\alpha_{j}+N_{j})+\sum_{k=1}^{r}N_{jk}\log(\alpha_{jk}+N_{jk})\Biggl{]}+\sum_{j=1}^{q}\Biggl{[}-\alpha_{j}\log(\alpha_{j}+N_{j})+\sum_{k=1}^{r}\alpha_{jk}\log(\alpha_{jk}+N_{jk})\Biggl{]}\\\
&+\frac{1}{2}\sum_{j=1}^{q}\Biggl{[}\log(\alpha_{j}+N_{j})-\sum_{k=1}^{r}\log(\alpha_{jk}+N_{jk})\Biggl{]}+const.\end{split}$
(S14)

In the last step of Eq. (S14), we used the facts that
$\sum_{k=1}^{r}{N_{jk}}=N_{j}$ and $\sum_{k=1}^{r}{\alpha_{jk}}=\alpha_{j}$,
and we applied these identities again to that equation to obtain the
following:

$\begin{split}&\lim_{N\rightarrow\infty}\log P(D_{o}|\mathbf{Z})=\\\
&\lim_{N\rightarrow\infty}\sum_{j=1}^{q}\sum_{k=1}^{r}\Biggl{[}N_{jk}\log(\frac{\alpha_{jk}+N_{jk}}{\alpha_{j}+N_{j}})+\alpha_{jk}\log(\frac{\alpha_{jk}+N_{jk}}{\alpha_{j}+N_{j}})\Biggl{]}+\frac{1}{2}\sum_{j=1}^{q}\Biggl{[}\log(\alpha_{j}+N_{j})-\sum_{k=1}^{r}\log(\alpha_{jk}+N_{jk})\Biggl{]}+const.\end{split}$
(S15)

Given that

$\lim_{N\rightarrow\infty}\frac{\alpha_{jk}+N_{jk}}{\alpha_{j}+N_{j}}=\frac{N_{jk}}{N_{j}}$

and

$\lim_{N\rightarrow\infty}\sum_{j=1}^{q}\sum_{k=1}^{r}\alpha_{jk}\log(\frac{\alpha_{jk}+N_{jk}}{\alpha_{j}+N_{j}})=const.,$

in the limit, Eq. (S15) becomes:

$\lim_{N\rightarrow\infty}\log
P(D_{o}|\mathbf{Z})=\lim_{N\rightarrow\infty}\sum_{j=1}^{q}\sum_{k=1}^{r}N_{jk}\log\frac{N_{jk}}{N_{j}}+\frac{1}{2}\sum_{j=1}^{q}\Biggl{[}\log(\alpha_{j}+N_{j})-\sum_{k=1}^{r}\log(\alpha_{jk}+N_{jk})\Biggl{]}+const.,$
(S16)

or equivalently:

$\begin{split}\lim_{N\rightarrow\infty}\log
P(D_{o}|\mathbf{Z})&=\lim_{N\rightarrow\infty}N\cdot\sum_{j=1}^{q}\sum_{k=1}^{r}\frac{N_{jk}}{N}\log\frac{N_{jk}}{N_{j}}+\frac{1}{2}\sum_{j=1}^{q}\Biggl{[}\log(\alpha_{j}+N_{j})-\sum_{k=1}^{r}\log(\alpha_{jk}+N_{jk})\Biggl{]}+const.\\\
&=\lim_{N\rightarrow\infty}-N\cdot
H(Y|\mathbf{Z})+\frac{1}{2}\sum_{j=1}^{q}\Biggl{[}\log(\alpha_{j}+N_{j})-\sum_{k=1}^{r}\log(\alpha_{jk}+N_{jk})\Biggl{]}+const.\end{split}$
(S17)

To simplify the second term in Eq. (S17), we divide the arguments in the log
terms by $N$ and equivalently add $\log N$ terms as follows:

$\begin{split}\lim_{N\rightarrow\infty}\frac{1}{2}\sum_{j=1}^{q}\Biggl{[}\log(\alpha_{j}+N_{j})&-\sum_{k=1}^{r}\log(\alpha_{jk}+N_{jk})\Biggl{]}=\lim_{N\rightarrow\infty}\frac{1}{2}\sum_{j=1}^{q}\Biggl{[}\log(\frac{\alpha_{j}+N_{j}}{N})+\log
N-\sum_{k=1}^{r}\log(\frac{\alpha_{jk}+N_{jk}}{N})+\log N\Biggl{]}\\\
&=\lim_{N\rightarrow\infty}\frac{1}{2}\sum_{j=1}^{q}\left(\log
N-\sum_{k=1}^{r}\log
N\right)+\frac{1}{2}\sum_{j=1}^{q}\Biggl{[}\log(\frac{\alpha_{j}+N_{j}}{N})-\sum_{k=1}^{r}\log(\frac{\alpha_{jk}+N_{jk}}{N})\Biggl{]}\\\
&=-\frac{q(r-1)}{2}\log N+const.\end{split}$ (S18)

Combining Equations (S17) and (S18), we obtain:

$\lim_{N\rightarrow\infty}\log
P(D_{o}|\mathbf{Z})=\lim_{N\rightarrow\infty}-N\cdot
H(Y|\mathbf{Z}^{\prime})-\frac{q\cdot(r-1)}{2}\log N+const.$ (S19)

∎

######  Theorem 4.1.

Given dataset $D_{o}$ that contains samples from a strictly positive
distribution $P$, which is a perfect map for a SMCM $\mathcal{G}$, the BD
score [Heckerman et al., 1995] will assign the highest score to the OMB of $Y$
in the large sample limit.

###### Proof.

Let $\mathbf{Z}^{\prime}$ be the OMB of $Y$ and
$\mathbf{Z}\subseteq\mathbf{V}$ be an arbitrary set. We want to show that:

$\lim_{N\rightarrow\infty}\frac{P(D_{o}|\mathbf{Z})}{P(D_{o}|\mathbf{Z}^{\prime})}=\begin{cases}1&\text{iff
$\mathbf{Z}$ is an OMB of $Y$}\\\ 0&\text{otherwise}\end{cases}\enspace,$
(S20)

Applying Lemma 9.5 we have:

$\lim_{N\rightarrow\infty}\log\frac{P(D_{o}|\mathbf{Z})}{P(D_{o}|\mathbf{Z}^{\prime})}=\lim_{N\rightarrow\infty}N\cdot\left[H(Y|\mathbf{Z}^{\prime})-H(Y|\mathbf{Z})\right]+\frac{(q^{\prime}-q)\cdot(r-1)}{2}\log
N.$ (S21)

where $q$ and $q^{\prime}$ are the number of possible parent instantiations of
$Y$ with $\mathbf{Z}$ and $\mathbf{Z}^{\prime}$ as the set of parents. There
are three possible cases:

Case 1: $\mathbf{Z}$ is a Markov blanket of $Y$ and its OMB.

Since both $\mathbf{Z}^{\prime}$ and $\mathbf{Z}$ are Markov blankets of $Y$,
$H(Y|\mathbf{Z})=H(Y|\mathbf{Z}^{\prime})$ by Lemma 9.3. Thus, the first term
in Eq. (S21) becomes $0$. Also, given that $\mathbf{Z}^{\prime}$ and
$\mathbf{Z}$ are OMBs, they have the same number of parameters $q^{\prime}=q$,
by which the second term in Eq. (S21) becomes $0$ in the limit as
$N\rightarrow\infty$, or equivalently Eq. (S20) approaches to 1.

Case 2: $\mathbf{Z}$ is a Markov blanket of $Y$ but not its OMB.

According to Lemma 9.3 $H(Y|\mathbf{Z})=H(Y|\mathbf{Z}^{\prime})$; therefore,
the first term in Eq. (S21) becomes $0$ and we obtain:

$\lim_{N\rightarrow\infty}\frac{P(D_{o}|\mathbf{Z})}{P(D_{o}|\mathbf{Z}^{\prime})}=\lim_{N\rightarrow\infty}\frac{(q^{\prime}-q)\cdot(r-1)}{2}\log
N.$ (S22)

Given that $\mathbf{Z}^{\prime}$ is the OMB with minimum number of variables,
and therefore, minimum number of parameters $q^{\prime}<q$. Thus, the term
$(q^{\prime}-q)$ becomes a negative constant. Also, the term $\frac{(r-1)}{2}$
is a positive constant. Consequently, Eq. (S22) goes to $-\infty$ in the limit
as $N\rightarrow\infty$, which implies that Eq. (S20) approaches to 0.

Case 3: $\mathbf{Z}$ is not a Markov blanket of $Y$.

The first term in Eq. (S21) is of $O(N)$ and dominates the second term, which
is $O(\log N)$. According to Lemma 9.4,
$H(Y|\mathbf{Z}^{\prime})<H(Y|\mathbf{Z})$; thus, the term
$H(Y|\mathbf{Z}^{\prime})-H(Y|\mathbf{Z})$ becomes a negative number. As a
result, Eq. (S21) becomes $-\infty$, which equivalently implies that Eq. (S20)
becomes 0. ∎