ON THE BOUNDARY BETWEEN QUALITATIVE AND QUANTITATIVE METHODS FOR
CAUSAL INFERENCE

arXiv:1711.04466v1 [math.ST] 13 Nov 2017

YUE WANG AND LINBO WANG
A BSTRACT. We consider how to quantify the causal effect from a random variable to a response variable. We
show that with multiple Markov boundaries, conditional mutual information (CMI) will produce 0, while causal
strength (CS) and part mutual information (PMI), which claim to behave better, are not well-defined, and have
other problems. The reason is that the quantitative causal inference with multiple Markov boundaries is an
ill-posed problem. We will give a criterion and some applicable algorithms to determine whether a distribution
has non-unique Markov boundaries.
KEY WORDS: causal inference, causal effect, Markov boundary

1. I NTRODUCTION
Inferring causal relationships is among the most important goals in artificial intelligence as humans tend
to think causally rather than probabilistically. A formal approach to represent causal relationships uses
causal directed acyclic graphs (DAG, Pearl, 2009), where random variables are represented as nodes and
causal relationships are represented as arrows. Besides qualitatively describing causal relationships via
causal DAGs, it is often desirable to obtain quantitative measures of the strength of arrows therein since they
provide more detailed information. There have been many measures proposed in the literature to quantify
the causal relationships between nodes in a causal DAG, such as the average treatment effect (ATE, Neyman,
1923; Rubin, 1974; Richardson and Robins, 2013), conditional mutual information (CMI, Dobrushin, 1963),
causal strength (CS, Janzing et al., 2013) and part mutual information (PMI, Zhao et al., 2016).
An interesting observation is that none of these measures is both identifiable and reasonable when the
causal system under consideration is degenerate. As a simple example, consider the confounder triangle
Z → X → Y with an edge Z → Y , where Z = X almost surely. For this example, the ATE of X on
Y is not identifiable, the conditional mutual information CMI(X, Y | Z) always equals zero even when
X influences Y strongly, while the causal strength and part mutual information for the arrow X → Y are
not well-defined. These problems should not be viewed as drawbacks of these existing measures, since the
causal effect of X on Y is not distinguishable from the causal effect of Z on Y . Instead, they reflect the
intrinsic difficulty of defining causal measures for this particular degenerate system.
In this paper, we generalize the observation above by providing a formal characterization of a degenerate
causal system, and argue that qualitative causal methods should be used in place of quantitative causal
methods when the causal system is degenerate. Our characterization of a degenerate causal system is based
on the notion of Markov boundary. The Markov boundary of a variable W in a variable set S is a minimal
subset of S, conditional on which all the remaining variables in S, excluding W , are rendered statistically
Date: November 15, 2017.
Yue Wang, Department of Applied Mathematics, University of Washington, Seattle, WA 98195, USA. Email address: yuewang@uw.edu.
Linbo Wang, Department of Biostatistics, Harvard University, Boston, MA 02115, USA.
1

2

YUE WANG AND LINBO WANG

independent of W (Statnikov et al., 2013). Under certain criteria for a “reasonable” quantitative measure of
causal effect, we show that if Y has multiple Markov boundaries in S, its possible parents, then in general
there is no measures for the causal effect of X on Y that is both reasonable and identifiable. In this case,
analysts may instead use qualitative causal methods and report all the causal explanations (i.e. Markov
boundaries) of the response variable.
We then propose practical approaches to determine the uniqueness of Markov boundary from data. Many
methods for Markov boundary discoveries exist in the literature, such as IAMB (Tsamardinos and Aliferis,
2003), KIAMB (Peña et al., 2007), Semi-Interleaved HITON-PC (Aliferis et al., 2010) and TIE* (Statnikov et al.,
2013). However, they often require strong assumptions or output all the Markov boundaries, which is not
necessary for our purpose. In this paper, we instead derive a necessary and sufficient condition for the
uniqueness of Markov boundary, based on which we develop algorithms that are computationally more
tractable. Moreover, the validity of our algorithms rely on fewer assumptions.
The rest of this article is organized as follows. In Section 2 we review backgrounds on causal effect
measures and Markov boundaries. We then propose several criteria for a reasonable quantitative causal effect
measure, and establish the equivalence between the uniqueness of Markov boundaries and the existence of a
reasonable quantitative causal effect measure in Section 3. In Section 4 we propose three algorithms to learn
the uniqueness of Markov boundary from observed data. The performance of these algorithms is evaluated
in Section 5. We finish with a discussion in Section 6.

2. BACKGROUND
2.1. Set-up. Consider a causal DAG G with vertices V and edges E. We say X is a parent of Y if the
path X → Y is present in G, and Y is a descendant of X is a path X → · · · → Y is present in G. A
variable is a descendant of itself, but not a parent of itself. For a variable W , we use Des(W ) to denote
the set consisting of all descendants of W , and Pa(W ) to denote the set consisting of all parents of W . We
assume the probability distribution P over variables V is Markov with respect to G in the sense that for
every W ∈ V , W is independent of V \ Des(W ) conditional on Pa(W ) (e.g. Spirtes et al., 2000).
Suppose X is a parent of Y and we are interested in the strength of the arrow X → Y . We assume that
we observe independent replications of V , and that Y has no descendants except for itself. In particular, we
exclude all variables that happen after Y from our analysis. If G is unknown, set L = V \ {X} so that L is
sufficient to control for confounding between X and Y . If G is known, set L = P a(Y ) \ {X}. We denote
the sample space of X, Y, L by X , Y, L, respectively.

2.2. Causal effect measures. We first review several causal effect measures commonly used in practice.
For simplicity, we only give the definitions in the discrete case. They can be easily extended to accommodate
continuous variables by replacing summation with integration.
Average treatment effect.
Under our setting, if X is binary, then the average treatment effect can be
defined as (e.g. Robins, 1986)
X XX
ATE(X → Y ) =
−(−1)x yf (y | x, l)f (l),
x∈X y∈Y l∈L

ON THE BOUNDARY BETWEEN QUALITATIVE AND QUANTITATIVE METHODS FOR CAUSAL INFERENCE

3

where f denotes the probability mass function. The average treatment effect describes the shift in the mean
of Y caused by switching X from 0 to 1, but fails to capture the influence of X on higher order moments of
the distribution of Y .
Conditional mutual information. The conditional mutual information (CMI) between X and Y conditional on L is defined as (Dobrushin, 1963)


X XX
f (x, y | l)
CMI(X, Y | L) =
f (x, y, l) log
.
f (x | l)f (y | l)
x∈X y∈Y l∈L

It can alternatively be expressed as
CMI(X, Y | L) = H(X, L) + H(Y, L) − H(L) − H(X, Y, L),

(1)

|=

P
where H(X) = − x∈X f (x) log f (x) is the (Shannon) entropy. It is easy to see from (1) that CMI is a
continuous function of the joint distribution of (X, Y, L). Furthermore, unlike the ATE, CMI captures all
conditional independences between variables in the sense that CMI(X, Y | L) = 0 if and only if X Y | L.
Causal strength. The causal strength of X on Y is defined as (Janzing et al., 2013)
X XX
f (y | x, l)
.
CS(X → Y ) =
f (x, y, l) log P
f (y | x′ , l)f (x′ )
x∈X y∈Y l∈L

Part mutual information.
as (Zhao et al., 2016)

x′ ∈X

The part mutual information between X and Y conditional on L is defined

PMI(X, Y | L) =

X XX

f (x, y, l) log

x∈X y∈Y l∈L

where f ∗ (x | l) =

P

y∈Y

f (x | y, l)f (y), f ∗ (y | l) =

PMI(X, Y | L) =

X

f (x, y, l) log P

X

f (x, y, l) log(

x,y,l

f (x, y | l)
∗
f (x | l)f ∗ (y | l)



,

f (y | x, l)f (x). One may also show that

f (x, y, l) log P

x,y,l

−

x∈X

X

x,y,l

+

P



f (y | x, l)
′
′
x′ f (y | x , l)f (x )
f (x | y, l)
′
′
y ′ f (x | y , l)f (y )

f (x, y | l)
).
f (x | l)f (y | l)

(2)

The first term on the right hand side of (2) equals CS(X → Y ), and the third term equals −CMI(X, Y | L).
Other measures. Other measures of causal effect between random variables include Pearson correlation
coefficient and partial correlation. For time series data, effect measures commonly used in practice include
Granger causality (Granger, 1969), conditional Granger causality (Geweke, 1984), partial Granger causality
(Guo et al., 2008), transfer entropy (Schreiber, 2000) and conditional transfer entropy (Kaiser and Schreiber,
2002; Frenzel and Pompe, 2007).
2.3. Markov blanket and Markov boundary.

|=

Definition 1. A Markov blanket M of a variable W within an observed variable set S (W ∈
/ S) is a subset
of S such that
(3)
W (S \ M ) | M .

4

YUE WANG AND LINBO WANG

