Kernel-Based Independence Tests for Causal Structure Learning on Functional Data

1. Introduction

Uncovering the causal relationships between measured variables, a discipline known as causal structure learning or causal discovery, is of great importance across various scientific fields, such as climatology [1], economics [2], or biology [3]. Doing so from passively collected (‘observational’) data enables the inference of causal interactions between variables without performing experiments or randomised control trials, which are often expensive, unethical, or impossible to conduct [4]. Causal structure learning is the inference, under a given set of assumptions, of directed and undirected edges in graphs representing the data generating process, where the nodes represent variables and the inferred edges capture causal (directed) or non-causal (undirected) relationships between them.

Much research in various areas collates functional data consisting of multiple series of measurements that are observed conjointly over a given continuum (e.g., time, space, or frequency), where each series is assumed to be a realisation of an underlying smooth process [5, §3]. By viewing the series of measurements as discretisations of functions, the observations are not required to be collected over regular meshes of points along the continuum. If the variables are measured over time as the underlying continuum, there is a long history of methods that have been developed to infer (time-based) causality between variables. Among those, the classic Granger-causality [6] declares that variable ‘X causes Y’ ($X\to Y$) if predicting the future of Y becomes more accurate with, as compared to without, access to the past of X, conditional on all other relevant variables. However, these methods assume that the observed time-series are stationary, the causal dependency of X on Y is linear, and no confounding variables are considered. More recently, Sugihara et al. [7] developed convergent cross mappings (ccms), a method that relaxes the assumption of linearity and finds causal relationships based on time-embeddings of the (stationary) time-series at each point, while also assuming no confounding variables are present. While useful in many situations, as we will show, Granger-causality and ccm can perform weakly when the time-series for X and Y are non-stationary and related through non-linear relationships, meaning that the influence of X on Y varies over time. Here we present a method that uses kernel-based independence tests to detect statistically significant causal relationships by extending constraint- and regression-based causal structure learning to functional data. The key advantage over Granger-causality and ccm is the consideration of confounders and the relaxation of their assumptions around linear relationships or stationarity in the data, which can lead to different causal relationships between variables. As a motivating example, consider the relationship between two variables: ‘corruption’ and ‘income inequality’, as measured by the World Governance Indicators (wgis) [8] and the World Bank [9], respectively. Using data for 48 African countries from 1996 to 2016, Sulemana and Kpienbaareh [2] investigated their cause-effect relationship and found that corruption ‘Granger-causes’ income inequality. Independently, we have also confirmed that applying ccm to the same data leads to the same conclusion. However, by considering the time-series data as realisations of functions over time, and thus avoiding linearity and stationarity assumptions, our proposed kernel-based approach suggests the reverse result, i.e., causal influence of income inequality on corruption appears as the more statistically likely direction. This is in agreement with other quantitative findings, which draw on different data sources [10,11,12]. We will return to this example in Section 4.2.2 where we analyse causal dependencies between all six wgis.

Methodologically, we extend the applicability of two popular paradigms in causal structure learning—constraint-based [13, § 5] and regression-based methods [14]—to functional data. Independence tests play a crucial role in uncovering causal relationships in both paradigms, and kernels provide a powerful framework for such tests by embedding probability distributions in reproducing kernel Hilbert space (rkhss) [15, § 2.2]. Until now, however, related methods for causal learning had only been applicable to univariate and multivariate data, but not to functional data. Here, we address this limitation by using recently derived kernels over functions [16] to embed probability distributions over functions in rkhss, thus widening the applicability of kernel-based independence tests to functional data settings. To test for conditional independence, we compute hscic [17] in a conditional permutation test (cpt) [18], and we propose a straightforward search to determine the optimised regularisation rate in hscic.

We structure our paper as follows. Section 2 provides a brief overview of prior literature on functional data analysis, kernel-based independence tests, and causal structure learning methods. Section 3 presents our methodological contributions, from the definition of a conditional independence test for functional data to causal structure learning on such data. We then empirically analyse the performance of our independence tests and causal structure learning algorithms on synthetic and real-world data in Section 4. We conclude with a discussion in Section 5.

We highlight the following contributions:

• In Section 3 we propose a conditional independence test for functional data that combines a novel test statistic based on hscic with cpt to generate samples under the null hypothesis. The algorithm also searches for the optimised regularisation strength ${\lambda }^{*}$ required to compute hscic, by pre-test permutations and calculating an evaluation rejection rate.
• In Section 4.1.2, we extend the historical functional linear model [19] to the multivariate case ($\{X_1,\dots ,X_i,\dots ,X_d\}\to Y$) for regression-based causal structure learning, and we show how a joint independence test can be used to verify candidate directed acyclic graph (dags) [20, § 5.2] that embed the causal structure of function-valued random variables. This model has been contributed to the Python package scikit-fda [21].
• On synthetic data, we show empirically that our bivariate, joint and conditional independence tests achieve high test power, and that our causal structure learning algorithms outperform previously proposed methods.
• Using a real-world data set (World Governance Indicators) we demonstrate how our method can yield insights into cause-effect relationships amongst socio-economic variables measured in every country worldwide.
• Implementations of our algorithms are made available at https://github.com/felix-laumann/causal-fda/ in an easily usable format that builds on top of scikit-fda and causaldag [22].

2. Background and Related Work

2.1. Functional Data Analysis

In functional data analysis [5], a variable X is described by a set of n samples (or realisations), $\mathbf{X}={\left\{x_i\left(s\right)\right\}}_{i=1}^n$, where each functional sample $x_i\left(s\right)$ corresponds to a series of observations over the continuum s, also called the functional dimension. Typical functional dimensions are time or space. In practical settings, the observations are taken at a set of S discrete values $s_1,\dots ,s_S$ of the continuum variable s. Examples of functional data sets include the vertical position of the lower lip over time when speaking out a given word [19], the muscle soreness over the duration of a tennis match [23], or the ultrafine particle concentration in air measured over the distance to the nearest motorway [24].

In applications, the functional samples are usually represented as linear combinations of a finite set of M basis functions ${\left\{{\varphi }_m\left(s\right)\right\}}_{m=1}^M$ (e.g., Fourier or monomial basis functions):

$$x_i\left(s\right)=\sum _{m=1}^M{c}_{i,m}{\varphi }_m\left(s\right),$$

(1)

where the coefficients $C_i=(c_{i,1},c_{i,2},\dots ,c_{i,M})$ characterise each sample. If the number of basis functions is equal to the number of observations, $M=S$, each observed value $x_i\left(s_k\right)$ can be fitted exactly by obtaining the coefficients $C_i$ using standard non-linear least squares fitting techniques (provided the ${\varphi }_m\left(s\right)$ are valid basis functions), and Equation (1) allows us to interpolate between any two observations. When the number of basis functions is smaller than the number of observations, $M<S$, as it is commonly the case in practice, the basis function expansion (1) provides a smoothed approximation to the set of observations, ${\widehat{x}}_i\left(s_k\right)$.

