A framework for causal discovery in non-intervenable systems

Causal discovery can be roughly divided into four different tasks: uncovering the causal network, building structural causal models, studying the influence of interventions, and counterfactual reasoning. This paper deals with the first task: uncovering the causal network with a goal to progress scientific knowledge of a system. The information source is time series of variables or processes of the system of interest.

Since the systems of interest are highly nonlinear, building structural equation models is a difficult task that we will not consider here. We refer to the literature on, e.g., relevant vector machines, Bayesian symbolic regression, and many other developments for interesting progress in that field, including using neural networks as structural equations; see, e.g., Refs. 1–5. We do not consider interventions because the systems of interest do not allow for interventions either because interventions are impossible and unethical or would change the dynamics of the system such that the intervention takes a completely different structure. Examples of application are systems where the internal dynamics is so complex that manipulating external forcings does not reveal much about the internal dynamics. One can think about systems such as the atmosphere or the ocean, biological food webs, complex chemical systems, the brain, etc. Interventions on internal variables of these systems typically push these systems off their attractor, resulting in causal inference on regions of state space that are not of interest. One could argue that a small enough intervention, small in magnitude, spatial and temporal extent, would still be useful in such systems. However, due to the strong feedback, it is hard to infer what the direct influence of the intervention is and what is related to reactions to the feedback. For this reason, we will call such systems non-intervenable.

In fact, excluding interventions means we cannot use Pearl and coworkers’ beautiful do-calculus, in which the “do” operator means a direct intervention. If a given causal network allows interventions, one can then use do-calculus to calculate the effect of the intervention without actually doing the intervention.⁶ While our framework potentially complements do-calculus, showing such is beyond the scope of this paper. Our intended applications do not allow either actual or theoretical interventions; therefore, we will use the term “non-intervenable,” or “observational causal inference,” to denote our path of study, restricting ourselves to use data only from observations, reliable simulations, or a combination. (See Ref. 7, for example, who explore a form of importance sampling to explore one data set to infer causal inference on another.)

Because of these limitations, we define causal relations between a target process $x$ and potential driver processes $y$ through two criteria: (1) the cause $y$ precedes the effect $x$⁠, and (2) a causal relation between processes or variables in a system exists if there is flow of information between them, hence information flow from $y$ to $x$⁠. Our goal is not to predict the future of $x$ from $y$⁠; that would be a next step. The goal is to increase understanding of a system by establishing how information flows through the system, where information is to be interpreted as reduction of uncertainty.

Precise mathematical descriptions of observational causal inference started with the seminal works of Wiener⁸ and Granger⁹ in the 1950s and 1960s. Their basic idea was to build a minimal structural equation model by defining a set of functions from observed variables and determine the regression coefficients of these driver functions, or driver processes, on a target process. A large regression coefficient suggests a large causal influence of that driver process on the target process. If the regression coefficient of a process is small, that process is not considered a cause for the target process. Pruning in this way leads to a minimal model, and this minimal model is then the causal model of the target process. In this framework, one has to define the potential driver processes directly, or nonlinear functions of them, beforehand, and the causal inference is in essence looking for linear cross correlations between linear or nonlinear functions.

Granger causality is based on the idea that a process is a driver of a target if it reduces the unexplained variance in the target process. The driver can be nonlinearly related to the target, but in that case, the functional form of the relation has to be specified. Methods based on information theory avoid the specification of this functional form by considering the reduction in unexplained entropy in the target process. The first example of this kind is transfer entropy,¹⁰ and many extensions are now available. The idea here is to define the causal strength of a process $y$ to a process $x$ as the conditional mutual information $I(xt;y<t|z<t)$ in which $z$ represents all other processes, and $y$ and $z$ are in the past of $xt$⁠. These methods pursue the identification of the causal network but are not useful to build an actual structural equation model because information theoretic measures such as (conditional) mutual information are invariant under single-variable nonlinear monotonic transformations. Hence, these methods cannot distinguish between a model in which a variable $x$ is present or say $exp⁡(x)$⁠.

Several of these methods rely on graphical representations, and algorithms typically start either from an empty graph and add strong relations or from a fully connected graph and prune weak relations until a minimal unidirectional acyclic graphical model is found that represents the causal network. An example of building up from an empty graph is the Peter and Clark (PC) algorithm,¹¹ and the so-called greedy equivalence search¹² is an example of a pruning algorithm. The strength of relations is determined via conditional independence tests (e.g., PC) or via scoring rules (Greedy equivalence Search), and the emphasis is more on establishment or removal of causal links, hence determining the causal structure rather than determining the actual causal strength (defined in whatever way). Recently, Sun et al.¹³ pose the problem as an information theoretic optimization method.

These methodologies have been extended to high dimensions and, in particular, applied to earth system processes by, e.g., Refs. 14–16. Many other formalisms have been proposed, and the excellent reviews of Runge et al.¹⁷ and Glymour et al.⁴ contain much of present-day efforts for systems in which interventions are not possible. As mentioned, our interest is in those kind of systems.

