A Bayesian Model for Bivariate Causal Inference

TL;DR

This paper introduces Bayesian Causal Inference (BCI), a novel Bayesian hierarchical model for determining causal direction between two variables from observational data, demonstrating competitive performance with existing methods.

Contribution

The paper presents BCI, a new Bayesian hierarchical model for bivariate causal inference that explicitly models data distributions and compares favorably to state-of-the-art methods.

Findings

BCI performs reliably on synthetic and real-world data.

BCI achieves accuracy comparable to existing methods.

BCI is robust in high noise, discretized, and sparse data scenarios.

Abstract

We address the problem of two-variable causal inference without intervention. This task is to infer an existing causal relation between two random variables, i.e. $X \rightarrow Y$ or $Y \rightarrow X$ , from purely observational data. As the option to modify a potential cause is not given in many situations only structural properties of the data can be used to solve this ill-posed problem. We briefly review a number of state-of-the-art methods for this, including very recent ones. A novel inference method is introduced, Bayesian Causal Inference (BCI), which assumes a generative Bayesian hierarchical model to pursue the strategy of Bayesian model selection. In the adopted model the distribution of the cause variable is given by a Poisson lognormal distribution, which allows to explicitly regard the discrete nature of datasets, correlations in the parameter spaces, as well as the variance of probability densities on logarithmic scales. We assume Fourier diagonal Field covariance operators. The model itself is restricted to use cases where a direct causal relation $X \rightarrow Y$ has to be decided against a relation $Y \rightarrow X$ , therefore we compare it other methods for this exact problem setting. The generative model assumed provides synthetic causal data for benchmarking our model in comparison to existing State-of-the-art models, namely LiNGAM , ANM-HSIC , ANM-MML , IGCI and CGNN . We explore how well the above methods perform in case of high noise settings, strongly discretized data and very sparse data. BCI performs generally reliable with synthetic data as well as with the real world TCEP benchmark set, with an accuracy comparable to state-of-the-art algorithms. We discuss directions for the future development of BCI .

Figures

Tables

Table 1. Table 1 : Accuracy for the synthetic data benchmark. All parameters for the forward model besides the mentioned one are kept to default values, namely n bins = 512 subscript 𝑛 bins 512 n_{\mathrm{bins}}=512 , N ~ = 300 ~ 𝑁 300 \tilde{N}=300 , ς 2 = 0.05 superscript 𝜍 2 0.05 \varsigma^{2}=0.05 .

Model	default	${\varsigma }^2$=0.2	${\varsigma }^2$=1	$n_{\mathrm{bins}}$=16	$n_{\mathrm{bins}}$=8	30 samples	10 samples
BCI	0.98	0.94	0.90	0.93	0.97	0.92	0.75
LiNGAM	0.30	0.31	0.40	0.23	0.21	0.44	0.45
ANM-HSIC	1.00	0.98	0.94	0.99	1.00	0.91	0.71
ANM-MML	1.00	0.99	0.99	1.00	1.00	0.98	0.69
IGCI	0.65	0.60	0.58	0.24	0.09	0.48	0.40
CGNN	0.72	0.75	0.77	0.57	0.22	0.46	0.39

Table 2. Table 2 : Accuracy for TCEP Benchmark

Model	TCEP	TCEP with 75 samples
BCI	0.64	0.60
ANM-HSIC	0.63	0.54
ANM-MML	0.58	0.56
IGCI	0.66	0.62
CGNN	0.70	0.69

Equations

$$X\to Y\iff \mathcal{P}(y|\mathrm{do}(x))\ne \mathcal{P}(y|\mathrm{do}(x^{\prime }))$$

$${\mathcal{O}}_{X\to Y}=\frac{\mathcal{P}(\mathbf{d}|X\to Y,M)}{\mathcal{P}(\mathbf{d}|Y\to X,M)},$$

$$\mathcal{P}(x|\beta )\propto e^{\beta (x)}$$

$${\mathbb{E}}_{\beta \sim \mathcal{P}(\beta )}[\beta (x)]=\mathbb{E}[\beta (x+t)]$$

$${\mathbb{E}}_{\beta \sim \mathcal{P}(\beta )}[\beta (x)\beta (y)]=\mathbb{E}[\beta (x+t)\beta (y+t)]$$

$$(\mathcal{F}B{\mathcal{F}}^{-1})(k,q)=2\pi \delta (k-q)P_{\beta }(k)$$

$$R_{jx}\equiv \delta (x-z_j)$$

$$\mathcal{P}(k_j|{\lambda }_j)=\frac{{\lambda }_j^{k_j}{e}^{-{\lambda }_j}}{k_j!}$$

$${\lambda }_j={\mathbb{E}}_{(k|\beta )}[k_j]=\rho e^{{\beta }_{z_j}}=\int dx{R}_{jx}{e}^{{\beta }_x}=\rho (R{e}^{\beta }{)}_j$$

$$\mathbf{\lambda }=\rho \mathbf{R}{e}^{\beta }=\rho e^{\beta (\mathbf{z})}$$

$$\mathcal{P}(\mathbf{k}|\mathbf{\lambda })=\prod _j\frac{{\lambda }_j^{k_j}{e}^{-{\lambda }_j}}{k_j!}=\prod _j\frac{(R_j{e}^{\beta }{)}^{k_j}{e}^{-R_j{e}^{\beta }}}{k_j!}=\frac{({\prod }_j(R_j{e}^{\beta }{)}^{k_j})e^{-1^{†}\mathbf{R}{e}^{\beta }}}{{\prod }_j{k}_j!}$$

$$\mathcal{P}(\mathbf{x}|\mathbf{k})=\frac{1}{N!}$$

$$\mathcal{P}(\mathbf{k}|\beta )=\frac{e^{{\sum }_j{k}_j\beta (z_j)}{e}^{-{\rho }^{†}{e}^{\beta (\mathbf{z})}}}{{\prod }_j{k}_j!}=\frac{e^{{\mathbf{k}}^{†}\beta (\mathbf{z})-{\mathbf{\rho }}^{†}{e}^{\beta (\mathbf{z})}}}{{\prod }_j{k}_j!}$$

$$\mathcal{P}(\mathbf{x}|P_{\beta },X\to Y)\approx \frac{1}{N!}\frac{e^{+{\mathbf{k}}^{†}{\mathbf{\beta }}_0-{\mathbf{\rho }}^{†}{e}^{{\mathbf{\beta }}_0}-\frac{1}{2}{\beta }_0^{†}{B}^{-1}{\beta }_0}}{{\left|\rho B\hat{e^{{\mathbf{\beta }}_0}}+\text{mathds}1\right|}^{\frac{1}{2}}{\prod }_j{k}_j!}$$