For the many applications where the continuum is time, historical functional linear model (hflms) [19] provide a comprehensive framework to map the relationship between two sets of functional samples. Let 0 and T be the initial and final time points for a set of samples $y_i\left(t\right)$. hflms describe dependencies that can vary over time using the function $\beta (s,t)$, which encapsulates the influence of $x\left(s\right),s\in [0,S]$ on another variable $y\left(t\right),t\in [0,T]$ at any two points in time, $s_k$ and $t_k$:

$$y_i\left(t\right)={\int }_{s_0\left(t\right)}^t{x}_i\left(s\right)\beta (s,t)ds,t\in [0,T]$$

(2)

where $s_0\left(t\right)$ is the maximum allowable lag for any influence of X on Y and $s_0\left(t\right)\le s\le t$. Typical choices for $\beta (s,t)$ are exponential decay and hyperbolic paraboloid (or “saddle”) functions. The continuum is not required to be time but can also be space, frequency or others (see Ramsay and Hooker [25] for an extensive collection of function-to-function models and applications).

2.2. Kernel Independence Tests

Let ${\mathcal{H}}_{k_1}$ and ${\mathcal{H}}_{k_2}$ be separable rkhss with kernels $k_1:{\mathbb{R}}^{\mathbf{X}}\times {\mathbb{R}}^{\mathbf{X}}\to \mathbb{R}$ and $k_2:{\mathbb{R}}^{\mathbf{Y}}\times {\mathbb{R}}^{\mathbf{Y}}\to \mathbb{R}$ such that the tensor product kernel $(k_1\otimes k_2)(x,y,x^{\prime },y^{\prime })=k_1(x,x^{\prime })·k_2(y,y^{\prime })$ implies $({\mathbb{R}}^{\mathbf{X}}\times {\mathbb{R}}^{\mathbf{Y}})\times ({\mathbb{R}}^{\mathbf{X}}\times {\mathbb{R}}^{\mathbf{Y}})\to \mathbb{R}$. If the kernel $k_1$ uniquely embeds a distribution ${\mathbb{P}}_{\mathbf{X}}$ in an rkhs by a mean embedding,

$${\mu }_{\mathbf{X}}:={\int }_{\mathcal{X}}{k}_1(\mathbf{x},·)d{\mathbb{P}}_{\mathbf{X}}$$

(3)

which captures any information about ${\mathbb{P}}_{\mathbf{X}}$, we call $k_1$ a characteristic kernel on ${\mathbb{R}}^{\mathbf{X}}\times {\mathbb{R}}^{\mathbf{Y}}$ [26]. Characteristic kernels have thus been extensively used in bivariate (${\mathbb{P}}_{XY}={\mathbb{P}}_X{\mathbb{P}}_Y$), joint (${\mathbb{P}}_{XYZ}={\mathbb{P}}_X{\mathbb{P}}_Y{\mathbb{P}}_Z$) and conditional (${\mathbb{P}}_{XY|Z}={\mathbb{P}}_{X|Z}{\mathbb{P}}_{Y|Z}$) independence tests [e.g., 20,27,28].

For the bivariate independence test, let ${\mathbb{P}}_{\mathbf{X}\mathbf{Y}}$ denote the joint distribution of $\mathbf{X}$ and $\mathbf{Y}$. Then hsic is defined as

$$\begin{array}{cc}\hfill \mathrm{HSIC}({\mathscr{H}}_{k_1},{\mathscr{H}}_{k_2},{\mathbb{P}}_{\mathbf{X}\mathbf{Y}})& :=\parallel {\mu }_{\mathbf{X}\mathbf{Y}}-{\mu }_{\mathbf{X}}\otimes {\mu }_{\mathbf{Y}}{\parallel }_{{\mathscr{H}}_{k_1}\otimes {\mathscr{H}}_{k_2}}^2\ge 0,\hfill \\ \hfill & \mathrm{with}\mathrm{equality}\mathrm{iff}{\mathbb{P}}_{\mathbf{X}\mathbf{Y}}={\mathbb{P}}_{\mathbf{X}}{\mathbb{P}}_{\mathbf{Y}}.\hfill \end{array}$$

(4)

We refer to Gretton et al. [27] for the definition of an estimator for finite samples that constitutes the test statistic in the bivariate independence test with null hypothesis $H_0:X⫫Y$. The test statistic is then computed on the original data $(\mathbf{X},\mathbf{Y})$ and statistically compared to random permutations ${\left\{(\mathbf{X},{\mathbf{Y}}_{p^{\prime }})\right\}}_{p^{\prime }=1}^P$ under the null hypothesis.

For distributions with more than two variables, let ${\mathbb{P}}_{{\mathbf{X}}_1,{\mathbf{X}}_2,\dots ,{\mathbf{X}}_d}$ denote the joint distribution on ${\mathbf{X}}_1,{\mathbf{X}}_2,\dots ,{\mathbf{X}}_d$. To test for joint independence we compute

$$\begin{array}{cc}\hfill d\text{-HSIC}({\mathscr{H}}_{k_1},{\mathscr{H}}_{k_2},\dots ,& {\mathscr{H}}_{k_d},{\mathbb{P}}_{{\mathbf{X}}_1,{\mathbf{X}}_2,\dots ,{\mathbf{X}}_d}):=\hfill \\ \hfill & \parallel {\mu }_{{\mathbf{X}}_1,{\mathbf{X}}_2,\dots ,{\mathbf{X}}_d}-{\mu }_{{\mathbf{X}}_1}\otimes {\mu }_{{\mathbf{X}}_2}\otimes \cdots \otimes {\mu }_{{\mathbf{X}}_d}{\parallel }_{{\mathscr{H}}_{k_1}\otimes {\mathscr{H}}_{k_2}\otimes \cdots \otimes {\mathscr{H}}_{k_d}}^2\ge 0\hfill \\ \hfill & \mathrm{with}\mathrm{equality}\mathrm{iff}{\mathbb{P}}_{{\mathbf{X}}_1,{\mathbf{X}}_2,\dots ,{\mathbf{X}}_d}={\mathbb{P}}_{{\mathbf{X}}_1}{\mathbb{P}}_{{\mathbf{X}}_2}\cdots {\mathbb{P}}_{{\mathbf{X}}_d}.\hfill \end{array}$$

(5)

Pfister et al. [20] derive a numerical estimator, which serves as the basis for a joint independence test on finite samples. Here, the distribution under the null hypothesis of joint independence is generated by randomly permuting all sample sets in the same way as $\mathbf{Y}$ is in the bivariate independence test.

Lastly, the conditional independence test relies on accurately sampling from the distribution under the null hypothesis $H_0:X⫫Y|Z$. At the core of the conditional permutation test [18] lies a randomisation procedure that generates permutations of $\mathbf{X}$, denoted ${\left\{{\mathbf{X}}_{p^{\prime }}\right\}}_{p^{\prime }=1}^P$, which are generated without altering the conditional distribution ${\mathbb{P}}_{\mathbf{X}|\mathbf{Z}}$, so

