On the Identifiability and Estimation of Functional Causal Models in
the Presence of Outcome-Dependent Selection

Kun Zhang⇤Jiji Zhang]kunz1@cmu.edu jijizhang@ln.edu.hk

Biwei Huang‡⇤

Bernhard Schölkopf‡

Clark Glymour⇤biweih@andrew.cmu.edu

bs@tuebingen.mpg.de

cg09@andrew.cmu.edu⇤

Department of philosophy, Carnegie Mellon University]
Lingnan University, Hong Kong‡
Max Planck Institute for Intelligent Systems, Germany

Abstract
We study the identifiability and estimationof functional causal models under selection
bias, with a focus on the situation where theselection depends solely on the e↵ect variable,
which is known as outcome-dependent selection. We address two questions of identifiability: the identifiability of the causal direction between two variables in the presence of
selection bias, and, given the causal direction,
the identifiability of the model with outcomedependent selection. Regarding the first, weshow that in the framework of post-nonlinear
causal models, once outcome-dependent selection is properly modeled, the causal direction between two variables is generically
identifiable; regarding the second, we identify
some mild conditions under which an additivenoise causal model with outcome-dependent
selection is to a large extent identifiable. We
also propose two methods for estimating an
additive noise model from data that are generated with outcome-dependent selection.

1

Introduction

Selection bias is an important issue in statistical inference. Ideally, samples should be drawn randomly
from the population of interest. In reality, however,
it is commonplace that the probability of including aunit in the sample depends on some attributes of theunit. Such selection bias, if not corrected, often distorts the results of statistical analysis. For example, itis well known that in a regression analysis, if there isselection on the dependent variable, the ordinary leastsquares estimation of the regression coefficients will bebiased and inconsistent (Heckman, 1979). The challenge is even bigger in causal inference; both the taskof learning causal structures from data and the taskof estimating causal mechanisms or parameters given
a causal structure are usually rendered more difficult
by the presence of selection bias.
In this paper, we are concerned with the approach tocausal inference based on (restricted) functional causal
models (Shimizu et al., 2006; Hoyer et al., 2009; Zhang& Hyvärinen, 2009), and aim to investigate the extent
to which selection bias can be handled within this approach. Specifically, we mainly focus on the outcomedependent selection bias, where the selection mechanism depends only on the e↵ect, and are interested in
the following two questions:• Is the causal direction between two random variables identifiable in the presence of selection bias?
• Is the causal mechanism as represented by a functional causal model identifiable in the presence of
selection bias?
These two questions have to do with the two main
aspects of causal inference, respectively. The former question is about the inference of causal structure. In the traditional conditional-independenceconstraint-based approach to learning causal structures (Spirtes et al., 2001; Pearl, 2000), some methods
have been developed to handle selection bias (Spirtes
et al., 1999; Zhang, 2008; Borboudakis & Tsamardinos, 2015). However, the structural information thatcan be learned via the constraint-based approach istypically limited to a Markov equivalence class. In
particular, the approach cannot distinguish cause frome↵ect with just two variables. In contrast, a distinctive virtue of the approach based on functional causal
models is that Markov equivalent causal structures canusually be distinguished. In particular, the direction
between two random variables is generically identifiable. Whether this virtue survives the challenge posedby selection bias is therefore worth investigating.
The latter question is related to the inference of causal
parameters (i.e., parameters or quantities that have

causal interpretations), including intervention e↵ects.
In addition to the work on various selection models
in econometrics and social science (Heckman, 1979;
Winship & Mare, 1992), recent literature has seen interesting work on the recoverability of causal parameters based on graphical models (Didelez et al., 2010;
Bareinboim & Pearl, 2012; Bareinboim et al., 2014;
Evans & Didelez, 2015). Much of this work, however,deals with linear models or discrete variables, whereas
we are concerned in this paper with continuous variables that may bear a nonlinear relationship.
We will proceed as follows. In Section 2, we introduce
the general setup and briefly discuss several types of
selection, before focusing our attention on the situation where the selection depends on the e↵ect variable,
known as outcome-dependent selection. In Section 3,
we show that in the framework of post-nonlinear causal
models, once outcome-dependent selection is properly
modeled, the causal direction between two variables isgenerically identifiable. In Section 4, we identify somemild conditions under which an additive noise causal
model with outcome-dependent selection is to a large
extent identifiable. We then propose, in Section 5, two
methods for estimating an additive noise model fromdata that are generated with outcome-dependent selection. Some experiments are reported in Section 6.

2

Outcome-Dependent Selection Bias

A common way to represent selection bias is to use abinary selection variable S encoding whether or nota unit is included in the sample. Suppose we are interested in the relationship between X and Y , where
X has a causal influence on Y . Let pXY denote thejoint distribution of X and Y in the population. The
selected sample follows pXY |S=1 instead of pXY . In
general, pXY |S=1 6= pXY , and that is how selection
may distort statistical and causal inference. However,di↵erent kinds of selection engender di↵erent levels ofdifficulty. In general, S may depend on any number ofsubstantive variables, as illustrated in Figure 1, whereX = (X1 , X2 ). 1
1

In this paper, we assume that we only know which variables the selection variable S depends on, but the selection
mechanism is unknown, i.e., the probability of S = 1 given
those variables is unknown. Notice that we do not haveaccess to the data points that were not selected. This isvery di↵erent from Heckman’s framework to correct thebias caused by a censored sample (Heckman, 1979), whichassumes access to an i.i.d. sample from the whole population, on which the Y values are observable only for the data
points that satisfy the selection criterion (implied by theselection equation), but other attributes of the “censored”
points are still available, enabling one to directly identify
the selection mechanism.