Using the notion of mutual information, condition (3) can be written as CMI(W, S | M ) = 0, or
equivalently, MI(W, S) = MI(W, M ). In this sense, the Markov blanket M contains all the information of
S on W .
The following proposition is a direct application of the weak union property of probability distributions
(Pearl, 1988).
Proposition 1. Any superset of a Markov blanket is still a Markov blanket.
Definition 2. A Markov boundary is a minimal Markov blanket. In other words, it is a Markov blanket such
that no proper subsets of which is also a Markov blanket.
Even though Markov boundaries are minimal, they need not be unique. As illustrated in our confounder
triangle example, both {X} and {Z} are Markov boundaries of Y . As another example, consider the DAG
in Figure 1. Nodes X, Y, Z take value in {0, 1, 2}, while W takes value in {0, 1}. One may verify that there
exist two Markov boundaries of Y : {X, W } and {Z, W }. See also Statnikov et al. (2013, Example 2).
X
Z

Y

0

0

0

1

1

1

2

2

2

W

0
1

F IGURE 1. A DAG for which variable Y has multiple Markov boundaries. Combinations of values connected with lines have positive joint probabilities. For example,
P(Z = 1, X = 0, Y = 0, W = 0) > 0, while P(Z = 1, X = 2, Y = 1, W = 1) = 0.
Unlike the confounder triangle example, the two Markov boundaries of Y in Figure 1 do not coincide
almost surely. However, they are still variation dependent since P(X = 2) > 0, P(Z = 0) > 0, but
P(X = 2, Z = 0) = 0. This variation dependence is in fact an essential property of the multiplicity of
Markov boundaries.
Lemma 1. Let M denote all Markov boundaries of Y in S. If X is contained in at least one, but not all
of the Markov boundaries of Y so that X ∈ ∪ M \ ∩ M , and L = S \ {X}, then X and L are
M ∈M

M ∈M

variation dependent in that there exist x ∈ X , l ∈ L such that f (x) > 0, f (l) > 0, but f (x, l) = 0.
Remark 1. Let M ′ ∈ {M ∈ M : X ∈
/ M } be a Markov boundary that does not contain X. We point out
′
that X and M may be variation independent. For example, suppose X, Z, W are independent Bernoulli
variables with mean 0.5, Y = X + Z mod 2 and U = W + Y mod 2, where a mod b denotes the
remainder of a divided by b. In this case Y has two Markov boundaries: {X, Z} and {W, U }, and X is
variation independent of {W, U }.
We now review some conditions under which the Markov boundary is unique.

ON THE BOUNDARY BETWEEN QUALITATIVE AND QUANTITATIVE METHODS FOR CAUSAL INFERENCE

5

Definition 3 (Faithfulness). Let I denote the collection of conditional independence relationships shared
by all probability distributions that are Markov with respect to G. A probability distribution is faithful to G
if and only if its conditional independence relationships are characterized by I (Peña et al., 2007).

|=

|=

|=

Definition 4 (Intersection). A probability distribution on V satisfies the intersection property if and only
if for any four subsets of variables P , Q, Z, W such that P Z | (Q, W ), P W | (Q, Z), we have
P (Z, W ) | Q.
Definition 5 (Strict positiveness). A probability distribution P on V is called strictly positive if and only if
for any two disjoint subsets of variables X and Y such that P(X = x) > 0, P(Y = y) > 0, we have
P(X = x, Y = y) > 0.
Proposition 2. If P (i) is faithful to G, or (ii) has intersection property, or (iii) is strictly positive, then any
variable Z ∈ V has a unique Markov boundary in V \ {Z}.

3. Q UANTIFYING

|=

|=

We refer readers to Pearl and Paz (1985) and Pearl (1988) for a proof of Proposition 2. None of the
three conditions in Proposition 2 is necessary. For example, suppose X, Y , Z, W ∈ {0, 1}, P(X = Z =
Y = W = 0) = 0.5, P(X = Z = 1, Y = W = 0) = 0.25, P(X = Z = Y = W = 1) = 0.25
and S = {X, Y, Z, W }. Since Y X | Z, Y Z | X, Y 6⊥⊥ X, Z, the joint distribution of (X, Z, Y, W )
does not have intersection property and is hence not faithful. Since P(X = 0) > 0, P(Z = 1) > 0,
P(X = 0, Z = 1) = 0, this distribution is not strictly positive. However, each variable in this example has
a unique Markov boundary within the other three variables.
There are many applications for which the Markov boundary of the response variable may not be unique.
For instance, in breast cancer studies, it may be the case that several gene sets have nearly the same effect
for survival prediction (Ein-Dor et al., 2004), such that each of the gene sets is a Markov boundary of the
survival indicator. We refer interested readers to Statnikov et al. (2013) for more real-world examples on
multiplicity of Markov boundaries.
CAUSAL EFFECTS WITH MULTIPLE

M ARKOV

BOUNDARIES

3.1. Problems of known causal effect measures. In this section, we show that the commonly used causal
effect measures may not be reasonable or identifiable when the response variable Y has multiple Markov
boundaries. On the one hand, for X ∈ ∪ M \ ∩ M , we have CMI(X, Y | L) = 0 even when X
M ∈M

M ∈M

influences Y strongly. This is not desirable, as pointed out by Janzing et al. (2013); Zhao et al. (2016). On
the other hand, Lemma 1 shows that X and L are variation dependent. It follows that f (y | x, l) may not
always be defined even when f (x) > 0 and f (l) > 0. It follows that ATE(X → Y ), CS(X → Y ) and
PMI(X, Y | L) are not defined.
A natural solution to the latter problem is to assign a value in these degenerate scenarios. However, Theorem 1 shows that the resulting measures are not continuous functions of the joint distribution of (X, Y, L). In
the following, for any given probability distribution P′ , we can calculate the corresponding causal quantities,
denoted as CS[P′ ](X → Y ), PMI[P′ ](X, Y | L) and so on.
Theorem 1. If X ∈

∪

M ∈M

M\

denoted as {P1 , P2 , · · · } and

∩

M , then there exist two sequences of distributions on (X, Y, L),

M ∈M
{P′1 , P′2 , · · · },

both of which converge to P under the total variation distance,

6

YUE WANG AND LINBO WANG

but limi→∞ CS[Pi ](X → Y ) 6= limi→∞ CS[P′i ](X → Y ). The same results apply to ATE(X → Y ) and
PMI(X, Y | L).
Theorem 1 may be proved using Lemma 1 and the following Lemma 2.
Lemma 2. Assume there exist x ∈ X , l ∈ L such that f (x) > 0, f (l) > 0, but f (x, l) = 0. Then there
exists an interval (g1 , g2 ) with positive length, such that for any g ∈ (g1 , g2 ), any δ > 0, there exists a
probability distribution P′ with d(P, P′ ) < δ, such that CS[P′ ](X → Y ) = g. The same results apply to
ATE(X → Y ) and PMI(X, Y | L).
Remark 2. Lemma 2 is similar in flavor to the Picard’s great theorem: If an analytic function f has an
essential singularity at a point w, then on any punctured neighborhood of w, f (z) takes on all possible
complex values, with at most a single exception. In this sense, CS and PMI are “essentially singular” at the
probability distribution that implies multiple Markov boundaries for Y .
To get intuition into Theorem 1, consider the following example: X, Z ∈ {1, 2, · · · , n}, Y ∈ {0, 1, 2, 3}.
The causal DAG is X → Y ← Z with an edge Z → X. If x = z, then P(X = x, Y = 0, Z = z) = η[1/n−
nǫ − (n − 1)δ], P(X = x, Y = 1, Z = z) = ηǫ, P(X = x, Y = 2, Z = z) = (1 − η)[1/n − nǫ − (n − 1)δ],
P(X = x, Y = 3, Z = z) = (1 − η)ǫ. If x 6= z, then P(X = x, Y = 2, Z = z) = (1 − η)δ,
P(X = x, Y = 3, Z = z) = (1 − η)ǫ.
Letting ǫ → 0, δ → 0, then
n
n
CS(X → Y ) → log
, PMI(X, Y | Z) → 2 log
.
δ
δ
1 + (n − 1) δ+ǫ
1 + (n − 1) δ+ǫ
By changing the value of δ/ǫ, PMI(X, Y | Z)/2 and CS(X → Y ) can take any value between 0 and
log n.
When δ = ǫ = 0, 0 < η < 1, Y has two Markov boundaries: {X} and {Z}. When δ, ǫ, η tend to 0,
H(Y ) → 0, but CS and PMI can still take large values, which is problematic.
We can see that near a distribution with multiple Markov boundaries, CS and PMI will change drastically.
Therefore calculating CS or PMI of a distribution which is close to another distribution with multiple Markov
boundaries might not be numerically feasible.
3.2. Criteria for reasonable causal effect measures. We propose the following criteria for reasonable
definitions of the strength of arrows in a causal graph.
C0. The strength of X → Y is identifiable from the joint distribution of Y and its possible parents S.
C1. If there is unique Markov boundary M of Y within S, and X ∈
/ M , then the strength of X → Y is
0.
C2. If there is unique Markov boundary M of Y within S, and X ∈ M , then the strength of X → Y is
at least CMI(X, Y | M \ {X}).
C3. The strength of X → Y is a continuous function of the joint distribution of Y and S, under total
variance distance.
Here we use “possible”, because we might have different levels of knowledge of the causal DAG. If we
know nothing about the causal DAG, then all the measured variables are possible parents. If we know the
full DAG, then we can abandon those variables which do not directly cause Y .