More general methods to generate causal models for non-intervenable systems have been developed since. For instance, convergent cross mapping¹⁸ tries to find the underlying dynamical system using Takens’ embeddings. The idea is that if a driver variable can be predicted from the time embedding of the target process, then that driver process is a cause of the target process. The reasoning is opposite to Granger causality, which tries to predict the target from the drivers. Sugihara et al.¹⁸ provide a very careful discussion of the connection of convergent cross mapping and Granger causality. Recently in Leng et al.,¹⁹ an extension called partial cross mapping was proposed that can distinguish direct (i.e., $y→x$⁠) from indirect (i.e., $y→z→x$⁠). Unfortunately, these CCM-based methods are less suited when the underlying process is strongly stochastic, or heavily corrupted by unknown processes (“noise”), because the embedding methodology is not robust to the presence of noise. Furthermore, each causal inference based on CCM looks for the behavior of binary connections, and it cannot infer how different causes work together to influence the target, which, as explained below, is a crucial motivation for the development of our methodology.

A recent surge of causal discovery methods originated in machine learning. The typical assumption is that information is available on interventions in a complex system, and large data sets are used to infer the average treatment effect, for instance, by including causal regularizers; see, e.g., Ref. 20 and references therein for many studies in these areas. Many studies search for causal features out of a finite set of features, while our focus is on processes where the causal feature set is continuous and hence infinitely large. The standard assumption is that the underlying processes can be represented on a standard Directed Acyclic Graph (DAG) as this is the underlying structure of a typical neural network, while, as argued below, graphs are too restrictive for the causal structures we are interested in. Machine learning has also been used to infer the influence of confounders on the causal net. The idea is that, although the confounder itself is not known, we do have proxy variables from which representations of the confounder influence can be estimated, for instance, via causal effect variational autoencoders.^21,22

One issue with the methods discussed above is that there are many examples in the real world where causes are nonlinear interactions between driver processes, and the above transfer-entropy-based and other methods cannot disentangle this properly; see, e.g., Refs. 14 and 23. A simple example is a transistor in which one process acts as a gate keeper for the connection between other processes. However, there is a more fundamental limitation of these methodologies.

All methods mentioned above can be represented on standard graphical networks, such as Bayesian networks or Markov-random fields, called causal graphs. Most theory and methodologies are based on so-called directed acyclic graphs (DAGs), in which each edge is an arrow, and there is at most one arrow between two vertices. These graphs represent the underlying joint probability of the system. However, these networks are designed to represent dyadic, so binary, interactions between the variables, while in many systems, the interactions are polyadic.^24,25 A simple example demonstrates the issue. Assume that driver $y=0$ if $z=0$ and $y$ can take on 0 or 1 with equal probability when $z=1$⁠. The target $x=yz+1−z−y+ηx$⁠, in which $ηx$ denotes random noise that is independent of $y$ and $z$⁠. All conditional mutual information terms of $x$ with any driver conditioned on the other driver are zero; therefore, a graph representation will consist of three nodes without edges. However, there is a nontrivial relation between the variables. One could argue that the problem is rather special in that it does not satisfy the faithfulness and causal Markov^6,26 conditions that are typically assumed in causal discovery, which together state that independency on the DAG means independency in the joint pdf and vice versa (see also the discussion in Ref. 27). However, the problem runs deeper. The issue is that many joint pdfs cannot be represented on a DAG, and hypergraphs are needed. Let us try to put this structure $x=yz+y+ηx$ on a graphical model, as in Fig. 1. The problem is that the combined influence of $y$ and $z$ cannot be represented by one edge, but edges need to be allowed to merge for which a hypergraph is needed. Many processes in nature are of the nonlinear form depicted in Fig. 1, for instance, advection terms in fluid dynamics, the interaction of radiation with matter, predator–prey biological systems, chemical reactions, etc. Standard graphical networks such as DAGs are not general enough to represent these processes. Our framework is not based on such graphical models and in fact is developed with this kind of polyadic interactions in mind.

A framework that can handle these polyadic interactions, at least in principle, is the interesting contribution by Williams and Beer.²⁸ They introduced a nonnegative decomposition of multivariate information, the so-called partial information decomposition (PID) that does allow for the inclusion of joint information; therefore, it does allow for polyadic interactions. The basic idea is that the total driver process information can be split into unique contributions $U$ from each driver, synergistic contributions $S$ with other driver processes, and redundant contributions $R$⁠. Redundant contributions are contributions to the target process that two or more driver processes have in common. These descriptions are rather vague, which allows for freedom, but also hampers applicability. For a system with only three processes, one target $x$ and two driver processes $y$ and $z$⁠, the mutual information between the target and the drivers is decomposed as

This system consists of three equations for the four unknown contributions and hence is underdetermined. The only general condition is that all four quantities have to be non-negative. We can eliminate the unique contribution $U$ by forming

and hence, the difference between $S$ and $R$ is defined in terms of mutual information but not each term individually. In information theory, this combination of mutual informations is minus the interaction information, and hence,

which can have any sign. Furthermore, from the basic PID equations and the conditional information relation $I(x;y|z)=I(x;y,z)−I(x;z)$⁠, we can derive

which shows that the conditional mutual information is interpreted as the sum of the unique and synergistic information in the PID framework.

Many definitions have been explored defining one of the variables in the PID framework and deducing the others from the framework, but all have their weaknesses. For instance, Barrett²⁹ showed that for dependent Gaussian source processes, three popular interpretations of PID,^28,30–32 all lead to the situation that the weakest source has zero unique contribution. However, this result is problematic. Suppose that we do have a weakest source process that has unique information on the target process, in the sense that it contains information on the target that none of the other driver processes have. All three PID interpretations mentioned above insist that unique contribution is zero, leading to a logical inconsistency.