W

X1

X2

Y

W

S

X1

X2

S

U

(a)
W

X1

(b)
X2

Y
S

(c)

Y

W

X1

X2

Y

S

(d)

Figure 1: Illustration of di↵erent situations with sample selection bias. (a) S depends on X = (X1 , X2 ) but
not on Y . (b) S depends on X and is also statisticallydependent on Y given X due to a confounder U . (c)S directly depends solely on Y (outcome-dependent
selection). (d) S depends on both X and Y .
Selection Bias on the Cause For the purpose of
causal inference, the least problematic kind of situation is depicted in Figure 1(a), in which S is independent of the e↵ect variable Y given the cause variable
X. It follows that pY |X,S=1 = pY |X . That is, theselection bias does not distort the conditional distribution of the e↵ect Y given the cause X or the structural equation model for the causal process. In sucha situation, causal inference can essentially proceed as
usual. However, if there is a (latent) confounder for
Y and S, as illustrated in Figure 1(b), S and Y are
not conditionally independent given X any more, thatis, pY |X,S=1 6= pY |X . Such a distortion may be corrected under rather restrictive assumptions; see, e.g.,
Heckman’s correction (Heckman, 1979).
Selection Bias on the E↵ect If the selection depends solely on the e↵ect, as depicted in Figure 1(c),then pY |X,S=1 6= pY |X , and the selection bias, if notcorrected, will mislead inference. Consider, for example, a standard assumption in functional causal modeling that the e↵ect Y is a function of the cause variable X and an noise variable E that is independentof X. Suppose this assumption holds in the population. With the outcome-dependent selection, X and
E are typically not independent in the selected sample, as they are typically not independent conditionalon S (which is a descendant of a collider between X
and E, i.e., Y ). Furthermore, even if one fits a regression model on selected sample, the estimated residual(which is usually di↵erent from the true noise term in
the causal process) is usually not independent from X;
we will get back to this issue in Section 4.1.
This kind of selection is known as outcome-dependent
selection bias (OSB) (Didelez et al., 2010; Bareinboim
et al., 2014), and will be our focus in this paper. We
will show that although outcome-dependent selection
seriously complicates analysis, it can be handled in
the identification and estimation of functional causal

models. Note that in the case of outcome-dependent
selection, X is independent of S given Y , and so we
can model the distribution of the observed sample as:
pX,Y,S=1
P (S = 1|X, Y )= pXY ·P (S = 1)
P (S = 1)
P (S = 1|Y )= pXY ·
= (y)pXY ,(1)
P (S = 1)

pXY , pXY |S=1 =where the nonnegative function (y) , P (S =1|Y )/P (S = 1) is a density ratio for biased samplingthat only depends on Y . We will adopt this representation of outcome-dependent selection in what follows.
Selection Bias on Both the Cause and the Effect An even more general situation is depicted in
Figure 1(d), where the selection depends on both X
and Y (and probably others). In such a situation,
the density ratio function
will depend on both X
and Y . The selected sample follows the distribution
pXY / pXY (x, y, w). Roughly speaking, the selection procedure is so flexible that without further
constraints on (x, y, w), we cannot see much information about the population pXY : if pXY is positive on ( 1, +1), the same pXY can be generated
from a large class of distributions pXY with a suitablychosen (x, y, w). Moreover, the causal direction is
generally not identifiable, for with a sufficiently flexible (x, y, w), either direction can be made compatible with whatever distribution. Interestingly, when
depends only on Y , as is the case under outcomedependent selection, the causal direction according toa restricted functional causal model is still generically
identifiable, without any substantial restriction on .
To this result we now turn.

3

Identifiability of Causal Direction

In this section we investigate whether it is possible tosuccessfully recover the causal direction between two
variables when the data are generated according to afunctional causal model, but with outcome-dependent
selection. Here we assume that both X and Y arescalar variables.
3.1

Identifiability Without Selection Bias

The traditional approaches to inferring causal structure from data, such as the constraint-based approach(Spirtes et al., 2001; Pearl, 2000) and the score-based
approach (Chickering, 2002; Heckerman et al., 1995)
cannot distinguish Markov equivalent causal structures without background knowledge. In particular,
with only two variables, those methods cannot distinguish cause from e↵ect. The more recent approachbased on restricted functional causal models is usuallymore powerful in this respect. In a functional causal
model, the e↵ect is taken to be a function of the directcauses together with an noise term that is independentof the direct causes (Pearl, 2000). When the class offunctions is constrained, the causal direction is usually
identifiable in that only one direction can satisfy themodel assumptions, such as the assumed independence
between the noise term and the direct causes. Available identifiability results include those on linear, nonGaussian, acyclic Model (LiNGAM) (Shimizu et al.,2006)), additive noise model (ANM) (Hoyer et al.,2009), and post-nonlinear (PNL) causal model (Zhang& Hyvärinen, 2009). In this section, we will establish amain result for the PNL causal model. The result also
applies to linear models and additive noise models, as
they are special cases of PNL models.
A PNL model for X ! Y is specified as follows:
Y = f2 (f1 (X) + E),(2)

where X and E are statistically independent, f1 is
a non-constant smooth function, f2 is an invertiblesmooth function, and f20 6= 0. This model is sufficiently
flexible to represent or approximate many causal processes in reality (Zhang & Hyvärinen, 2009).
Similarly, for the reverse direction Y ! X, a PNL
model would take the following form:
X = g2 (g1 (Y ) + Ẽ),(3)

