Identifying confounders using additive noise models

Dominik Janzing
MPI for Biol. Cybernetics
Spemannstr. 38
72076 Tübingen
Germany

Jonas Peters
MPI for Biol. Cybernetics
Spemannstr. 38
72076 Tübingen
Germany

Abstract
We propose a method for inferring the existence of a latent common cause (“confounder”) of two observed random variables.
The method assumes that the two effects ofthe confounder are (possibly nonlinear) functions of the confounder plus independent, additive noise. We discuss under which conditions the model is identifiable (up to an arbitrary reparameterization of the confounder)from the joint distribution of the effects. We
state and prove a theoretical result that provides evidence for the conjecture that themodel is generically identifiable under suitable technical conditions. In addition, we
propose a practical method to estimate theconfounder from a finite i.i.d. sample of theeffects and illustrate that the method works
well on both simulated and real-world data.

1

Introduction

Since the pioneering work on causal inference methods(described for example in [Pearl, 2000] and [Spirtes
et al., 1993]), much work has been done under theassumption that all relevant variables have been observed. An interesting, and possibly even more important, question is how to proceed if not all the relevant
variables have been observed. In that case, dependencies between observed variables may also be explained by confounders — for instance, if a dependence
between the incidence of storks and the birthrate istraced back to a common cause influencing both variables. In general, the difficulty not only lies in the factthat the values of the latent variables have not been observed, but also that the causal structure is unknown.
In other words, it is in general not clear whether andhow many confounders are needed to explain the data
and which observed variables are directly caused by

Joris Mooij
MPI for Biol. Cybernetics
Spemannstr. 38
72076 Tübingen
Germany

Bernhard Schölkopf
MPI for Biol. Cybernetics
Spemannstr. 38
72076 Tübingen
Germany

which confounder. Under the assumption of linear
relationships between variables, confounders may be
identified by means of Independent Component Analysis, as shown recently by Hoyer et al. [2008], if the
distributions are non-Gaussian. Other results for the
linear case are presented in Silva et al. [2006]. In thispaper, we will not assume linear relationships, but try
to tackle the more general, nonlinear case.
In the case of two variables without confounder, Hoyer
et al. [2009] have argued that the causal inference task
(surprisingly) becomes easier in case of nonlinear functional relationships. They have described a method toinfer whether X → Y (“X causes Y ”) or Y → X fromthe joint distribution P (X, Y ) of two real-valued random variables X and Y . They consider models where
Y is a function f of X up to an additive noise term,i.e.,
Y = f (X) + N ,(1)
where N is an unobserved noise variable that is statistically independent of X. They show in their paperthat generic choices of functions f , distributions of X
and distributions of N induce joint distributions on
X, Y that do not admit such an additive noise model inthe inverse direction, i.e., from Y to X. If P (X, Y ) admits an additive model in one and only one direction,
this direction can be interpreted as the true causal
direction. We believe that the situation with a confounder between the two variables is similar in that
respect: nonlinear functional relationships enlarge theclass of models for which the causal structure is identifiable.
We now state explicitly which assumptions we make inthe rest of this paper. First of all, we focus on the caseof only two observed and one latent continuous randomvariables, all with values in R. We assume that there
is no feedback, or in other words, the true causal structure is described by a DAG (directed acyclic graph).
Also, we assume that selection effects are absent, thatis, the data samples are drawn i.i.d. from the probability distribution described by the model.

TXu(t), v(t)

y

Y



x

Figure 1: Directed acyclic graph and a scatter plotcorresponding to a CAN model for two observed variables X and Y that are influenced by an unobserved
confounder variable T .
Definition 1 Let X, Y and T be random variables
taking values in R. We define a model for Confounders with Additive Noise (CAN) by
X = u(T ) + NX
Y = v(T ) + NY

with NX , NY , T
jointly independent.

(2)

where u, v : R → R are continuously differentiable
functions and NX , NY are real-valued random variables.
The random variables NX and NY describe additive“noise”, of which one may think of as the net effectof all other causes which are not shared by X and Y .
This model can be represented graphically by the DAG
shown in Figure 1.
Definition 2 We call two CAN models equivalent if
they induce the same distributions of NX , NY and thesame joint distribution of (u(T ), v(T )).
This definition removes the ambiguity arising from unobservable reparameterizations of T . We further adopt
the convention E(NX ) = E(NY ) = 0.
The method we propose here enables one to distinguish between (i) X → Y , (ii) Y → X, and (iii)X ← T → Y for the class of models defined in (2),and (iv) to detect that no CAN model fits the data(which includes, for example, generic instances of thecase that X causes Y and in addition T confounds
both X and Y ). If NX = 0 a.s. (“almost surely”) andu is invertible, the model reduces to the model in (1)by setting f := v ◦ u−1 . Given that we have observed ajoint density on R2 of two variables X, Y that admits aunique CAN model, the method we propose identifies
this CAN model and therefore enables us to do causal
inference by employing the following decision rule: weinfer X → Y whenever NX is zero a.s. and u invertible,
infer Y → X whenever NY is zero a.s. and v invertible,
and infer X ← T → Y if neither of the alternatives
hold, which corresponds in spirit with Reichenbach’s
principle of common cause [Reichenbach, 1956]. Notethat the case of NX = NY = 0 a.s. and u and v invertible implies a deterministic relation between X and Y ,which we exclude here.