Besides this, it is not clear if a unique contribution between a driver and target can be well defined in the first place. It can when the system can be decomposed on a graphical network, and in that case, conditioning out other processes as in transfer entropy logically provides the unique contribution. However, when a cause consists of nonlinear interactions between drivers, which is reality for many systems in the natural world, conditioning is insufficient to define a unique contribution. As mentioned above, conditioning can open gateways that are otherwise closed. The unique contribution of process $y$ on $x$⁠, $U(x;y|z)$⁠, is supposed to mean something like “when driver $z$ is not present.” However, $z$ is always present. We would need to find measures that exclude all influence of $z$⁠, but a general way to do that does not exist (unless we allow for interventions). These problems, to us, seem to point to serious issues with present-day interpretations of the PID formalism. That does not mean that the formalism is not useful, just that more work on its interpretation is needed.

In this paper, we provide a new causal discovery framework for non-intervenable systems that is unique in several ways. It does not rely on causal graphs, and hence, its application is not restricted by the issues discussed above. It is based on the notion of certainty instead of entropy. Certainty increases monotonically with the amount of information we have about a target process. Another reason to introduce this new concept is that it is always nonnegative for unbounded variables, while differential entropy can be negative, obscuring clear interpretations. The total mutual information of all driver processes with the target process is interpreted as the increase of certainty compared to having only the time series of the target process. We decompose the total mutual information in direct contributions from each driver to the target and joint contributions between two processes, between three processes, etc. We normalize the contribution of each process and define direct, joint, and total causal strengths from one process to another. By normalizing each contribution, different studies can be compared, and the certainty from the original time series of the target process, the so-called self-certainty, can be reinterpreted as the contribution from unknown processes. Hence, we can quantify the contribution of unknown processes (such as “noise”) and show that including new processes can only decrease this contribution from unknown processes. This is different from confounder influence estimation in machine learning, in which proxies for the confounder influence have to be present, while we need no extra information on the confounders.

This paper is organized as follows. In Sec. II, the basic ingredients of the new framework are introduced, followed by an example of how to decompose the mutual information when three processes are involved. Then, we show in Sec. IV the general theory of the decomposition of the total mutual information, discuss confounders in Sec. V, and apply the framework to several examples in Sec. VI. The paper is finalized by a discussion and concluding section.

The problem we want to solve is to identify the relative influence of a set of nonlinearly interacting random processes $yi$⁠, $i=1,2,…,N$ on a target random process $x$⁠. The subindex $i$ denotes a separate process, not a time index. Each of these processes themselves is a time series. Relative influence is defined as the extent to which process $yi$ increases our knowledge about process $x$⁠. Hence, we want to decompose our predictive knowledge about $x$ in its contributions from all processes $yi$ (which can include the past of any driver or even the past of $x$ itself), written symbolically as

where $(yi→x)$ contains all direct and joint contributions to $x$ in which $yi$ is involved. We will show in Secs. III and IV that this decomposition can indeed meaningfully be made with a proper definition of the arrow, even when the interaction of drivers with each other to influence the target is inseparable.

The full framework consists of a set of tools and operations as defined below. The framework defines a causal web that consists of driver nodes and a target node and directed edges (arrows) pointing into the target. An arrow can directly connect a driver node to the target or it can merge with one or more arrows from other driver nodes. Figure 1 is the simplest example, but more complex examples are provided later.

Each arrow that points to the target has a causal strength attached to it that describes the relative strength of the causal link between the one node it originates from or the relative causal strength due to the interaction of several nodes that the arrow originates from. The causal strength is normalized to enhance understanding and comparison with other systems. Subsections II A–II C describe how these relative causal strengths are calculated and how they should be interpreted. One special directed edge points to the target from unresolved processes. This can be noise or other physical processes that could be of interest, but we forgot to include in the driver set, the so-called confounders. The relative causal strength attached to this arrow allows for confounder detection, as explained in the chapter on confounders.

The causal strength calculations are based on mutual information and conditional mutual information among drivers conditioned on other drivers, as detailed below. In principle, the strength of every possible interaction between drivers toward the target needs to be calculated to obtain a full description of the causal structure of the system. Our method is minimally causally sufficient in accordance with definitions in Ref. 33. It is sufficient because for every target variable, the set of mlinks represents all possible causal influences due to combination of two or more drivers. Also, it is minimal as any subset of these would not be causally sufficient in general.

The time-lagged mutual information $I(x;y1:N)$ between a target process $x$ and a possible driver process $y$⁠, or a whole range of driver processes $y1:N$⁠, is defined via the Shannon entropy $H(..)$ as

where we assume a positive time lag between process $x$ and driver processes $y1:N$ for a causal link.

This lagged mutual information denotes the reduction in entropy of process $x$ when we condition on the processes $y1:N$⁠. We want to interpret the entropy in terms of information as in Shannon’s entropy, but we are interested in the case that each process lives on an unbounded domain. The differential entropy is defined as

where $p(x)$⁠, the probability density function (pdf) of a process $x$⁠, cannot be used because it can be negative, and it is not invariant under nonlinear transformations of the variables. Alternatively, we use the relative entropy, relative to a reference process with probability density $q(x)$ as