$${\mathbb{P}}_{{\mathbf{X}}_{p^{\prime }},\mathbf{Y},\mathbf{Z}}={\mathbb{P}}_{\mathbf{X},\mathbf{Y},\mathbf{Z}}$$

(6)

under $H_0$ while breaking any dependence between X and Y. The null distribution can therefore be generated by repeating this procedure multiple times, and we can decide whether $H_0$ should be rejected by comparing a test statistic on the original data against its results on the generated null distribution.

The existing literature on kernel-based independence tests is extensive, see e.g., Berrett et al. [18] for a relevant review, but only a small part of those tests investigates independence among functional variables [29,30]. There have been particularly strong efforts in developing conditional independence tests and understanding their applicable settings. The kernel conditional independence test (kcit) [28], for example, gave promising results in use with univariate data but increasingly suffered when the number of conditional variables was large. In contrast, the kernel conditional independence permutation test (kcipt) [31] repurposed the well-established kernel two-sample test [32] to a conditional independence setting which delivered stable results for multiple conditional variables. However, Lee and Honavar [33] pointed out that as the number of permutations increases, while its power increases, its calibratedness decreases. This issue was overcome by their proposed self-discrepancy conditional independence test (sdcit), which is based on a modified unbiased estimate of the maximum mean discrepancy (mmd).

2.3. Causal Structure Learning

The aim of causal structure learning, or causal discovery, is to infer the qualitative causal relationships among a set of observed variables, typically in the form of a causal diagram or dag. Once learnt, such a causal structure can then be used to construct a causal model such as a causal Bayesian network or a structural causal model (scm) [34]. Causal models are endowed with a notion of manipulation and, unlike a statistical model, do not just describe a single distribution, but many distributions indexed by different interventions and counterfactuals. They can be used for causal reasoning, that is, to answer causal questions such as computing the average causal effect of a treatment on a given outcome variable. Such questions are of interest across many disciplines, and causal discovery is thus a highly topical area. We refer to Glymour et al. [4], Mooij et al. [35], Peters et al. [36], Schölkopf and von Kügelgen [37], Squires and Uhler [38], Vowels et al. [39] for comprehensive surveys and accounts of the main research concepts. In particular, we focus here on causal discovery methods for causally sufficient systems, for which there are no unobserved confounders influencing two or more of the observed variables. Existing causal discovery methods can roughly be categorised into three families:

• Score-based approaches assign a score, such as a penalised likelihood, to each candidate graph and then pick the highest scoring graph(s). A common drawback of score-based approaches is the need for a combinatorial enumeration of all dags in the optimisation, although greedy approaches have been proposed to alleviate such issues [40].
• Constraint-based methods start by characterising the set of conditional independences in the observed data [13]. They then determine the graph(s) consistent with the detected conditional independences by using a graphical criterion called d-separation, as well as the causal Markov and faithfulness assumptions, which establish a one-to-one connection between d-separation and conditional independence (see Appendix D for definitions). When only observational i.i.d. data are available, this yields a so-called Markov equivalence class, possibly containing multiple candidate graphs. For example, the graphs $X\to Y\to Z$, $X\leftarrow Y\leftarrow Z$, and $X\leftarrow Y\to Z$ are Markov equivalent, as they all imply $X⫫Z|Y$ and no other conditional independence relations.
• Regression-based approaches directly fit the structural equations $X_i:=f_i({\mathbf{PA}}_i,U_i)$ of an underlying scm for each $X_i$, where ${\mathbf{PA}}_i$ denote the parents of $X_i$ in the causal dag and $U=(U_1,...,U_n)$ are jointly independent exogenous noise variables. Provided that the function class of the $f_i$ is sufficiently restricted, e.g., by considering only linear relationships [41] or additive noise models [42], the true causal graph is identified as the unique choice of parents for each i such that the resulting residuals ${\widehat{U}}_i=X_i-{\widehat{f}}_i\left({\mathbf{PA}}_i\right)$ are jointly independent.

As can be seen from these definitions, conditional, bivariate, and joint independence tests are an integral part of constraint- and regression-based causal discovery methods. Our main focus in the present work is therefore to extend the applicability of these causal discovery frameworks to functional data by generalising the underlying independence tests to such domains.

3. Methods

In all three of our independence tests (bivariate, joint, conditional), we employ kernels over functions, also known as squared-exponential T (se-t) kernels [16]. Let $\mathcal{X}$ and $\mathcal{Y}$ be real, separable Hilbert spaces with norms ${\parallel ·\parallel }_{\mathcal{X}}$ and ${\parallel ·\parallel }_{\mathcal{Y}}$, respectively. Then, for $T:\mathcal{X}\to \mathcal{Y}$, se-t kernels are defined as

$$k_T(x,y)=e^{-\frac{1}{2{\sigma }^2}{\parallel T\left(x\right)-T\left(y\right)\parallel }_{\mathcal{Y}}^2}.$$

(7)

where ${\sigma }^2$ is commonly defined as the median heuristic, ${\sigma }^2={\mathrm{Median}\{\parallel T\left(a\right)-T\left(b\right)\parallel }_{\mathcal{Y}}^2:a,b\in {\left\{x_i\right\}}_{i=1}^n\cup {\left\{y_i\right\}}_{i=1}^n,a\ne b\}$. Replacing any characteristic kernel by the se-t kernels for bivariate independence tests (based on hsic) and for joint independence tests (based on d-variable Hilbert-Schmidt independence criterion ($d$-hsic)) is straightforward and does not require further theoretical investigation besides evaluating numerically the validity and power of the tests (see Section 4). However, the application of se-t kernels to conditional independence tests needs further theoretical results, as we discuss now.

3.1. Conditional Independence Test on Functional Data

We consider the conditional independence test, which generates samples under the null hypothesis based on the cpt, and uses the sum of hscics over all samples $z\in \mathbf{Z}$ as its test statistic. The cpt defines a permutation procedure that preserves the dependence of Z on both X and Y while resampling data for X that eliminates any potential dependence between X and Y. This procedure results in samples according to the null hypothesis, $X⫫Y|Z$. We use this procedure whenever permutations are required as part of the conditional independence test. Given the computation of the hscic is based on a kernel ridge regression, it requires to set a regularisation strength $\lambda$, which has been manually chosen previously [17]. Generally, our aim is to have type-I error rates close to the allowable false positive rate $\alpha$. However, choosing $\lambda$ inappropriately may result in an invalid test (type-I error rates exceed $\alpha$ if $\lambda$ is chosen too large), or in a deflated test power (type-I error rates are well below $\alpha$ and type-II error rates are high if $\lambda$ is chosen too small). Thus, we must define an algorithm that conducts a search over a range of potentially suitable values for $\lambda$ and assesses each candidate value by, what we will call, an evaluation rejection rate.

