Causal Inference on Multivariate and Mixed-Type Data

TL;DR

This paper introduces a new information-theoretic method using MDL and a greedy algorithm called Crack to accurately infer causal directions in multivariate and mixed-type data, addressing a complex problem in causal inference.

Contribution

It proposes a novel MDL-based scoring method for mixed data types and a fast greedy algorithm, Crack, for causal inference, extending previous approaches to multivariate and mixed data.

Findings

Crack reliably infers causal direction with high accuracy.

Method works well on both univariate and multivariate cause-effect pairs.

Effective on data with mixed data types.

Abstract

Given data over the joint distribution of two random variables $X$ and $Y$, we consider the problem of inferring the most likely causal direction between $X$ and $Y$. In particular, we consider the general case where both $X$ and $Y$ may be univariate or multivariate, and of the same or mixed data types. We take an information theoretic approach, based on Kolmogorov complexity, from which it follows that first describing the data over cause and then that of effect given cause is shorter than the reverse direction. The ideal score is not computable, but can be approximated through the Minimum Description Length (MDL) principle. Based on MDL, we propose two scores, one for when both $X$ and $Y$ are of the same single data type, and one for when they are mixed-type. We model dependencies between $X$ and $Y$ using classification and regression trees. As inferring the optimal model is NP-hard, we propose Crack, a fast greedy algorithm to determine the most likely causal direction directly from the data. Empirical evaluation on a wide range of data shows that Crack reliably, and with high accuracy, infers the correct causal direction on both univariate and multivariate cause-effect pairs over both single and mixed-type data.

Tables

Table 1. Table 1 : Comparison of LTR , Ergo , Origo and Crack on eleven multivariate data sets. We write (n/a) whenever a method is not applicable on the pair.

					Decisions
Dataset	$m$	$k$	$l$	LTR	Ergo	Origo	Crack
Climate	$\mathrm{10\hspace{0.17em}226}$	4	4	✓	✓	–	–
Ozone	$989$	1	3	(n/a)	✓	✓	✓
Car	$392$	3	2	–	✓	✓	✓
Radiation	$72$	16	16	–	–	–	✓
Symptoms	$120$	6	2	✓	✓	–	✓
Brightness	$\mathrm{1\hspace{0.17em}000}$	9	1	(n/a)	(n/a)	–	✓
Chemnitz	$\mathrm{1\hspace{0.17em}440}$	3	7	✓	✓	✓	✓
Precip.	$\mathrm{4\hspace{0.17em}748}$	3	12	✓	–	–	✓
Stock 7	$\mathrm{2\hspace{0.17em}394}$	4	3	–	✓	–	✓
Stock 9	$\mathrm{2\hspace{0.17em}394}$	4	5	–	✓	–	✓
Haberman	$306$	3	1	✓	✓	–	–
Iris flower	$150$	4	1	(n/a)	(n/a)	–	✓
Canis	$\mathrm{2\hspace{0.17em}183}$	4	2	(n/a)	(n/a)	✓	✓
Lepus	$\mathrm{2\hspace{0.17em}183}$	4	3	(n/a)	(n/a)	✓	✓
Martes	$\mathrm{2\hspace{0.17em}183}$	4	2	(n/a)	(n/a)	✓	✓
Mammals	$\mathrm{2\hspace{0.17em}183}$	4	7	(n/a)	(n/a)	✓	✓
Octet	$82$	1	10	(n/a)	✓	✓	✓
Accuracy				$0.56$	$0.82$	$0.47$	$0.88$

Equations