The relative entropy is positive for any choice of $q(x)$⁠, as long as its support is equal to or larger than that of $p(x)$⁠. This density $q(x)$ will provide an offset relative to $p(x)$⁠, the pdf of the process of interest. Although this offset density cancels in Eq. (6), it influences the size of our causal strengths between processes.

Many choices can be made for this reference density. Ideally, it has as little structure as possible, such as a uniform density. However, a uniform density does not exist on an unbounded domain. In this paper, we will mainly use the Cauchy or Lorentz pdf, given by

defined by width parameter $γ$ and mean $μx$⁠. A logical choice for the mean of the pdf is the sample mean. We choose $γ$ such that the reference density has the same entropy as the density with maximum entropy based on the mean and variance of the original process, so is the Gaussian. Since the entropy of the Cauchy distribution is $log⁡(4πγ)$ and that of the corresponding Gaussian $(1/2)log⁡(2πeσx2)$⁠, we can identify $γ=(e/8π)σx$⁠, where $e$ is the base of the natural logarithm. With this choice, the reference density can be interpreted as a maximum entropy pdf in the sense that it has the same maximum entropy as a Gaussian, but on top of that, it has infinite variance.

Other choices can be used too, e.g., a Gaussian with sample mean and variance or a uniform pdf with boundaries defined by the sample minimum and maximum values. We will discuss the influence of the reference density in Sec. VII-D.

In Sec. II A, we introduced the relative entropy as an important quantity in our framework. We now introduce a related quantity called certainty, defined as

This quantity is a relative entropy with a reference density that still has to be determined. We also introduce the unconditioned version, called the self-certainty in this context, as

which is also a relative entropy. The reference density is an important quantity for the size of the noise term and hence is related to confounder detection. It should be chosen as uninformative as possible, so uniform for discrete variables (which is not the focus of this paper), and a wide distribution, wider than the target, for continuous variables.

We have $0≤W(x)≤∞$⁠, with boundaries attained when $p(x)=q(x)$ or a delta Dirac function, respectively. This is in contrast to entropy, which is a measure of uncertainty. Hence, $W$ can be seen as a measure of certainty: given a wide reference density, the narrower the pdf of $x$⁠, the more we know about $x$⁠, and indeed the higher our certainty about $x$⁠. Similarly, for the conditional variant, we have $0≤W(x)≤W(x|y1:N)≤∞$⁠, as can easily be verified.

The certainty and the self-certainty are related through the information theoretic relation

as follows directly from the definition of the terms. Expression (12), which only contains non-negative terms, will be the basis for our causal inference. The term $W(x)$ denotes the amount of self-certainty we have on process $x$⁠. The mutual information term $I(x;y1:N)$ is the increase in certainty on $x$⁠, due to knowledge of $y1:N$⁠. $W(x|y1:N)$ denotes the information we have on process $x$ when we condition on processes $y1:N$⁠; therefore, we know what these processes $y1:N$ are doing.

Section II C will introduce normalization, which will allow for a more direct interpretation of the terms in the theory and will make different experiments comparable.

We can calculate the quantities above, but they would have little direct meaning. What does a mutual information of, say, 2.6 mean? Some meaning can be extracted if we compare what this value would mean for a standard process, such as a Gaussian, but if the process is far from Gaussian, e.g., multimodal, this explains very little. Since our quantity of interest is the relative contribution to the certainty in $x$ brought by each process, we normalize (12) by the certainty conditioned on all these processes, $W(x|y1:N)$⁠,

Using normalization by $W(x|y1:N)$⁠, we find as the relative influence of all processes $y1:N$ on process $x$ or the causal strength of processes $y1:N$ toward process $x$,

and hence,

This last equation is the same as (13) showing the contributions to $x$ by processes $y1:N$ and its self-certainty. The importance of the normalization is that now, we can compare different studies on causal discovery. Instead of having to infer if mutual information of, say, 2.6 is large or not, we know immediately if a causal strength of, say, 1/2, is large, this means that process contributes $50%$ to explaining the target process.

There is, however, another reason for introducing normalization. To understand the framework further, we assume that the underlying equation that governs process $x$ can be written as

for some function $g(..)$⁠, in which $η$ denotes all processes not included in $y1:N$⁠; therefore, all unresolved or unknown processes are typically considered noise. This assumption is completely general. The process $η$ is included because any real-world time series will always contain unknown or unresolved processes as well as observation noise; therefore, process $η$ does play a role in reality. If we would know the process $η$⁠, we could calculate $I(x;y1:n,η)$ and the result would be $∞$⁠. In that case, $W(x)$ would be insignificant compared to the mutual information. This suggests that the ratio between the self-certainty and the mutual information of the known processes $y1:N$ gives us a measure of how close we are in taking all relevant processes for $x$ into account. This ratio contains the same information as the ratio between $W(x)$ and $I(x;y1:n)+W(x)=W(x|y1:N)$⁠. This, then, suggests that the smaller $W(x)/W(x|y1:N)$⁠, the more complete the processes $y1:N$ are in the causal description of $x$⁠.