In practical applications, however, we propose to prefer the causal hypothesis X → Y already if the variance of NX is small compared to the variance of NY(after we have normalized both X and Y to variance1). To justify this, consider the case that X causes Y
and the joint distribution admits a model (1), but by aslight measurement error, we observe X̃ instead of X,
differing by a small additive noise term. Then P (X̃, Y )admits a proper CAN model because X is the latentcommon cause, but we infer X̃ → Y because, from acoarse-grained point of view, we should not distinguishbetween the quantity X and the measurement resultX̃ if both variables almost coincide.
Finding the precise conditions under which the identification of CAN models is unique up to equivalence,is a non-trivial problem: If u and v are linear and
NX , NY , T are Gaussian, one obtains a whole familyof models inducing the same bivariate Gaussian joint
distribution. Other examples where the model is notuniquely identifiable are given in Hoyer et al. [2009]:
any joint density which admits additive noise models
from X to Y and also from Y to X is a special case ofa non-identifiable CAN model.
The remainder of this paper is organized as follows. In
the next section, we provide theoretical motivation forour belief that in the generic case, CAN models are
uniquely identifiable. A practical algorithm for thetask is proposed in Section 3. It builds on a combination of nonlinear dimensionality reduction and kernel
dependence measures. Section 4 provides empirical results on synthetic and real world data, and Section 5
concludes the paper.

2

Theoretical motivation

Below, we prove a partial identifiability result for the
special case that both u, v are invertible, where we
consider the following limit: first, let the variances of
the noise terms NX and NY be small compared to thecurvature of the graph (u(t), v(t)); second, we assume
that the curvature is non-vanishing nevertheless (ruling out the linear case), and third, that the density on
the graph (u(t), v(t)) changes slowly compared to thevariance of the noise.
Because of the assumption that u is invertible, we canreparameterize T such that the CAN model (2) simplifies to
X = T + NX and Y = v(T ) + NY ,i.e., the joint probability density is given byp(x, y, t) = q(x − t)r(y − v(t))s(t)where q, r, s are the densities of NX , NY and T respectively. We will further assume that s is differentiable

y
y2
y1
y0

*
*
* ** v(x)
* *
* **
* **
** * *** ** **** * *** *** * * ** * ** * * *
*
***
w(y0)

because NX and NY are centralized. We can thus approximately determine the function v from observing
all the conditional expectations Ey (X). After a shiftby the estimated function w(y), the conditional distribution of X, given y, is approximately given by aconvolution of NX with the scaled noise βy NY :
X − Ey (X) | y ∼ NX + βy NY
xFigure 2: Sample from a distribution p(x, y) obtainedby adding small noise in x- and y-direction to a distribution supported by a graph (t, v(t)). We show thatfor fixed y, X is approximately an affine function in
the noise variables NX and NY .
with bounded derivative and is nonzero at an infinite
number of points. We also assume that the inverse
function w := v −1 of v is two times differentiable andw00 is bounded. We further assume that the density of
NX is several times differentiable, that all its moments
are finite and that E(|NX |k )/k! → 0. This implies
that the characteristic function of NX , and therefore
its distribution, can be uniquely expressed in terms of
its moments. We assume likewise for NY .
We will show that under these assumptions, one canestimate the function v from the conditional expectation E(X|Y = y) and that all moments of NX andNY can be approximated from higher order conditional
moments E(X k |Y = y), with vanishing approximation
error in the limit.
Defining ry (t) := p(y, t) = r(y − v(t))s(t), the conditional distribution of T given Y = y is given byry (t)/p(y). For fixed y, we can locally (i.e., for t ≈t0 := w(y)) approximatery (t) ≈ r − v 0 (t0 )(t − t0 ) s(t0 )if s is almost constant around t0 and if v has small
curvature at t0 . Hence,
T | y ∼ w(y) + βy NY

We cannot observe T directly, but we can observe thenoisy version X = T + NX of it (see also Figure 2).
From (3), we conclude that the conditional distribution of X given Y = y is approximately given by(4)