where Y and Ẽ are independent, g1 is non-constantand smooth, g2 is invertible and smooth, and g20 6= 0.
As shown in (Zhang & Hyvärinen, 2009), (2) and (3)can generate the same distribution of X and Y only
for very special configurations of the functions and distributions. In generic cases, if data are generated according to a model of form (2), there is no model ofform (3) that generates the same distribution. Hencethe causal direction is generically identifiable.
3.2

Identifiability of Causal Direction in
PNL-OSB

We now show that the generic identifiability of causal
direction based on PNL models still holds even if weallow the possibilty of outcome-dependent selection.
Suppose the data distribution is generated by a PNL
causal model from X to Y in the form of (2), denotedby F! , followed by an outcome-dependent selection
with an density ratio (y), as in (1). Call (F! , (y))
a PNL-OSB model, and let p!
XY denote the joint density of X and Y resulting from (F! , (y)). We areinterested in whether there is a PNL-OSB model inthe reverse direction that can generate the same data
distribution. That is, consider (F , v(x)), where F

is a PNL causal model from Y to X in the form of(3), and v(x) is an density ratio function that depends
on X. Let pXY denote the joint density of X and Y
resulting from (F , v(x)). When is it the case thatp!
XY = pXY ?
To simplify the presentation, we define random variables T , g2 1 (X), Z , f2 1 (Y ), and function h ,f1 g2 . That is, h(t) = f1 (g2 (t)) = f1 (x). Similarly, h1 , g1 f2 is a function of Z. Moreover,
we let ⌘1 (t) , log pT (t) = log pX (x) + log |g20 (t)|, and⌘2 (e) , log pE (e).
Note that T and E are independent (for X and E areassumed to be independent), and Z and Ẽ are independent (for Y and Ẽ are assumed to be independent).
It follows thatF!
0p!
XY = (y)pXY = (y)pXE /|f2 | =0 0
f2 (z)pT pE /|f2 g2 |,0
0 0
pXY = v(x)pF
XY = v(x)pY Ẽ /|g2 | = vg2 (t)pZ Ẽ /|f2 g2 |,

where

f2 =

f2 , and vg2 = v g2 .

Now suppose

p!
XY

= pXY

(4)

This implies
pZ Ẽ =

f2 (z)

vg2 (t)pT pE ,

or equivalently
log pZ Ẽ = log= log

f2 (z)

log vg2 (t) + log pT + log pE

˜1 (t) + ⌘2 (e),f2 (z) + ⌘where ⌘˜1 (t) , log pT log vg2 (t) = ⌘1 (t)Since Z and Ẽ are independent, we have(5)
log vg2 (t).

@ 2 log pZ Ẽ
⌘ 0.
(6)
@z@ẽ
(5) and (6) entail very strong constraints on the distribution of E, as stated in the following theorem.
Theorem 1 Suppose that the densities of E and
T and the functions f1 , f2 , g1 , g2 , and v(x) are
third-order di↵erentiable and that pE is positive on( 1, +1). The condition (4) implies that for every
point of (X, Y ) satisfying ⌘200 h0 6= 0:⇣ ⌘ 0 ⌘ 000⌘⌘˜100 h00
⌘2000 0 00
2 2
00
0 00
⌘˜1000=2⌘·hh· h ⌘˜1
2
00
h0
⌘2
⌘200
⇣
h002 ⌘+ ⌘20 · h000,
(7)h0
and h1 depends on ⌘˜1 , ⌘2 , and h in the following way:1
⌘˜100 + ⌘200 h02 ⌘20 h00=.
h01⌘200 h0(8)Further assume that ⌘200 h0 6= 0 almost everywhere.
Then in order for (7) to hold, pE and h must satisfy
one of the five conditions listed in Table 1.

Table 1: All situations in which the causal direction
implied by the PNL-OSB model may be unidentifiable.
1
2
3
4

pE
Gaussianlog-mix-lin-exp
log-mix-lin-exp
generalized mixtureof two exponentials

h = f 1 g2
linearlinearh strictly monotonic,
and h0 ! 0, as t1 !
+1 or as t1 ! 1Same as above

All proofs are given in the Supplementary material.
In the five situations given in Table 1, the causal direction may not be identifiable according to the PNLOSB model, and the involved distribution pE is very
specific. For the definition of distributions of the form
log-mix-lin-exp or generalized mixture of two
exponentials, see (Zhang & Hyvärinen, 2009). As aconsequence, generally speaking, the causal direction
implied by PNL-OSB is identifiable.
This identifiability result regarding the causal direction implied by PNL-OSB is similar to the original result on PNL, which was given in (Zhang & Hyvärinen,2009). The di↵erence is that ⌘1 (t) = log pT (t) in theoriginal identifiability result on PNL is replaced by⌘˜1 (t) = log vpgT (t)(t) . Recall that vg2 (t) can be any valid
2
density ratio; if pT (t) is positive on ( 1, +1), one
can always adjust vg2 (t) so that vpgT (t)meets the con2 (t)straint on ⌘1 in (Zhang & Hyvärinen, 2009). That is,in our result any pT (t) that is positive on ( 1, +1)
is allowed. Therefore, our non-identifiable situations(Table 1) do not contain any constraints on pT , but
still have very strong constraints on PE and h = f1 g2 .

4

Identifiability of ANM-OSB Model