$$\mathcal{H}(\mathbf{x}|P_{\beta },X\to Y)\approx {\mathcal{H}}_0+\frac{1}{2}\log |\rho B\hat{e^{{\mathbf{\beta }}_0}}+\widehat{\text{mathds}1}|+\log (\prod _j{k}_j!)-{\mathbf{k}}^{†}{\mathbf{\beta }}_0+{\mathbf{\rho }}^{†}{e}^{{\mathbf{\beta }}_0}+\frac{1}{2}{\beta }_0^{†}{B}^{-1}{\beta }_0$$

$${\mathbb{R}}^{[0,1]}\ni f\sim \mathcal{N}(0|f,F)$$

$$(\mathcal{F}F{\mathcal{F}}^{-1})(k,q)=2\pi \delta (k-q)P_f(k),$$

$$\mathbf{y}$$

$$\mathbf{ϵ}$$

$$\mathcal{E}$$

$$\mathcal{P}(\mathbf{y}|\mathbf{x},P_f,\varsigma ,X\to Y)$$

$$=(2\pi {)}^{-\frac{N}{2}}{\left|\tilde{F}+\mathcal{E}\right|}^{-\frac{1}{2}}{e}^{-\frac{1}{2}{\mathbf{y}}^{†}(\tilde{F}+\mathcal{E}{)}^{-1}\mathbf{y}}$$

$$O_{X\to Y}$$

$$=\exp [\mathcal{H}(\mathbf{d}|P_{\beta },P_f,\varsigma ,Y\to X)-\mathcal{H}(\mathbf{d}|P_{\beta },P_f,\varsigma ,X\to Y)]$$

$$\mathcal{H}(\mathbf{d}|P_{\beta },P_f,\varsigma ,X\to Y)$$

$$={\mathcal{H}}_0+\log (\prod _j{k}_j!)+\frac{1}{2}\log |\rho B\hat{e^{{\mathbf{\beta }}_0}}+\text{mathds}1|-{\mathbf{k}}^{†}{\mathbf{\beta }}_0+$$

$$+{\mathbf{\rho }}^{†}{e}^{{\mathbf{\beta }}_0}+\frac{1}{2}{\beta }_0^{†}{B}^{-1}{\beta }_0+\frac{1}{2}{\mathbf{y}}^{†}(\tilde{F}+\mathcal{E}{)}^{-1}\mathbf{y}+\frac{1}{2}\left|\tilde{F}+\mathcal{E}\right|$$

$$P_{\beta }=P_f\propto \frac{1}{q^4+1},$$

$${\varsigma }^2=0.01,$$

$$r=512,$$

$$\rho =1.$$

$$\mathcal{P}(\mathbf{x}|P_{\beta },X\to Y)$$

$$=\frac{1}{N!}|2\pi B{|}^{-\frac{1}{2}}\int \mathcal{D}[\beta ]\frac{e^{{\mathbf{k}}^{†}\beta (\mathbf{z})-{\mathbf{\rho }}^{†}{e}^{\beta (\mathbf{z})}}}{{\prod }_j{k}_j!}{e}^{-\frac{1}{2}{\beta }^{†}{B}^{-1}\beta }=$$

$$=\frac{|2\pi B{|}^{-\frac{1}{2}}}{N!{\prod }_j{k}_j!}\int \mathcal{D}[\beta ]e^{-\gamma [\beta ]}$$

$$\gamma [\beta ]\equiv -{\mathbf{k}}^{†}\beta (\mathbf{z})+{\mathbf{\rho }}^{†}{e}^{\beta (\mathbf{z})}+\frac{1}{2}{\beta }^{†}{B}^{-1}\beta .$$

$${\partial }_{\beta }\gamma [\beta ]=-{\mathbf{k}}^{†}+\rho (e^{\beta (\mathbf{z})}{)}^{†}+{\beta }^{†}{B}^{-1}$$

$${\partial }_{\beta }{\partial }_{\beta }\gamma [\beta ]=\widehat{\rho e^{\beta (\mathbf{z})}}+B^{-1}.$$

$${\mathbf{k}}_u^{†}\equiv \sum _{j=1}^{n_{\mathrm{bins}}}{k}_j({\tilde{R}}_j{)}_u$$

$$(\widehat{\rho e^{\beta (\mathbf{z})}}{)}_{uv}=\rho \sum _{j=1}^{n_{\mathrm{bins}}}({\tilde{R}}_j{)}_u({\tilde{R}}_j{)}_v{e}^{\beta (u)}.$$

Full text

, which assumes a generative Bayesian hierarchical model to pursue the strategy of Bayesian model selection. In the adopted model the distribution of the cause variable is given by a Poisson lognormal distribution, which allows to explicitly regard the discrete nature of datasets, correlations in the parameter spaces, as well as the variance of probability densities on logarithmic scales. We assume Fourier diagonal Field covariance operators. The model itself is restricted to use cases where a direct causal relation $X\rightarrow Y$ has to be decided against a relation $Y\rightarrow X$, therefore we compare it other methods for this exact problem setting. The generative model assumed provides synthetic causal data for benchmarking our model in comparison to existing State-of-the-art models, namely LiNGAM, ANM-HSIC, ANM-MML, IGCI and CGNN. We explore how well the above methods perform in case of high noise settings, strongly discretized data and very sparse data. BCI performs generally reliable with synthetic data as well as with the real world TCEP benchmark set, with an accuracy comparable to state-of-the-art algorithms. We discuss directions for the future development of BCI.

\pubvolume

xx \issuenum1 \articlenumber5

\historyReceived: date; Accepted: date; Published: date

\TitleA Bayesian Model for Bivariate Causal Inference \AuthorMaximilian Kurthen1,∗\orcidB and Torsten Enßlin1\orcidA \AuthorNamesMaximilian Kurthen and Torsten Enßlin

\corresCorrespondence: MKurthen@live.de

1 Introduction

1.1 Motivation and Significance of the Topic