Using Ey (. . . ) as shorthand for the conditional expectation E(. . . | Y = y), we conclude that
Ey (X) ≈ w(y) ,Lemma 1 Let Z and W be independent random variables for which all moments E(Z k ) and E(W k ) fork = 1, . . . , n exist. Then all these moments can be reconstructed by observing the nth moments of Z + cj W
for n + 1 different values c0 , . . . , cn .
Proof: The nth moments of Z + cj W are given by
 
n
 X
n
k n
E (Z + cj W ) =cj
E(Z n−k )E(W k ) .
kk=0

The (n + 1) ×
 (n + 1) matrix M given by the entries
mjk := ckj nk is invertible because it coincides with
the Vandermonde matrix up to multiplying the kth
column with nk . Hence we can compute all productsE(Z n−k )E(W k ) for k = 0, . . . , n by matrix inversion
and obtain in particular E(Z n ) and E(W n ). Taking
into account only a subset of points, we can obviously
compute lower moments, too. 
We now choose n + 1 values y0 , . . . , yn such that we
obtain n + 1 different values βy0 , . . . , βyn and apply
Lemma 1 to identify the first n moments of NX and
NY . For higher moments, the approximation will typically cause larger errors.
We now focus on the error made by the approximations
above. We introduce the error termsn (y) := Ey (T − w(y))n − E(βyn NYn ) ,and hence obtain(3)

1
approximately, where βy := − v0 (t= −w0 (y).
0)

X | y ∼ w(y) + βy NY + NX

Since we can observe this convolution for different values βy , we can compute the moments of NX and NY
using the following Lemma:(5)

Ey (X) = Ey (T ) = w(y) + 1 (y) .
Some calculations then yield the following exact relation between the moments of NX and NY :
 
n
X
n
n−k
βyk
E(NX
)E(NYk )(6)
kk=0= Ey ((X − Ey (X))n )
 
n
X
k n
+
1 (y)Ey ((X − Ey (X))n−k )
k
k=1
 
n
X
n
n−k
−
k (y)
E(NX
)
k
k=0

(7)
(8)
(9)

Defining vectors b, c and d by letting their j’th component be the value of (7)–(9) for y = yj , respectively,
and defining the vector q and the matrix M by
 
n
n−k
k
k
qk := E(NX )E(NY ),Mjk := βyj,kwe can write (6)–(9) for y = y0 , . . . , yn in matrix notation as
Mq = b + c + dand therefore
q = M−1 b + M−1 (c + d).
−1

The above approximation yielded q = M b.
The remaining terms are bounded from above by
M−1 (kck + kdk ). The following lemma (which
is proved in the Appendix) shows that the errors
k (y) are small. This then implies that the error inq ≈ M−1 b is also of order O(ks0 k∞ + kw00 k∞ ).
Lemma 2 For any n and y such that s(w(y)) 6= 0:
n (y) = O(ks0 k∞ + kw00 k∞ ),where the terms involve moments of NY and positive
powers of βy and of s(w(y)), which are all bounded forgiven y and n.
We consider now a sequence of distributions p` (x, y)obtained by scaling the graph (t, v(t)) up to the larger
graph (`t, `v(t)), while keeping the distributions of thenoise NX and NY fixed. Then we consider the conditional distributions of X, given Y = y for y =`y0 , . . . , `yn and can determine the moments of NX
and NY up to an error that converges to zero for` → ∞, as the following Theorem shows.
Theorem 1 Define a sequence of joint densitiesp` (x, y) byZp` (x, y) := q(x − `t)r(y − `v(t))s(t)dt .
(10)
Let, as above, y0 , . . . , yn be chosen such that allkβy0 , . . . , βyn are different. Then every moment E(NX)
kand E(NY ) can be computed from the conditional moments E`yj (X m ) for j = 0, . . . , n and m = 1, . . . , n
up to an error that vanishes for ` → ∞ under the assumptions made above.
Proof: We rewrite eq. (10) asZp` (x, y) := q(x − t̃` )r(y − ṽ` (t̃` ))s̃` (t̃` )dt̃` .
witht̃` := `t,ṽ` (t̃` ) := `v(t̃` /`),s̃` (t̃` ) := s(t̃` /`)/` .

w̃` is then defined as ṽ`−1 , and one checks easily thatks̃0` k∞ and kw̃`00 k∞ tend to zero with O(1/`2 ). Further,note that the βyj are invariant with respect to thescaling: β̃yj ,` = βyj . Hence, by Lemma 2, all k (yj )converge to zero.
To summarize, we have sketched a proof that the
CAN model becomes identifiable in a particular limit.
We expect that a stronger statement holds (undersuitable technical conditions), but postpone the nontrivial problem of finding the right technical conditions
and the corresponding proof to future work.
The analysis above also shows that it should be possible to identify a confounder by estimating the variances of the noises NX and NY and comparing theirsizes (as discussed in the previous section). Indeed, inorder to estimate the variance of the noise variable NX ,one observes the conditional expectation E(X | Y = y)for three different values of y; if the conditional expectations are sufficiently different, one can apply Lemma
1 to calculate the moments of NX and NY up to second order. Assume, we observe a dependence between
X and Y (which are normalized to have variance 1).
In case one of the noise variances is much smaller thanthe other (say NX  NY ), this would indicate a direct causal influence (X → Y in that case), but if bothnoise variances are large, this indicates the existenceof a confounder.

3

Identifying Confounders from Data

In this section we propose an algorithm (ICAN) thatis able to identify a confounder in CAN models. While
we only addressed the low noise regime in the previous
theoretical section, the practical method we propose
here should work even for strong noise, although inthat case more data points are needed.
Assume that X, Y are distributed according to the
CAN model (2). We write s(t) := u(t), v(t) for
the “true” curve of the confounder in R2 . A scatterplot of the samples (X, Y ) (right panel of Figure 1,
for example) suggests a simplistic method for detecting the confounder: for every curve s : [0, 1] → R2
project the data points (Xk , Yk ) onto this curve s,such that the Euclidean distance is minimized: T̂k =argmint∈[0,1] k(Xk , Yk )−s(t)k2 . From a set of all possible paths S P
now choose the ŝ that minimizes the globaln`2 distance k=1 k(Xk , Yk ) − s(T̂k )k2 (dimensionality
reduction) and propose T̂k to be the confounder for
Xk and Yk . This results in the estimated residuals(N̂X,k , N̂Y,k ) = (Xk , Yk ) − ŝ(T̂k ). If the hypothesesT̂ ⊥⊥ N̂X , T̂ ⊥⊥ N̂Y , N̂X ⊥⊥ N̂Y cannot be rejected, propose that there is a confounder whose values are given
by T̂k .