To clarify this further, assume we discover a new important process $w$⁠. Because $W(x|y1:N,w)=W(x|y1:N)+I(x;w|y1:N)$ and $I(x;w|y1:N)≥0$ because it is a bivariate mutual information, we have $W(x|y1:N,w)≥W(x|y1:N)$⁠. Since $W(x)$ does not change by incorporating $w$⁠, the ratio $W(x)/W(x|y1:N,w)$ will be smaller than $W(x)/W(x|y1:N)$⁠. This means that the more relevant driver processes we include, the smaller the ratio between the self-certainty and the certainty. Hence, we can attribute this ratio to unmodeled processes. We thus find that the normalization by $W(x|y1:N)$ changes the interpretation of the $W(x)$ term from self-information to the causal strength of unmodeled processes, and hence, we identify $cs(x;x)=cs(x;η)$⁠. We consider this a very useful property that other frameworks lack.³⁴ McGill develops this decomposition in the discrete setting and does notice that what we call self-certainty is related to what he calls noise. However, using this as a measure on the accuracy of the causal discovery is new. In Sec. V, we discuss the important case when the missed processes contain important information on the causal structure, the so-called confounders.

Now that we have defined the general framework a method to quantify the individual contributions $(yi→x)$ is developed. As an example of how we determine individual contributions to the target process $x$⁠, we first study the case of target process $x$ and two driver or source processes $y$ and $z$⁠. For completeness, we note that each of these processes $y$ or $z$ could be process $x$ itself but lagged in time. Equation (12) for three processes reads

and our task now is to decompose $I(x;y,z)$ into the contributions from $y$ and $z$⁠.

The influence of each process on $x$ can be divided in two contributions: a contribution when we fix the other process, which we will call the 1link contribution, and a correction to that. That correction by $y$ and $z$ together; therefore, a 2link contribution is often only partially taken into account. For instance, the situation depicted in Fig. 1, where the red arrow denotes a nonlinear interaction between the two drivers, is often ignored in the literature that bases the causal structure on DAGs (e.g., Refs. 6, 14, and 25) because that structure cannot be represented on such a graph. While ignoring this contribution might be useful for some systems, we will show in the examples that the present generalization is necessary for a full description of the causal network.

The conditional 1link contribution of process $y$ is found by conditioning on all other processes, so on process $z$ in this case. This means that we study the influence of $y$ on $x$ when the influence of $z$ has already been taken into account because it is given. This 1link contribution can be quantified by the conditional mutual information of $y$ to $x$ given process $z$⁠,

That the 1link can be considered the direct contribution of process $y$ on $x$ can be seen from the conditioned version of Eq. (17),

therefore, the increase in certainty of $x$ when $y$ becomes available, given that we know the influence of process $z$⁠. Similarly, for the 1link contribution from process $z$⁠, we find

The correction term of both contributions has to be related to the combined influence of $y$ and $z$ on $x$⁠. Since the full contributions of $y$ and $z$ should add up to $I(x;y,z)$⁠, as shown in Sec. II-B, the correction term has to be

This is half the total contribution of both processes, minus their conditional 1link contributions. If this term is positive, it can be interpreted as the contribution of the combination of $y$ and $z$ not contained in the conditional 1link contributions from $y$ to $x$ and from $z$ to $x$⁠, which can be termed the “synergy.” A more direct phrasing would be that this term provides a measure on how $y$ and $z$ enhance each other’s influence on $x$⁠. On the other hand, when it is negative, it can be seen as the “redundant” information in the conditional information. A more direct phrasing would be that the two drivers hinder each other’s influence on $x$⁠. Since this contribution is purely combined, i.e., it only acts when both $y$ and $z$ are active, the symmetry between $y$ and $z$ in this term demands that it must be divided equally between the two processes. Hence, the total contribution from $y$ to $x$ becomes

Using the standard relation $I(x;y,z)=I(x;z|y)+I(x;y)$⁠, we find

The quantity between the parentheses is known as the interaction information, defined as

Interaction information measures the influence of a variable $z$ on the amount of information shared between $x$ and $y$⁠, but it can do this in a non-intuitive way. For instance, when $y$ and $z$ are enhancing each other’s influence on $x$⁠, conditioning on $z$⁠, therefore, fixing $z$ can diminish this enhancement; therefore, $I(x;y)>I(x;y|z)$ and the interaction information is positive. On the other hand, $z$ can open a pathway between $y$ and $x$ that is not present without $z$⁠. In that case, one would expect $I(x;y)<I(x;y|z)$⁠; therefore, the 2link is negative. The fact that the 2link can be negative shows that one cannot identify the 1link with “unique” information from $y$ to $x$ in a PID framework interpretation because non-unique information should still be positive. For completeness, the total contribution from $z$ is

Now, we find the causal strength of $y$ to $x$ as

and similarly for $z$⁠. The unmodeled or noise relative contribution to $x$ is given by

leading to the total causal strength toward $x$ as

As mentioned above, a large portion of the previous literature on causal inference using standard graphs has systematically ignored the corrections to the “pure” 1link contributions. They thus missed potentially important parts of the causal network. It is true that the order of importance of processes $y$ and $z$ for $x$ will not change when the 2link is included as that term is the same for $z$ and $y$⁠. However, the ratio of the contributions will change. Furthermore, when more processes are present, 2links (and higher-order links) can change the order of importance compared to the 1link order and hence can lead to a completely different interpretation of the causal structure of the system. We will see examples of this later.

We can make the link to the PID framework by using (3) and (4) and decomposing our total contribution from $y$ to $x$ as