Causal Inference regards the problem of drawing conclusions about how some entity we can observe does - or does not - influence or is being influenced by another entity. Having knowledge about such law-like causal relations enables us to predict what will happen ( $\widehat{=}$ the effect) if we know how the circumstances ( $\widehat{=}$ the cause) do change. For example, one can draw the conclusion that a street will be wet (the effect) whenever it rains (the cause). Knowing that it will rain, or indeed observing the rainfall itself, enables one to predict that the street will be wet. Less trivial examples can be found in the fields of epidemiology (identifying some bacteria as the cause of a disease) or economics (knowing how taxes will influence the GDP of a country). Under ideal conditions the system under investigation can be manipulated. Such interventions might allow to set causal variables to specific values which allows to study their effects statistically. In many applications, like in astronomy, geology, global economics and others, this is hardly possible. E.g. in the area of astrophysics, observed properties of galaxies or galaxy clusters include redshift, size, spectral features, fluxes at different wavelength for line of sights through the Milky Way. There an identification of causal directions has to rest on the hints the causal relation imprints onto the available data only. Restricting on the decision between a direct causal relation $X\rightarrow Y$ vs. $Y\rightarrow X$ ignores important possibilities (such as hidden confounders), however is still important as decision within of a subset of possible causal relations and has been giving rise to own benchmark datasets Mooij et al. (2016).

Especially within the fields of data science and machine learning specific tasks from causal inference have been attracting much interest recently. Hernán et al. (2018) propose that causal inference stands as a third main task of data science besides description and prediction. Pearl (2018) claims that the task of causal inference will be the next “big problem” for Machine Learning. Such a specific problem is the two variable causal inference, also addressed as the cause-effect problem by Peters et al. (2017). Given purely observational data from two random variables, $X$ and $Y$, which are directly causally related, the challenge is to infer the correct causal direction. Interestingly, this is an incorporation of a fundamental asymmetry between cause and effect which does always hold and can be exploited to tackle such an inference problem. Given two random variables, $X$ and $Y$ which are related causally, $X\rightarrow Y$ (“$X$ causes $Y$”), there exists a fundamental independence between the distribution of the cause $\mathcal{P}(X)$ and the mechanism which relates the cause $X$ to the effect $Y$. This independence however does not hold in the reverse direction. Most of the proposed methods for the inference of such a causal direction make use of this asymmetry in some way, either by considering the independence directly (Daniusis et al., 2010; Mooij et al., 2016), or by taking into account the algorithmic complexity for the description of the factorization $\mathcal{P}(X)\mathcal{P}(Y|X)$ and comparing it to the complexity of the reverse factorization $\mathcal{P}(Y)\mathcal{P}(X|Y)$.

The point we want to make here is, that from the perspective of Bayesian probability theory, causal inference looks like a classical hypothesis testing problem. For this, the probability $\mathcal{P}(\bm{d}|X\rightarrow Y,\mathcal{M})$ is to be compared to $\mathcal{P}(\bm{d}|Y\rightarrow X,\mathcal{M})$. Here, $\bm{d}$ is the data, $X\rightarrow Y$ is the causal direction and $\mathcal{M}$ is a (meta-) model of how causal relations between $X$ and $Y$ are typically realized. In the following we will adopt a specif choice for the model $\mathcal{M}$. This choice and those of our numerical approximations could and should be criticized and eventually improved. The hypothesis we are following here, however, is that, given $\mathcal{M}$, Bayesian inference is everything that is needed to judge causal directions. Different algorithms for causal inference might just result from different assumptions and different numerical approximations made.

1.2 Structure of the Work

The rest of the paper will be structured as following. In Section 2 we will briefly outline and specify the problem setting. We also will review existing methods here, namely Additive Noise Models, Information Geometric Causal Inference and Learning Methods. Section 3 will describe our inference model which is based on a hierarchical Bayesian model. In Section 4 we will accompany the theoretical framework with experimental results. To that end we outline the “forward model” which allows to sample causally related data in 4.1. We describe a specific algorithm for the inference model in 4.2, which is then tested on various benchmark data (4.3). The performance is evaluated and compared to state-of-the-art methods mentioned in Section 2. We conclude in Section 5 by assessing that our model generally can show competitive classification accuracy and propose possibilities to further advance the model.

2 Problem Setting and Related Work

Here and in the following we assume two random variables, $X$ and $Y$, which map onto measurable spaces $\mathcal{X}$ and $\mathcal{Y}$. Our problem, the two-variable causal inference, is to determine if $X$ causes $Y$ or $Y$ causes $X$, given only observations from these random variables. The possibilities that they are unrelated ($X\perp\!\!\!\perp Y$ ), connected by a confounder Z ( $X\leftarrow Z\rightarrow Y$ ) or interrelated ($X\leftrightarrow Y$ ) are ignored here for the sake of clarity.

2.1 Problem Setting

Regarding the definition of causality we refer to the do-calculus introduced by Pearl (2000). Informally, the intervention $\mathrm{do}(X=x)$ can be described as setting the random variable $X$ to attain the value $x$. This allows for a very compact definition of a causal relation $X\rightarrow Y$ (“$X$ causes $Y$”) via

[TABLE]

for some $x,x^{\prime}$ being realizations of $X$ and $y$ being a realization of $Y$ (Mooij et al., 2016). 222 In general the probabilities P(y—X=x) and P(y—do(X=x)) are different. While the conditional probability P(y—X=x) corresponds to just observing the value x for the variable X, the do-probability P(y—do(X=x)) corresponds to a direct manipulation of the system, only modifying x, without changing any other variable directly.

We want to focus on the case of two observed variables, where either $X\rightarrow Y$ or $Y\rightarrow X$ holds. Our focus is on the specific problem to decide, in a case where two variables $X$ and $Y$ are observed, whether $X\rightarrow Y$ holds or $Y\rightarrow X$ We suppose to have access to a finite number of samples from the two variables, i.e. samples $\bm{x}=(x_{1},...,x_{N})$ from $X$ and $\bm{y}=(y_{1},...,y_{N})$ from $Y$. Our task is to decide the true causal direction using only these samples:

Problem \thetheorem.

**Prediction of causal direction for two variables

Input: A finite number of sample data $\bm{d}\equiv(\bm{x},\bm{y})$, where $\bm{x}=(x_{1},...,x_{N}),\bm{y}=(y_{1},...,y_{N})$

Output: A predicted causal direction $\mathcal{D}\in\{X\rightarrow Y,Y\rightarrow X\}$**

2.2 Related Work