This idea turns out to be too naive: even if the data
have been generated according to the model (2), theprocedure results in dependent residuals. As an example, consider a data set simulated from the following
model:
X = 4 · ϕ−0.1 (T ) + 4 · ϕ1.1 (T ) + NX
Y = 1 · ϕ−0.1 (T ) − 1 · ϕ1.1 (T ) + NY
where ϕµ is the probability density of a
N (µ, 0.12 ) distributed random variable and
NX , NY ∼ U ([−0.1, 0.1]) and T ∼ U ([0, 1]) are
jointly independent. We now minimize 
the global `2distance over the set of functions S = s : si (t) =αi · ϕ−0.1 (t) + βi · ϕ1.1 (t); i = 1, 2 . Since there are
only four parameters to fit, the problem is numericallysolvable and gives the following optimal solution:
α1 = 3.9216, β1 = 4.0112, α2 = 0.9776, β2 = −0.9911.
The `2 projections T̂ result in a lower global `2distance (6.92) than the true values of T (11.87).

estimated NX

1
0.8
0.6

0

−0.1

0.4

estimated NY

0.2
Y

0.1

0

0

0.5
estimated T

1

0

0.5
estimated T

1

0.1
0

−0.1

−0.2

estimated NY

−0.4
−0.6
−0.8

0.1
0

−0.1

−1

3.4

3.8
X

4.2

0
0.1
estimated NX

3. Estimated residuals and confounder values clearly
depend on each other. Also independence tests likethe Hilbert-Schmidt Independence Criterion (see below) reject the hypotheses of independence: p-valuesof 5 × 10−2 , 7 × 10−5 and 1 × 10−4 are computed corresponding to the right plots in Figure 3 from top tobottom. These dependencies often occur when the
projections onto the curve ŝ are chosen to minimize
the global `2 distance, which can be seen as follows:1
in our example ∂s∂t is small for T ≈ 0.5 or Y ≈ 0.
Since the points are projected onto the curve orthogonally, the projection results in very small residuals N̂Y .
This introduces a dependency between N̂Y and T̂ . Dependency between the residuals N̂Y and N̂X can arise
from regions, where ŝ is approximately linear, like inthe bottom right part of Figure 3: positive residuals
NY lead to positive residuals NX and vice versa.
Summarizing, projecting the pairs (X, Y ) onto thepath (ŝ(t)) by minimizing the `2 distance to the path is
the wrong approach for our purpose. Instead, the data
(X, Y ) should be projected in a way that minimizes thedependence between residuals and confounder (N̂X , T̂and N̂Y , T̂ ) and between the residuals itself (N̂X , N̂Y ).
Let DEP(W, Z) denote any non-negative dependence
measure between random variables W and Z, which iszero if and only if W and Z are independent (we latersuggest to use the Hilbert-Schmidt Independence Criterion). In the example above we can use the red curveas an initial guess, but choosing the projections byminimizing the sum of the three dependence measuresinstead of `2 distances. In our example this indeed
leads to residuals that fulfill the independence constraints (p-values of 1.00, 0.65, 0.80). For the generalcase, we propose Algorithm 1 as a method for identifying the hidden confounder T given an i.i.d. sampleof (X, Y ).
If a CAN model can be found we interpret the outcomeN̂X
of our algorithm as X → Y if Var
 1 and û inVarN̂