ON THE BOUNDARY BETWEEN QUALITATIVE AND QUANTITATIVE METHODS FOR CAUSAL INFERENCE

7

C0 is similar to the “Locality” postulate P2 in Janzing et al. (2013). If Z is not a parent of Y , then it
impacts Y through some parents of Y . In other words, the knowledge of all possible parents should suffice
to determine the strength of X → Y . C1 is reasonable since Markov boundary contains all the information
on Y from S. It is natural to say that X has no causal effect on Y if X ∈
/ M . On the other hand, if
X is contained in the unique Markov boundary M , then it contains non-trivial information of Y since
MI(Y, S \ {X}) < MI(Y, S). In this case, it is reasonable to assign a positive value to X → Y . In C2, the
concrete form of the lower bound is not important. We only require that it is positive, and relies only on the
joint distribution of M and Y . Finally, if C3 fails, then a small observational error might lead to significant
error in the estimated strength of arrow X → Y .
3.3. An impossibility result. In this section, we show that with multiple Markov boundaries, reasonable
causal effect measures are impossible, since the above criteria are incompatible.
Theorem 2. Assume there is no extra knowledge of the causal DAG. Assume X is contained in some Markov
boundary, but not all the Markov boundaries. Then we cannot define the strength of X → Y satisfying
criteria C0, C1, C2 and C3.
To provide intuition, we first show Theorem 2 for a simple distribution. Consider X, Y, Z which only
take 0, 1, 2. P(X = 0, Y = 0, Z = 0) = 1/2 − ǫ, P(X = 1, Y = 2, Z = 1) = 1/4, P(X = 1, Y =
0, Z = 1) = 1/4, P(X = 2, Y = 2, Z = 0) = ǫ, where ǫ is a small positive number. Here {X}
is the unique Markov boundary of Y , so from criterion C2, the strength of X → Y should be at least
CMI(X, Y | ∅) = MI(X, Y ) ≈ 0.216. Now consider a perturbation of the above example: P(X =
0, Y = 0, Z = 0) = 1/2 − ǫ, P(X = 1, Y = 2, Z = 1) = 1/4, P(X = 1, Y = 0, Z = 1) = 1/4,
P(X = 0, Y = 2, Z = 2) = ǫ. Here we just switch X and Z in the joint distribution. Now {Z} is the
unique Markov boundary of Y , so from criterion C1, the strength of X → Y should be 0. We can see that
with an arbitrarily small perturbation, the strength of X → Y has to change drastically. Thus we cannot give
a satisfactory value when ǫ = 0, while keeping the strength continuous. Therefore criterion C3 is violated.
Near ǫ = 0, such quantity cannot be uniformly continuous, so no matter how we increase the accuracy of
measuring joint distribution, we might still have significant errors.
The key idea above is to show that under some arbitrarily small perturbation, there is only one Markov
boundary M which contains X. Let the magnitude of perturbation tend to 0, then the strength of X → Y
in the original distribution should be at least CMI(X, Y | M \ {X}). However, we can also add arbitrarily
small perturbation such that there is only one Markov boundary, and it does not contain X. When we let
the magnitude of perturbation tend to 0, we can see that the strength of X → Y in the original distribution
should be 0. Such dilemma finishes the proof. First we need the description of perturbations and some
lemmas.
Definition 6. Given a variables Zi . Let ξi be an independent Bernoulli variable with mean ǫ. If ξi = 0, then
Ziǫ = Zi . If ξi = 1, then Ziǫ = Ui , where Ui is a random variable that is independent of Zi and ξi . When we
add ǫ-noise on multiple variables Zi , Zj , · · · , those variables ξi , ξj , · · · , Ui , Uj , · · · are independent, and
they are independent with X, Y, Z1 , · · · , Zk .
We can see that substituting X by X ǫ does not affect the joint distribution of Y, Z1 , · · · , Zk . If we add
ǫ-noise on m variables, then the total variance distance of joint distributions before and after substituting

8

YUE WANG AND LINBO WANG

is at most mǫ. Therefore, when we let ǫ → 0, the perturbed joint distribution will converge to the original
distribution.
Lemma 3 (Strict Data Processing Inequality). S1 is a group of variables without X, Y . If we add ǫnoise on X to get X ǫ , then CMI(X ǫ , Y | S1 ) ≤ CMI(X, Y | S1 ), and the equality holds if and only if
CMI(X, Y | S1 ) = 0.
This is a special case of “data processing inequality” in information theory (Cover and Thomas, 2012).
The intuition is that transmitting data through a noisy channel cannot increase information, namely “garbage
in, garbage out”. However, the original data processing inequality only indicates that the sufficient and
necessary condition under which the equality holds is CMI(X, Y | X ǫ , S1 ) = 0, which is rather abstract
and difficult to verify, since it depends on the form of ǫ noise. Here we propose a full proof, since we need
to show that adding non-zero ǫ-noise will always decrease (instead of just “not increase”) the information of
Y , unless there is no information to lose. In other words, we propose the concrete form of the sufficient and
necessary condition under which the equality holds. That is why the proof is much lengthier than that of the
original data processing inequality.
Expand CMI(X ǫ , Y | S1 ) ≤ CMI(X, Y | S1 ), then we have MI({X ǫ } ∪ S1 , Y ) ≤ MI({X} ∪ S1 , Y ).
Especially, the total information MI(S, Y ) will not increase if we add ǫ-noise on possible parents of Y .
Lemma 4. Assume Y has multiple Markov boundaries within S = {X, Z1 , · · · , Zk }. Choose one Markov
boundary: M0 . If we add ǫ noise on each variable in S \M0 , then in the new distribution, M0 is the unique
Markov boundary.
Proof of Theorem 2. Lemma 4 shows that if X ∈ M0 , X ∈
/ M1 , then we can add ǫ-noise on S \ M0 , such
that M0 is the only Markov boundary. Therefore, in the perturbed distribution, the strength of X → Y is at
least CMI(X, Y | M0 \ {X}), required by criterion C2. Let ǫ → 0, the joint distribution converges to the
original joint distribution, and the continuous criterion C3 shows that the strength of X → Y in the original
distribution is also no less than CMI(X, Y | M0 \ {X}). Similarly, we can add ǫ-noise on S \ M1 , such
that M1 is the only Markov boundary. Therefore, in the perturbed distribution, the strength of X → Y is
0, required by criterion C1. When we let ǫ → 0, criterion C3 requires that the strength of X → Y in the
original distribution is also 0. This is the contradiction we need to finish the proof of Theorem 2.

In Theorem 2, we assume there is no extra knowledge of the causal DAG. If we know partial or full causal
DAG, then we know some variables do not cause the response variable Y directly. After discarding these
variables, we have similar results.
Theorem 3. Assume we have partial or full knowledge of the causal DAG. Discard all variables that have
been known not causing Y directly. For the remaining variables, namely all the possible parents of Y ,
assume X is contained in some Markov boundary, but not all Markov boundaries. Then we cannot define
the strength of X → Y satisfying criteria C0, C1, C2 and C3.
The proof is the same with Theorem 2. The reason of discarding variables that do not cause Y directly
is to guarantee the joint distribution with ǫ-noise is consistent with the causal DAG. Also it is required by
criterion C0. Consider an example with DAG Z → X → Y , where X = Z. Here Z does not cause Y
directly. If we do not discard Z, then there are two Markov boundaries for Y : {X} and {Z}. However, in

ON THE BOUNDARY BETWEEN QUALITATIVE AND QUANTITATIVE METHODS FOR CAUSAL INFERENCE

9

|=

this case the strength of X → Y is well-defined. If we try to utilize the proof of Theorem 2, in one step we
need to perturb the joint distribution by adding ǫ-noise on X. Then we have CMI(Z, Y | X ǫ ) > 0, but the
causal DAG requires Z Y | X ǫ . Therefore adding ǫ-noise on X is prohibited, and the proof is invalid.
In sum, when there are multiple Markov boundaries, quantifying the causal strength is an ill-posed problem. We can only qualitatively describe the causal relationship. That is why the uniqueness of Markov
boundary is important.
4. A LGORITHMS

FOR TESTING THE UNIQUENESS OF

M ARKOV

BOUNDARY