The search proceeds by iterating over a range of values ${\left\{{\lambda }_l\right\}}_{l=1}^L$ to find the optimised value ${\lambda }^{*}$, as follows. For each ${\lambda }_l$, we start by producing one permutation of the samples $\mathbf{X}$, which we denote as ${\mathbf{X}}_b$, and we compute its corresponding evaluation test statistic given by the sum of hscics over all samples $z\in \mathbf{Z}$. We then apply the usual strategy for permutation-based statistical tests: we produce an additional set of P permuted sample sets of $\mathbf{X}$, which we denote ${\left\{{\mathbf{X}}_{\pi }\right\}}_{\pi =1}^P$, and for each ${\mathbf{X}}_{\pi }$, we determine the sum of hscics over $z\in \mathbf{Z}$ to generate a distribution over statistics under the null hypothesis, which we call the evaluation null statistics. We then compute the percentile where the evaluation test statistic on ${\mathbf{X}}_b$ falls within the distribution of evaluation null statistics on the permutation set ${\left\{{\mathbf{X}}_{\pi }\right\}}_{\pi =1}^P$. This results in an evaluation p-value which is compared to the allowable false positive rate $\alpha$ to establish if the hypothesis $H_0:{\mathbf{X}}_b\cancel{⫫}\mathbf{Y}|\mathbf{Z}$ can be rejected. Given that both the evaluation test statistic and the evaluation null statistics are computed on conditionally independent samples, we repeat this procedure for $b=1,\dots ,B\ge 100$ times to estimate an evaluation rejection rate for each value of ${\lambda }_l$. Having completed this procedure over all values ${\left\{{\lambda }_l\right\}}_{l=1}^L$, we select the ${\lambda }_l$ that produces an evaluation rejection rate closest to $\alpha$ as the optimised regularisation strength, ${\lambda }^{*}$. Finally, we apply a cpt-based conditional independence test using the optimised ${\lambda }^{*}$ to test the null hypothesis $H_0:X\cancel{⫫}Y|Z$. This entire procedure is summarised in Algorithm 1.

Algorithm 1 Search for ${\lambda }^{*}$ with subsequent conditional independence test

Require:  Samples ($\mathbf{X},\mathbf{Y},\mathbf{Z}$), Range ${\left\{{\lambda }_l\right\}}_{l=1}^L$, Significance level $\alpha$, Permutation iterations P, Rejection iterations B Initialise:  $\text{evaluation}\text{rejection}\text{rate}\left[{\lambda }^{*}\right]=0$ 1: for $1\le l\le L$ do                                                                                                                                                                                                                                  ▹Start: Search for ${\lambda }^{*}$ 2:     for $1\le b\le B$ do 3:         Permute $\mathbf{X}$ by cpt, call them ${\mathbf{X}}_b$ 4:         $\mathrm{evaluation}\mathrm{test}\mathrm{statistic}\left[b\right]\leftarrow {\sum }_{z\in \mathbf{Z}}\mathrm{HSCIC}({\mathbf{X}}_b,\mathbf{Y},\mathbf{Z},{\lambda }_l)$ 5:         for $1\le \pi \le P$ do 6:            Permute ${\mathbf{X}}_b$ by cpt, call them ${\mathbf{X}}_{b,\pi }$ 7:            $\mathrm{evaluation}\mathrm{null}\mathrm{statistics}[b,\pi ]\leftarrow {\sum }_{z\in \mathbf{Z}}hscic({\mathbf{X}}_{b,\pi },\mathbf{Y},\mathbf{Z},{\lambda }_l)$ 8:         $\mathrm{evaluation}p\text{-value}\leftarrow 1-$ 9:                     $\mathrm{percentile}\left(\mathrm{evaluation}\mathrm{test}\mathrm{statistic}\right[b],\mathrm{evaluation}\mathrm{null}\mathrm{statistics}[b,:\left]\right)$ 10:         if $\mathrm{evaluation}\mathrm{p}-\mathrm{value}\ge \alpha$ then 11:            Fail to reject $H_0:{\mathbf{X}}_b⫫\mathbf{Y}|\mathbf{Z}$ 12:            $\mathrm{rejects}\left[b\right]=0$ 13:         else 14:            Reject $H_0:{\mathbf{X}}_b\cancel{⫫}\mathbf{Y}|\mathbf{Z}$ 15:            $\mathrm{rejects}\left[b\right]=1$ 16:     $\mathrm{evaluation}\mathrm{rejection}\mathrm{rate}\left[{\lambda }_l\right]\leftarrow {\mathrm{mean}}_{1\le b\le B}\left(\mathrm{rejects}\right)$ 17:     if $|\mathrm{evaluation}\mathrm{rejection}\mathrm{rate}\left[{\lambda }^{*}\right]-\alpha |\ge |\mathrm{evaluation}\mathrm{rejection}\mathrm{rate}\left[{\lambda }_l\right]-\alpha |$ then 18:         ${\lambda }^{*}\leftarrow {\lambda }_l$ 19: return ${\lambda }^{*}$                                                                                                                                                                                                                                         ▹End: Search for ${\lambda }^{*}$       20: $\mathrm{test}\mathrm{statistic}\leftarrow {\sum }_{z\in \mathbf{Z}}\text{HSCIC}(\mathbf{X},\mathbf{Y},\mathbf{Z},{\lambda }^{*})$                                                                                                                                                    ▹Start: Conditional independence test 21: for $1\le p^{\prime }\le P$  do 22:     Permute $\mathbf{X}$ by cpt, call them ${\mathbf{X}}_{p^{\prime }}$ 23:     $\mathrm{null}\mathrm{statistics}\left[p^{\prime }\right]\leftarrow {\sum }_{z\in \mathbf{Z}}\text{HSCIC}({\mathbf{X}}_{p^{\prime }},\mathbf{Y},\mathbf{Z},{\lambda }^{*})$ 24: $p\text{-value}\leftarrow 1-\mathrm{percentile}(\mathrm{test}\mathrm{statistic},\mathrm{null}\mathrm{statistics})$ 25: if  $p\text{-value}\ge \alpha$  then 26:     Fail to reject $H_0:X⫫Y|Z$ 27: else 28:     Reject $H_0:X\cancel{⫫}Y|Z$                                                                                                                                                                                                                                   ▹End: Conditional independence test

The following Theorem 1 guarantees the consistency of the conditional independence test in Algorithm 1 with respect to the regularisation parameter $\lambda$.

Theorem 1.