Y

Figure 3: Left: a scatter plot of the data, true paths and projections (black and solid), estimated path ŝand projections (red and dashed). Right: residuals
plotted against each other and estimated confounder.
Figure 3 shows the true function s (black line), a scatter plot of (X, Y ) (black circles) and the computed
curve ŝ that minimizes the global `2 distance (dashed
red line). Additionally, for some data points projections onto s and ŝ are shown, too: black crosses correspond to the “true projections” (i.e., the points without noise) onto s and red crosses correspond to projections onto the estimated function ŝ minimizing the`2 distance. The latter result in the proposed residuals, which are shown together with the estimatedvalues of the confounder on the right side of FigureN̂X
vertible and as Y → X if Var
 1 and v̂ invertible.
VarN̂Y
There is no mathematical rule that tells whether one
should identify a variable X and its (possibly noisy)
measurement X̃ or consider them as separate variables
instead. Thus we cannot be more explicit about thethreshold for the factor between VarN̂X and VarN̂Y
that tells us when to accept X → Y or Y → X or
X ←T →Y.

To implement the method we still need an algorithm
for the initial dimensionality reduction, a dependence
criterion DEP, a way to minimize it and an algorithm for non-linear regression. Surely, many different
choices are possible. We will now briefly justify our
choices for the implementation.

Algorithm 1 Identifying Confounders using Additive
Noise Models (ICAN)
1: Input: (X1 , Y1 ), . . . , (Xn , Yn ) (normalized)
2: Initialization:
3: Fit a curve ŝ to the data
P that minimizes `2 dis-

tance: ŝ := argmins∈Sn
k=1 dist s, (Xk , Yk ) .

4: repeat
5:
Projection:
6:T̂ := argminT DEP(N̂X , N̂Y ) + DEP(N̂X , T ) +

7:
8:
9:
10:

DEP(N̂Y , T ) with (N̂X,k , N̂Y,k ) = (Xk , Yk ) −ŝ(Tk )
if N̂X ⊥⊥ N̂Y and N̂X ⊥⊥ T̂ and N̂Y ⊥⊥ T̂ thenOutput: (T̂1 , . . . , T̂n ), û = ŝ1 , v̂ = ŝ2 , andVarN̂X.
VarN̂Y
Break.
end ifRegression:
Estimate ŝ by regression (X, Y ) = ŝ(T̂ )+ N̂. Setû = ŝ1 , v̂ = ŝ2 .
13: until K iterations
11:
12:

14: Output: Data cannot be fitted by a CAN model.

they should be able both to deal with continuous data
and to detect non-linear dependencies. Since there
is no canonical way of discretizing continuous variables, methods that work for discrete data (like a χ2test) are not suitable for our purpose. In our method
we chose the Hilbert-Schmidt Independence Criterion(HSIC) [Gretton et al., 2005, 2008]. It can be defined as the distance between the joint distribution
and the product of the marginal distribution represented in a Reproducing Kernel Hilbert Space. For
specific choices of the kernel (e.g., a Gaussian kernel)
it has been shown that HSIC is zero if and only if
the two distributions are independent. Furthermore
the distribution of HSIC under the hypothesis of independence can be approximated by a Gamma distribution [Kankainen, 1995]. Thus we can constructa statistical test for the null hypothesis of independence. In our experiments we used Gaussian kernels
and chose their kernel sizes to be the median distances
between the points [Schölkopf and Smola, 2002]: e.g.
2σ 2 = median{kXk − Xl k2 : k < l}. We will use theterm HSIC for the value of the Hilbert-Schmidt normand pHSIC for the corresponding p-value. For a small
p-value (< 0.05, say) the hypothesis of independence
is rejected.

Initial Dimensionality Reduction
It is difficult to solve the optimization problem (line 3of the algorithm) for a big function class S. Our approach thus separates the problem into two parts: westart with an initial guess for the projection values T̂k(this is chosen using an implementation of the Isomap
algorithm [Tenenbaum et al., 2000] by van der Maaten[2007]) and then iterate between two steps: In the firststep we keep the projection values T̂k fixed and choose
a new function ŝ = (û, v̂), where û and v̂ are chosen
by regression from X on T̂ and Y on T̂ , respectively.
To this end we used Gaussian Process Regression [Rasmussen and Williams, 2006], using the implementationof Rasmussen and Williams [2007], with hyperparameters set by maximizing marginal likelihood. In the second step the curve is fixed and each data point (Xk , Yk )is projected to the nearest point of the curve: Tk is
chosen such that kŝ(Tk ) − (Xk , Yk )k`2 is minimized.
We then perform the first step again. A similar iterative procedure for dimensionality reduction has been
proposed by Hastie and Stuetzle [1989].