Given the causal direction, a further important question is whether the causal mechanism, represented by
the functional causal model, and the selection procedure, represented by (y), can be recovered from data.
For simplicity of the derivation and presentation, weshall consider the ANM for the causal mechanism (nota PNL one in this section):Y = f AN (X) + E,

(9)where E ?
? X. Here we further assume that f AN issmooth. The observed data are generated by applying
the selection bias on Y , i.e., they were drawn from the
distributionF
pXY = (y)pF(10)
X pY |X ,where pF
Y |X is specified by the causal model (9) andF
pX denotes the distribution of X before applying theSuppose that the observed data are generated from an
ANM-OSB (F1 , 1 (y)). We are interested in whether
another ANM-OSB (F2 , 2 (y)) can generate the samedata distribution. Suppose it does. The observed data
distribution is then1
pXY = pF
XY

Let

r (y)

(1)
pX pE1 (Y

2 (y)

,
(1)

1 (y)

.

F2
1 (y) = pXY

1

(11)

1
Bear in mind that pF
XY2(X)) and pF
XYf
If (11) holds, we haver

2 (y).

(1)(2)= pX pE2 (Yf(2)(y)pX (x)pE1 (e1 ) = pX (x)pE2 (e2 ).

(2)

=

(X)).

log

(1)
(2)
r (y)+log pX +log pE1 (e1 ) = log pX +log pE2 (e2 ).

(13)
Now let us see whether it is possible for (13) tohold and, if yes, what constraints the functions r (y),(1)
log pX , and log pE1 (e1 ) must satisfy. Denote by(2)
JAN , log pX + log pE2 (e2 ). As seen from the RHS,
(13) implies@ 2 JAN
⌘ 0.
(14)
@x@e2
(1)

Let l (y) = log r (y), ⌘X1 (x) , log pX and ⌘E1 (e1 ) ,log pE1 (e1 ). By solving (14), we can establish the relationship between the two ANM-OSB models.
4.1

General Results

Interestingly, as stated in the following theorem, if
the noise E1 is non-Gaussian, then f (2) (x) must be ashifted version of f (1) (x); in other words, the underlying function f AN is identifiable up to a constant. Furthermore, if E1 is non-Gaussian, the selection weight(y) can be recovered up to a factor which is an exponential function of y, i.e., 2 (y) / 2 (y) · ec2 y , wherec2 is a constant; accordingly, pE2 / pE1 · e c2 e1 .
Theorem 2 Let Assumptions A1 and A2 hold true:(1)
A1 . pX and pE1 are positive on ( 1, +1).
000
A2 . ⌘ E(e1 )f (1) (x) = 0 only at finite points.2
1
Then if (11) is true, the following statements hold.
a) If E1 is not Gaussian, then f (2) (x) = f (1) (x)+c1 ,and 2 (y) = 1 (y) r (y), where r (y) = ec2 y+d1 =20This excludes the special case where f (1) ⌘ 0, i.e.,
where X and Y are independent; in this case clearly theselection procedure is not identifiable.

(1)

ed1 · ec2 f (x) · ec2 e1 . Accordingly, pX / pX ·e c2 f1 (x) , and pE2 / pE1 · e c2 e1 . Here c1 , c2 , andd1 are constants, and d1 guarantees that 2 (y) is
a valid density ratio.
That is, f (2) (x) is equal to f (1) (x) (up to a constant), and with proper scaling, 2 (y) equals 1 (y)times an exponential function of y.
b) If E1 is Gaussian, then 2 (y) = 1 (y) r (y),ab 2where r (y) = e 2 y +c4 y+d4 , and f (2) (x) =1(1)(x)+d3 . Here a, b, c4 , d3 , and d4 are con1+b f
stants, d4 guarantees that 2 (y) is a valid density
ratio, and a 6= 0.
That is, with proper scaling, 2 (y) equals 1 (y)times a Gaussian function of y (which includes
the exponential function of y as a special case by
setting b = 0).

(12)

Taking the logarithm of both sides gives

(2)

(1)

selection procedure. Note that generally speaking, pF
X
is not identical to pX . Call the model (F, (y)) an
ANM-OSB model.

An interesting implication of Theorem 2 is that generally speaking, fitting an ordinary ANM on the datathat were generated by an ANM-OSB will not produce an independent error term. That is, under
mild assumptions, if one sets 2 (y) ⌘ 1, (F2 , 2 (y))
cannot produce the same distribution over (X, Y ) as
(F1 , 1 (y)) does, as is given in the following corollary.
Corollary 3 Let Assumptions A1 and A2 hold. Then
under either of the following conditions, there does notexist an ANM, specified by (9), to generate the same
distribution over (X, Y ) as (F1 , 1 (y)) does.
a) E1 is not Gaussian and0
ec y for any c0 .

1 (y) is not proportional tob) E1 is Gaussian and 1 (y) is not proportional to0 2
0ea y +c y for any a0 and c0 (i.e., 1 (y) is not proportional to an exponential function of any polynomial of
y of degree 1 or 2).
4.2

When the Noise is Gaussian