4.1. An assumption-free algorithm. We have seen that when Markov boundary is not unique, causal
quantities are problematic. The next question is to determine the condition under which Markov boundary
is unique.
Our first algorithm is based on the following necessary and sufficient condition for the uniqueness of
Markov boundary.
Definition 7 (Essential variable). A variable W ∈ S is called an essential variable with respect to Y , if
CMI(Y, S | S \ {W }) > 0. It means that W itself contains some irreplaceable information of Y . Denote
the set of all essential variables by E.
Lemma 5. E is the intersection of all Markov boundaries of Y within S.
Proof. Assume there exists a Markov boundary M such that W ∈ E, W ∈
/ M . Since CMI(Y, S | M ) =
CMI(Y, S | S) = 0, we also have CMI(Y, S | S \ {W }) = 0 (Proposition 1), which is a contradiction.
If W ∈
/ E, then CMI(Y, S | S \ {W }) = 0, S \ {W } is a Markov blanket, which means it contains a
Markov boundary. This Markov boundary does not contain W .

Theorem 4. The Markov boundary of Y is unique if and only if CMI(Y, S | E) = 0. In this case, the unique
Markov boundary is E.
Proof. If Markov boundary is unique, then E is just the Markov boundary, therefore CMI(Y, S | E) = 0.
If CMI(Y, S | E) = 0, then E is a Markov blanket, which means it should contain a Markov boundary.
But E should be contained in every Markov boundary, therefore E itself is a Markov boundary. E as a
Markov boundary cannot be a proper subset of another Markov boundary, thus the only Markov boundary
is E.

|=

Theorem 4 is the foundation of Algorithm 1. We construct E, then check whether Y S | E.
To perform the condition independence tests in Algorithm 1, one may use G-test (Neapolitan, 2004) or
Pearson’s χ2 -test if variables in S are discrete. For continuous variables, we can use KCI-test (Zhang et al.,
2011), maximal nonlinear conditional correlation (Huang, 2010) or a weighted Hellinger distance based test
(Su and White, 2008).
The conditional independence relationships (4) concerns all observed variables, so that Algorithm 1 may
be highly data inefficient (Peña et al., 2007) in that to achieve the same accuracy, it needs much more observations, compared with other algorithms in the next section. We now consider algorithms with better data
efficiency.
4.2. Data efficient algorithms. Alternatives for detecting the uniqueness of Markov boundaries usually
depend on an algorithm that generates one Markov boundary. Therefore we first consider such algorithms.

10

YUE WANG AND LINBO WANG

Algorithm 1 An assumption-free algorithm for determining the uniqueness of Markov boundary

(4)

|=

|=

(1) Input
Observations of S = {X1 , . . . , Xk } and Y
(2) Set E = ∅
(3) For i = 1, · · · , k,
Test
Xi Y | S \ {Xi }
If Xi 6⊥⊥ Y | S \ {Xi }
E = E ∪ {Xi }
(4) If Y S | E output: Y has a unique Markov boundary
else
output: Y has multiple Markov boundaries

4.2.1. Algorithm for producing a single Markov boundary. There are many algorithms to produce Markov
boundary, such as IAMB (Tsamardinos and Aliferis, 2003), KIAMB (Peña et al., 2007), Semi-Interleaved
HITON-PC (Aliferis et al., 2010). There is a summary of different algorithms in Statnikov et al. (2013).
Remark 3. There are also other approaches for finding Markov boundaries, such as MBOR (de Morais and Aussem,
2010), BLCD (Mani and Cooper, 2004), PCMB (Peña et al., 2007) and GLL-PC (Aliferis et al., 2010).
However, these algorithms require faithfulness of the observed data distribution so that they cannot be
used to detect uniqueness of Markov boundary.
IAMB and KIAMB can produce one Markov boundary if the distribution has composition property.
KIAMB has a parameter K ∈ [0, 1]. When K = 1, it coincides with IAMB. When K < 1, it is stochastic.
When K = 0, every Markov boundary has a positive probability to be produced in one execution. When
K < 1, if there are multiple Markov boundaries, then at least two of them have positive probabilities to be
produced in one execution (Peña et al., 2007).
Semi-Interleaved HITON-PC requires that there exists a Markov boundary whose members are all dependent and conditional dependent with the response variable, unless some variables have exact the same
information (Aliferis et al., 2010; Statnikov et al., 2013). Under this condition, multiple Markov boundaries
is also possible.
According to Peña et al. (2007), KIAMB outperforms IAMB considerably often. The complexity of
Semi-Interleaved HITON-PC grows exponentially with the size of the Markov boundary, therefore it is much
slower than KIAMB (Statnikov et al., 2013). Thus we utilize KIAMB in the implementation of Algorithm
3 below, where a general Markov boundary producing algorithm is needed.
The validity of all these algorithms relies on extra assumptions. Instead, we propose an assumption-free
algorithm to produce one Markov boundary.
Here the “Association” is the negative p-value of the conditional independence test.
Proof of correctness of Algorithm 2. Obviously the output M is a Markov blanket. In the last step, we
have checked that X0 6⊥⊥ Y | M \ {X0 }. For Xi ∈ M , since Association(Xi , Y | M \ {Xi }) ≥
Association(X0 , Y | M \ {X0 }), we also have Xi 6⊥⊥ Y | M \ {Xi }. Therefore the output of Algorithm 2
is a Markov boundary.

In Algorithm 2, there are some conditional independence tests which concern all or almost all variables,
which is data inefficient. This is the price we pay for requiring no condition. As a result, these tests may
have lower power detecting conditional dependencies. Therefore, we need to avoid executing such tests on

ON THE BOUNDARY BETWEEN QUALITATIVE AND QUANTITATIVE METHODS FOR CAUSAL INFERENCE

11

Algorithm 2 An assumption-free algorithm for producing one Markov boundary

|=

(1) Input
Observations of S = {X1 , . . . , Xk } and Y
(2) Set M0 = S
(3) Repeat
Set X0 = arg minX∈M0 Association(X, Y | M0 \ {Xi })
If X0 Y | M0 \ {X0 } Set M0 = M0 \ {X0 }
Until X0 6⊥⊥ Y | M0 \ {X0 }
(4) Output M0 is a Markov boundary
almost all variables if the true result is not independent. Our solution is to feed these tests the variables with
the least associations, namely with the highest possibilities to be independent, so as to reduce the error rate.
4.2.2. Algorithm for testing the uniqueness of Markov boundary. TIE* is an algorithm to produce all
Markov boundaries (Statnikov et al., 2013), which needs to utilize another algorithm which can correctly
produce one Markov boundary (for example KIAMB). To test the uniqueness of Markov boundary, we can
execute TIE* until it produces two Markov boundaries, or it terminates with the unique Markov boundary.
This is Algorithm 3.
Algorithm 3 A general algorithm for determining the uniqueness of Markov boundary

|=

(1) Input
Observations of S = {X1 , . . . , Xk } and Y
Algorithm A which can correctly produce one
Markov boundary
(2) Set M0 = {X1 , X2 , · · · , Xm } to be the result of Algorithm A on S
(3) For i = 1, · · · , m,
Set Mi to be the result of Algorithm A on S \ {Xi }
Output Y has multiple Markov boundaries
If Y M0 | Mi
Terminate
(4) Output Y has a unique Markov boundary

Proof of correctness of Algorithm 3. Markov boundary is not unique if and only if there exists variable Xi ∈
M0 which is not essential, namely MI(Y, S \{Xi }) = MI(Y, S). Since MI(Y, S \{Xi }) = MI(Y, Mi ) and
MI(Y, M0 ) = MI(Y, S), it is equivalent to MI(Y, Mi ) = MI(Y, M0 ), namely CMI(Y, M0 | Mi ) = 0. 

|=

Specifically, in Algorithm 3, if we set A to be Algorithm 2, then we acquire an assumption free algorithm
to test the uniqueness of Markov boundnary, which is denoted as Algorithm 3-AF (AF means assumption
free). Algorithm 3 with A set as KIAMB is denoted as Algorithm 3-KI.
In Algorithm 3, we first produce one Markov boundary M0 = {X1 , X2 , · · · , Xm }. Then for each
Xi ∈ M0 , we produce a Markov boundary within S \ {Xi }, Mi , and check if Mi , which is supposed to
contain all the information in S \ {Xi }, has the same information with S. In fact, for Xi ∈ M0 , we can
directly test whether Y Xi | S \ {Xi }. This is Algorithm 4.
Proof of correctness of Algorithm 4. For a Markov boundary M0 , it is the unique Markov boundary if and
only if it coincides with E. Therefore, we only need to check whether there exists a variable Xi ∈ M0
which is not essential, namely CMI(Y, Xi | S, M \ {Xi }) = 0.


12

YUE WANG AND LINBO WANG

Algorithm 4 An assumption-free algorithm for determining the uniqueness of Markov boundary

|=

(1) Input
Observations of S = {X1 , . . . , Xk } and Y
(2) Set M0 = {X1 , X2 , · · · , Xm } to be the result of Algorithm 2 on S
(3) For i = 1, · · · , m,
If Y Xi | S \ {Xi }
Output Y has multiple Markov boundaries
Terminate
(4) Output Y has a unique Markov boundary

5. S IMULATION

RESULTS AND ANALYSIS

In this section, we compare different algorithms for testing the uniqueness of Markov boundary. We
will first implement these algorithms on some randomly generated data, and then conduct some theoretical
analyses.

|=

|=