For the projection step we now minimize
HSIC(N̂X , N̂Y ) + HSIC(N̂X , T̂ ) + HSIC(N̂Y , T̂ )with respect to T̂ . Note that at this part of thealgorithm the function ŝ (and thus û and v̂) remain
fixed and the residuals are computed according to
NX = X − û(T̂ ) and NY = Y − v̂(T̂ ). We useda standard optimization algorithm for this task
(fminsearch in MatLab) initializing it with thevalues of T̂ obtained in the previous iteration. Insteadof the sum of the three dependence criterion the
maximum can be used, too. This is theoreticallypossible, but complicates the optimization problem
since it introduces non-differentiability.

This initial step of the algorithm is used for stabilization. Although the true curve s may differ from the`2 minimizer ŝ, the difference is not expected to bevery large. Minimizing dependence criteria from the
beginning often results in very bad fits.

Non-linear Regression
Here, again, we used Gaussian process regression forboth variables separately. Since the confounder values
T̂ are fixed we can fit X = û(T̂ ) + N̂X and Y = v̂(T̂ ) +N̂Y to obtain ŝ = (û, v̂).

Dependence Criterion and its Minimization
There are various choices of dependence criteria thatcan be used for the algorithm. Notice, however, thatIn the experiments this step was mostly not necessary:
whenever the algorithm was able to find a solution

It should be mentioned that sometimes (not for alldata sets though) a regularization for the T values
may be needed. Even for dependent noise very positive (or negative) values of T result in large residuals,
which may be regarded as independent. In our implementation we used a heuristic and just performed
5000 iterations of the minimization algorithm, which
proved to work well in practice.

In this section we show that our method is able todetect confounders both in simulated and in real data.
4.1

Simulated data

Data set 1.
We show on a simulated data set that our algorithm
finds the confounder if the data comes from the model
assumed in (2). We simulated 200 data points froma curve whose components u and v consist of a random linear combination of Gaussian bumps each. The
noise is uniformly distributed on [−0.035, 0.035]. Notethat in contrast to the example given in Section 3
we are now doing the regression using Gaussian processes. The algorithm finds a curve and corresponding
projections of the data points onto this curve, suchthat N̂X , N̂Y and T̂ are pairwise independent, which
can be seen from the p-values pHSIC (N̂X , N̂Y ) = 0.94,
pHSIC (N̂X , T̂ ) = 0.78 and pHSIC (N̂Y , T̂ ) = 0.23. The
top panel of Figure 4 shows the data and both true
(black) and estimated (red) curve. The bottom panel
shows estimated and true values of the confounder.
Recall that the confounder can be estimated only up
to an arbitrary reparameterization (e.g. t 7→ −t).

Data set 2.
This data set is produced in the same way as data
set 1, but this time using an invertible v and unequal
scaled noises. We sampled NX ∼ U ([−0.008, 0.008])
and NY ∼ U ([−0.0015, 0]). We argued above thatfor finite sample sizes this case should rather be regarded as Y → X and not as an example with ahidden common cause. The algorithm again identifies a curve and projections, such that the independence constraints are satisfied (pHSIC (N̂X , N̂Y ) = 1.00,
pHSIC (N̂X , T̂ ) = 1.00 and pHSIC (N̂Y , T̂ ) = 1.00, see
Figure 5), and it is important to note that the different
scales of the variances are reduced, but still noticeableN̂X )≈ 5). In such a case we indeed interpret the( Var(
Var(N̂Y )outcome of our algorithm as “Y causes X”.
0.5
Y

Experiments

0−0.5

estimated NX

4ship between X and Y . It is obvious that the datacannot be explained by X = g(Y ) + N with a noise N
that is independent of Y . It turns out that also themodel corresponding to the other direction X → Y
can be rejected since a regression of Y onto X leads todependent residuals (pHSIC (X, Y − fˆ(X)) = 0.0015).

0.6

0.62

estimated NY

Y

−4
−4.2
0

0.5

1

X

1.5

0.66

0.68

0.7

0
−0.01

0

0.2

0.4
0.6
estimated T

0.8

1

0

0.2

0.4
0.6
estimated T

0.8

1

−2
0
2
estimated NX

4

6
−3
x 10

0.01
0
−0.01

−3.6
−3.8

X

0.01

estimated NY

with independent residuals, it did so in the first or second iteration after optimizing the projections according to the dependence measures. We still think that
this step can be useful for difficult data sets, where thecurve that minimizes the `2 distance is very different
from the ground truth.