When pE1 is Gaussian, as stated in b) of Theorem 2,
the function f AN is not identifiable any more: f (2) (x)and f (1) (x)can di↵er by an affine transformation (notsimply a shift). Accordingly, 2 (y) can di↵er from
1 (y) in a Gaussian function of y. Compared to thecase where E1 is non-Gaussian, the Gaussian case suffers from more indeterminacies because the productof two Gaussian functions is still a Gaussian function. In particular, since r 1 (y) and pE1 (y f (1) (x))
1(or pF
Y |X (y|x)) are both Gaussian functions in y, theirproduct is still a Gaussian function in y. Accordingly,(12) will hold true, by setting pE2 to another appropriate Gaussian density; in other words, in this case two
additive noise models F1 and F2 can both generatethe same observed data, by applying the bias selection
procedures 1 (y) and 2 (y) = 1 (y) r (y), respectively.
More specifically, we can derive the function f (2)(x)and noise distribution pE2 (e2 ) for the model F2 . As
1shown above, f (2)(x) = 1+b
f (1) (x) + d3 . Eq. 12,
combined with (25) and (26), implies that pE2 (e2 ) /a(1+b) 2
a 2
ab 2
a 2
e 2 e2 + 2 e2 = e 2 e2 , while pE1 (e1 ) / e 2 e1 .
Figure 2 gives an illustration of this result. Noticethat the identifiability results imply some constraints
on r (y) = 2 (y)/ 1 (y), not on 1 (y), so without lossof generality, we set 1 (y) ⌘ 1, leading to 2 (y) =r (y). The circles denote the data points generated byapplying the density ratio 1 (y) ⌘ 1 on additive noisemodel F1 : the dash line shows the nonlinear function
f (1) (x), and the red solid line shows the shape of pE1 .
In contrast, additive noise model F2 uses nonlinear
function f (2) (x), which is di↵erent from f (1) (x), andits noise variance is slightly larger than that in F1 .
The crosses denote the data points generated by F2 .
Although F1 and F2 are not the same in this case,2
applying the density ratio function r (y) on pF
XY gives
F1
1
the same joint distribution as pF
XY , i.e., 1 (y)pXY =F1F2F2
pXY = r (y)pXY = 2 (y)pXY .

the ANM model, including the function, the noise distribution, and the density ratio, can be fully identifiable. Below is an example showing that this is the case
if we know that pE is symmetric and non-Gaussian.
Corollary 4 Let the assumptions made in Theorem 2
A1 and A2 hold. Suppose E1 is not Gaussian. Then If
both pE1 and pE2 are symmetric about the origin, thenf (2) (x) = f (1) (x), E1 = E2 , pE1 (e1 ) = pE2 (e2 ), andr (y) ⌘ 1, i.e., 1 (y) = 2 (y).

5

Estimation of ANM-OSB

Eq. 10 gives the distribution for the observed data. In
theory, we can then estimate the parameters involvedin (y), pF
Y |X , as well as pX , by maximum likelihood.
However, when using maximum likelihood, we have toguarantee that the quantity on the right hand side of(10) is a valid density. This constraint is notoriously
difficult to enforce in the optimization procedure. Below we propose two methods to estimate the underlying additive noise model and (y); one is maximum
likelihood with the above constraint enforced approximately, and the other makes use of the score matching
technique.
5.1

Maximum Likelihood Estimation with a
Sample Approximation ConstraintF
To estimate the involved functions (y), pF
X , and pY |X ,we can maximize the data likelihood:

L=
=

n
X

k=1
n
X
k=1

Figure 2: Illustration of the non-identifiability of theadditive noise model, especially f AN , when the noise
is Gaussian. Red circles denote data points generated
by the ANM F1 , or by the ANM-OSB (F1 , 1 (y) ⌘1). The gray crosses denote data generated by the
ANM F2 . The two ANM-OSB models, (F2 , r (y)) and(F1 , 1), produce the same distribution of (X, Y ).
4.3

With Further Constraints

Not surprisingly, if we have more knowledge about thenoise distribution pE or the density ratio function (y),log pXY (xk , yk )⇥⇤Flog (yk ) + log pF
X (xk ) + log pY |X (yk |xk ) . (15)

According to the theory shown in Section 4, the soF
lution to y, pF
X , and pY |X su↵ers from some indeterminacies, e.g., the solution to y may di↵er from thetrue one by an exponential transformation. To find thesolution for which the biased selection procedure is as
weak as possible, we regularize the likelihood function
with the constraint that log (y) is close to 0. That is,
we maximize
Lr = Lrn
Xp(log (yk ))2 + r,

(16)

k=1

where the regularization parameter r was set to 10 3
in our experiments, and r is a small positive number
and was set to 0.02.
Now we have two issues to consider. One is how toparameterize the involved functions. The other is how
to enforce that pXY , specified in (10), corresponds toa valid density. More specifically, the constraint isF(y)pF
X pY |X > 0, or equivalently (y) > 0, and
Z
F
(y)pF
X pY |X dxdy = 1.

(17)
(18)

Without constraint (18), the scale of pXY will go toinfinity during the process of maximizing (15).
Parameterization The additive noise model for thedata-generating process, (9), implies that pF
Y |X =pE (y f AN (x)). We parameterize (y) as the exponential transformation of a nonlinear function represented by MLP’s (with the tanh activation function);
this automatically guarantees the nonnegativity constraint of (y), as required in (17). Furthermore, we
represent pF
X with a mixture of Gaussians, the nonlinear function f AN with MLP’s (with the tanh activation function), and pE with another mixture of Gaussians.
Enforcing pXY to Be a Valid Density We present
a sample-average approximation scheme to approximately enforce the condition that the right hand side of(10) corresponds to a valid distribution, or more specifically, to impose the constraint (18). Notice that thegiven data points {xk , yk }nk=1 were drawn from pXY .
As a matter of fact, we have
Z
Z
F
(y)pF
X pY |X
F
o
(y)pF
p
=
p
dxdy (19)X Y |X
XY
poXY
n
F
F
1 X (yk )pX (xk )pY |X (yk |xk )⇡(20)
npoXY (xk , yk )
1
⇡
n