This suggests that the total contribution of $y$ to $x$ is a unique contribution and half the sum of the synergy and redundancy, all as defined in the PID framework. This makes sense if we invoke the symmetry argument that synergy and redundancy should only be included half for the $y$ contribution, with the other half for the $z$ contribution, but remember that none of the terms are defined uniquely in the PID framework. Given the difficulty in defining a unique contribution, our framework makes perhaps more sense than the PID framework. Our decomposition is based on the number and identity of the “active” (as opposed to conditioned on) variables in the mutual information, as explained further in Sec. IV.

When $N$ processes $yi$⁠, $i=1,2,…,N$ influence process $x$⁠, we can generalize the above as follows. To find the total contribution of each process $yi$⁠, we first quantify how much each of them contributes to $I(x;y1:N)$ on top of what all others contribute. Then, we quantify how much each process contributes in combination only with one other process. This is followed by how much each process contributes in combination only with two other processes, etc., until we reach how much each process contributes in combination only with all other processes. The word “only” is important as we have to avoid double counting. This leads to a decomposition of the total contribution of process $yi$ to $W(x|y1:N)$ as

Factors such as $1/2$ appear because each $2link$ process $yi,yj$ appears both in the contribution from $yi$ and in the contribution from $yj$⁠. Hence, this contribution needs to be distributed between these two process contributions. Since they both serve in equal capacity to this term, each process contributes 1/2 of this term. A similar argument holds for all higher-link terms in this decomposition. The decomposition of the total mutual information in terms of all the links was also noted by McGill,³⁴ who developed the theory for discrete systems, and Brown,³⁵ who worked in arbitrary spaces but neither did decompose this further into the contributions from each driver separately and hence neither did find our interpretation.

Each $mlink$ contains conditional mutual informations of the form $I(x;yi,z|w)$⁠, in which $z$ is a (m$−$1) subset of $y≠i$⁠, and $w$ contains those processes that are not process $yi$ and not in $z$⁠. This conditional mutual information contains all possible interactions between the active variables $yi$ and all variables in $z$⁠, including lower order links. To make sure this term only contains pure $mlinks$⁠, we need to subtract all links of lower order; therefore, $(m−1)links$⁠, $(m−2)links$⁠, etc., all the way to the conditional $1links$⁠, contained in the original $mlink$ set, to avoid double counting. An expression for the mlink of driver $yi$ can be written as follows. First, define $Ii$ as a set of $m−1$ non-overlapping indices from ${1,…,N}$ that do not contain $i$⁠, and $mi$ the set ${i,Ii}$⁠. Then, define $mi−1$ as a set of $m−1$ non-overlapping indices from the set ${mi}$ and $mi−2$ as a set of $m−2$ non-overlapping indices from the set ${mi−1}$⁠, etc. With these definitions, the mlink of driver $i$ in an N-driver system can be written recursively as

As an example, when three processes influence $x$ (⁠$N=3$⁠), we find for each $i$

in which the $2links$ are given by

and for the $1links$

Note the structure of this decomposition. Every term $I(x;yi,yj|yk)$ contains two 1link contributions that need to be subtracted to define the 2link. This is similar to what we did in Sec. III, but now with the extra conditioning on $yk$⁠. For the 3link, we have that every term $I(x;y1,y2,y3)$ contains both 2links and 1links that need to be subtracted. Indeed, all three possible 2links and all three possible 1links are subtracted.

Let us now evaluate this completely for the first $y$ process; therefore, $i=1$⁠. There is only one 1link, namely, $I(x;y1|y2,y3)$⁠. There are two 2links, namely, $y1,y2$ and $y1,y3$⁠. From each of them, we need to subtract the two 1links; therefore, in total, we need to subtract four 1links. Finally, there is only one 3link $y1,y2,y3$⁠, which will be decomposed as above. Hence, we find

Because of the symmetry of the $3links$ term, it is the same for all processes $yi$⁠. However, both the $1link$ and the $2links$ terms are dependent on the driver process under study. In general, for a system with $N$ drivers, all links smaller than the $Nlink$ will have links that are driver-process specific.

By grouping similar terms, it is possible to simplify the expression above further as

which now defines the total contribution of $y1$ to target $x$⁠, decomposed in its mlink contributions.

Adding all contributions from $y1$ to $y3$ together, we can show that indeed

as expected. It is straightforward to extend the decomposition for $N>3$⁠.

In terms of computational expense, the number of terms grows rapidly with the number of processes. However, two features of the theory keep the work manageable. First, the scheme is recursive, and second, the contributions from the different terms contain many terms that are the same. In fact, for $N=3$⁠, we need to calculate three terms of the form $I(x;yi|yj,yk)$⁠, three terms of the form $I(x;yi,yj|yk)$⁠, the term $I(x;yi,yj,yk,yl)$⁠, and $W(x^)$ (or $W(x^|y1:N)$ but that is more expensive to calculate), so 8 terms in total. It is easy to show that the number of terms to be calculated is equal to

This growth with the number of driver processes $N$ is exponential, but all mutual information calculations are independent and can be performed in parallel.

It is important to mention upfront that we do not need to calculate probability density functions in high-dimensional spaces but instead can use the time series directly in our calculation of the mutual informations by using the k-nearest-neighbor algorithm of Ref. 36. Before we discuss how the new framework deals with a few well-chosen systems that illustrate its strengths and weaknesses, we say a few words on how the system deals with confounders.