5.1. Simulation setup. We implement Algorithm 1, Algorithm 3 with A set as KIAMB (denoted as Algorithm 3-KI), Algorithm 3 with A set as Algorithm 2 (denoted as Algorithm 3-AF) and Algorithm 4. We use
R 3.3.2 and an R package “infotheo” (Meyer, 2008). The significance level α is set to 0.05. In KIAMB,
the parameter K is set as 0.8, which is adopted in Peña et al. (2007); Statnikov et al. (2013). We generate
various random data with different mechanisms, and count the frequency of success for each algorithm.
Since all algorithms are theoretically correct if the needed assumptions are satisfied, when the number of
observations is large enough, these algorithms can almost always produce correct results. To evaluate these
algorithms, we need to consider dataset with limited number of observations. We will also consider a case
where the needed assumption of KIAMB is violated. In each case, we run each algorithm 1000 times for
each number of observations.
In all simulations, we consider 1 response variable Y and 10 possible parents of Y , S = {X1 , · · · , X10 }.
All variables are Bernoulli with mean 0.5.
We consider three cases with different mechanisms. In each case, we execute all algorithms with multiple
numbers of observations: 300, 500, 800, 1000, 2000, 2400, 4000, 5000, 6000, 10000, 14000, 17000, 20000,
23000, 30000, 40000, 50000, 60000.
Case A, X1 , X2 , · · · , X10 are independent. Y = X1 with probability 0.8, Y = X2 with probability 0.1,
Y = X3 with probability 0.1. There is unique Markov boundary: {X1 , X2 , X3 }. The composition property
is satisfied.
Case B, X1 , X2 , X4 , X5 , · · · , X10 are independent. X3 = X4 . Y = X1 with probability 0.8, Y = X4
with probability 0.1, Y = X5 with probability 0.1. There are multiple Markov boundaries: {X1 , X3 /X4 , X5 }.
The composition property is satisfied.
Case C, define Z = X1 + X2 mod 2. X1 , X2 , · · · , X8 are independent. Y = Z with probability
0.8, Y = X3 with probability 0.1, Y = X4 with probability 0.1. X9 = Z with probability 0.95, X9 =
1 − Z with probability 0.05. X10 = X9 . There is unique Markov boundary {X1 , X2 , X3 , X4 }. We have
X1 Y | (X3 , · · · , X10 ), X2 Y | (X3 , · · · , X10 ), but (X1 , X2 ) 6⊥⊥ Y | (X3 , · · · , X10 ), therefore the
composition property, which is required by KIAMB, is violated. Theoretically, in this case, Algorithm 3KI cannot recognize the information of X1 and X2 , and will claim there are multiple Markov boundaries:

ON THE BOUNDARY BETWEEN QUALITATIVE AND QUANTITATIVE METHODS FOR CAUSAL INFERENCE

13

{X3 , X4 , X9 /X10 }. Even if we increase the number of observations, Algorithm 3-KI still cannot provide
correct results.
5.2. Simulation results and analysis. See Figures 2,3,4 for success rates and average time costs per execution of different algorithms in Cases A, B, C. Detailed results can be found in Table 1, 2, 3 in Appendix.
2

1

10
Alg 1
Alg 3-KI
Alg 3-AF
Alg 4

0.8

10

1

10

0

0.6

Time

Success rate

Success rates in Case A

0.4
0.2

Time costs in Case A
Alg 1
Alg 3-KI
Alg 3-AF
Alg 4

10

-1

10

-2

0
10

2

10

3

10

4

10

5

10

2

Number of observations

10

3

10

4

10

5

Number of observations

F IGURE 2. Simulation results of Case A

2

1

10
Alg 1
Alg 3-KI
Alg 3-AF
Alg 4

0.8

10

1

0.6

Time

Success rate

Success rates in Case B

0.4

Time costs in Case B
Alg 1
Alg 3-KI
Alg 3-AF
Alg 4

100

10-1

0.2
0
102

103

104

Number of observations

105

10-2 2
10

103

104

105

Number of observations

F IGURE 3. Simulation results of Case B
In all situations, The time costs are ascending for Algorithm 1, 4, 3-KI, 3-AF. Time cost increases with
the number of observations. Compared to Algorithm 1 and 4, Algorithm 3-AF and 3-KI calculates much
more CMI, which significantly increases its time cost. Besides, when we execute multiple tests, we might
need to make false discover rate adjustment (Benjamini and Yekutieli, 2001), but in most algorithms, the
number of tests cannot be predicted, which might be problematic.
For success rate, Algorithm 3-AF is always satisfactory. Algorithm 3-KI performs slightly better than
Algorithm 3-AF if the composition property is satisfied (Case A, Case B), but performs disastrously if the
composition property is violated (Case C). In Case C, Algorithm 3-KI still has a very small probability to
produce correct answer, since it is stochastic. Algorithm 1 tends to claim that Markov boundary is unique,

14

YUE WANG AND LINBO WANG

2

1

10
Alg 1
Alg 3-KI
Alg 3-AF
Alg 4

0.8

10

1

0.6

Time

Success rate

Success rates in Case C

0.4
0.2

Time costs in Case C
Alg 1
Alg 3-KI
Alg 3-AF
Alg 4

100

10

-1

0
102

103

104

105

Number of observations

10-2 2
10

103

104

105

Number of observations

F IGURE 4. Simulation results of Case C

|=

|=

while Algorithm 4 tends to claim that Markov boundary is not unique, regardless of whether the Markov
boundary is unique or not.
Notice that in the last step of Algorithm 1, Y S | E means Markov boundary is unique, while in
the last step of Algorithm 4, Y Xi | S \ {Xi } means Markov boundary is not unique. For such tests
which concern all variables, when the number of observations is not large enough, they cannot reject the
conditional independence hypothesis, even if variables are dependent. Therefore, Algorithm 1 tends to
claim Markov boundary is unique (false negative), Algorithm 4 tends to claim Markov boundary is not
unique (false positive), whether Markov boundary is really unique or not. In Algorithm 3-AF, generally the
last step test does not concern many variables. Although in the first few steps some tests concern almost all
variables, these tests do not directly determine the output, and we have argued that our implementation can
significantly lower the error rate. In Algorithm 3-KI, generally there is no test with many variables. This
explains why Algorithm 3-KI is slightly better than Algorithm 3-AF under composition property, and why
on average Algorithm 3-AF and 3-KI are better than Algorithm 1 and 4.
Wrong conditional independence tests are necessary but not sufficient for wrong final result (except Algorithm 3-KI in Case C). In some cases, we can see that when we increase the number of observations, the
success rate of an algorithm might first decrease and then increase. The reason is that, the success rate of
each test is positively related to the number of observations, but the success rate of the whole algorithm has
complicated relationship with wrong tests. Initially this algorithm has several tests with wrong results, but
these wrong tests happen to make the correct output. When we increase the number of observations, some
of these tests become correct, but the remaining wrong tests lead to wrong output. Finally when the number
of observations is large enough, all tests are correct, and the output is correct.
On average, Algorithm 3-KI performs the best when the composition property is satisfied. However, if we
set A to be any algorithm which requires additional assumption, such as KIAMB, we need to guarantee that
theoretically algorithm A produces a correct Markov boundary in each execution, otherwise Algorithm 3
might produce wrong results, just like Case C. To make sure algorithm A can produce correct result, we can
prove the distribution satisfies the condition required by algorithm A. However, generally such verification
process is slower and less efficient than Algorithm 1, 4 or 3-AF. If we do not verify the requirement of
algorithm A, then we need to verify the result of A is a Markov boundary after each execution. If the

ON THE BOUNDARY BETWEEN QUALITATIVE AND QUANTITATIVE METHODS FOR CAUSAL INFERENCE

15

Markov boundary produced by A is fake, then we need to apply Algorithm 2 to produce a correct Markov
boundary, which is cumbersome.
In sum, if we know that the condition required by A is satisfied with no or very small cost, then utilizing A in Algorithm 3 is better. Otherwise, Algorithm 3-AF is better. Also, we can directly utilize A in
Algorithm 3, to obtain a less reliable result. Since multiple Markov boundaries are problematic, we prefer
false positive (regard unique as non-unique) rather than false negative (regard non-unique as unique). In this
sense, Algorithm 4 has the lowest false negative rate. If time cost is much more concerned than accuracy,
then Algorithm 1 is preferred.
As expected, all these algorithms are faster than executing full TIE* algorithm.
6. D ISCUSSION
We show that when there are multiple Markov boundaries, we cannot concordantly define the causal
effect from a random variable to a response variable. Especially, CS and PMI are ill-defined with multiple
Markov boundaries. Thus we propose a sufficient and necessary condition under which Markov boundary
is unique. We also propose several applicable algorithms to check whether Markov boundary is unique.
ACKNOWLEDGEMENTS
The authors would like to thank Siqi He, Yifei Liu, Daniel Malinsky, Hong Qian, Thomas Richardson,
Weili Wang, Daxin Xu, Qingyuan Zhao, Xiao-Hua Zhou for helpful advice and discussions.
A PPENDIX
A. Proof of Lemma 1.

|=