k=1
n
X

F
(yk )pF
X (xk )pY |X (yk |xk )p̂oXY (xk , yk )

k=1

,

(21)

where poXY denotes the data distribution of (X, Y ),and p̂oXY (xk , yk ) denotes its estimate at point (xk , yk ).
Here the expression in (20) is an empirical estimate of(19) on the sample drawn from the distribution poXY ;furthermore, (20) replaces the density poXY (xk , yk )with its empirical estimate p̂oXY (xk , yk ). As a consequence, the constraint (18) can be (approximately)achieved by enforcingnFF
1 X (yk )pX (xk )pY |X (yk |xk )= 1.
np̂oXY (xk , yk )

(22)

k=1

In our experiments, we used kernel density estimation with a Gaussian kernel for p̂oXY (xk , yk ); for eachdimension, we set the kernel width to the median distance between points in the sample, as in (Gretton
et al., 2007).
Under the parameterization given in (27) and with theabove approach to guarantee that poXY is (approximately) a valid density, one can then maximize thelikelihood function given in (15) to estimate the function f AN , the noise distribution, and (y).
5.2With Score Matching

Alternatively, we can estimate the parameters by scorematching (Hyvärinen, 2005), i.e., by minimizing the
expected squared distance between the gradient of thelog-density given by the model and the gradient of thelog-density of the observed data. This procedure aims
to match the shape of the density given by the modeland that of the empirical density of the observed data,and is invariant to the scaling factor of the model density. As a clear advantage, in the optimization procedure one does not need to guarantee that pXY is avalid density.
Given any model density pZ (z; ✓) of a m-dimensionalrandom vector Z, the score function is the gradient of the log-density w.r.t.
the data vector, i.e.,(z; ✓) = ( 1 (z; ✓), ..., m (z; ✓))| =pZ (z;✓)pZ (z;✓) |, ..., @ log@z) . Note that the score( @ log@z
1mfunction is invariant to scale transformations inpZ (z), i.e., it is invariant to the normalization constant for a valid density. One can then estimate
model parameters by minimize the expected squared
distance between the model score function (·; ✓)andR the data score function 2Z (·; ✓), i.e., minimize
1
Z (z)|| dz. It has been shown
2 z2Rm pZ (z)|| (z; ✓)in (Hyvärinen, 2005) that minimizing the above
squared distance is equivalent to minimizing
J SM (✓) =Z

pZ (z)
z2Rm

m
X
⇥
i=1

˜i (z; ✓) + 1
2⇤2i (z; ✓) dz,i (z;✓)
where ˜i (z; ✓) = @ @z
. The sample version ofiSMJ(✓) over the sample z1 , ..., zn is

1JˆSM (✓) =
n

n X
m
X
⇥

k=1 i=1

˜i (zk ; ✓) + 1
2⇤

2i (zk ; ✓) .

(23)

In particular, here we have 1 = X and 2 = Y ;noting that pF
f AN (x)), we can write
Y |X = pE (ydown the involved derivatives involved in (28), andthen minimize the regularized score function (with thesame regularization term as in Eq. 16) to estimate theinvolved parameter.

6

Experiments

Simulations The simulated data are generated byapplying the biased selection procedure on the data
generated by a additive noise model with function
f AN , i.e., by (9) and (10). As shown in Section 4,
the function f AN is identifiable up to some shift when

1

β(y)

0.8
0.6
β (1)(y)β (2)(y)β (3)(y)

0.4
0.2
0
-4

-2

0
y

2

AN
fˆML
vs. fˆGPfˆAN vs. fˆGP
SM

10

-1

10 -2

10 -3

4

10 -2
M SE(fˆGP , f AN )

10 -1

AN
fˆML
vs. fˆGPfˆAN vs. fˆGP
SM

10

-1

10 -2

10 -3

10 -4
10 -4

10

-3

10

-2

M SE(fˆGP , f AN )

10

-1

4

2

2

f AN

10

10 -1

-2
AN
fˆML
vs. fˆGPANfˆSM
vs. fˆGP

-1

0

AN
fˆML

fAN
0

-2
1

10 -3

fˆGP

Selected points

-1

0

X

10 -2

10 -4
10 -4

0

Generated by ANM

0log(βy )0ANAN
M SE(fˆML, f AN ) (’+’) & M SE(fˆSM, f AN ) (’O’)

ANAN
M SE(fˆML, f AN ) (’+’) & M SE(fˆSM, f AN ) (’O’)(a) (y) used in simulations. (b) Results in settings 1⇠4.
10

4
10 -3

Y

1.2

10 0

ing method. The estimation performance is evaluatedby the mean square error (MSE) between the estimate
Pn
and the true function, n1 i=1 (fˆAN (xi ) f AN (xi ))2 .
We compare the estimates produced by our methods
with that estimated by Gaussian process (GP) regression on the given sample, denoted by fˆGP . Figure 3(bANˆANd) compares the estimation quality of fˆM
L and fSM
ˆagainst fGP ; note that they are plotted on a log scale.

Y

AN
AN
M SE(fˆML, f AN ) (’+’) & M SE(fˆSM, f AN ) (’O’)

the noise is non-Gaussian. We shall study the estimation quality of the regression function f AN underdi↵erent settings.

1

X

6

2

4

0

2

-2

0

-4

log(p̂X )log(p̂E )

10

-3

10

-2

10

-1

M SE(fˆGP , f AN )(c) Results in settings 5⇠8. (d) Results in settings 9⇠12.

Figure 3: Simulation settings and results. (a) shows
the density ratio functions for OSB, 1 (y), 2 (y), and
3 (y), which are used in settings 1-4, 5-8, and 9-12,
respectively. (b), (c), and (d) show the pairwise MSEof the estimated function for the proposed methods
against GP regression on the given sample, in settings
1-4, 5-8, and 9-12, respectively. The dashed line marks
the threshold where the proposed methods and GP
regression on the given sample perform equally well.
We consider three settings for OSB, by setting (y)(see Eq. 10) to di↵erent functions. As shown in Figure 3(a), (1) (y) is a piecewise linear function, (2) (y)is a (scaled) Gaussian function with mean -1 and standard deviation 2, and (3) (y) corresponding to a hardbiased selection procedure: it drops all data points
corresponding to the 10% largest values of Y .
We use two ways to generate the distributions of X andE; one is the uniform distribution, and the other the
mixture of three Gaussians with random coefficients.
The function f AN is a mixture of the linear, tanh, andcubic function with random coefficients (the coefficientfor the cubic function is constrained to be small to
avoid extreme values in Y ).
In total there are 2 ⇥ 2 ⇥ 3 simulation settings. For
each setting we repeat the experiment with 15 random
replications. We use the methods proposed in SecAN
tion 5 to recover this function. Denote by fˆM
L the estimate given by the (approximate) maximum likelihood
ANmethod, and by fˆSMthat given by the score match-

-2
-5

0

5

-6
-2

-1

0

1

y

Figure 4: Results of a typical run estimated by the
maximum likelihood approach. The four subfigures
show the data produced by the ANM, the selected
sample and the estimates of f AN , the estimate of (y),and the estimates of pX and pE , respectively.
As one can see from Figure 3(b-d), the proposed methods may converge to unwanted solutions, as shown by
the few points above the dashed lines. However, inmost cases the proposed method provides a better estimate of the function f AN . As suggested by (Demšar,2006), we use the Wilcoxon signed ranks test to check
whether M SE(fˆAN , f AN ) is significantly better thanM SE(fˆGP , f AN ) under all the three settings for (y).
It is a nonparametric test to detect shifts in populations given a number of paired samples. Under the null
hypothesis the distribution of di↵erences between thetwo populations is symmetric about 0. We find that
under all the three sets of settings, for the score matching method, the null hypothesis is always rejected at
the 0.01 level (the p-values are 6 ⇥ 10 7 , 2 ⇥ 10 6 , and8 ⇥ 10 5 , respectively). For the maximum likelihood
method, the null is rejected in settings 1-4 and 5-9 (thep-values are 2 ⇥ 10 7 , and 8 ⇥ 10 5 , respectively); we
fail to reject the null in settings 5-8 (with the p-value0.22): this seems to be caused by local optima, which
will be discussed below. This means that the proposed
method outperforms the method that fits GP regression on the observed data, in terms of the estimation
quality of the true function.
Figure 4 gives the result of a typical run (with 3 (y))
produced by maximum likelihood. Interestingly, one

2

1.8

1

1.6

0

β̂(y)

Y (after PNL correction)

can see that compared to the true (t), which is 3 (y)in Figure 3(a), ˆ(y) contains an additional factor of
the form ec2 y with some constant c2 . The estimates ofpX and pE are also skewed accordingly. This verifies
the statements given in Theorem 2(a).

-1
-2
-3

1.4
1.2

Data
fˆGP
fˆAN

1

2

0.8
-4

ML

As seen from Figure 3, both algorithms may get stuck
in local optima. Let us have a closer look at the results given by maximum likelihood. We found thatfor each simulation setting, in all runs where the function f AN was not accurately recovered, or more specifAN
AN
ically, where M SE(fˆM) > M SE(fˆGP , f AN ),L, fthe corresponding likelihood values are among the lowest across all 15 runs. That is, the attained likelihood
value suggests whether the algorithm converges to a
local optimum. Therefore, in practice one may run
the algorithms multiple times with random initializations and choose the one which gives the best modelfit (e.g., the highest likelihood). However, this is not
the case for the score matching-based approach: we
did not find that the unwanted solutions always correspond to relatively large score distances. The precise
reason for this phenomenon is still under investigation.
Hence, below we only report the results given by the(approximate) maximum likelihood approach.
Experiments on Real Data We went through thecause-e↵ect pairs (http://webdav.tuebingen.mpg.
de/cause-effect/) to find data sets which are likely
to su↵er the OSB issue according to commonsense orbackground knowledge. We selected Pairs 25, 40, and
41. Here to save space, we only report the results on
Pair 25; it is about the relationship between the age(X) and the concrete compressive strength (Y ) of different samples of concrete.
The empirical distribution of the data in Pair 25 suggests that it is very likely for the e↵ect to su↵er froma PNL distortion. We use a rough way to take into
account both the PNL distortion in the causal process
and the OSB. We first fit the PNL causal model (Zhang& Hyvärinen, 2009) on the data and correct the datawith the estimated PNL transformation on Y . We
then fit the ANM-OSB procedure on the corrected
data. To avoid local optima, we run the (approximate) maximum likelihood algorithm presented in Section 5.1 five times with random initializations andchoose the one with the highest likelihood. Figure 7
shows the result on Pair 25. As seen from ˆ(y), it
seems for some reason, the samples whose compressive
strength is very high were not selected. The estimatedGP
function fˆM
L seems to address this issue.

7

Conclusion and Discussions

As we have shown, in the presence of outcomedependent selection, the causal direction is still generi-

-4
-2

-1

0

1

3

-2

0

X

y

(a) Data & estimated functions.

(b) ˆ(y).

2Figure 5: Results on pair 25 of the cause-e↵ect pairs.
(a) The scatterplot of the data (after correcting thenonlinear distortion in Y with the PNL causal model),the nonlinear regression function fˆGP on the data, andANthe estimated function fˆM
L by the maximum likelihood approach. (b) The estimated density ratio (y).
cally identifiable if the causal relationship can be modeled by a post-nonlinear causal model. Moreover, inthe case of an additive noise model, the causal mechanism as represented by the function in the modelis identifiable up to a constant if the noise term is
non-Gaussian (and completely identified if the noise
term follows a symmetric, non-Gaussian distribution).
However, due to the selection bias, the estimation requires more care than standard methods for fitting
such models, and we developed two estimation procedures in this paper.
This is a first step towards a better understanding of
the bearing of selection bias on the identifiability andestimation of functional causal models. There are several interesting problems for future work. First, theidentifiability result on additive noise models can begeneralized to post-nonlinear models, but it will takemore work to put the more general result in a sufficiently simple form. Second, our positive results here
are confined to outcome-dependent selection. Thanks
to this restriction, our results do not rely on any substantial assumption on the selection mechanism. For
more complex structures of selection, such as when
the selection depends on both cause and e↵ect, identifiability will require more specifications of the selection model. Third, our result on the identification of
causal direction is confined to the two-variable casewithout latent confounders; how to handle selection
in multi-variable structural learning, with or without
latent confounders, remains an open problem.
Acknowledgement: Research reported in this publication was supported by the National Institutes of
Health (NIH) under Award Number U54HG008540.
The content is solely the responsibility of the authorsand does not necessarily represent the official views of
the NIH. The research of J. Zhang was supported inpart by the Research Grants Council of Hong Kong
under the General Research Fund LU342213.

References
Bareinboim, E. and Pearl, J. Controlling selection bias
in causal inference. In Proceedings of AAAI, pp.
100–108, 2012.
Bareinboim, E., Tian, J., and Pearl, J. Recoveringfrom selection bias in causal and statistical inference. In Proceedings of AAAI, 2014.
Borboudakis, G. and Tsamardinos, I. Bayesian network learning with discrete case-control data. In
Proceedings of the 31st Conference on Uncertainty
in Artificial Intelligence, 2015.
Castillo, E. Functional Equations and Modelling in
Science and Engineering. CRC Press, 1992.
Chickering, D. M. Optimal structure identification
with greedy search. Journal of Machine Learning
Research, 3:507–554, 2002.
Demšar, J. Statistical comparisons of classifiers overmultiple data sets. Journal of Machine Learning
Research, 7:1–30, 2006.
Didelez, V., Kreiner, S., and Keiding, N. Graphicalmodels for inference under outcome-dependent sampling. Statistical Science, 25:368–387, 2010.
Evans, R. J. and Didelez, V. Recovering from selection
bias using marginal structure in discrete models. In
Proceedings of the 31st Conference on Uncertainty
in Artificial Intelligence, 2015.
Gretton, A., Borgwardt, K., Rasch, M., Schölkopf, B.,
and Smola, A. A kernel method for the two-sampleproblem. In NIPS 19, pp. 513–520, Cambridge, MA,
2007. MIT Press.
Heckerman, D., Geiger, D., and Chickering, D. M.
Learning bayesian networks: The combination ofknowledge and statistical data. Machine Learning,
20:197–243, 1995.
Heckman, J. Sample selection bias as a specification
error. Econometrica, 47:153–161, 1979.
Hoyer, P.O., Janzing, D., Mooji, J., Peters, J., andSchölkopf, B. Nonlinear causal discovery with additive noise models. In Advances in Neural Information Processing Systems 21, Vancouver, B.C.,
Canada, 2009.
Hyvärinen, A. Estimation of non-normalized statistical models by score matching. Journal of Machine
Learning Research, 6:695–708, 2005.
Kagan, A. M., Linnik, Y. V., and Rao, C. R. Characterization Problems in Mathematical Statistics. Wiley, New York, 1973.
Pearl, J. Causality: Models, Reasoning, and Inference.
Cambridge University Press, Cambridge, 2000.
Shimizu, S., Hoyer, P.O., Hyvärinen, A., and Kerminen, A.J. A linear non-Gaussian acyclic model forcausal discovery. Journal of Machine Learning Research, 7:2003–2030, 2006.
Spirtes, P., Meek, C., and Richardson, T. An algo-rithm for causal inference in the presence of latentvariables and selection bias. In Glymour, C. andCooper, G. (eds.), Computation, Causation, andDiscovery, pp. 211–252. MIT Press, 1999.
Spirtes, P., Glymour, C., and Scheines, R. Causation,
Prediction, and Search. MIT Press, Cambridge,
MA, 2nd edition, 2001.
Winship, C. and Mare, R. D. Models for sample selection bias. Annual Review of Sociology, 18:327–350,
1992.
Zhang, J. On the completeness of orientation rules
for causal discovery in the presence of latent confounders and selection bias. Artificial Intelligence,
172:1873–1896, 2008.
Zhang, K. and Hyvärinen, A. On the identifiabilityof the post-nonlinear causal model. In Proceedingsof the 25th Conference on Uncertainty in Artificial
Intelligence, Montreal, Canada, 2009.