Missing processes are processes that are missed when potential drivers are identified and that if included would have a strong influence on the causal strength of one or several other drivers toward the target process. An example of such a process is a confounder, which is a missed process that drives both one or more identified drivers and the target process. Let us see how these effects are represented in the framework. Assume, for ease of notation, that the system contains three processes, a target $x$⁠, a known process $y$⁠, and a missing process $z$⁠.

In the new framework, there is a way to infer that a missing process is present via the causal strength of the noise term $cs(x;η)=W(x)/W(x|y)$ in this case. If this term is larger than expected from observational noise, it is likely that we missed an important driver. Practically, the effect from observational noise can be estimated by perturbing the target process series with three times the observational noise (so that the effective perturbation from the true process is twice as large) and recalculate the causal noise term and do this several times. If the noise term in these perturbed experiments remains the same, a missing process will be present. On the other hand, if the causal noise term is sensitive to these perturbations, it is unlikely that an important missing process is present.

We illustrate the idea in the following simple experiment. The underlying true system is described by

in which $εx∼N(0,10−4)$⁠, $εy,εz∼N(0,10−2)$⁠. Using 20 000 time steps, we calculate the causal strength of the noise without taking process $z$ into account as $cs(x;η)noz=0.301$⁠. We add a realization of noise to the target and recalculate the causal strength; repeating this ten times gives $cs(x;η)noz=0.308±0.005$⁠. The difference between the causal strengths before and after the noise perturbations is about 2%. This insensitivity to the observational noise suggests that there is a much larger contributor to the causal strength of the noise, and hence, we can assume that an important driver process has been missed.

After realizing that $z$ is an important process, we recalculate the causal strength of the noise as $cs(x;η)=0.086$⁠. Repeating the noise perturbation experiments for this case leads to $cs(x,η)=0.124$⁠, a change of about $50%$⁠. The larger difference between these two suggests that no further missing processes are present.

If the missing process is a confounder, it is a strong driver of both an identified driver and the target. The now standard way to define the presence of a confounder is via do-calculus,

in which $do(y)$ means that $y$ is given a certain value and all processes that influence $y$ are blocked. However, we cannot use this idea here because we are especially also interested in systems in which such an intervention would change the whole dynamics, resulting in a different system we are not interested in.

The definition used here is that $I(xn+1;y−∞:n|z−∞:n)=0$⁠, where $n$ is a time index and in which the conditioning is on the whole past of $z$⁠, and we consider the causal strength of the whole past of $y$ on $xn+1$⁠. On a DAG, this would mean that there is no direct arrow from $y$ to $x$ that does not go via $z$⁠: either $y$ influences $x$ via $z$⁠, i.e., $yn−1→zn→xn+1$ , or $z$ is a common driver of both $y$ and $x$⁠, i.e., $yn+1←zn→xn+1$⁠. In the first case, when $yn−1→zn→xn+1$⁠, we would find a nonzero causal strength between $yn−1$ and $xn+1$⁠, which is correct because $y$ does drive $x$⁠, via $z$⁠. The framework will give the correct answer in this case. One could argue that this is not a true confounder case because of this relation.

The second case, in which $yn+1←zn→xn+1$⁠, is more interesting. If we do not know about the existence of $z$⁠, the causal strength $cs(xn+1;yn)$ will be nonzero if any of the three processes has memory, meaning that its past is a driver of its present. This is because $xn$ and $yn$ are driven by $zn−1$⁠, and memory in either $x$⁠, $y$⁠, or $z$ will result in a connection between $yn$ and $xn+1$⁠. Without having knowledge of $z$⁠, this is a correct answer.

The question is what happens when $z$ comes to our attention. Remember that the confounding nature of $z$ on the $(x,y)$ relation is defined as $I(xn+1;y−∞:n|z−∞:n)=0$⁠. In this case, the total causal strength of the whole past of $y$ would be

Hence, this nonzero causal strength shows that the framework suggests that $y$ does influence $x$⁠, while in fact, $z$ is driving both $y$ and $x$ separately. However, if $z$ is driving $y$⁠, then $y$ has information from $z$⁠; therefore, $y$ has information that drives $x$⁠. The framework does recognize the information $y$ has on driving $x$⁠. Since the framework does not only calculate causal strength but also all mlink contributions including $I(xn+1;y−∞:n|z−∞:n)=0$⁠, it recognizes the confounding nature of $z$⁠.

This short discussion demonstrates how missing processes are identified and how subsequently confounders are identified in the causal framework. We have a way of estimating the influence on the target of missing processes relative to the observational noise in the system, allowing for identification of their presence, and we identify true confounders via 1links.

Since our framework is specifically developed with nonlinear systems in mind, we only compare to other methods that allow for nonlinear interactions. In Table I, a comparison of capabilities is provided. As mentioned in Sec. I, transfer entropy as introduced by Schreiber¹⁰ and further developed for high-dimensional systems by, e.g., Ref. 25, does only a partial decomposition, similar to a 1link in the new framework. The PID framework has no unique definition of unique, synergistic, and redundant information. As we will see in the first example in Sec. VII, the simplest three-variable system already leads to interpretation problems that suggest that this decomposition is not that useful. Only the new framework provides a full decomposition of the causal strength over its different contributions, has a unique normalization that does not depend on application, is complete in the sense that it allows for confounder detection, and generates complete causal webs that are much richer than DAGs because they allow for merging of directed edges. For two different causal web styles, see Figs. 2 and 8.