Let ${\mathcal{H}}_{k_1}$ and ${\mathcal{H}}_{k_2}$ be separable RKHSs with kernels $k_1:{\mathbb{R}}^{\mathbf{X}}\times {\mathbb{R}}^{\mathbf{X}}\to \mathbb{R}$ and $k_2:{\mathbb{R}}^{\mathbf{Y}}\times {\mathbb{R}}^{\mathbf{Y}}\to \mathbb{R}$ such that the tensor product kernel $k_1\otimes k_2:({\mathbb{R}}^{\mathbf{X}}\times {\mathbb{R}}^{\mathbf{Y}})\times ({\mathbb{R}}^{\mathbf{X}}\times {\mathbb{R}}^{\mathbf{Y}})\to \mathbb{R}$ is a characteristic kernel on ${\mathbb{R}}^{\mathbf{X}}\times {\mathbb{R}}^{\mathbf{Y}}$ [43]. If the regularisation parameter $\lambda ={\lambda }_n$ decays as $n\to \infty$ at a slower rate than $n^{-1/2}$, then the test based on the test statistic in Algorithm 1 (lines 19–27) is consistent.

Proof. 

See Appendix B.    □

3.2. Causal Structure Learning on Functional Data

To infer the existence of (directed) edges in an underlying causal graph G based on the joint data distribution $\mathbb{P}$ over a set of observed variables, we must assume that G and $\mathbb{P}$ are intrinsically linked. Spirtes et al. [13, § 2.3.3] define the faithfulness assumption, from which it follows that if two random variables are (conditionally) independent in the observed distribution $\mathbb{P}$, then they are d-separated in the underlying causal graph G [34, § 1.2.3]. Using this fact, constraint-based causal structure learning methods [13, § 5] take advantage of bivariate and conditional independence tests to infer whether two nodes are d-separated in the underlying graph. These methods yield completed partially directed acyclic graphs (cpdags), which are graphs with undirected and/or directed edges. In contrast, regression-based methods [14], utilise the joint independence test to discover dags that have all edges oriented. We now describe both of these approaches in more detail.

3.2.1. Constraint-Based Causal Structure Learning

Constraint-based causal structure learning relies on performing conditional independence tests for each pair of variables, X and Y, conditioned on any possible subset of the remaining variables, that is, any subset within the power set of the $d\setminus \{X,Y\}$ remaining variables $Z_1,Z_2,\dots ,Z_{*}$. Therefore, for d variables we need to carry out $2^{d-2}$ conditional independence tests for every pair of variables, which results in $\left(\genfrac{}{}{0pt}{}{d}{2}\right)2^{d-2}$ tests in total.

Conditional dependences found in the data can then be used to delete and orient edges in the graph G, as follows. We start with a complete graph with d nodes, where each node corresponds to one of the variables. The edge connecting X and Y is deleted if there exists a subset ${\mathcal{Z}}_1$ of the $d-2$ remaining variables such that $Z_{*}\notin {\mathcal{Z}}_1$ and $X⫫Y|{\mathcal{Z}}_1$. If we then find another subset ${\mathcal{Z}}_2$, such that $Z_{*}\in {\mathcal{Z}}_2$, ${\mathcal{Z}}_1={\mathcal{Z}}_2\setminus \left\{Z_{*}\right\}$, and for which $X\cancel{⫫}Y|{\mathcal{Z}}_2$, we can orient the edges $X-Z$ and $Y-Z$ to form colliders (or v-structures), $X\to Z_{*}\leftarrow Y$.

Based on these oriented edges, ‘Meek’s orientation rules’ defined by Meek [44] [see e.g., 36, § 7.2.1] can be applied to direct additional edges based on certain graphical compositions, as shown in Appendix F. Briefly, these rules follow from the avoidance of cycles, which would violate the acyclicity of a dag, and colliders, which violate the found conditional independence. Algorithms that implement these criteria are the SGS algorithm and the more efficient PC algorithm [13], which we summarise in Appendix F.

3.2.2. Regression-Based Causal Structure Learning

To carry out our regression-based causal learning, we choose additive noise models (anms), a special class of scms where the noise terms are assumed to be independent and additive. This assumption guarantees that the causal graph can be identified if the function $f_i$ is non-linear with Gaussian additive noise [42,45], a typical setting in functional data. Our model over the set of random variables ${\left\{X\right\}}_{i=1}^d$ is given by:

$${\mathbf{X}}_i:=f_i\left({\mathbf{PA}}_i\right)+U_i,$$

(8)

where the additive noise terms $U_i$ are jointly independent, i.e., any noise term $U_i$ is independent of any intersection of the other noise terms,

$$\mathbb{P}\left({\cap }_{j=1}^k{U}_{i_j}\right)=\prod _{j=1}^k\mathbb{P}\left(U_{i_j}\right),\forall k\le d.$$

Based on these assumptions, we follow resit [Regression with Subsequent Independence Testing; 14], an approach to causal discovery that can be briefly described as follows. If we only have two variables X and Y, we start by assuming that the samples in $\mathbf{Y}$ are the effect and the samples in $\mathbf{X}$ are the cause. Therefore we write $\mathbf{Y}$ as a function ${\widehat{f}}_{\mathbf{Y}}$ of $\mathbf{X}$ plus some noise, and we test whether the residual $r_{\mathbf{Y}}=\mathbf{Y}-{\widehat{f}}_{\mathbf{Y}}\left(\mathbf{X}\right)$ is independent of $\mathbf{X}$. We then exchange the roles of assumed effect and assumed cause to obtain the residual $r_{\mathbf{X}}=\mathbf{X}-{\widehat{f}}_{\mathbf{X}}\left(\mathbf{Y}\right)$, which is tested for independence from $\mathbf{Y}$. To overcome issues with finite samples and test power, Peters et al. [14] used the p-value of the independence tests as a measure of strength of independence, which we follow in our experiments in Section 4. Alternatively, one may determine the causal direction for this bivariate case by comparing both cause-effect directions with respect to a score defined by Bühlmann et al. [46].

For d variables, the joint independence test replaces the bivariate independence test. Firstly, we consider the set of candidate causal graphs, $\mathsf{\Gamma }=\{G_1,\dots ,G_{\gamma },\dots ,G_{\mathsf{\Gamma }}\}$, which contains every potential dag with d nodes. For each candidate dag, $G_{\gamma }$, we regress each descendant ${\mathbf{X}}_i$ on its parents ${\mathbf{PA}}_i$:

$${\mathbf{X}}_i:=\sum _{k\in {\mathbf{PA}}_i}{\widehat{f}}_{i,k}\left({\mathbf{X}}_k\right)+U_i,i\in 1,\dots ,d,$$

(9)

and we compute the residuals $r_{{\mathbf{X}}_i}={\mathbf{X}}_i-{\sum }_{k\in {\mathbf{PA}}_i}{\widehat{f}}_{i,k}\left({\mathbf{X}}_k\right)$. We then apply the joint independence test to all d residuals. The candidate dag $G_{\gamma }$ is accepted as the ‘true’ causal graph if the null hypothesis of joint independence amongst the residuals is not rejected. This process is repeated over all candidate causal graphs in the set $\mathsf{\Gamma }$. Again, because in finite samples this procedure may not always lead to only one candidate dag being accepted, the one with the highest p-value is chosen [14].

4. Experiments