Proof. Consider two Markov boundaries M1 , M2 within {X} ∪ L. Let M1 = {X} ∪ Z 1 , X ∈
/ M2 ,
2
3
1
2
′
′ },
M2 \ M1 = Z , L ∪ {X} \ (M1 ∪ M2 ) = Z , where Z = {Z1 , · · · , Zn }, Z = {Z1 , · · · , Zm
Z 3 = {Z1′′ , · · · , Zl′′ }. Therefore L = Z 1 ∪ Z 2 ∪ Z 3 .
Fix z01 ∈ Z 1 such that f (z01 ) > 0. Assume that for xi ∈ X , f (xi , z01 ) > 0 is true for i ∈ {1, 2, · · · , p}.
Assume that for zj2 ∈ Z 2 , f (z01 , zj2 ) > 0 is true for j ∈ {1, 2, · · · , q}. Consider any y ∈ Y.
To obtain contradiction, we assume that f (xi , z01 , zj2 ) > 0, ∀i ∈ {1, 2, · · · , p}, ∀j ∈ {1, 2, · · · , q}.
Since X Y | (Z 1 , Z 2 ) (Proposition 1), ∀i, r ∈ {1, 2, · · · , p}, ∀j ∈ {1, 2, · · · , q},
f (y | xi , z01 , zj2 ) = f (y | xr , z01 , zj2 ).
Y | (X, Z 1 ), ∀r ∈ {1, 2, · · · , p}, ∀j, s ∈ {1, 2, · · · , q},

|=

Since Z 2

f (y | xr , z01 , zj2 ) = f (y | xr , z01 , zs2 ).
All the conditions have positive probabilities, so the conditional probabilities are well-defined.
Then we have
f (y | xi , z01 , zj2 ) = f (y | xr , z01 , zs2 ),
∀i, r ∈ {1, 2, · · · , p}, ∀j, s ∈ {1, 2, · · · , q}.
Since this is true for any possible values of X and Z 2 when Z 1 = z01 , we know that
f (y | xi , z01 , zj2 ) = f (y | z01 ).

16

YUE WANG AND LINBO WANG

Therefore, ∀z11 ∈ Z 1 with f (z11 ) > 0, ∀y ∈ Y, ∀i, j,
f (xi , zj2 , y | z11 ) = f (xi , zj2 | z11 )f (y | z11 )
|=

|=

is valid.
This implies that (X, Z 2 ) Y | Z 1 , therefore X Y | Z 1 , MI(Y, Z 1 ) = MI(Y, (X, Z 1 )). Since
M1 = {X} ∪ Z 1 , MI(Y, (X, Z 1 )) = MI(Y, L). Thus MI(Y, Z 1 ) = MI(Y, {X} ∪ L), implying that
Z 1 is a Markov blanket, which is a contradiction. So there exists x ∈ X , z01 ∈ Z 1 , z12 ∈ Z 2 such that
f (x, z01 ) > 0 (implies f (x) > 0), f (z01 , z12 ) > 0, but f (x, z01 , z12 ) = 0. Choose z13 ∈ Z 3 such that
f (z01 , z12 , z13 ), and let l = (z01 , z12 , z13 ), then f (x) > 0, f (l) > 0, but f (x, l) = 0.


B. Proof of Lemma 2.
Proof. In the following we will assume there is only one pair of (x, l) such that f (x) > 0, f (l) > 0,
f (x, l) = 0. If there are multiple pairs, we can treat them one by one.
We construct a family of probability distributions Pηi with mass functions fiη based on P. For (x′ , l′ ) 6=
P
(x, l), fiη (x′ , y, l′ ) = (1 − η)f (x′ , y, l′ ). fiη (x, l) = η > 0, fiη (yj | x, l) = αji , where αji ≥ 0, j αji = 1.
Then for each i, CS[Pηi ](X → Y ) can be defined, and when η → 0, fiη converges to f . The total variance
distance between f and fiη is η.
When η → 0,
X XX η
fiη (y ′ | l′ , x′ )
CS[Pηi ](X → Y ) =
fi (x′ , y ′ , l′ ) log P
η ′
′ ′′ η ′′
x′′ ∈X fi (y | l , x )fi (x )
′
′
′
x ∈X y ∈Y l 6=l

+

X X

x′ ∈X y ′ ∈Y

fiη (x′ , y ′ , l) log P

fiη (y ′ | l, x′ )
η ′
′′ η ′′
x′′ ∈X fi (y | l, x )fi (x )

f (y ′ | l′ , x′ )
′
′ ′′
′′
x′′ ∈X f (y | l , x )f (x )
x′ ∈X y ′ ∈Y l′ 6=l
X X
X
X
+
f (x′ , y ′ , l) log f (y ′ | x′ , l) −
f (yj , l) log[f (x)αji +
f (x′ )f (yj | x′ , l)].
→

X XX

f (x′ , y ′ , l′ ) log P
j

x′ 6=x y ′ ∈Y

x′ 6=x

P
For different i, when we let η → 0, the only different terms are − j f (yj , l) log[f (x)αji + x′ 6=x f (x′ )f (yj |
x′ , l)]. We will show that the above term is not a constant with {αji }. Therefore we can find two groups of
{αji } for i = 1, 2 such that g1 = limη→0 CS[Pη1 ](X → Y ) < limη→0 CS[Pη2 ](X → Y ) = g2 .
If there is only one y1 such that f (y1 , l) > 0, then
X
X
−
f (yj , l) log[f (x)αji +
f (x′ )f (yj | x′ , l)]
P

j

= −f (y1 , l) log[f (x)α1i +

x′ 6=x

X

f (x′ )f (y1 | x′ , l)].

x′ 6=x

It is not a constant when we change α1i .
If there are at least two values y1 , y2 of Y , such that f (y1 , l) > 0, f (y2 , l) > 0, then we can change α1i
while keeping α1i + α2i = d, and leave other αji fixed.

ON THE BOUNDARY BETWEEN QUALITATIVE AND QUANTITATIVE METHODS FOR CAUSAL INFERENCE

17

P
P
Set f (y1 , l) = a1 , f (y2 , l) = a2 , f (x) = c, x′ 6=x f (x′ )f (y1 | x′ , l) = b1 , x′ 6=x f (x′ )f (y2 | x′ , l) =
P
P
b2 . All these terms are positive. Then in − j f (yj , l) log[f (x)αji + x′ 6=x f (x′ )f (yj | x′ , l)], terms
containing α1i and α2i are
−a1 log(cα1i + b1 ) − a2 log[c(d − α1i ) + b2 ].
Its derivative with respect to α1i is
−

a1 c
a2 c
+
.
cα1i + b1 c(d − α1i ) + b2

If the derivative always equal 0 in an interval, then we should have
a1
cα1i + b1
≡
,
a2
c(d − α1i ) + b2
which is incorrect.
Now we have two groups of {αji } for i = 1, 2 such that g1 = limη→0 CS[Pη1 ](X → Y ) < limη→0
CS[Pη2 ](X → Y ) = g2 . Then for any g ∈ (g1 , g2 ), any δ > 0, we can find η < δ small enough such that
CS[Pη1 ](X → Y ) < g, CS[Pη2 ](X → Y ) > g. Then we change {αj1 } continuously to {αj2 }. During this
process CS is always defined, and there exists {α3 } such that CS[Pη3 ](X → Y ) = g.
This shows that CS(X → Y ) is essentially ill-defined.
Since CS(X → Y ) and PMI(X, Y | L) have the same non-zero terms containing f (· | x, l), the same
argument shows that PMI(X, Y | L) is not well-defined.
The proof of ATE(X → Y ) is similar, but much more straightforward. We can change the values of
f (yi | x, l) for different i, while keeping their sum. Since ATE(X → Y ) is linear with each f (y | x, l), and
the corresponding coefficients are different, the conclusion is obvious.


C. Proof of Lemma 3 when X is discrete. The proofs for discrete and continuous X are different, therefore we state them separately. Whether Y is discrete or continuous does not matter, therefore we assume Y
is discrete/continuous when X is discrete/continuous. We impose some restrictions to simplify the proofs.
If Zi is discrete, then Ui is an arbitrary discrete random variable which takes all and only the values of Zi
with positive probabilities. If Zi is continuous, then Ui is continuous, and its density function is always
positive.
P
Proof. CMI(X, Y | S1 ) = s1 P(S1 = s1 )CMI(X, Y | S1 = s1 ). For a fixed s1 , assume X can take
values 1, 2, · · · , r ′ , UX can take 1, 2, · · · , r ′ , · · · , r. Y can take values 1, 2, · · · , t with positive probabilities.
P
P
Denote P(X = i, Y = j | S1 = s1 ) by pij . Define p−j = i pij , pi− = j pij . With ǫ-noise, pǫ−j = p−j ,
pǫij = (1 − ǫ)pij + ǫqi p−j , pǫi− = (1 − ǫ)pi− + ǫqi . Then we have
′

CMI(X, Y | S1 = s1 ) =

r
t X
X
j=1 i=1

pij log

pij
,
pi− p−j

r
t X
X
(1 − ǫ)pij + ǫqi p−j
.
[(1 − ǫ)pij + ǫqi p−j ] log
CMI(X , Y | S1 = s1 ) =
[(1 − ǫ)pi− + ǫqi ]p−j
ǫ

j=1 i=1

CMI(X, Y | S1 = s1 ) − CMI(X ǫ , Y | S1 = s1 ) =

18

YUE WANG AND LINBO WANG
r n
t X
X
(1 − ǫ + qi ǫ)pij log

X
pkj
pij
+
qi ǫpkj log
pi− p−j
pk− p−j
j=1 i=1
k6=i
(1 − ǫ)pij + ǫqi p−j o
.
−[(1 − ǫ)pij + ǫqi p−j ] log
p−j [(1 − ǫ)pi− + ǫqi ]
For fixed i, j and k = 1, · · · , r, set
pkj
akij =
,
pk− p−j
ǫqi pk−
bkij =
for k 6= i,
(1 − ǫ)pi− + ǫqi
(1 − ǫ + qi ǫ)pi−
biij =
,
(1 − ǫ)pi− + ǫqi
cij = p−j [(1 − ǫ)pi− + ǫqi ].
Here we know that p−j > 0, (1 − ǫ)pi− + ǫqi > 0. If pk− = 0, namely k = r ′ + 1, · · · , r, then we stipulate
akij = 1. In this case, the corresponding bkij = 0.
Then we have
CMI(X, Y | S1 = s1 ) − CMI(X ǫ , Y | S1 = s1 )
r
r
r
r
t X
X
X
X
X
akij bkij )] ≥ 0.
akij bkij ) log(
bkij akij log akij − (
cij [
=
j=1 i=1

k=1

k=1

k=1

Pr

The last step is Jensen’s inequality, since akij ≥ 0, bkij ≥ 0, k=1 bkij = 1, cij > 0, f (x) = x log x is
strictly convex down when x ≥ 0 (stipulate 0 log 0 = 0).
The equality holds if and only if for each j, a1j = a2j = · · · = ar′ j , which means pij /pi− are equal for
P′
P′
all i ≤ r ′ . Since ri=1 pi− (pij /pi− ) = p−j , ri=1 pi− = 1, we have pij /pi− = p−j for each i, j such that
pi− > 0 and p−j > 0. This is equivalent with that X and Y are independent conditioned on S1 = s1 .
CMI(X, Y | S1 ) = 0 if and only if X and Y are independent conditioned on any possible value of S1 .
Therefore, CMI(X ǫ , Y | S1 ) ≤ CMI(X, Y | S1 ), and the equality holds if and only if CMI(X, Y | S1 ) =
0.

D. Proof of Lemma 3 when X is continuous.
Z ∞
CMI(X, Y | S1 = s1 )h(s1 )ds1 ,
CMI(X, Y | S1 ) =
−∞

where h(s1 ) is the probability density function of S1 . For a fixed s1 , denote the joint probability denR∞
sity function of X, Y conditioned on S1 = s1 by p(x, y). Define p1 (x) = −∞ p(x, y)dy, p2 (y) =
R∞
density function of X, Y conditioned on S1 = s1 is
−∞ p(x, y)dx. With ǫ-noise, the joint probability
R∞
R∞
(1 − ǫ)p(x, y) + ǫq(x)p2 (y). Notice that −∞ q(x)dx = 1, −∞ [(1 − ǫ)p(x, y) + ǫq(x)p2 (y)]dx = p2 (y).
Then we have
CMI(X, Y | S1 = s1 ) − CMI(X ǫ , Y | S1 = s1 )
Z ∞Z ∞
p(x, y)
=
p(x, y) log
dxdy
p1 (x)p2 (y)
−∞ −∞
Z ∞Z ∞
(1 − ǫ)p(x, y) + ǫq(x)p2 (y)
dxdy
−
[(1 − ǫ)p(x, y) + ǫq(x)p2 (y)] log
[(1 − ǫ)p1 (x) + ǫq(x)]p2 (y)
−∞ −∞

∞

∞

i
hZ ∞
p(x, y)
p(x0 , y)
+ q(x)ǫ
dx0
p(x0 , y) log
p1 (x)p2 (y)
p1 (x0 )p2 (y)
−∞ −∞
−∞
(1 − ǫ)p(x, y) + ǫq(x)p2 (y) o
dxdy
−[(1 − ǫ)p(x, y) + ǫq(x)p2 (y)] log
[(1 − ǫ)p1 (x) + ǫq(x)]p2 (y)
For fixed x, y, we can define a probability measure on R, µx,y (x0 ), which is a mixture of discrete and
(1−ǫ)p1 (x)
to take x. For
continuous type measures. For the discrete component, it has probability (1−ǫ)p
1 (x)+ǫq(x)
=

Z

Z

n

(1 − ǫ)p(x, y) log

the continuous component, the probability density function at x0 is

q(x)ǫp1 (x0 )
(1−ǫ)p1 (x)+ǫq(x) .

Define Fx,y (x0 ) =

p(x0 ,y)
p1 (x0 )p2 (y) .

If p1 (x0 ) = 0 or p2 (y) = 0, stipulate Fx,y (x0 ) = 1.
Now we have
CMI(X, Y | S1 = s1 ) − CMI(X ǫ , Y | S1 = s1 )
Z ∞Z ∞
nZ ∞
=
[(1 − ǫ)p1 (x) + ǫq(x)]p2 (y)
Fx,y (x0 ) log Fx,y (x0 )dµx,y (x0 )
−∞ −∞
−∞
hZ ∞
i
hZ ∞
io
−
Fx,y (x0 )dµx,y (x0 ) log
Fx,y (x0 )dµx,y (x0 ) dxdy ≥ 0.
−∞

−∞

The last step is the probabilistic form of Jensen’s inequality, since Fx,y (x0 ) is non-negative and integrable
with probability measure µx,y (x0 ), [(1 − ǫ)p1 (x) + ǫq(x)]p2 (y) > 0 if p2 (y) > 0, and f (x) = x log x is
strictly convex down when x ≥ 0 (stipulate 0 log 0 = 0).
The equality holds if and only if for p1 (x0 ) > 0 and p2 (y) > 0, Fx,y (x0 ) is a constant with x0 , which
R∞
R∞
means p(x0 , y)/p1 (x0 ) is a constant. Since −∞ p1 (x0 )p(x0 , y)/p1 (x0 )dx0 = p2 (y), −∞ p1 (x0 ) = 1, we
have p(x0 , y)/p1 (x0 ) = p2 (y) for each x0 , y such that p1 (x0 ) > 0 and p2 (y) > 0. This is equivalent with
that X and Y are independent conditioned on S1 = s1 .
CMI(X, Y | S1 ) = 0 if and only if X and Y are independent conditioned on any possible value of S1 ,
except a zero-measure set. Therefore, CMI(X ǫ , Y | S1 ) ≤ CMI(X, Y | S1 ), and the equality holds if and
only if CMI(X, Y | S1 ) = 0.
E. Proof of Lemma 4.
Proof. Remember that a Markov boundary M is a minimal subset of S such that MI(M , Y ) = MI(S, Y ).
Denote S with ǫ-noise on Zi ∈
/ M0 by S ǫ . Since MI(M0 , Y ) = MI(S, Y ), MI(M0 , Y ) ≤ MI(S ǫ , Y ),
MI(S ǫ , Y ) ≤ MI(S, Y ), we have MI(S ǫ , Y ) = MI(S, Y ). Therefore, M0 is still a Markov boundary
after adding ǫ-noise. Assume in the new distribution, there is another Markov boundary, then it contains a
variable with ǫ-noise: Ziǫ . Denote this Markov boundary by {Ziǫ } ∪ S1 . Therefore, CMI(Ziǫ , Y | S1 ) > 0.
However, from Lemma 3, this implies CMI(Ziǫ , Y | S1 ) < CMI(Zi , Y | S1 ), namely MI({Ziǫ } ∪ S1 , Y ) <
MI({Zi } ∪ S1 , Y ). But MI({Ziǫ } ∪ S1 , Y ) = MI(S ǫ , Y ) = MI(S, Y ) ≥ MI({Zi } ∪ S1 , Y ), which is a
contradiction.

F. Detailed simulation results. See Tables 1, 2, 3 for detailed performances of different algorithms in
different cases.
R EFERENCES
Aliferis, C., Statnikov, A., Tsamardinos, I., Mani, S., and Koutsoukos, X. (2010). Local causal and Markov
blanket induction for causal discovery and feature selection for classification part i: Algorithms and
empirical evaluation. J. Mach. Learn. Res., 11(Jan):171–234.

20

YUE WANG AND LINBO WANG

TABLE 1. Performances of Algorithms 1, 3-KI, 3-AF, 4 in Case A
Number of
observations
300
500
800
1000
2000
2400
4000
5000
6000
10000
14000
17000
20000
23000
30000
40000
50000
60000

Alg 1
1
1
1
1
1
1
1
1
1
1
1
0.951
1
1
1
1
1
1

Success rate
Alg 3-KI Alg 3-AF
0.957
0.938
0.989
0.981
0.997
0.998
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1

Alg 4
0
0
0
0
0
0
0
0
0
0
0
0.015
0.936
1
1
1
1
1

Alg 1
0.065
0.076
0.092
0.099
0.146
0.156
0.263
0.267
0.343
0.482
0.596
0.803
0.868
1.020
1.285
1.757
2.256
2.890

Time cost
Alg 3-KI Alg 3-AF
0.285
0.453
0.547
0.738
0.665
0.876
0.688
1.271
1.287
1.843
1.350
2.097
2.170
3.332
2.534
3.934
2.905
4.535
4.486
6.794
5.999
9.007
7.670
11.458
9.068
13.412
10.185
15.048
12.787
19.937
20.142
27.387
27.461
33.680
29.637
42.677

Alg 4
0.244
0.283
0.393
0.372
0.597
0.666
0.984
1.200
1.316
1.959
2.629
3.279
4.043
4.547
5.716
7.936
10.769
13.166

TABLE 2. Performances of Algorithms 1, 3-KI, 3-AF, 4 in Case B
Number of
observations
300
500
800
1000
2000
2400
4000
5000
6000
10000
14000
17000
20000
23000
30000
40000
50000
60000

Alg 1
0
0
0
0
0
0
0
0
0
0
0
0
0
0.002
0
0.691
1
1

Success rate
Alg 3-KI Alg 3-AF
0.326
0.282
0.750
0.680
0.990
0.944
1
0.970
1
0.999
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1

Alg 4
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1
1

Alg 1
0.061
0.071
0.082
0.087
0.130
0.146
0.227
0.264
0.292
0.425
0.574
0.742
0.811
1.013
1.270
1.696
2.188
3.075

Time cost
Alg 3-KI Alg 3-AF
0.175
0.461
0.396
0.798
0.482
0.912
0.583
0.964
1.021
1.414
1.125
1.575
1.828
2.552
2.106
2.920
2.538
3.443
3.696
5.050
4.933
6.837
6.246
8.535
7.370
9.855
8.206
11.275
10.788
14.260
15.350
20.513
19.669
25.121
24.790
32.176

Alg 4
0.261
0.287
0.352
0.360
0.549
0.605
0.987
1.144
1.295
1.905
2.566
3.251
3.765
4.352
5.559
7.966
9.894
13.322

ON THE BOUNDARY BETWEEN QUALITATIVE AND QUANTITATIVE METHODS FOR CAUSAL INFERENCE

21

TABLE 3. Performances of Algorithms 1, 3-KI, 3-AF, 4 in Case C
Number of
observations
300
500
800
1000
2000
2400
4000
5000
6000
10000
14000
17000
20000
23000
30000
40000
50000
60000

Alg 1
1
1
1
1
1
0.356
0
0.008
0.694
1
1
1
1
1
1
1
1
1

Success rate
Alg 3-KI Alg 3-AF
0.004
0.542
0.003
0.921
0.002
0.999
0.005
1
0.001
1
0.001
1
0.001
1
0.001
1
0.001
1
0.002
1
0.001
1
0
1
0.001
1
0.002
1
0.001
1
0.005
1
0
1
0.001
1

Alg 4
0
0
0
0
0
0
0
0
0
0
0
0
0
0.752
1
1
1
1

Alg 1
0.054
0.062
0.074
0.081
0.119
0.127
0.202
0.236
0.275
0.434
0.566
0.701
0.844
0.912
1.195
1.694
2.183
2.788

Time cost
Alg 3-KI Alg 3-AF
0.162
0.444
0.186
0.722
0.520
1.280
0.799
1.407
1.241
2.145
1.383
2.438
2.246
3.824
2.697
4.478
3.118
5.338
4.774
7.887
6.397
10.525
8.147
13.392
9.450
15.485
10.858
17.535
13.658
22.432
20.002
32.030
28.318
39.914
30.695
48.107

Alg 4
0.245
0.276
0.299
0.335
0.507
0.547
0.857
1.019
1.194
1.829
2.414
3.054
3.476
4.039
5.429
7.731
9.628
12.182

Benjamini, Y. and Yekutieli, D. (2001). The control of the false discovery rate in multiple testing under
dependency. Ann. Stat., 29(4):1165–1188.
Cover, T. M. and Thomas, J. A. (2012). Elements of Information Theory. John Wiley & Sons.
de Morais, S. R. and Aussem, A. (2010). A novel Markov boundary based feature subset selection algorithm.
Neurocomput., 73(4):578–584.
Dobrushin, R. L. (1963). General formulation of Shannon’s main theorem in information theory. Amer.
Math. Soc. Trans., 33:323–438.
Ein-Dor, L., Kela, I., Getz, G., Givol, D., and Domany, E. (2004). Outcome signature genes in breast cancer:
is there a unique set? Bioinformatics, 21(2):171–178.
Frenzel, S. and Pompe, B. (2007). Partial mutual information for coupling analysis of multivariate time
series. Phys. Rev. Lett., 99(20):204101.
Geweke, J. F. (1984). Measures of conditional linear dependence and feedback between time series. J. Am.
Stat. Assoc., 79(388):907–915.
Granger, C. W. J. (1969). Investigating causal relations by econometric models and cross-spectral methods.
Econometrica, pages 424–438.
Guo, S., Seth, A. K., Kendrick, K. M., Zhou, C., and Feng, J. (2008). Partial Granger causality–eliminating
exogenous inputs and latent variables. J. Neurosci. Meth., 172(1):79–93.
Huang, T.-M. (2010). Testing conditional independence using maximal nonlinear conditional correlation.
Ann. Stat., 38(4):2047–2091.
Janzing, D., Balduzzi, D., Grosse-Wentrup, M., and Schölkopf, B. (2013). Quantifying causal influences.
Ann. Stat., 41(5):2324–2358.
Kaiser, A. and Schreiber, T. (2002). Information transfer in continuous processes. Physica D, 166(1):43–62.

22

YUE WANG AND LINBO WANG

Mani, S. and Cooper, G. (2004). Causal discovery using a bayesian local causal discovery algorithm.
Medinfo, 11(Pt 1):731–735.
Meyer, P. E. (2008). Information-theoretic variable selection and network inference from microarray data.
PhD thesis, Université Libre de Bruxelles.
Neapolitan, R. (2004). Learning Bayesian Networks. Pearson Prentice Hall Upper Saddle River, NJ.
Neyman, J. (1923). Sur les applications de la thar des probabilities aux experiences Agaricales: Essay des
principle. English translation of excerpts by Dabrowska, D. and Speed, T. (1990). Statistical Science,
5:463–472.
Pearl, J. (1988). Probabilistic Inference in Intelligent Systems. Morgan Kaufmann, San Mateo.
Pearl, J. (2009). Causality. Cambridge University Press.
Pearl, J. and Paz, A. (1985). Graphoids: A Graph-based Logic for Reasoning about Relevance Relations.
University of California (Los Angeles). Computer Science Department.
Peña, J., Nilsson, R., Björkegren, J., and Tegnér, J. (2007). Towards scalable and data efficient learning of
Markov boundaries. Int. J. Approx. Reason., 45(2):211–232.
Richardson, T. S. and Robins, J. M. (2013). Single world intervention graphs (swigs): A unification of the
counterfactual and graphical approaches to causality. Center for the Statistics and the Social Sciences,
University of Washington Series. Working Paper, 128(30):2013.
Robins, J. (1986). A new approach to causal inference in mortality studies with a sustained exposure periodapplication to control of the healthy worker survivor effect. Mathematical Modelling, 7(9-12):1393–
1512.
Rubin, D. B. (1974). Estimating causal effects of treatments in randomized and nonrandomized studies.
Journal of Educational Psychology, 66(5):688.
Schreiber, T. (2000). Measuring information transfer. Phys. Rev. Lett., 85(2):461.
Spirtes, P., Glymour, C., and Scheines, R. (2000). Causation, Prediction, and Search. MIT press, 2nd
edition.
Statnikov, A., Lytkin, N., Lemeire, J., and Aliferis, C. (2013). Algorithms for discovery of multiple Markov
boundaries. J. Mach. Learn. Res., 14(Feb):499–566.
Su, L. and White, H. (2008). A nonparametric hellinger metric test for conditional independence. Econ.
Theory, 24(04):829–864.
Tsamardinos, I. and Aliferis, C. (2003). Towards principled feature selection: relevancy, filters and wrappers. In Proceedings of the Ninth International Workshop on Artificial Intelligence and Statistics.
Zhang, K., Peters, J., Janzing, D., and Schölkopf, B. (2011). Kernel-based conditional independence test
and application in causal discovery. In Proceedings of the Twenty-Seventh Conference on Uncertainty in
Artificial Intelligence, pages 804–813. AUAI Press.
Zhao, J., Zhou, Y., Zhang, X., and Chen, L. (2016). Part mutual information for quantifying direct associations in networks. Proc. Natl. Acad. Sci., 113(18):5130–5135.