Approaches to causal inference from purely observational data are often divided into three groups (Spirtes and Zhang, 2016; Mitrovic et al., 2018), namely constraint-based, score-based and asymmetry-based methods. Sometimes this categorization is extended by considering learning methods as a fourth, separate group. Constraint-based and score-based methods are using conditioning on external variables. In a two-variable case there are no external variables so they are of little interest here.

Asymmetry-based methods exploit an inherent asymmetry between cause and effect. This asymmetry can be framed in different terms. We will follow a similar overview here as (Peters et al., 2017), where we also refer to for a more detailed discussion. One way is to use the concept of algorithmic complexity - given a true direction $X\rightarrow Y$, the factorization of the joint probability into $\mathcal{P}(X,Y)=\mathcal{P}(X)\mathcal{P}(Y|X)$ will be less complex than the reverse factorization $\mathcal{P}(Y)\mathcal{P}(X|Y)$. Such an approach is often used by Additive Noise Models (ANMs). This family of inference models assume additive noise, i.e. in the case $X\rightarrow Y$, $Y$ is determined by some function $f$, mapping $X$ to $Y$, and some collective noise variable $E_{Y}$, i.e. $Y=f(X)+E_{Y}$, where $X$ is independent of $E_{Y}$.

An early model called LiNGAM (Shimizu et al., 2006) uses Independent Component Analysis on the data belonging to the variables. This model however makes the assumptions of linear relations and non-Gaussian noise.