In practice, functional data analysis requires samples of continuous functions evaluated over a mesh of discrete observations. In our synthetic data, we consider the interval $s=[0,1]$ for our functional dimension and draw 100 equally spaced points over s to evaluate the functions. For real-world data, we map the space in which the data live (e.g., the years 1996 to 2020) to the interval $s=[0,1]$ and interpolate over the discrete measurement points. We then use this interpolation to evaluate the functional samples on 100 equally spaced points. Unless otherwise mentioned, henceforth we use se-t kernels with $T=I$ where I is the identity matrix and ${\sigma }^2$ as the median heuristic, as previously defined in Section 3.

4.1. Evaluation of Independence Tests for Functional Data

Before applying our proposed independence tests to causal structure learning problems, we use appropriate synthetic data to evaluate their type-I error rate (test size) and type-II error rate (test power). Type-I errors are made when the samples are independent (achieved by setting $a=0$ in our experiments below) but the test concludes that they are dependent, i.e., the test falsely rejects the null hypothesis $H_0$. In contrast, type-II errors appear when the samples are dependent but the test fails to reject $H_0$ (achieved by setting $a>0$). Specifically, we compute the error rates on 200 independent trials, which correspond to different random realisations of the data sets, i.e., synthetically generated data sets with random Fourier basis coefficients, random coefficients of the $\beta$-functions and random additive noise. We set the test significance level at $\alpha =0.05$ and approximate the null hypothesis of independence by 1000 permutations using the cpt scheme described above. Note that although our synthetic data are designed to replicate typical behaviour in time-series, the independence tests are not limited to time as the functional dimension, and can be used more generally.

4.1.1. Bivariate Independence Test

We consider n functional samples ${\left\{x_i\left(s\right)\right\}}_{i=1}^n$ of a random variable X defined over the interval $s=[0,1]$. To generate our samples, we sum $M=3$ Fourier basis functions ${\varphi }_{m,T}\left(s\right)$ (definitions in Appendix A.1) with period $T=0.1$ and coefficients $c_{i,m}$ randomly drawn from a standard normal distribution:

$$x_i\left(s\right)=\sum _{m=1}^M{c}_{i,m}{\varphi }_{m,T}\left(s\right).$$

(10)

To mimic real-world measurements, we draw each sample over a distinct, irregular mesh of evaluation points. We then interpolate them with splines, and generate realisations over a regular mesh of points. Multivariate random noise ${ϵ}_i^X\left(s\right)\sim \mathcal{N}(0,I)$, where I is the identity matrix, is then added. The samples ${\left\{y_i\left(t\right)\right\}}_{i=1}^n$ of random variable Y are defined as a hflm[19] over $t=[0,1]$ by:

$$y_i\left(t\right)=a{\int }_0^t{x}_i\left(s\right)\beta (s,t)ds,$$

(11)

where $a\in [0,1]$. The function $\beta (s,t)$ maps the dependence from $x_i$ at any point s to $y_i$ at any point t, and is defined here as a hyperbolic paraboloid function:

$$\begin{array}{c}\hfill \beta (s,t)=8{(s-c_1)}^2-8{(t-c_2)}^2,\end{array}$$

(12)

with coefficients $c_1$ and $c_2$ drawn from uniform distributions (see Appendix A for lower and upper bounds). Afterwards, the samples $y_i\left(t\right)$ are evaluated over a regular mesh of 100 points $t\in [0,1]$ and random noise ${ϵ}_i^Y\left(t\right)\sim \mathcal{N}(0,I)$ is added. Clearly, for $a=0$ our samples are independent, as the samples $y_i\left(t\right)={ϵ}_i^Y\left(t\right)$ are just random noise. As a is increased, the dependence between the samples $x_i\left(s\right)$ and $y_i\left(t\right)$ becomes easier to detect. Figure 1 shows that our test has no difficulties in staying within the allowable false-positive rate $\alpha$ and quickly detects the dependence as soon as $a>0$ even with very low number of samples.

4.1.2. Joint Independence Test

To produce the synthetic data for joint independence tests, we first generate random dags with d nodes using an Erdös-Rényi model with density 0.5 [20， § 5.2]. For each dag, we first use Equation (10) to produce function-valued random samples for the variables $X^i$ without parents. We then generate samples for other variables using a historical functional linear model:

$$x_i^j\left(t\right)=a\sum _{p\in {\mathbf{PA}}^j}{\int }_0^t{x}_i^p\left(s\right){\beta }^p(s,t)ds,$$

(13)

where ${\mathbf{PA}}^j$ are the parents of node $X^j$, and the function ${\beta }^p(s,t)$ is given by Equation (12) with random coefficients $c_1$ and $c_2$ independently generated for each descendant-parent pair indexed by p. After being evaluated at a regular mesh of 100 points within $t\in [0,1]$, random noise ${ϵ}_i^j\left(t\right)\sim \mathcal{N}(0,I)$ is added. Note that, again, an increase in the factor $a\in [0,1]$ should make the dependence structure of the dag easier to detect. Figure 1 shows the test size where $a=0$, resulting in independent variables, and test power where $a>0$. We evaluate the joint independence test for $d=4$ variables over various values of a and a range of sample sizes.

4.1.3. Conditional Independence Test

To evaluate whether X is independent of Y given $\mathcal{Z}$, where $\mathcal{Z}$ may be any subset within the power set of the $d\setminus \{X,Y\}$ remaining variables, we generate data samples according to:

$$\begin{array}{cc}\hfill z_i^j\left(u\right)& =\sum _{m=1}^M{c}_{i,m}^j{\varphi }_{m,T}\left(u\right),j=1,\dots ,\left|\mathcal{Z}\right|\hfill \\ \hfill x_i\left(s\right)& =\sum _{p\in {\mathbf{PA}}^X}{\int }_0^s{z}_i^p\left(u\right){\beta }^p(u,s)du\hfill \\ \hfill y_i\left(t\right)& =\sum _{p\in {\mathbf{PA}}^Y}{\int }_0^t{z}_i^p\left(u\right){\beta }^p(u,t)du+a^{\prime }{\int }_0^t{x}_i\left(s\right){\beta }^Y(s,t)ds,\hfill \end{array}$$

(14)

where $\left|\mathcal{Z}\right|$ is the cardinality of the set $\mathcal{Z}$, $\beta$’s are given by Equation (12), and noise terms ${ϵ}_i^{Z^j}\left(u\right)$, ${ϵ}_i^X\left(s\right)$ and ${ϵ}_i^Y\left(t\right)$, are added to the discrete observation values.