0.01
0

−0.01
−8

−6

−4

estimated T

1

0.5

0
0

0.2

0.4

T

0.6

0.8

1

Figure 4: Data set 1. Top: true (black) and estimated (red) curve. Bottom: The estimated valuesof the confounder are plotted against the true values.
Apart from the arbitrary reparameterization t 7→ −t
the method inferred confounder values close to the true
ones.

Figure 5: Data set 2. Top: true (black) and estimated(red) curve. Others: Scatter plots of the fitted residuals against each other and against estimated values forthe confounder. The fact that the noise NX has beensampled with a higher variance than NY can also be
detected in the fitted residuals.
Since the variances of NX and NY differ significantly
and the sample size is small, we can (as expected)
even fit a direct causal relationship between X andY : Assuming the model
X = g(Y ) + N

In this example the empirical joint distribution of(X, Y ) does not allow a simple direct causal relation-(11)

and fitting the function ĝ by Gaussian Process re-

Y

7

6
0

X

0.5

1

0.2

estimated T

−1.5

−1

−0.5
0
pressure KABE

0.5

1

1.5

150
100
50
0
0

50

100

150

1.5

0.4

Real data

ASOS data.
The Automated Surface Observations Systems(ASOS) consists of several stations that automaticallycollect and transmit weather data every minute. We
used 150 values for air pressure that were collected bystations KABE and KABI in January 2000 [NCDC,
2009]. We expect the time to be a confounder. As in
the other experiments a projection minimizing the `2
distance would not be sufficient: after the initialization
step we obtain p-values, which reject independence(pHSIC (N̂X , N̂Y ) = 0.00, pHSIC (N̂X , T̂ ) = 0.00,
pHSIC (N̂Y , T̂ ) = 0.02). After the projection step
minimizing the sum of HSICs the residuals are
regarded as independent: pHSIC (N̂X , N̂Y ) = 1.00,

estimated NY

Figure 6: Data set 3. To check whether our method
does not always find a confounder we simulated a data
set where the noise clearly depends on T . Indeed thealgorithm does not find an independent solution andstops after K = 5 iterations. Top: true (black) and estimated (red) curve. Bottom: the estimated residuals
clearly depend on the estimated confounder.
4.2−2

Figure 7: ASOS data. Top: scatter plot of the data,
together with the estimated path ŝ (note that it is notinterpolating between the data points). Bottom: ordering of the estimated confounder values against thetrue ordering. The true ordering is almost completely
recovered.

0

0

0
−1

true order of T

estimated NX

−0.5

1

−2
−2.5

estimated NY

estimated NY−1

2
pressure KABI

Data set 3.
We also simulated a data set for which the noise terms
NX and NY clearly depend on T . Figure 6 shows ascatter plot of the data set, the outcome curve of thealgorithm after K = 5 iterations (top) and a scatter plot between the estimated residuals N̂Y and confounder values T̂ (bottom). The method did not finda curve and corresponding projections for which theresiduals were independent (pHSIC (N̂Y , T̂ ) = 0.00, forexample), and thus results in “Data cannot be fitted
by a CAN model”. This makes sense since the data
set was not simulated according to model (2).

pHSIC (N̂X , T̂ ) = 1.00, pHSIC (N̂Y , T̂ ) = 0.16. Figures
7 and 8 show the results. The confounder has beensuccessfully identified.

estimated order of T

gression, for example, results in independent residuals:
pHSIC (Y, X − ĝ(Y )) = 0.97. Thus we regard the model(11) and thus Y → X to be true. This does not contradict the identifiability conjecture because the dependencies introduced by setting the noise N̂Y mistakenly
to zero are not detectable at this sample size.

0.1
0
−0.1

0

0.2

0.4
0.6
estimated T

0.8

1

0

0.2

0.4
0.6
estimated T

0.8

1

0.02

0.04

0
−0.2

0
−0.2
−0.06

−0.04

−0.02

0
estimated NX

Figure 8: ASOS data. Residuals plotted against eachother and against the estimated confounder. The hypothesis of independence is not rejected, which means
the method identified the confounder.

5

Conclusion and Future Work

We have proposed a method to identify the confounderof a model where two observable variables are functions of an unobserved confounder plus additive noise.
We have provided a theoretical motivation for the
method, and showed that the algorithm works on both
simulated and real world data sets.
Our initial results are encouraging, and our theoreticalmotivation provides some insight into why the problem
is solvable in the first place.

A complete identifiability result in the style of Hoyer
et al. [2009], however, would clearly be desirable, alongwith further experimental evidence.
References
A. Gretton, O. Bousquet, A. Smola, and
B. Schölkopf.
Measuring statistical dependencewith Hilbert-Schmidt norms. In ALT, pages 63–78.
Springer-Verlag, 2005.
A. Gretton, K. Fukumizu, C.-H. Teo, L. Song,B. Schölkopf, and A. Smola. A kernel statistical testof independence. In Advances in Neural Information
Processing Systems 20, pages 585–592, Cambridge,
MA, 2008. MIT Press.
T. Hastie and W. Stuetzle. Principal curves. Journalof the American Statistical Association, 84:502–516,
1989.
P. O. Hoyer, S. Shimizu, A. J. Kerminen, and
M. Palviainen. Estimation of causal effects using
linear non-Gaussian causal models with hidden variables. Int. J. Approx. Reasoning, 49(2):362–378,
2008. ISSN 0888-613X. doi: http://dx.doi.org/10.
1016/j.ijar.2008.02.006.
P. O. Hoyer, D. Janzing, J. M. Mooij, J. Peters, and
B. Schölkopf. Nonlinear causal discovery with additive noise models. In D. Koller, D. Schuurmans,
Y. Bengio, and L. Bottou, editors, Advances in Neural Information Processing Systems 21 (NIPS*2008),
pages 689–696. 2009.
A. Kankainen. Consistent testing of total independence based on the empirical characteristic function.
PhD Thesis, University of Jyväskylä, 1995.
NCDC. ASOS 1-minute data (DSI 6406/page 2data), 2009. URL http://www.ncdc.noaa.gov/oa/
climate/climatedata.html.
J. Pearl. Causality: Models, Reasoning, and Inference. Cambridge University Press, 2000.
C. E. Rasmussen and C. Williams.
GPML
code.
http://www.gaussianprocess.org/gpml/
code, 2007.
C. E. Rasmussen and C. Williams. Gaussian Processes for Machine Learning. MIT Press, 2006.
H. Reichenbach. The direction of time. Universityof Los Angeles Press, Berkeley, 1956.
B. Schölkopf and A. Smola. Learning with Kernels.
MIT Press, 2002.

R. Silva, R. Scheines, C. Glymour, and P. Spirtes.
Learning the structure of linear latent variable models. Journal of Machine Learning Research, 7:191–
246, 2006.
P. Spirtes, C. Glymour, and R. Scheines. Causation,
Prediction, and Search. Springer-Verlag, 1993. (2nded. MIT Press 2000).
J. B. Tenenbaum, V. de Silva, and J. C. Langford. A
global geometric framework for nonlinear dimensionality reduction. Science, 290:2319–2323, 2000.
L. J. P. van der Maaten. An introduction to dimensionality reduction using matlab.
Technical
Report MICC-IKAT 07-07, Maastricht University,
Maastricht, 2007.

Appendix: proof of Lemma 2
Fix y and set ty := w(y), sy := s(w(y)), βy := −w0 (y).
We compute:Z
p(y)Ey ((T − ty )n ) = (t − ty )n r(y − v(t))s(t) dt.
We want to make the substitution ỹ = y − v(t), i.e.,t = w(y − ỹ). Application of Taylor’s theorem yields:w(y − ỹ) − ty = βy ỹ + ỹ 2 || ≤ kw00 k∞ /2,0−w (y − ỹ) = βy + η ỹ

|η| ≤ kw00 k∞ ,s(w(y − ỹ)) = sy + ζ(βy ỹ + ỹ 2 )

|ζ| ≤ ks0 k∞ .

where we suppress the dependencies of , η, ζ on ỹ inthe notation. Therefore:Z(t − ty )n r(y − v(t))s(t) dtZ= (βy ỹ + ỹ 2 )n r(ỹ) sy + ζ(βy ỹ + ỹ 2 ) (βy + η ỹ) dỹ.
The special case n = 0 yields
Z
p(y) = r(ỹ) sy + ζ(βy ỹ + ỹ 2 ) (βy + η ỹ) dỹ
= sy βy + O(ks0 k∞ + kw00 k∞ ).
For arbitrary n, we obtainZ(βy ỹ + ỹ 2 )n r(ỹ) sy + ζ(βy ỹ + ỹ 2 ) (βy + η ỹ) dỹ= sy βy E (βy NY )n + O(ks0 k∞ + kw00 k∞ ).
The error terms contain moments of NY , which are all
finite by assumption, and (positive) powers of sy andβy , which are also finite. Thus we conclude thatn (y) = Ey ((T −ty )n )−E(βyn NYn ) = O(ks0 k∞ +kw00 k∞ ).