A more common approach is to use some kind of regression (e.g. Gaussian process regression) to get an estimate on the function $f$ and measure how well the model such obtained fits the data. The latter is done by measuring independence between the cause variable and the regression residuum (ANM-HSIC, Hoyer et al., 2009; Mooij et al., 2016) or by employing a Bayesian model selection (introduced besides other methods as (ANM-MML, Stegle et al., 2010)

Another way of framing the asymmetry mentioned above is to state that the mechanism relating cause and effect should be independent of the cause. This formulation is employed by the concept of IGCI (Information Geometric Causal Inference, Daniusis et al., 2010).

The recent advances in the field of deep learning are represented in an approach called CGNN (Causal Generative Neural Networks, Goudet et al., 2017). The authors use Generative Neural Networks to model the distribution of one variable given samples from the other variable. As Neural Networks are able to approximate nearly arbitrary functions, the direction where such a artificial modelling is closer to the real distributions (inferred from the samples) is preferred.

Finally, KCDC (Kernel Conditional Deviance for Causal Inference, Mitrovic et al., 2018) uses the thought of asymmetry in the algorithmic complexity directly on the conditional distributions $\mathcal{P}(X|Y=y),\mathcal{P}(Y|X=x)$. The model measures the variance of the conditional mean embeddings of the above distributions and prefers the direction with the less varying embedding.

A somewhat related approach employed by Friedman and Koller (2003) uses a Bayesian model in combination with MCMC sampling in order to reconstruct Bayesian networks. Heckerman et al. (2006) use a Bayesian approach specifically in the domain of causal model discovery and compare it to constraint-based approaches. However, an approach which considers multi-variable causal models is fundamentally different from a 2-variable scenario (as conditioning on a third variable is not possible in the latter case).

3 A Bayesian Inference Model

Our contribution incorporates the concept of Bayesian model selection and builds on the formalism of Information Field Theory (IFT, the information theory for fields Enßlin et al., 2009). Bayesian model selection compares two competing models, in our case $X\rightarrow Y$ and $Y\rightarrow X$, and asks for the ratio of the marginalized likelihoods,

[TABLE]

the so called Bayes Factor. Here $M$ denotes the hyperparameters which are assumed to be the same for both models and are yet to be specified.

In the setting of the present causal inference problem a similar approach has already been used by Stegle et al. (2010). This approach however does use a Gaussian mixture model for the distribution of the cause variable while we model the logarithmic cause distribution as a more flexible Gaussian random field (or Gaussian process) Rasmussen and Williams (2006) and explicitly consider an additional discretization via introducing counts in bins (using a Poissonian statistic on top ) . Gaussian random fields have, in principle, an infinite number of degrees of freedom, making them an interesting choice to model our distribution of the cause variable and the function relating cause and effect. The formalism of IFT borrows computational methods from quantum field theory for computations with such random fields.

Throughout the following we will consider $X\rightarrow Y$ as the true underlying direction which we derive our formalism for. The derivation for $Y\rightarrow X$ will follow analogously by switching the variables. We will begin with deriving in 3.1 the distribution of the cause variable, $\mathcal{P}(X|X\rightarrow Y,M)$. In 3.2 we continue by considering the conditional distribution $\mathcal{P}(Y|X,X\rightarrow Y,M)$. Combining those results, we compute then the full Bayes factor in section 3.3.

3.1 Distribution of the Cause Variable

Without imposing any constraints, we reduce our problem to the interval $[0,1]$ by assuming that $\mathcal{X}=\mathcal{Y}=[0,1]$. This can always be ensured by rescaling the data. We make the assumption that in principle the cause variable $X$ follows a lognormal distribution.

[TABLE]

with $\beta\in\mathbb{R}^{[0,1]}$ 333 Throughout the paper we will use the set theory notation for functions, i.e. for a function $f:X\rightarrow Y$ we write $f\in Y^{X}$, which allows a concise statement of domain and codomain without imposing any further restrictions. , being some signal field which follows a zero-centered normal distribution, $\beta\sim\mathcal{N}(\beta|0,B)$.

Here we write $B$ for the covariance operator $\mathbb{E}_{\beta\sim\mathcal{P}(\beta)}[\beta(x_{0})\beta(x_{1})]=B(x_{0},x_{1})$.

We postulate statistical homogeneity for the covariance, that is

[TABLE]

i.e. first and second moments should be independent on the absolute location. The Wiener-Khintchine Theorem now states that the covariance has a spectral decomposition, i.e. it is diagonal in Fourier space, under this condition (see e.g. Chatfield, 2016). Denoting the Fourier transform by $\mathcal{F}$, i.e. in the one dimensional case, $\mathcal{F}[f](q)=\int\mathop{}\!\mathrm{d}x\,e^{-iqx}f(x)$. Therefore, the covariance can be completely specified by a one dimensional function:

[TABLE]

Here, $P_{\beta}(k)$ is called the power spectrum.

Building on these considerations we now regard the problem of discretization. Measurement data itself is usually not purely continuous but can only be given in a somewhat discretized way (e.g. by the measurement device itself or by precision restrictions imposed from storing the data). Another problem is that many numerical approaches to inference tasks, such as Gaussian Process regression, use finite bases as approximations in order to efficiently obtain results (Mooij et al., 2016; Stegle et al., 2010). Here, we aim to directly confront these problems by imposing a formalism where the discretization is inherent.

So instead of taking a direct approach with the above formulation, we use a Poissonian approach and consider an equidistant grid $\{z_{1},...,z_{n_{\mathrm{bins}}}\}$ in the $[0,1]$ interval. This is equivalent to defining bins, where the $z_{j}$ are the midpoints of the bins. We now take the measurement counts, $k_{i}$ which gives the number of $x$-measurements within the $i$-th bin. For these measurement counts we now take a Poisson lognormal distribution as an Ansatz, that is, we assume that the measurement counts for the bins are Poisson distributed, where the means follow a lognormal distribution. We argue that this is indeed a justified approach here, as in a discretized scenario we have to deal with count-like integer data (where a Poisson distribution is natural). The lognormal distribution of the Poisson parameters is in our eyes well-justified here. On the one hand, it takes into account the non-negativity of the Poisson parameter. On the other hand, only proposing a normal distribution of the log allows for a uncertainty in the order of magnitude while permitting for spatial correlation in this log density.

We can model this discretization by applying a response operator $R:\mathbb{R}^{[0,1]}\rightarrow\mathbb{R}^{n_{\mathrm{bins}}}$ to the lognormal field. This is done in the most direct way via employing a Dirac delta distribution

[TABLE]

In order to allow for a more compact notation we will use an index notation from now on, e.g. $f_{x}=f(x)$ for some function $f$ or $O_{xy}=O(x,y)$ for some operator $O$. Whenever the indices are suppressed, an integration (in the continuous case) or dot product (in the discrete case) is understood, e.g. $(Of)_{x}\equiv O_{xy}f_{y}=\int\mathop{}\!\mathrm{d}yO_{xy}f_{y}=\int\mathop{}\!\mathrm{d}yO(x,y)f(y)$ In the following we will use bold characters for finite dimensional vectors, e.g. $\bm{\lambda}\equiv(\lambda_{1},...,\lambda_{n_{\mathrm{bins}}})^{T}$. By inserting such a finite dimensional vector in the argument of a function, e.g. $\beta(\bm{x})$ we refer to a vector consisting of the function evaluated at each entry of $\bm{x}$, that is $\beta(\bm{z})\equiv(\beta(z_{1}),...,\beta(z_{n_{\mathrm{bins}}})$. Later on we will use the notation $\widehat{\cdot}$ which raises some vector to a diagonal matrix ($\widehat{\bm{x}}_{ij}\equiv\delta_{ij}x_{i}$ (no summation implicated)). We will use this notation analogously for fields, e.g. ($\widehat{\beta}_{uv}\equiv\delta(u-v)\beta(u)$). Writing $\bm{1}^{\dagger}\bm{R}$ denotes the dot product of the vector $\bm{R}$ with a vector of ones and hence corresponds to the summation of the entries of $\bm{R}$ ($\bm{1}^{\dagger}\bm{R}=\sum_{j}R_{j}$). Now we can state the probability distribution for $k_{j}$, the measurement count in bin $j$:

[TABLE]

Therefore, considering the whole vector $\bm{\lambda}$ of bin means and the vector $\bm{k}$ of bin counts at once:

[TABLE]

The last equation follows from the consideration that given the counts $(k_{1},...,k_{n_{\mathrm{bins}}})$ for the bins, only the positions of the observations $(x_{1},...,x_{N})$ is fixed, but the ordering is not. The $N$ observations can be ordered in $N!$ ways.

Now considering the whole vector of bin counts $\bm{k}$ at once, we get

[TABLE]

A marginalization in $\beta$ involving a Laplace approximation around the most probable $\beta=\beta_{0}$ leads to (see Appendix for a detailed derivation):

[TABLE]

where $\mathcal{H}(\cdot)\equiv-\log(\mathcal{P}(\cdot))$ is called the information Hamiltonian in IFT and $\mathcal{H}_{0}$ collects all terms which do not depend on the data $\bm{d}$.

3.2 Functional Relation of Cause and Effect

Similarly to $\beta$, we suppose a Gaussian distribution for the function $f$, relating $Y$ to $X$:

[TABLE]

Proposing a Fourier diagonal covariance $F$ once more, determined by a power spectrum $P_{f}$,

[TABLE]

we assume additive Gaussian noise, using the notation $f(\bm{x})\equiv(f(x_{1}),...,f(x_{N}))^{T}$ and $\bm{\epsilon}\equiv(\epsilon_{1},...,\epsilon_{N})^{T}$. So, in essence we are proposing a Gaussian Process Regression with a Fourier diagonal covariance. We have

[TABLE]

that is, each independent noise sample is drawn from a zero-mean Gaussian distribution with given variance $\varsigma^{2}$.

Knowing the noise $\bm{\epsilon}$, the cause $\bm{x}$ and the causal mechanism $f$ completely determines $\bm{y}$ via Equation 6. Therefore, $\mathcal{P}(\bm{y}|\bm{x},f,\bm{\epsilon},X\rightarrow Y)=\delta(\bm{y}-f(\bm{x})-\bm{\epsilon})$. We can now state the conditional distribution for the effect variable measurements $\bm{y}$, given the cause variable measurements $\bm{x}$. Marginalizing out the dependence on the relating function $f$ and the noise $\bm{\epsilon}$ we get:

[TABLE]

In the equation above, $\tilde{F}$ denotes a the $N\times N$-matrix with entries $\tilde{F}_{ij}=F(x_{i},x_{j})$ 444 This type of matrix, i.e. the evaluation of covariance or kernel at certain positions, is sometimes called a Gram matrix. . Again, we give a detailed computation in the appendix.

3.3 Computing the Bayes factor

Now we are able to calculate the full likelihood of the data $\bm{d}=(\bm{x},\bm{y})$ given our assumptions $P_{\beta},P_{f},\varsigma$ for the direction $X\rightarrow Y$ and vice versa $Y\rightarrow X$. As we are only interested in the ratio of the probabilities and not in the absolute probabilities itself, it suffices to calculate the Bayes factor:

[TABLE]

Here we used again the information Hamiltonian $\mathcal{H}(\cdot)\equiv-\log\mathcal{P}(\cdot)$

Making use of Equations 4 and 7 we get, using the calculus for conditional distributions on the Hamiltonians, $\mathcal{H}(A,B)=\mathcal{H}(A|B)+\mathcal{H}(B)$,

[TABLE]

In this formula we suppressed the dependence of $\tilde{F},\beta_{0}$ on $\bm{x}$ (for the latter, the dependence is not explicit, but rather implicit as $\beta_{0}$ is determined by the minimum of the $\bm{x}$-dependent functional $\gamma$). We omit stating $\mathcal{H}(\bm{d}|P_{\beta},P_{f},\varsigma,Y\rightarrow X)$ explicitly as the expression is just given by taking Equation 8 and switching $\bm{x}$ and $\bm{y}$ or $X$ and $Y$, respectively.

4 Implementation and Benchmarks

We can use our model in a forward direction to generate synthetic data with a certain underlying causal direction. We describe this process in 4.1. In 4.2 we give an outline on the numerical implementation of the inference algorithm. This algorithm is tested on and compared on benchmark data. To that end we use synthetic data and real world data. We describe the specific datasets and give the results in 4.3.

4.1 Sampling Causal Data via a Forward Model

To estimate the performance of our algorithm and compare it with other existing approaches a benchmark dataset is of interest to us. Such benchmark data is usually either real world data or synthetically produced. While we will use the TCEP benchmark set of Mooij et al. (2016) in 4.3.2, we also want to use our outlined formalism to generate artificial data representing causal structures. Based on our derivation for cause and effect we implement a forward model to generate data $\bm{d}$ as following.

Algorithm 1**.**

Sampling of causal data via forward model

**Input: Power spectra $P_{\beta},P_{f}$, noise variance $\varsigma^{2}$, number of bins $n_{\mathrm{bins}}$, desired number222 As we draw the number of samples from Poisson distribution in each bin, we do not deterministically control the total number of samples

of samples $\tilde{N}$

Output: $N$ samples $(d_{i})=(x_{i},y_{i})$ generated from a causal relation of either $X\rightarrow Y$ or $Y\rightarrow X$**

Draw a sample field $\beta\in\mathbb{R}^{[0,1]}$ from the distribution $\mathcal{N}(\beta|0,B)$ 2. 2.

Set an equally spaced grid with $n_{\mathrm{bins}}$ points in the interval $[0,1]$: $\bm{z}=(z_{1},...,z_{n_{\mathrm{bins}}}),z_{i}=\frac{i-0.5}{n_{\mathrm{bins}}}$ 3. 3.

*Calculate the vector of Poisson means $\bm{\lambda}=(\lambda_{1},...\lambda_{n_{\mathrm{bins}}})$ with $\lambda_{i}\propto e^{\beta(z_{i})}$ * 4. 4.

At each grid point $i\in\{1,...,n_{\mathrm{bins}}\}$, draw a sample $k_{i}$ from a Poisson distribution with mean $\lambda_{i}$: $k_{i}\sim\mathcal{P}_{\lambda_{i}}(k_{i})$ 5. 5.

Set $N=\sum_{i=1}^{n_{\mathrm{bins}}}k_{i}$, 6. 6.

For each $i\in\{1,...,n_{\mathrm{bins}}\}$ add $k_{i}$ times the element $z_{i}$ to the set of measured $x_{j}$. Construct the vector $\bm{x}=(...,\underbrace{z_{i},z_{i},z_{i}}_{k_{i}\text{ times}},...)$ 7. 7.

Draw a sample field $f\in\mathbb{R}^{[0,1]}$ from the distribution $\mathcal{N}(f|0,F)$. Rescale $f$ s.th. $f\in[0,1]^{[0,1]}$. 8. 8.

Draw a multivariate noise sample $\bm{\epsilon}\in\mathbb{R}^{N}$ from a normal distribution with zero mean and variance $\varsigma^{2}$, $\bm{\epsilon}\sim\mathcal{N}(\bm{\epsilon}|0,\varsigma^{2})$ 9. 9.

Generate the effect data $\bm{y}$ by applying $f$ to $\bm{x}$ and adding $\bm{\epsilon}$: $\bm{y}=f(\bm{x})+\bm{\epsilon}$ 10. 10.

With probability $\frac{1}{2}$ return $\bm{d}=(\bm{x}^{T},\bm{y}^{T})$, otherwise return $\bm{d}=(\bm{y}^{T},\bm{x}^{T})$,

Comparing the samples for different power spectra (see Fig. 2), we decide to sample data with power spectra $P(q)=\frac{1000}{q^{4}+1}$ and $P(q)=\frac{1000}{q^{6}+1}$ , as these seem to resemble “natural” mechanisms, see Fig. 2.

4.2 Implementation of the Bayesian Causal Inference Model

Based on our derivation in Section 3 we propose a specific algorithm to decide the causal direction of a given dataset and therefore give detailed answer for Problem 2.1. Basically the task comes down to find the minimum $\beta_{0}$ for the saddle point approximation and calculate the terms given in Equation 8:

Algorithm 2**.**

2-variable causal inference

**Input: Finite sample data $\bm{d}\equiv(\bm{x},\bm{y})\in\mathbb{R}^{N\times 2}$, Hyperparameters $P_{\beta},P_{f},\varsigma^{2},r$

Output: Predicted causal direction $\mathcal{D}_{X\rightarrow Y}\in\{X\rightarrow Y,Y\rightarrow X\}$

**

Rescale the data to the $[0,1]$ interval. I.e. $\min\{x_{1},...,x_{N}\}=\min\{y_{1},...,y_{N}\}=0$ and $\max\{x_{1},...,x_{N}\}=\max\{y_{1},...,y_{N}\}=1$ 2. 2.

Define an equally spaced grid of $(z_{1},...,z_{n_{\mathrm{bins}}})$ in the interval $[0,1]$ 3. 3.

Calculate matrices $\bm{B},\bm{F}$ representing the covariance operators $B$ and $F$ evaluated at the positions of the grid, i.e. $\bm{B}_{ij}=B(z_{i},z_{j})$ 4. 4.

Find the $\beta_{0}\in\mathbb{R}^{[0,1]}$ for which $\gamma$, as defined in Appendix A.1 ( Equation 14), becomes minimal 5. 5.

Calculate the $\bm{d}$-dependent terms of the information Hamiltonian in Equation 8 (i.e. all terms except $\mathcal{H}_{0}$) 6. 6.

Repeat steps 4 and 5 with $\bm{y}$ and $\bm{x}$ switched 7. 7.

Calculate the Bayes factor $\mathcal{O}_{X\rightarrow Y}$ 8. 8.

If $\mathcal{O}_{X\rightarrow Y}>1$, return $X\rightarrow Y$, else return $Y\rightarrow X$

We provide an implementation of Algorithm 2 in Python333 https://gitlab.mpcdf.mpg.de/ift/bayesian_causal_inference . We approximate the operators $B,F$ as matrices $\in\mathbb{R}^{n_{\mathrm{bins}}\times n_{\mathrm{bins}}}$, which allows us to explicitly numerically compute the determinants and the inverse. As the most critical part we consider the minimization of $\beta$, i.e. step 4 in Algorithm 2. As we are however able to analytically give the curvature $\Gamma_{\beta}$ and the gradient $\partial_{\beta}\gamma$ of the energy $\gamma$ to minimize, we can use a Newton-scheme here. We derive satisfying results (see Fig. LABEL:fig:demonstration_beta_min ) using the Newton-CG algorithm (Nocedal and Wright, 2006), provided by the SciPy-Library (Jones et al., 2001). After testing our algorithm on different benchmark data, we choose the default hyperparameters as

[TABLE]

While fixing the power spectra might seem somewhat arbitrary, we remark that this corresponds to fixing a kernel e.g. as a squared exponential kernel, which is done in many publications (e.g. Mitrovic et al. (2018); Goudet et al. (2017)). Future extensions of our method might learn $P_{\beta}$ and $P_{f}$ if the data is rich enough.

4.3 Benchmark Results

We compare our outlined model, in the following called BCI (Bayesian Causal Inference), to a number of state-of-the-art approaches. The selection of the considered methods is influenced by the ones in recent publications, e.g. Mitrovic et al. (2018); Goudet et al. (2017). Namely, we include the LiNGAM algorithm, acknowledging it as one of the oldest models in this field and a standard reference in many publications. We also use the ANM algorithm (Mooij et al., 2016) with HSIC and Gaussian Process Regression (ANM-HSIC) as well as the ANM-MML approach (Stegle et al., 2010). The latter uses a Bayesian Model Selection, arguably the closest to the algorithm proposed in this publication, at least to our best knowledge. We further include the IGCI algorithm, as it differs fundamentally in its formulation from the ANM algorithms and has shown strong results in recent publications (Mooij et al., 2016; Goudet et al., 2017; Mitrovic et al., 2018). We employ the IGCI algorithm with entropy estimation for scoring and a Gaussian distribution as reference distribution.

Finally, CGNN (Goudet et al., 2017) represents the rather novel influence of deep learning methods. We use the implementation provided by the authors, with itself uses Python with the Tensorflow (Abadi et al., 2015) library. The most critical hyper-parameter here is, as the authors themselves mention, the number of hidden neurons which we set to a value of $n_{h}=30$, as this is the default in the given implementation and delivers generally good results. We use 32 runs each, as recommended by the authors of the algorithm.

A comparison with the KCDC algorithm would be interesting, unfortunately the authors did not provide any computational implementation so far (October 2019). We compare the mentioned algorithms to BCI on basis of synthetic and real world data. For the synthetic data we use our outlined forward model as outlined in 4.1 with varying parameters. For the real world data we use the well-known TCEP dataset (Tuebingen Cause Effect Pairs, Mooij et al., 2016).

4.3.1 Results for synthetic benchmark data

We generate our synthetic data adopting the power spectra $P(q)=\frac{1}{q^{4}+1}$ for both, $P_{\beta}$, $P_{f}$. We further set $n_{\mathrm{bins}}$=512, $\tilde{N}$=300 and $\varsigma^{2}$=0.05 as default settings. We provide the results of the benchmarks in Table 1. While BCI achieves almost perfect results (98%), the assessed ANM algorithms provide a perfect performance here. As a first variation, we explore the influence of high and very high noise on the performance of the inference models. Therefore we set the parameter $\varsigma^{2}$=0.2 for high noise and $\varsigma^{2}$=1 for very high noise in Algorithm 1, while keeping the other parameters set to the default values. While our BCI algorithm is affected but still performs reliably with an accuracy of $\geq 90\%$, especially the ANM algorithms perform remarkably robust in presence of the noise. This is likely due to the fact that the distribution of the true cause $\mathcal{P}(X)$ is not influenced by the high noise and this distribution is assessed in its own by those.

As our model uses a Poissonian approach, which explicitly considers discretization effects of data measurement, it is of interest how the performance behaves when using a strong discretization. We emulate such a situation by employing our forward model with a very low number of bins. Again we keep all parameters to default values and set $n_{\mathrm{bins}}$=16 and $n_{\mathrm{bins}}$=8 for synthetic data with high and very high discretization. Again the ANM models turn out to be robust again discretization. CGNN and IGCI perform significantly worse here. In the case of IGCI this can be explained by the entropy estimation, which simply removes non-unique samples. Our BCI algorithm is able to achieve over $90\%$ accuracy here.

We explore another challenge for inference algorithms by strongly reducing the number of samples. While we sampled about 300 observations with our other forward models so far, here we reduce the number of observed samples to 30 and 10 samples. In this case BCI performs very well compared to the other models, in fact it is able to outperform them in the case of just 10 samples being given. We note that of course BCI does have the advantage that it “knows” the hyperparameters of the underlying forward model. Yet we consider the results as encouraging, the advantage will be removed in the confrontation with real world data.

4.3.2 Results for Real World Benchmark Data

The most widely used benchmark set with real world data is the TCEP (Mooij et al., 2016). We use the 102 2-variable datasets from the collection with weights as proposed by the maintainers. As some of the contained datasets include a high number of samples (up to 11000) we randomly subsample large datasets to 500 samples each in order to keep computation time maintainable. We did not include the LiNGAM algorithm here, as we experienced computational problems with obtaining results here for certain datasets (namely pair0098). Goudet et al. (2017) report the accuracy of LiNGAM on the TCEP dataset to be around $40\%$. BCI performs generally comparable to established approaches as ANM and IGCI. CGNN performs best with an accuracy about $70\%$ here, a bit lower than the one reported by Goudet et al. (2017) of around $80\%$. The reason for this is arguably to be found in the fact that we set all hyperparameters to fixed values, while Goudet et al. (2017) used a leave-one-out-approach to find the best setting for the hyperparameter $n_{h}$.

Motivated by the generally strong performance of our approach in the case of sparse synthetic data, we also explore a situation where real world data is only sparsely available. To that end, we subsample all TCEP datasets down to 75 randomly chosen samples kept for each one. To circumvent the chance of an influence of the subsampling procedure we average the results over 20 different subsamplings. The results are as well given in Table 2. The loss in accuracy of our model is rather small.

5 Discussion

The problem of purely bivariate causal discovery is a rather restricted one as it ignores the possibility of causal independence or hidden confounders. However it is a fundamental one and still remains as a part of the decision within subsets of possible causal structures. It also might be a valuable contribution to real world problems such as in astrophysics, where one might be interested in discovering the main causal direction within the multitude of discovered variables. The Bayesian Causal Inference method introduced builds on the formalism of information field theory. In this regard, we employed the concept of Bayesian model selection and made the assumption of additive noise, i.e. $\bm{x}=f(\bm{y})+\bm{\epsilon}$. In contrast to other methods which do so, such as ANM-MML, we do not model the cause distribution by a Gaussian mixture model but by a Poisson Lognormal statistic.

We could show that our model is able to provide classification accuracy in the present causal inference task that is comparable to the one of other methods in the field. Our method is of course restricted to a bivariate problem, i.e. determining the direction of a causal relation between two one-dimensional variables. One difference from our model to existing ones is arguably to be found in the choice of the covariance operators. While our method uses, at its heart, Gaussian Process Regression, most other publications which do so use squared exponential kernels. We, however, choose a covariance which is governed by a $\frac{1}{q^{4}+1}$ power spectrum. This permits detecting more structure at small scales than methods using a squared exponential kernel.

As a certain weak point of BCI we consider the approximation of the uncomputable path integrals via the Laplace approximation. A thorough investigation of error bounds (e.g. Majerski, 2015) is yet to be carried out. As an alternative, one can think about sampling-based approaches to approximate the integrals. A recent publication (Caldwell et al., 2018) introduced a harmonic-mean based sampling approach to approximate moderate dimensional integrals. Adopting such a technique to our very high dimensional case might be promising to improve BCI. Furthermore, the novel Metric Gaussian Variational Inference method (MGVI Knollmüller and Enßlin, 2019) might allow to go beyond the saddle point approximation method used here. At this point we also want to mention that our numerical implementation will have difficulties to scale well to very large (lots of data points) problems. This is however also an issue for competing models and a subsampling procedure can always be used to decrease the scale of the problem.

Another interesting perspective is provided by deeper hierarchical models. While the outlined method took the power spectra and the noise variance as fixed hyperparameters it would also be possible to infer these as well in an extension of the method. MGVI has already permitted the inference of rather complex hierarchical Bayeisan models Hutschenreuter and Enßlin (2019). Yet, the implementation of our model with fixed noise variance and power spectra was able to deliver competitive results with regard to state-of-the-art methods in the benchmarks. In particular, our method seems to be slightly superior in the low sample regime, probably due to the more appropriate Poisson statistic used. We consider this as an encouraging result for a first work in the context of information field theory-based causal inference.

\appendixtitles

yes

Appendix A. Explicit Derivations

We give explicit derivations for the obtained results

A.1 Saddle Point Approximation for the Derivation of the Cause Likelihood

Marginalizing $\beta$ we get

[TABLE]

where

[TABLE]

We approach this integration by a saddle point approximation. In the following we will denote the functional derivative by $\partial$, i.e. $\partial_{f_{z}}\equiv\frac{\delta}{\delta f(z)}$.

Taking the first and second order functional derivative of $\gamma$ w.r.t. $\beta$ we get

[TABLE]

The above derivatives are still defined in the space of functions $\mathbb{R}^{[0,1]}$, that is

[TABLE]

The latter expression therefore represents a diagonal operator with $e^{\beta(\bm{z})}$ as diagonal entries.

Let $\beta_{0}$ denote the function that minimizes the functional $\gamma$, i.e.

[TABLE]

We expand the functional $\gamma$ up to second order around $\beta_{0}$,

[TABLE]

where we dropped higher order terms of $\beta$, used that the gradient at $\beta=\beta_{0}$ vanishes and evaluated the remaining Gaussian integral.

Plugging the result of Equation 15 into Equation 13 and using

[TABLE]

we get

[TABLE]

A.2 Explicit Derivation of the Effect Likelihood

Here we will derive the explicit expression of the likelihood of the effect, which can be carried out analytically. We start with the expression for the likelihood of the effect data $\bm{y}$, given the cause data $\bm{x}$, the power spectrum $P_{f}$, the noise variance $\varsigma^{2}$ and the causal direction $X\rightarrow Y$. That is, we simply marginalize over the possible functions $f$ and noise $\bm{\epsilon}$.

[TABLE]

Above we just used the equation $\bm{y}=f(\bm{x})+\bm{\epsilon}$ and used the distributions for $f$ and $\bm{\epsilon}$. We will now use the Fourier representation of the delta distribution, specifically $\delta(x)=\int\frac{\mathop{}\!\mathrm{d}q}{2\pi}e^{iqx}$.

[TABLE]

Once more we employ a vector of response operators, mapping $\mathbb{R}^{\mathbb{R}}$ to $\mathbb{R}^{N}$,

[TABLE]

This allows to represent the evaluation $f(\bm{x})=\bm{R}^{\dagger}f$, i.e. as a linear dot-product. Using the well known result for Gaussian integrals with linear terms (see e.g. Greiner et al., 2013),

[TABLE]

We are able to analytically do the path integral over $f$,

[TABLE]

Now, we do the integration over the noise variable, $\bm{\epsilon}$, by using the equivalent of Equation 16 for the vector-valued case:

[TABLE]

In the following we will write $\mathbb{R}^{N\times N}\ni\tilde{F}=\bm{R}^{\dagger}F\bm{R}$, with entries $\tilde{F}_{ij}=F(x_{i},x_{j})$. The integration over the Fourier modes $\bm{q}$, again via the multivariate equivalent of Equation 16, will give the preliminary result:

[TABLE]

\authorcontributions

conceptualization, M.K. and T.E.; methodology,; M.K. and T.E.; software, M.K.; validation, M.K. and T.E.; formal analysis, M.K.; investigation, M.K. and T.E.; resources, T.E.; data curation, M.K..; writing–original draft preparation, M.K.; writing–review and editing, T.E.; visualization, M.K.; supervision, T.E.; project administration, T.E.

\funding

This research received no external funding

\conflictsofinterest

The authors declare no conflict of interest.

\externalbibliography

yes