We then apply Algorithm 1 to compute the statistics for our conditional independence test. Firstly, we find that the optimised regularisation strength ${\lambda }^{*}$ is robust and reproducible for different random realisations with the same sample size and model parameters. Consequently, we search for one optimised ${\lambda }^{*}$ (line 1–18 in Algorithm 1) for each sample size, and fix this value for all 200 independent trials. Figure 2 summarises the results for ${\lambda }^{*}$ for $n\in \{100,200,300\}$ samples after conducting a grid search over a range of 18 possible values ${10}^{-4}\le \lambda \le {10}^{-1}$, with $P=1000$ permutations, and $B=100$ rejection iterations. We find that the optimised ${\lambda }^{*}$ exhibits a saturation as the number of conditional variables (the dimension of the conditional set) is increased. Given that we perform a kernel ridge regression, a larger number of samples n should result in lower requirements for regularisation—which aligns with our observations of a decreasing ${\lambda }^{*}$ over the increase in number of samples n. Therefore, Algorithm 1 optimises the regularisation parameter ${\lambda }^{*}$ such that the evaluation rejection rate of the test is closest to the allowable false positive rate $\alpha =0.05$.

Figure 3a–d show the results of the conditional independence test for increasing dimension d of the conditional set (i.e., number of conditional variables). We find that the test power is well preserved as d increases through a concomitant increase of ${\lambda }^{*}$ that partly ameliorates the “curse of dimensionality” [47]. Furthermore, the values of the test power for $a^{\prime }=0$ correspond to the type-I error rates.

4.2. Causal Structure Learning

We use the three independence tests evaluated numerically in Section 4.1 to learn different causal structures among function-valued random variables. We start with regression-based methods for the bivariate case $X\to Y$ and extend to the multivariate case through regression-, constraint-based, and a combination of constraint- and regression-based causal structure learning approaches.

4.2.1. Synthetic Data

Regression-based causal discovery. The synthetic data for a bivariate system $X\to Y$ is generated according to Equations (10)–(11) with $a=1$, and we measure the performance of the method by the structural Hemming distance (shd), whose definition can be found in Appendix E. As seen by the low value of the shd over 200 independent trials, Figure 4a shows that for this case our regression-based algorithm using resit outperforms two well-known algorithms for causal discovery: Granger-causality, which assumes linearity in the relationships [6], and ccm, which allows for non-linearity through a time-series embedding [7].

We also evaluate how the regression-based approach performs in a system of three variables where data are generated according to random dags with three nodes that entail historical non-linear dependence. The above methods (Granger, ccm) are commonly applied to sets with two variables only, hence we compare our method on a system with three variables to PC algorithm with momentary conditional independence (pcmci), an algorithm where each time point is represented as a node in a graph [48] from which we extract a directed graph for the variables. Figure 4b also shows improved performance of our regression-based algorithm, with a reduced shd. Details of the comparing methods are described in Appendix G.

Constraint-based causal discovery. The accuracy of our constraint-based approach is evaluated by computing three metrics: normalised shd ($sh{d}_{\mathrm{norm}}$), precision, and recall (see Appendix E for definitions). Figure 5 shows the accuracy of our algorithm when applied to $d\in \{3,4,5,6\}$ variables.

Combined approach. To be able to scale up our causal discovery techniques more efficiently to larger graphs, it is also possible to combine constraint- and regression-based causal learning, yielding dags. In this “combined” approach we start with a constraint-based step through which we learn cpdags by applying conditional independence tests to any two variables conditioned on any subset of the remaining $d-2$ variables and orient edges according to the Algorithm in Appendix F, which finds v-structures, applies Meek’s orientation rules and returns the Markov equivalence class. Often, the Markov equivalence class entails undirected edges, but we can then take a second step using a regression-based approach, under a different set of assumptions, to orient edges by applying resit. This two-step process finally yields a dag, i.e., every edge in the graph is oriented, in a more scalable manner than applying regression-based causal discovery directly. We then measure the accuracy of our approach computing the normalised shd (A4) as above. Figure 6 shows that our results for $d=3,4,5,6$ variables compare favourably to pcmci.

4.2.2. Real-World Data

To showcase the application of our methods to real data, we return now to the motivating example on socio-economic indicators from the Introduction (Section 1). Sulemana and Kpienbaareh [2] tested for the causal direction between corruption and income inequality, as measured by the Control of Corruption index [8] and Gini coefficient [49], respectively. They applied Granger-causality with Wald-${\chi }^2$ tests, which statistically test against the null hypothesis that an assumed explanatory variable does not have an influence on other considered variables. Both variables, corruption and income inequality, were tested against being the explanatory variable (i.e., the cause) in their bivariate relationship. Their findings showed that income inequality as the explanatory variable results in a p-value of $0.079$, whereas corruption is more likely to be the cause with a p-value of $0.276$. Analysing the same data using ccm, which does not assume linearity, we also find the same cause-effect direction: $70.9\%$ of the sample pairs $(x_i,y_i)$ validate corruption as the cause, against $29.7\%$ of pairs which detect the opposite direction. However, using our proposed regression-based kernel approach, which does not assume linearity or stationarity, we find the opposite result, i.e., the causal dependence of corruption on income inequality is statistically more likely (p $=0.0859$) than the reverse direction (p $=0.0130$).

Going beyond pairwise interactions, we have also applied our causal structure learning method to the set of World Governance Indicators (wgis) [8]. The wgis are a collection of six variables, which have been measured in 182 countries over 25 years from 1996 to 2020 (see Appendix H). We view the time-series of the countries as independent functional samples of each variable (so that $n=182$ and the functional dimension is time), and we apply our methods as described above to establish the causal structure amongst the six wgis. In this case, we use the “combined” approach by successively applying our constraint- and regression-based causal structure learning methods. We first learn an undirected causal skeleton from the wgis (Figure 7a), and we find a triangle between Government Effectiveness (GE), Rule of Law (RL) and Control of Corruption (CC), and a separate link between Voice and Accountability (VA) and Regulatory Quality (RQ). We then orient these edges (Figure 7b) and find that Government Effectiveness (GE) causes both Rule of Law (RL) and Control of Corruption (CC), and Rule of Law (RL) causes Control of Corruption (CC). We also find that Voice and Accountability (VA) causes Regulatory Quality (RQ).

5. Discussion and Conclusion

We present a causal structure learning framework on functional data that utilises kernel-based independence tests to extend the applicability of the widely used regression- and constraint-based approaches. The foundation of the framework originates from the functional data analysis literature by interpreting any discrete-measured observations of a random variable as finite-dimensional realisations of functions.

Using synthetic data, we demonstrate that our regression-based approach significantly outperforms existing methods such as Granger-causality, ccm and pcmci when learning causal relationships between two and three variables. Further, we show that our conditional independence test, which is the cornerstone of the constraint-based causal discovery approach, achieves type-I error rates close to the acceptable false-positive rate $\alpha$ and high type-II error rates, even when the number of variables in the conditional set increases. Shah et al. [47] rightly state that any conditional independence test can suffer from an undesirable test size or low test power in finite samples, and our method is obviously no exception. However, Figure 3 demonstrates the counterbalance of the “curse of dimensionality” through the optimised regularisation strength ${\lambda }^{*}$. Indeed, while with larger numbers of conditional variables, we would generally expect the test power to diminish, the optimisation of ${\lambda }^{*}$ offsets this reduction, resulting in no significant decrease in test power with a growing number of conditional variables. Although our proposed method is computationally expensive and significantly benefits from high sample sizes, a suitable regularisation strength could be chosen for large conditional sets in principle.