Our framework has some parallels with the framework that Runge developed in Ref. 14. His framework aims to answer the question how strong the indirect causal influence is of a process on a target process, where the direct causal influence is defined via a transfer entropy. Specifically, the paper concentrates on the specific influence of a process $y$ that is a few time steps in the past of the target process $x$ and where $y$ influences other processes $z$ that in their turn influence $x$⁠. The interaction information from $y$ via $z$ is defined as the mutual information of all paths between $y$ and $x$ minus the mutual information of all paths between $y$ and $x$ conditioned on process $z$⁠. The paper restricts the analysis to causal systems that can be represented by a DAG, while our framework is more general than that because we explicitly take nonlinear interactions between processes into account, which cannot be represented on a such a graph.

Several examples are discussed to illustrate the behavior of the new framework. We start with linear models with Gaussian noise, then discuss nonlinear models without interactions between the terms, followed by models with nonlinear interactions and finally the Lorenz 1963 model.

All information theoretic quantities were calculated using the k-nearest-neighbor algorithm of Ref. 36, where the number of nearest neighbors is set to $4$⁠, with little sensitivity to the actual number. To increase the numerical accuracy, the target and all drivers are transformed via the CDF of a Gaussian to a truncated Gaussian, specifically cutting 4.25% of the full distribution off each tail. The reference density for the certainty calculation is the Lorentz–Cauchy density. In Sec. VII-D, the influence of different reference densities is investigated in detail.

The following models are special in that their temporal memory is strongly limited, allowing us to concentrate on local-in-time relations. Furthermore, the models are simple enough so that they can be represented on a standard graph except for model 6. Table II shows the first three models that we used to generate time series, on which we then test the causal discovery framework. We generated 100 time series from each model of length 50 000 steps and calculated the mutual information and conditional mutual informations as needed. The results of the experiments are presented in Table III.

Model 1 is perhaps the most simple model one can think of and can be represented on a graph as $y→x←z$⁠. It is linear and has no memory; therefore, interpretation of the terms should be straightforward. The conditional mutual information, the 1link, is larger than the mutual information between $y$ and $x$ and between $z$ and $x$⁠. This means that the interaction information is negative, and the reason is that without conditioning, the variable $z$ acts as noise in the mutual information calculation of $y$ and $x$ and similarly for $y$⁠. The causal strength of $y$ to $x$ is 1.6 times larger than that of $z$ to $x$ (0.56/0.35), with a small contribution for the noise process. If only the 1links would be taken into account, the ratio of the $y$ contribution to the $z$ contribution is much lower, 1.1, due to the omission of the 2link contributions. In this simple model, the conditional mutual information is equal to the transfer entropy, and it is interesting to see how transfer entropy suggests very similar contributions from $y$ and $z$⁠, while the causal strength of the former is expected to be much larger considering the actual model, as correctly indicated by our framework.

It is also interesting to connect these results to the PID framework. The form of model 1 suggests that there is no synergy and no redundancy in this system because $y$ and $z$ are completely independent when driving $x$⁠. However, $I(x;y|z)>I(x;y)$⁠, and hence, $R(x;y,z)−S(x;y,z)<0$ in the PID framework so that the synergy has to be nonzero, as all contributions are non-negative in the PID framework. The only way to keep consistency in the PID framework is to introduce the influence of unaccounted-for processes (“noise”), as we do in our discussion above on the 2links, but PID does not include “noise” terms. The words “unique,” “synergy,” and “redundancy” are difficult to define even in this simple system. If we make the straightforward choice $U(x;y|z)=I(x;y|z)$ we find, see Eq. (4) in Sec. I, that $S(x;y,z)=0$⁠, which leads to $R(x;y,z)<0$⁠, which is inconsistent with the non-negativeness of the variables in the framework. If, instead, we choose $R(x;y,z)=0$⁠, the unique contribution of $y$ on $x$ becomes the mutual information $I(x;y)$⁠, while we know that this mutual information also contains influence from $z$⁠.

Model 2 is a system in which $z$ acts as a gateway for the information flow between $y$ and $x$⁠; graphically, $y→z→x$⁠. The only nonzero contributions in model 2 are those between $xn+1$⁠, $zn$⁠, and $yn−1$ since the system variables have no memory. This leads to positive interaction information because conditioning on $z$ in $I(x;y|z)$ destroys the connection between $y$ and $x$⁠. The mutual information between $y$ and $x$ is nonzero without this conditioning, showing that there is information flow from $y$ to $x$ in this system because knowing $y$ does provide information about $z$ and hence information about $x$⁠; the 1links are just unable to pick this up. As expected, the causal strength of $z$ to $x$ is much higher than that of $y$ to $x$⁠. The reason for this small $y$ contribution is that the noise $ηz$ is of the same order of the signal $y$⁠, making the $I(x;y)$ dominated by noise. This can be seen clearly when we substitute the expression for $zn$ in model 2: $xn+1=yn−1+ηzn+ηxn+1$⁠. Indeed, lowering the noise in $z$ does make $I(x;y)$ much larger (not shown). It is also interesting that if we remove $y$ from the causal calculations, the underlying model remains model 2, but if we only consider $x$ and $z$⁠, the total causal strength from $z$ to $x$ increases to $0.90$⁠. Hence, in this case, $z$ takes up the causal strength of $y$⁠, which is exactly what the framework should do.