$$|\mathbb{D}(a)|=\{\begin{array}{ll}\#\text{values}& \text{if type(a) is nominal}\\ \frac{\max (a)-\min (a)}{res(a)}+1& \text{if type(a) is numeric},\end{array}$$

$$K(x)=\min \{|p|\mid p\in \{0,1{\}}^{\ast },\mathcal{U}(p)=x\}.$$

$$L(D,M)=L(M)+L(D\mid M),$$

$$K(P(X))+K(P(Y\mid X))\le K(P(Y))+K(P(X\mid Y)).$$

$${\Delta }_{X\to Y}^{\ast }=\frac{K(P(X))+K(P(Y\mid X))}{K(P(X))+K(P(Y))},$$

$$L(D,M_{X\to Y})=L(X,M_X)+L(Y,M_{Y\mid X}\mid X),$$

$${\Delta }_{X\to Y}=\frac{L(X,M_X)+L(Y,M_{Y\mid X}\mid X)}{L(X,M_X)+L(Y,M_Y)}.$$

$$L(X,M_X)-L(X,M_{X|Y}|Y)<L(Y,M_Y)-L(Y,M_{Y|X}|X).$$

$$\frac{L(X,M_X)-L(X,M_{X|Y}|Y)}{L(X,M_X)}=1-\frac{L(X,M_{X|Y}|Y)}{L(X,M_X)}.$$

$$\frac{L(X,M_{X\mid Y}\mid Y)}{L(X,M_X)}>\frac{L(Y,M_{Y\mid X}\mid X)}{L(Y,M_Y)}.$$

$${\delta }_{X\to Y}=\frac{1}{|Y|}\sum _{Y_i\in Y}\frac{L(Y_i,M_{Y_i\mid X}\mid X)}{L(Y_i,M_{Y_i})}$$

$$L(T)=|T|+\sum _{v\in int(T)}(1+L(v))+\sum _{l\in lvs(T)}L(l).$$

$$L_{1split}(v)=1+\log |A|+\{\begin{array}{l}\log |\mathbb{D}(a_j)|\text{ if ai is categorical,}\\ \log |\mathbb{D}(a_j)-1|\text{ else.}\end{array}$$

$$L_{\text{k}split}(v)=1+\log |A|+\{\begin{array}{l}0\text{ if ai is categorical,}\\ L_{\mathbb{N}}(k)\text{ numeric,}\end{array}$$

$$L_{reg}(v)=\log |A|+\sum _{\varphi \in \Phi (v)}\left(1+L_{\mathbb{N}}(s)+L_{\mathbb{N}}(\lfloor \varphi \cdot 10^s\rfloor )\right),$$

$$L_{nom}(l)=$$

$$-|l|\sum _{c\in \mathbb{D}(a_i)}\mathrm{Pr}(a_i=c\mid l)\log \mathrm{Pr}(a_i=c\mid l).$$

$$L_{num}(l\mid \sigma ,\mu )=$$

$$L_{num}(l\mid \min ,\max )=$$

$$L_{num}(l)$$

$$+\min \{L_{num}(l\mid \sigma ,\mu ),L_{num}(l\mid \min ,\max )\}.$$

Full text

The ideal score is not computable, but can be approximated through the Minimum Description Length (MDL) principle. Based on MDL, we propose two scores, one for when both $X$ and $Y$ are of the same single data type, and one for when they are mixed-type. We model dependencies between $X$ and $Y$ using classification and regression trees. As inferring the optimal model is NP-hard, we propose Crack, a fast greedy algorithm to determine the most likely causal direction directly from the data.

Empirical evaluation on a wide range of data shows that Crack reliably, and with high accuracy, infers the correct causal direction on both univariate and multivariate cause-effect pairs over both single and mixed-type data.

1 Introduction

Telling cause from effect is one of the core problems in science. It is often difficult, expensive, or impossible to obtain data through randomized trials, and hence we often have to infer causality from, what is called, observational data [21]. We consider the setting where, given data over the joint distribution of two random variables $X$ and $Y$, we have to infer the causal direction between $X$ and $Y$. In other words, our task is to identify whether it is more likely that $X$ causes $Y$, or vice versa, that $Y$ causes $X$, or that the two are merely correlated.

In practice, $X$ and $Y$ do not have to be of the same type. The altitude of a location (real-valued), for example, determines whether it is a good habitat (binary) for a mountain hare. In fact, neither $X$ nor $Y$ have to be univariate. Whether or not a location is a good habitat for an animal, is not just caused by a single aspect, but by a combination of conditions, which not necessarily are of the same type. We are therefore interested in the general case where $X$ and $Y$ may be of any cardinality, and may be single or mixed-type.

To the best of our knowledge there exists no method for this general setting. Causal inference based on conditional independence tests, for example, requires three variables, and cannot decide between ${X\rightarrow Y}$ and ${Y\rightarrow X}$ [21]. All existing methods that consider two variables are only defined for single-type pairs. Additive Noise Models (ANMs), for example, have only been proposed for univariate pairs of real-valued [24] or discrete variables [23], and similarly so for methods based on the independence of $P(X)$ and $P(Y\mid X)$ [28, 16]. Trace-based methods require both $X$ and $Y$ to be strictly multivariate real-valued [9, 2], and whereas Ergo [33] also works for univariate pairs, these again have to be real-valued. We refer the reader to Sec. 6 for a more detailed overview of related work.

Our approach is based on algorithmic information theory. That is, we follow the postulate that if ${X\rightarrow Y}$, it will be easier—in terms of Kolmogorov complexity—to first describe $X$, and then describe $Y$ given $X$, than vice-versa [11, 33, 1]. Kolmogorov complexity is not computable, but can be approximated through the Minimum Description Length (MDL) principle [26, 5], which we use to instantiate this framework. In addition, we develop a causal indicator that is able to handle multivariate and mixed-type data.

To this end, we define an MDL score for coding forests, a model class where a model consists of classification and regression trees. By allowing dependencies from $X$ to $Y$, or vice versa, we can measure the difference in complexity between ${X\rightarrow Y}$ and ${Y\rightarrow X}$. Discovering a single optimal decision tree is already NP-hard [20], and hence we cannot efficiently discover the coding forest that describes the data most succinctly. We therefore propose Crack, an efficient greedy algorithm for discovering good models directly from data.

Through extensive empirical evaluation on synthetic, benchmark, and real-world data, we show that Crack performs very well in practice. It performs on par with existing methods for univariate single-type pairs, is the first to handle pairs of mixed data type, and outperforms the state of the art on multivariate pairs with a large margin. It is also very fast, taking less than 4 seconds over any pair in our experiments.

2 Preliminaries

First, we introduce notation and give brief primers to Kolmogorov complexity and the MDL principle.

2.1 Notation

In this work we consider data $D$ over the joint distribution of random variables $X$ and $Y$. Such data $D$ contains $n$ records over a set $A$ of $|A|=|X|+|Y|=m$ attributes, $a_{1},\dots,a_{m}\in A$. An attribute $a$ has a type $\textit{type}(a)$ where $\textit{type}(a)\in\{\text{{binary}, {categorical}, {numeric}}\}$. We will refer to binary and categorical attributes as nominal attributes. The size of the domain of an attribute $a$ is defined as

[TABLE]

where $\mathit{res}(a)$ is the resolution at which the data over attribute $a$ was recorded. For example, a resolution of 1 means that we consider integers, of $0.01$ means that $a$ was recorded with a precision of up to a hundredth.

We will consider decision and regression trees. In general, a tree $T$ consist of $|T|$ nodes. We identify internal nodes as $\mathit{v}\in\mathit{int}(T)$, and leaf nodes as $\mathit{l}\in\mathit{lvs}(T)$. A leaf node $l$ contains $|l|$ data points.

All logarithms are to base 2, and we use $0\log 0=0$.

2.2 Kolmogorov Complexity, a brief primer

The Kolmogorov complexity of a finite binary string $x$ is the length of the shortest binary program $p^{*}$ for a universal Turing machine $\mathcal{U}$ that generates $x$, and then halts [13, 15]. Formally, we have

[TABLE]

Simply put, $p^{*}$ is the most succinct algorithmic description of $x$, and the Kolmogorov complexity of $x$ is the length of its ultimate lossless compression. Conditional Kolmogorov complexity, $K(x\mid y)\leq K(x)$, is then the length of the shortest binary program $p^{*}$ that generates $x$, and halts, given $y$ as input. For more details see [15].

2.3 MDL, a brief primer

The Minimum Description Length (MDL) principle [26, 5] is a practical variant of Kolmogorov Complexity. Intuitively, instead of all programs, it considers only those programs that we know that output $x$ and halt. Formally, given a model class $\mathcal{M}$, MDL identifies the best model $M\in\mathcal{M}$ for data $D$ as the one minimizing

[TABLE]

where $L(M)$ is the length in bits of the description of $M$, and $L(D\mid M)$ is the length in bits of the description of data $D$ given $M$. This is known as two-part MDL. There also exists one-part, or refined MDL, where we encode data and model together. Refined MDL is superior in that it avoids arbitrary choices in the description language $L$, but in practice only computable for certain model classes. Note that in either case we are only concerned with code lengths — our goal is to measure the complexity of a dataset under a model class, not to actually compress it [5].

3 Causal Inference by Compression

We pursue the goal of causal inference by compression. Below we give a short introduction to the key concepts.

3.1 Causal Inference by Complexity

The problem we consider is to infer, given data over two correlated variables $X$ and $Y$, whether $X$ caused $Y$, whether $Y$ caused $X$, or whether $X$ and $Y$ are only correlated. As is common, we assume causal sufficiency. That is, we assume there exists no hidden confounding variable $Z$ that is the common cause of both $X$ and $Y$.

The Algorithmic Markov condition, as recently postulated by Janzing and Schölkopf [11], states that factorizing the joint distribution over $\mathit{cause}$ and $\mathit{effect}$ into $P(\mathit{cause})$ and $P(\mathit{effect}\mid\mathit{cause})$, will lead to simpler—in terms of Kolmogorov complexity—models than factorizing it into $P(\mathit{effect})$ and $P(\mathit{cause}\mid\mathit{effect})$. Formally, they postulate that if $X$ causes $Y$,

[TABLE]

While in general the symmetry of information, $K(x)+K(y\mid x)=K(y)+K(x\mid y)$, holds up to an additive constant [15], Janzing and Schölkopf [11] showed it does not hold when $X$ causes $Y$, or vice versa. Based on this, Budhathoki & Vreeken [1] proposed

[TABLE]

as a causal indicator that uses this asymmetry to infer that $X\rightarrow Y$ as the most likely causal direction if $\Delta_{{X\rightarrow Y}}^{*}<\Delta_{{Y\rightarrow X}}^{*}$, and vice versa. The normalisation has no function during inference, but does help to interpret the confidence of the indicator.

Both scores assume access to the true distribution $P(\cdot)$, whereas in practice we only have access to empirical data. Moreover, following from the halting problem, Kolmogorov complexity is not computable. We can approximate it, however, via MDL [15, 5], which also allows us to directly work with empirical distributions.

3.2 Causal Inference by MDL

For causal inference by MDL, we will need to approximate both $K(P(X))$ and $K(P(Y\mid X))$. For the former, we need to consider the model classes $\mathcal{M}_{X}$ and $\mathcal{M}_{Y}$, while for the latter we need to consider class $\mathcal{M}_{Y\mid X}$ of models $M_{Y\mid X}$ that describe the data of $Y$ dependent the data of $X$.

That is, we are after the causal model $M_{{X\rightarrow Y}}=(M_{X},M_{Y\mid X})$ from the class $\mathcal{M}_{X\rightarrow Y}=\mathcal{M}_{X}\times\mathcal{M}_{Y\mid X}$ that best describes the data $Y$ by exploiting as much as possible structure of $X$ to save bits. By MDL, we identify the optimal model $M_{X\rightarrow Y}\in\mathcal{M}_{X\rightarrow Y}$ for data $D$ over $X$ and $Y$ as the one minimizing

[TABLE]

where the encoded length of the data of $X$ under a given model is encoded using two-part MDL, similarly so for $Y$, if we consider the inverse direction.

To identify the most likely causal direction between $X$ and $Y$ by MDL we can now simply rewrite Eq. (3.3)

[TABLE]

Similar to the original score, we infer that $X$ is a likely cause of $Y$ if $\Delta_{{X\rightarrow Y}}<\Delta_{{Y\rightarrow X}}$, $Y$ is a likely cause of $X$ if $\Delta_{{Y\rightarrow X}}<\Delta_{{X\rightarrow Y}}$, and that $X$ and $Y$ are only correlated or might have a common cause if $\Delta_{{X\rightarrow Y}}=\Delta_{{Y\rightarrow X}}$.

3.3 Normalized Causal Indicator

Although $\Delta$ has nice theoretical properties, it has a mayor drawback. It assumes that a bit gain in the description of the data over one attribute has the same importance as one bit gain in the description of the data over another attribute. This does not hold true if these attributes have different intrinsic complexities, such as when their domain sizes strongly differ. For example, a continuous valued attribute is very likely to have a much higher intrinsic complexity than a binary attribute. This means that gaining $k$ bits from an attribute with a large domain is not comparable to gaining $k$ bits from an attribute with a small domain. Since the $\Delta$ indicator compares the absolute difference in bits, it does not account for differences w.r.t. the intrinsic complexity. Hence, $\Delta$ is highly likely to be a bad choice when $X$ and $Y$ are of different, or of mixed-type data.

We therefore propose an alternative indicator for causal inference on mixed-type data. Instead of taking the absolute difference between the conditioned and unconditioned score, we instead consider relative differences w.r.t. the marginal. We can derive the Normalized Causal Indicator ($\mathit{NCI}$) starting from the numerator of the $\Delta$ indicator. By subtracting the conditional costs on both sides, we have

[TABLE]

Since the aim of the $\mathit{NCI}$ is to measure the relative gain, we divide by the costs of the unconditioned data

[TABLE]

After this step, we can conclude that for the relative gain it holds, if ${X\rightarrow Y}$

[TABLE]

This score can be understood as an instantiation of the Ergo indicator proposed by Vreeken [33]. From the derivation, we can easily see that the difference between the score of both indicators depends only on the normalization factor and hence both are based on the Algorithmic Markov condition. It turns out, however, that the Ergo indicator is also biased. Although it balances the gain between $X$ and $Y$, we need a score that does not impose prior assumptions to the individual attributes of $X$ and $Y$. With the Ergo indicator, it could happen that a single $X_{i}\in X$ dominates the whole score for $X$. To account for this, we assume independence among the variables within $X$ and $Y$, meaning that the domain of two individual attributes within $X$ or $Y$ is allowed to differ. Hence, we formulate the $\mathit{NCI}$, which we from now on denote by $\delta$, from $X$ to $Y$ as

[TABLE]

and analogously $\delta_{{Y\rightarrow X}}$. To avoid bias towards dimensionality, we normalize by the number of attributes. As above, we infer ${X\rightarrow Y}$ if $\delta_{{X\rightarrow Y}}<\delta_{{Y\rightarrow X}}$ and vice versa.

In practice, we expect that $\Delta$ performs well on data where $X$ and $Y$ are of the same type, especially when $|X|=|Y|$ and the domain sizes of their attributes are balanced. For unbalanced domains, dimensionality, and especially for mixed-type data, we expect $\delta$ to perform much better. The experiments indeed confirm this.

4 MDL for Tree Models

To use the above defined causal indicators in practice, we need to define a casual model class $\mathcal{M}_{X\rightarrow Y}$, how to encode a model $M\in\mathcal{M}$ in bits, and how to encode a dataset $D$ using a model $M$. As models we consider tree models, or, coding forests. A coding forest $M$ contains per attribute $a_{i}\in A$ one coding tree $T_{i}$. A coding tree $T_{i}$ encodes the values of $a_{i}$ in its leaves, splitting or regressing the data of $a_{i}$ on attribute $a_{j}$ ($i\neq j$) in its internal nodes to encode the data of $a_{i}$ more succinctly.

We encode the data over attribute $a_{i}$ with the corresponding coding tree $T_{i}$. The encoded length of data $D$ and $M$ then is $L(D,M)=\sum_{a_{i}\in A}L(T_{i})$, which corresponds to the sum of costs of the individual trees.

To ensure lossless decoding, there needs to exist an order on the trees $T\in M$ such that we can transmit these one by one. In other words, in a valid tree model there are no cyclic dependencies between the trees $T\in M$, and a valid model can hence be represented by a DAG. Let $\mathcal{M}(D)$ be the set of all valid tree models for $D$, that is, $M\in\mathcal{M}(D)$ is a set of $|A|$ trees such that the data types of the leafs in $T_{i}$ corresponds to the data type of attribute $a_{i}$, and its dependency graph is acyclic.

We write $\mathcal{M}_{X}(X)$ and $\mathcal{M}_{Y}(Y)$ to denote the subset of valid coding forests for $X$ and $Y$, where we do not allow dependencies. To describe the possible set of models where we allow attributes of $X$ to only depend on attributes of $Y$ we write $\mathcal{M}_{X\mid Y}(X)$ and do so accordingly for $Y$ depending only on $X$. If an attribute does not have any incoming dependencies, its tree is a stump. Fig. 1 shows the DAG for a toy data set, and an example tree for $Y_{2}$. From the DAG, the set of purple edges would be a valid model in $\mathcal{M}_{Y\mid X}(Y)$, whereas the orange edges are a valid model from $\mathcal{M}_{X\mid Y}(X)$.

Cost of a Tree

The encoded cost of a tree consists of two parts. First, we transmit the topology of the tree. From the root node on we indicate with one bit per node whether it is a leaf or an internal node, and if the latter, one further bit to identify whether it is a split or regression node. Formally we have that

[TABLE]

Next, we explain how we encode internal nodes and then specify the encoding for leaf nodes.

Cost of a Single Split

The length of a split node $\mathit{v}$ is

[TABLE]

whereas we first need one bit to indicate this is a single-split node, then identify in $\log|A|$ bits on which attribute $a_{j}$ we split, and third the split condition.

The split condition can be any value in the domain for categorical, and can lie in between two consecutive values of a numeric attribute ($|\mathbb{D}(a_{j})-1|$ choices). For binary we only have one option, resulting in zero cost.

Costs of a Multiway split

A multiway split is only possible for categorical and real valued data. As there are exponentially many multiway splits, we consider only a subset. The costs for a multiway split are

[TABLE]

where the first two terms are similar to above. For categorical data, we only consider splitting on all values, and hence have no further cost. For numeric data, we only split non-deterministic cases, i.e. if there exist duplicate values. To do so, we split on every such value that occurs at least $k$ times, and one residual split for all remaining data points. To encode such a split, we transmit $k$ using $L_{\mathbb{N}}(k)$ bits, where $L_{\mathbb{N}}$ is the MDL optimal encoding for integers $z\geq 1$ [27].

Cost of Regressing

For a regression node we also first encode the target attribute, and then the parameters of the regression, i.e.

[TABLE]

where $\Phi(\mathit{v})$ denotes the set of parameters for the regression. For linear regression, it consists of $\alpha$ and $\beta$, while for quadratic regression it further contains $\gamma$. To describe each parameter $\phi\in\Phi$ we transmit it up to a user defined precision, e.g. $0.001$, we first encode the corresponding number of significant digits $s$, e.g. $3$, and then the shifted parameter value.

Next, we describe how to encode the data in a leaf $l$. As we consider both nominal and numeric attributes, we need to define $L_{\mathit{nom}}(l)$ for nominal and $L_{\mathit{num}}(l)$ for numeric data.

Cost of a Nominal Leaf

To encode the data in a leaf of a nominal attribute, we can use refined MDL [14]. That is, we encode minimax optimally, without having to make design choices [5]. In particular, we encode the data of a nominal leaf using the normalized maximum likelihood (NML) distribution as

[TABLE]

Kontkanen & Myllymäki [14] derived a recursive formula to calculate this in linear time.

Cost of a Numerical Leaf

For numeric data existing refined MDL encodings have high computational complexity [14]. Hence, we encode the data in numeric leaves using two-part MDL, using point models with Gaussian or uniform noise. A split or a regression on an attribute aims to reduce the variance or the domain in the leaf. We encode the costs of a numeric leaf as

[TABLE]

given empirical mean $\mu$ and variance $\sigma$ or as uniform given $\min$ and $\max$ as

[TABLE]

We encode the data as Gaussian if this costs fewer bits than encoding it as uniform. To indicate this decision, we use one bit and encode the minimum of both plus the corresponding parameters. As we consider empirical data, we can safely assume that all parameters lie in the domain of the given attribute. Since we do not have any preference on the parameter values, the encoded costs of a numeric leaf $\mathit{l}$ are

[TABLE]

Putting it all together, we now know how to compute $L(D,M)$, by which we can formally define the Minimal Coding Forest problem.

Minimal Coding Forest Problem * Given a data set $D$ over a set of attributes $A=\{a_{1},\ldots,a_{m}\}$, and $\mathcal{M}$ a valid model class for $A$. Find the smallest model $M\in\mathcal{M}$ such that $L(D,M)$ is minimal. *

From the fact that both inferring optimal decision trees and structure learning of Bayesian networks—to which our tree-models reduce for nominal-only data and splitting on all values—are NP-hard [20], it trivially follows that the Minimal Coding Forest problem is also NP-hard. Hence, we resort to heuristics.

5 The Crack Algorithm

Knowing the score $L(D,M)$ and the problem, we can now introduce the Crack algorithm, which stands for classification and regression based packing of data. Crack is an efficient greedy heuristic for discovering a coding forest $M$ from given model class $\mathcal{M}$ with low $L(D,M)$. It builds upon the well-known ID3 algorithm [25]. In the next section we explain the main aspects of the algorithm.

Greedy algorithm

We give the pseudocode of Crack as Algorithm 1. Before running the algorithm, we set the resolution per attribute, which is $1$ for nominal data (line 1). For numeric data, we calculate the differences between adjacent values, and to reduce sensitivity to outliers take the $k^{\mathit{th}}$ smallest difference as resolution. In general, setting $k$ to $0.1n$ works well in practice.

Crack starts with an empty model consisting of only trivial trees, i.e. leaf nodes containing all records, per attribute (line 1). The given model class $\mathcal{M}$ implicitly defines a graph $\mathcal{G}$ of dependencies between attributes that we are allowed to consider (line 1). To make sure the returned model is valid, we need to maintain a graph representing its dependencies (lines 1–1). We iteratively discover that refinement of the current model that maximizes compression. To find the best refinement, we consider every attribute (line 1), and every legal additional split or regression of its corresponding tree (line 1). A refinement is only legal when the dependency is allowed by the model family (line 1), the dependency graph remains acyclic, and we do not split or regress twice on the same attribute (line 1). We keep track of the best found refinement.

The key subroutine of Crack is RefineLeaf, in which we discover the optimal refinement of a leaf $l$ in tree $T_{i}$. That is, it finds the optimal split of $l$ over all candidate attributes $a_{j}$ such that we minimize the encoded length. In case both $a_{i}$ and $a_{j}$ are numeric, RefineLeaf also considers the best linear and quadratic regression and decides for the variant with the best compression—choosing to split in case of a tie. In the interest of efficiency, we do not allow splitting or regressing multiple times on the same candidate.

Since we use a greedy heuristic to construct the coding trees, we have a worst case runtime of $O(2^{m}n)$, where $m$ is the number of attributes and $n$ is the number of rows. Although the worst case runtime is exponential, in practice, Crack takes only a few seconds.

Causal Inference with Crack

To compute our causal indicators we have to run Crack twice on $D$. First with model class $\mathcal{M}_{X\mid Y}(X)$ to obtain $M_{X\mid Y}(X)$ and second with $\mathcal{M}_{Y\mid X}(Y)$, to obtain $M_{Y\mid X}(Y)$. To estimate $\mathcal{M}_{X}(X)$, we assume a uniform prior $L(X\mid M_{X})=-n\sum_{a_{i}\in X}\log res(a_{i})$ and similarly for $M_{Y}(Y)$. We can use these scores to calculate both the $\delta$ score and the $\Delta$ score. We will refer to Crack using the $\delta$ indicator as $\textsc{Crack}_{\delta}$, and Crack with the $\Delta$ indicator as $\textsc{Crack}_{\Delta}$.

6 Related Work

Causal inference on observational data is a challenging problem, and has recently attracted a lot of attention [21, 11, 29, 1]. Most existing proposals, however, are highly specific in the type of causal dependencies and type of variables they can consider.

Clasical constrained-based approaches, such as conditional independence tests, require three observed random variables [30, 21], cannot distinguish Markov equivalent causal DAGs [32] and hence cannot decide between ${X\rightarrow Y}$ and ${Y\rightarrow X}$. Recent approaches use properties of the joint distribution to break the symmetry.

Additive Noise Models (ANMs) [29], for example, assume that the effect is a function of the cause and cause-independent additive noise. ANMs exist for univariate real-valued [29, 8, 34, 24] and discrete data [22]. A related approach considers the asymmetry in the joint distribution of $\mathit{cause}$ and $\mathit{effect}$ for causal inference. The linear trace method (LTR) [9] and the kernelized trace method (KTR) [2] aim to find a structure matrix $A$ and the covariance matrix $\Sigma_{X}$ to express $Y$ as $AX$. Both methods are restricted to multivariate continuous data. Sgouritsa et al. [28] show that the marginal distribution $P(\mathit{cause})$ of the cause is independent of the conditional distribution $P(\mathit{effect}\mid\mathit{cause})$ of the effect. They proposed Cure, using unsupervised reverse regression on univariate continuous pairs. Liu et al [16] use distance correlation to identify the weakest dependency between univariate pairs of discrete data.

The algorithmic information-theoretic approach views causality in terms of Kolmogorov complexity. The key idea is that if $X$ causes $Y$, the shortest description of the joint distribution $P(X,Y)$ is given by the separate descriptions of the distributions $P(X)$ and $P(Y\mid X)$ [11], and justifies additive noise model based causal inference [12]. However, as Kolmogorov complexity is not computable [15], causal inference using algorithmic information theory requires practical implementations, or notions of independence. For instance, the information-geometric approach [10] defines independence via orthogonality in information space for univariate continuous pairs. Vreeken [33] instantiates it with the cumulative entropy to infer the causal direction in continuous univariate and multivariate data. Mooij instantiates the first practical compression-based approach [18] using the Minimum Message Length. Budhathoki and Vreeken approximate $K(X)$ and $K(Y\mid X)$ through MDL, and propose Origo, a decision tree based approach for causal inference on multivariate binary data [1]. Marx and Vreeken[17] proposed Slope, an MDL based method employing local and global regression for univariate numeric data.

In contrast to all methods above, Crack can consider pairs of any cardinality, univariate or multivariate, and of same, different, or even mixed-type data.

7 Experiments

In this section, we evaluate Crack empirically. We implemented Crack in C++, and provide the source code including the synthetic data generator along with the tested datasets for research purposes.​111http://eda.mmci.uni-saarland.de/crack/ The experiments concerning Crack were executed single-threaded. All tested data sets could be processed within seconds; over all pairs the longest runtime for Crack was $3.8$ seconds.

We compare Crack to Cure [28], IGCI [10], LTR [9], Origo [1], Ergo [33] and Slope [17] using their publicly available implementations and recommended parameter settings.

7.1 Synthetic data

The aim of our experiments on synthetic data is to show the advantages of either score. In particular, we expect $\textsc{Crack}_{\Delta}$ to perform well on nominal data and numeric data with balanced domain sizes and dimensions. On the other hand, $\textsc{Crack}_{\delta}$ should have an advantage when it comes to numeric data with varying domain sizes and mixed-type data.

We generate synthetic data with assumed ground truth ${X\rightarrow Y}$ with $|X|=k$ and $|Y|=l$, each having $n=5\,000$ rows, in the following way. First, we randomly assign the type for each attribute in $X$. For nominal data, we randomly draw the number of classes between two (binary) and five and distribute the classes uniformly. Numeric data is generated following a normal distribution taken to the power of $q$ by keeping the sign, leading to a sub-Gaussian ($q<1.0$) or super-Gaussian ($q>1.0$) distribution.​222We use super- and sub-Gaussians to ensure identifiability.

To create data with the true causal direction ${X\rightarrow Y}$, we introduce dependencies from $X$ to $Y$, where we distinguish between splits and refinements. We call the probability threshold to create a dependency $\varphi\in[0,1]$. For each $j\in\{1,\dots,l\}$, we throw a biased coin based on $\varphi$ for each $X_{i}\in X$ that determines if we model a dependency from $X_{i}$ to $Y_{j}$. A split means that we find a category (nominal) or a split-point (numeric) on $X_{i}$ to split $Y_{j}$ into two groups, for which we model its distribution independently. As refinement, we either do a multiway split or model $Y_{j}$ as a linear or quadratic function of $X_{i}$ plus independent Gaussian noise.

Accuracy

First, we compare the accuracies of $\textsc{Crack}_{\delta}$ and $\textsc{Crack}_{\Delta}$ with regard to single-type and mixed-type data. To do so, we generate $200$ synthetic data sets with $|X|=|Y|=3$ for each dependency level where $\varphi\in\{0.0,0.1,\dots 1.0\}$. Figure 2 shows the results for numeric, nominal and mixed-type data. At $\varphi=1.0$ both approaches reach nearly $100\%$ accuracy on single-type data. For single-type data, the accuracy of both methods increases with the dependency. At $\varphi=0$, both approaches correctly do not decide instead of taking wrong decisions. As expected $\textsc{Crack}_{\delta}$ strongly outperforms $\textsc{Crack}_{\Delta}$ on mixed-type data, reaching near $100\%$ accuracy, whereas $\textsc{Crack}_{\Delta}$ reaches only $72\%$. On nominal data, $\textsc{Crack}_{\Delta}$ picks up the correct signal faster than $\textsc{Crack}_{\delta}$.

Dimensionality

Next, we check how sensitive both scores are to dimensionality, whereas we discriminate between asymmetric $k\neq l$ and symmetric $k=l$. We evaluated $200$ data sets per dimensionality. For the symmetric case, both methods are near to $100\%$ on single-type data, whereas only $\textsc{Crack}_{\delta}$ also reaches this target on mixed-type data, as can be seen in the appendix.​11footnotemark: 1 We now discuss the more interesting case for asymmetric pairs in detail.

To test asymmetric pairs, we keep the dimension of one variable at three, $k=3$, while we increase the dimension of the second variable $l$ from one to eleven. To avoid bias, we assigned the dimension $k$ to $X$ and $l$ to $Y$ and swap the dimensions in every other test. We show the results in Figure 3. As expected, we observe that $\textsc{Crack}_{\delta}$ has much fewer difficulties with the asymmetric data sets than $\textsc{Crack}_{\Delta}$. From $l=3$ onwards, $\textsc{Crack}_{\delta}$ is close to $100\%$. On nominal data, $\textsc{Crack}_{\Delta}$ performs near perfect and also has the clear advantage for $l=1$.

7.2 Real world data

Based on the evaluation on synthetic data, we test our approach on univariate benchmark data and multivariate data consisting of known test sets and new causal pairs with known ground truth that we present in the current paper.

Univariate benchmark

To evaluate Crack on univariate data, we apply it to the well-known Tuebingen benchmark data set that consists of $100$ univariate pairs.​333https://webdav.tuebingen.mpg.de/cause-effect/ The pairs mainly consist of numeric data and a few categoric instances. Therefore, we apply $\textsc{Crack}_{\Delta}$. We compare to the state of the art methods that are applicable to multivariate and univariate data, Origo [1] and Ergo [33], and methods specialized for univariate pairs, Cure [28], IGCI [10] and Slope [17]. For each approach, we sort the results by their confidence. According to this order, we calculate for each position $k$ the percentage of correct inferences up to this point, called the decision rate. We weigh the decisions as specified by the benchmark, plot the results in Fig. 4 and show the $95\%$ confidence interval of a fair coin flip as a grey area. Except to Crack and Slope, all methods are insignificant w.r.t. the fair coin flip. In particular, Crack has an accuracy of over $90\%$ for the first $41\%$ of its decisions and reaches $77.2\%$ overall. Regarding the whole decision rate, Crack is nearly on par with Slope, which is as far as we know, the current state of the art for univariate continuous data.

Multivariate data

To test $\textsc{Crack}_{\delta}$ on multivariate mixed-type and single-type data, we collected 17 data sets. The information of the dimensionality for each data set is listed in Table 1. The first six data sets belong to the Tuebingen benchmark data set [19] and the next four were published by Janzing et al. [9]. Further, we extracted cause-effect pairs form the Haberman [6], Iris [3], Mammals [7] and Octet [4, 31] data sets. Those are described in more detail in the appendix.

We compare $\textsc{Crack}_{\delta}$ with LTR, Ergo and Origo. Ergo and LTR do not consider categoric data, and are hence not applicable on all data sets. In addition, LTR is only applicable to strictly multivariate data sets. $\textsc{Crack}_{\delta}$ is applicable to all data sets, infers $15/17$ causal directions correctly, by which it has an overall accuracy of $88.2\%$. Importantly, the two wrong decisions have low confidences compared to the correct inferences.

8 Conclusion

We considered the problem of inferring the causal direction from the joint distribution of two univariate or multivariate random variables $X$ and $Y$ consisting of single-, or mixed-type data. We point out weaknesses of known causal indicators and propose the Normalized Causal Indicator for mixed-type data and data with highly unbalanced domains. Further, we propose a practical encoding based on classification and regression trees to instantiate these causal indicators and provide a fast greedy heuristic to compute good solutions.

In the experiments we evaluate the advantages of the NCI and the common indicator and give advice on when to use them. On real world benchmark data, we are on par with the state of the art for univariate continuous data and beat the state of the art on multivariate data with a wide margin.

For future work, we aim to investigate in the application of Crack for causal discovery, meaning that we would like to infer causal networks. In addition, we only selected a subset of possible refinements to exploit dependencies from candidates. This choice could be expanded by considering more complex functions, finding combinations of categories for splitting. However, unless specific care is taken many of such extensions will likely have repercussions on the runtime of our algorithm, which is why besides being out of scope here, we leave this for future work.

Acknowledgements

The authors wish to thank Kailash Budhathoki for insightful discussions. Alexander Marx is supported by the International Max Planck Research School for Computer Science (IMPRS-CS). Both authors are supported by the Cluster of Excellence “Multimodal Computing and Interaction” within the Excellence Initiative of the German Federal Government.