Furthermore, we demonstrate that constraint- and regression-based causal discovery methods can be combined to learn dags with large number of nodes, which would otherwise be computationally very expensive when relying on regression-based methods to yield dags. By comparing Figure 5a and Figure 6 we see however that it does not necessarily result in a lower shd. After learning the Markov equivalence class through the constraint-based approach, we apply regression with subsequent independence test (resit) to orient undirected edges in the “combined” approach. Here, mistakes in the orientation of undirected edges add 2 to the shd whereas an undirected edge only adds 1 to the shd.

The presented work contributes to opening the field of causal discovery to functional data. Beyond the results presented here, we believe that more research needs to be conducted to (i) increase the efficiency of the independence tests, meaning smaller sample sets can achieve higher test power; (ii) learn about upper bounds of the regularisation strength ${\lambda }^{*}$ with respect to the size of the sample and conditional sets; and (iii) make the computations in the conditional independence test less expensive.

Appendix A.1. Fourier Basis Functions

Throughout our paper, we use the following Fourier basis functions:

$$\begin{array}{cc}\hfill {\varphi }_0& =\frac{1}{\sqrt{2}}\hfill \\ \hfill {\varphi }_{2m-1,T}\left(s\right)& =\frac{sin\left(\frac{2\pi m}{T}s\right)}{\sqrt{\frac{T}{2}}}\hfill \\ \hfill {\varphi }_{2m,T}\left(s\right)& =\frac{cos\left(\frac{2\pi m}{T}s\right)}{\sqrt{\frac{T}{2}}}\hfill \end{array}$$

Appendix A.2. Lower and Upper Bounds of β-Functions

Table A1. Lower and upper bounds of uniform distribution that define $c_1$ and $c_2$ in $\beta$.

	Lower Bound	Upper Bound
$c_1$	0.25	0.75
$c_2$	0.25	0.75

Table A1. Lower and upper bounds of uniform distribution that define $c_1$ and $c_2$ in $\beta$.

	Lower Bound	Upper Bound
$c_1$	0.25	0.75
$c_2$	0.25	0.75

Figure A1. Function $\beta (s,t)$ generated with random coefficients $c_1$ and $c_2$ as specified in Table A1.

Figure A1. Function $\beta (s,t)$ generated with random coefficients $c_1$ and $c_2$ as specified in Table A1.

Appendix D.1. D-Separation

If $X,Y$ and Z are three pairwise disjoint (sets of) nodes in a dagG, then Z d-separates X from Y in G if Z blocks every possible path between (any node in) X and (any node in) Y in G [34, § 1.2.3]. We then write $X{\perp }_{G}Y|Z$. The symbol ${⫫}_G$ is commonly used to denote d-separation.

Appendix D.2. Faithfulness Assumption

If two random variables are (conditionally) independent in the observed distribution $\mathbb{P}$, then they are d-separated in the underlying causal graph G [Definition 6.33 [36].

Appendix D.3. Causal Markov Assumption

The distribution ${\mathbb{P}}_{XYZ}$ is Markov with respect to a dagG if $X{⫫}_{G}Y|Z\Rightarrow X⫫Y|Z$ for all disjoint (sets of) nodes $X,Y,Z$, where ${⫫}_G$ denotes d-separation [Definition 6.21 [36].

Appendix E.1. shd

For two graphs, $G_1$ and $G_2$, shd is the sum of the element-wise absolute differences between the adjacency matrices $A_1$ and $A_2$ of $G_1$ and $G_2$, respectively.

$$\mathrm{SHD}=\sum _i|A_1^i-A_2^i|,$$

(A3)

where we iterate over each element i in the adjacency matrices. Thus, when a learnt causal graph $G_1$ including an edge, e.g., $X\to Y$ is oriented opposite to the corresponding edge in the true graph $G_2$, meaning $X\leftarrow Y$ is true, we add 2 to the shd score (“double penalty”). Note that others [e.g., 51] only penalise by 1 for a learnt edge being directed opposite to the true edge.

Appendix E.2. Normalised shd

We normalise the shd by the number of possible directed edges in a graph of size d, thus allowing comparison across graphs with different number of nodes.

$${\mathrm{SHD}}_{\mathrm{norm}}=\frac{shd}{d(d-1)}.$$

(A4)

Appendix E.3. Precision

We define the total number of edges in a learnt graph $G_1$ as the sum of truly and falsely learnt edges, $e=e_t+e_f$. The truly learnt edges are the edges that correspond to the zero elements in the matrix containing the element-wise differences between the learnt $A_1$ and the true $A_2$. In contrast, the falsely learnt edges are the non-zero elements. The total number of truly and falsely learnt edges are $e_t$ and $e_f$, respectively. Precision is the proportion of correctly oriented edges $e_t$ in a learnt graph $G_1$ in comparison to a true graph $G_2$ out of all learnt edges e in $G_1$.

$$\mathrm{Precision}=\frac{e_t}{e}$$

(A5)

Appendix E.4. Recall

We separate between oriented and oriented edges, $e_o$ and $e_u$, respectively, which can again be summed up to the total number of edges, $e=e_o+e_u$. Recall is the fraction of edges in the cpdag that are oriented:

$$\mathrm{Recall}=\frac{e_o}{e}.$$

(A6)

Appendix G.1. Granger-Causality

We implement Granger-causality through software provided by Seabold and Perktold [52] and define the optimal lag as the first local minimum in the mutual information function over measurement points. Given that Granger-causality is applied to a pair of samples, $(x_i,y_i)$, we iterate over every pair and determine the percentage of pairs that result in $X\to Y$ against $X\leftarrow Y$.

Appendix G.2. ccm

We implement ccm through software provided by 53 and define the optimal lag as for Granger-causality. The embedding dimensionality is determined by a false nearest neighbour method [53]. As with Granger-causality, ccm is also applied to a pair of samples, $(x_i,y_i)$. We take the same summarising approach and iterate over every pair and determine the percentage of of pairs that result in $X\to Y$ against $X\leftarrow Y$.

Appendix G.3. pcmci

Runge et al. [48] proposed pcmci, a method that produces time-series graphs where each node represents a variable at a certain state in time. When applying the method in our experiments, we estimate the maximum possible lag as the average of the first local minimum of mutual information of all considered data series. Then, we test for the presence of an edge by applying distance correlation-based independence tests [54] between the residuals of X and Y that remain after regressing out the influence of the remaining nodes in the graph through Gaussian processes. The distance correlation coefficients of the significant edges between two variables are summed up to measure the prevalence of a causal direction of edges connecting these two variables. The causal direction with the greater sum is then considered to be the directed edge, and assessed against the edge between the same two variables in the true dag.