In this article, we use the definition of ‘causality’ that is common in statistics and intuitive to scientists: $X$ has a causal effect on $Y$ (‘$X$ causes $Y$’ or ‘$X$ is a causer; $Y$ is a causee’ or ‘$X$ is a cause; $Y$ is an effect’) if some externally applied perturbation of $X$ can result in a perturbation in $Y$ (Figure 1A). We say that $X$ and $Y$ are causally related if $X$ causes $Y$, $Y$ causes $X$, or some other variable causes both. Otherwise, $X$ and $Y$ are causally unrelated. Additionally, one can talk about direct versus indirect causality (Figure 1B; see legend for definitions). A surprising result from past several decades of causality research is that there are in fact some conditions under which directional causal structures can be correctly inferred (‘identified’) from purely observational data (Spirtes and Zhang, 2016; Peters et al., 2017; Hitchcock, 2020b) (e.g. Appendix 2—figure 2, last row). However, empirical time series often do not contain enough information for easy causal identifiability (Spirtes and Zhang, 2016; Glymour et al., 2019).

Figure 1

Download asset Open asset

The adage ‘correlation is not causality’ is well-known to the point of being cliché (Sugihara et al., 2012; Coenen and Weitz, 2018; Carr et al., 2019; Mainali et al., 2019). Yet, to dismiss correlative evidence altogether seems too extreme. To make use of correlative evidence without being reckless, it helps to distinguish between the terms ‘correlation’ and ‘dependence’. When applied to ecological time series, the term ‘correlation’ is often used to describe some statistic that quantifies the similarity between two observed time series (Weiss et al., 2016; Coenen and Weitz, 2018). Examples include Pearson’s correlation coefficient and local similarity (Ruan et al., 2006). In contrast, statistical dependence is a hypothesis about the probability distributions that produced those time series, and has close connections to causality.

Dependence has a precise definition in statistics, and is most easily described for two binary events. For instance, if the incidence of vision loss is higher among diabetics than among the general population, then vision loss and diabetes are statistically dependent. In general, events $A$ and $B$ are dependent if across many independent trials (e.g. patients), the probability that $A$ occurs given that $B$ has occurred (e.g. incidence of vision loss among diabetics only) is different from the background probability that $A$ occurs (e.g. background incidence of vision loss). If $A$ and $B$ are not dependent, then they are called independent. The concept of dependence is readily generalized from binary events to numerical variables, and also to vectors such as time series (Appendix 1).

Dependence is connected to causation by the widely accepted ‘Common Cause Principle’: if two variables are dependent, then they are causally related (i.e. one causes the other, or both share a common cause; Peters et al., 2017; Runge et al., 2019a; Hitchcock, 2020b; Hitchcock and Rédei, 2020a). Note however that if one mistakenly introduces selection bias, then two independent variables can appear to be dependent (Appendix 2—figure 3). The closely related property of conditional dependence (i.e. whether two variables are dependent after statistically controlling for certain other variables; Appendix 1) can be even more causally informative. In fact, when conditional dependence (and conditional independence) relationships are known, it is sometimes possible to infer most or all of the direct causal relationships at play, even without manipulative experiments or temporal information. Many of the algorithms that accomplish this rely on two technical but often reasonable assumptions: the ‘causal Markov condition’, which allows one to infer causal information from conditional dependence, and the ‘causal faithfulness condition’, which allows one to infer causal information from conditional independence (Appendix 2; Peters et al., 2017; Glymour et al., 2019; Hitchcock, 2020b).

In sum, whereas a correlation is a statistical description of data, statistical dependence is a hypothesis about the relationship between the underlying probability distributions. Dependence is in turn linked to causality. Below, we discuss tests that use correlation to detect dependence in time series.

Despite its scientific usefulness, dependence between time series can be treacherous to test for. This is because time series are often autocorrelated (e.g. what occurs today influences what occurs tomorrow), so that a single pair of time series contains information from only a single trial. If one has many trials that are independent and free of systematic differences (e.g. $\ge 20$ as in some laboratory microcosm experiments), the task is relatively easy: One can test whether the abundances of species $X$ and $Y$ are statistically dependent by comparing the correlation between $X$ and $Y$ abundance series from the same trial with those between $X$ and $Y$ abundance series from different trials (Appendix 1—figure 4; see also Moulder et al., 2018). However, a large trial number is generally a luxury and often only one trial is available. In such cases, attempting to discern whether two time series are statistically dependent is like attempting to divine whether diabetes and vision loss are dependent with only a single patient (i.e. we have an ‘$n$-of-one problem’). As one possible remedy, there are parametric tests using the Pearson correlation coefficient that account for autocorrelation. In these tests, one estimates the correlation coefficient between time series, and evaluates its statistical significance using the variance of the null distribution (Afyouni et al., 2019). However, the calculation of this variance relies on estimates of the autocorrelation at each lag for both time series, which can be highly uncertain (Pyper and Peterman, 1998; Ebisuzaki, 1997). Furthermore, after estimating the variance, one must also assume the shape of the null distribution before a p-value can be assigned to the correlation.

Alternatively, the $n$-of-one problem is often addressed by a technique called surrogate data testing. Specifically, one computes some measure of correlation between two time series $X$ and $Y$. Next, one uses a computer to simulate replicates of $Y$ that might have been obtained if $X$ and $Y$ were independent (see below). Each simulated replicate is called a ‘surrogate’ $Y$. Finally, one computes the correlation between $X$ and each surrogate $Y$. A p-value (representing evidence against the null hypothesis that $X$ and $Y$ are independent) is then determined by counting how many of the surrogate ${\displaystyle Y}$s produce a correlation at least as strong as the real $Y$. For example, if we produced 19 surrogates and found the real correlation to be stronger than all 19 surrogate correlations, then we would write down a p-value of $1/(1+19)=0.05$. Ideally, if two time series are independent, then we should register a a p-value of 0.05 (or less) in only 5% of cases.

Several procedures can be used to produce surrogate time series, each corresponding to an assumption about how the original time series was generated (Lancaster et al., 2018). One popular procedure is to simply shuffle the values of a time series (Ruan et al., 2006; Eiler et al., 2012; Shade et al., 2013; Cyriaque et al., 2020). This procedure, often called permutation, assumes that all possible orderings of the time points in the series are equally likely. This assumption is commonly violated in time series due to autocorrelation, and thus the test is often invalid. For example, for independent time series in Figure 2B–C, this test returns ${\displaystyle p<0.05}$ at rates of $30\sim 92\%$, much higher than 5%. Nevertheless, permutation testing has appeared in many applied works, perhaps because it has been the default option in some popular software packages. Another procedure for generating surrogates is called phase randomization. It first uses the Fourier transform to represent a time series as a sum of sine waves, then randomly shifts each of the component sine waves in time, and finally sums the phase-shifted components (Ebisuzaki, 1997; Schreiber and Schmitz, 2000; Andrzejak et al., 2003; Appendix 3—figure 1). This procedure is considered appropriate when the original time series is obtained from a linear, Gaussian, and stationary process (Andrzejak et al., 2003; Lancaster et al., 2018), where ‘linear’ means that future values depend linearly on past values, ‘Gaussian’ means that any subsequence follows a multivariate Gaussian distribution, and ‘stationary’ means that this distribution does not change over time. See Chan, 1997 for a discussion of exact requirements. Indeed, this test performed well (with a false positive rate of 4%) when time series satisfied its assumptions (Figure 2C), and poorly when the stationarity assumption was violated (with a false positive rate of 21%; Figure 2B). Other surrogate data procedures include time shifting (Andrzejak et al., 2003), the block bootstrap (Papana et al., 2017), and the twin method (Thiel et al., 2006). Some surrogate data tests have been shown to perform reasonably well even when the exact theoretical requirements are unmet or unknown (Thiel et al., 2006; Papana et al., 2017), but a more comprehensive benchmarking effort is needed to map out each method’s valid domain in practice.

Figure 2

Download asset Open asset

In sum, surrogate data allow a researcher to use an observed correlation statistic to test for dependence under some assumption about the data-generating process. Dependence indicates the presence of a causal relationship, and conditional dependence can sometimes even indicate the direction (Hitchcock, 2020b; Glymour et al., 2019; Heinze-Deml et al., 2018; Appendix 2—figure 2). Below we consider Granger causality and state space reconstruction, two approaches that can be used to directly infer the direction of causality from time series.

In simple language, $X$ is said to Granger-cause $Y$ if a collection of time series containing all historical measurements predicts $Y$’s future behavior better than a similar collection that excludes the history of $X$. An important consequence of this definition is that Granger causality excludes indirect causes, as illustrated in Figure 3A. In practice, whether a causal relationship is direct or indirect depends on which variables are observed. For instance, in Figure 3A, if $Y$ were not observed, then $X$ would “directly” cause (and Granger-cause) $Z$.

Figure 3

Download asset Open asset

Granger causality has many related but nonequivalent quantitative incarnations in the literature, including several that were proposed by Granger himself (Granger, 1969; Granger, 1980). Box 1 presents two definitions: one based on a linear regression which we call ‘linear Granger causality’ (Gibbons et al., 2017; Ai et al., 2019; Barraquand et al., 2020; Mainali et al., 2019) and another more general definition which we call ‘general Granger causality’ (also sometimes called nonlinear Granger causality; Granger, 1980; Diks and Panchenko, 2006; Bekiros and Diks, 2008; Vicente et al., 2011; Roux et al., 2013; Papana et al., 2017). See theorem 10.3 of Peters et al., 2017 for a discussion of the theoretical relationship between general Granger causality and (true) causality.

We discuss four important instances where Granger causality can fail as an indicator of direct causality (Figure 3B). These pathologies can be understood intuitively and can apply to both linear and general Granger causality. First, if a system has deterministic dynamics (see Appendix 3), then Granger causality may fail to detect causal relations (Figure 3Bi). More generally, if dynamics have a low degree of randomness, Granger causality signals can be very weak (e.g. knowing $X$’s past improves predictions of $Y$’s future only slightly; Janzing et al., 2013; Peters et al., 2017). Moreover, as we will discuss later, this limitation has motivated other methods that take a primarily deterministic view (Sugihara et al., 2012). Second, Granger causality may erroneously assign a direct causal relation between a pair of variables that have an unobserved common cause (Figure 3Bii). Third, recording data at a frequency below that of the original process by ‘subsampling’ (e.g. taking weekly measurements of a daily process) or by ‘temporal aggregation’ (e.g. taking weekly averages of a daily process) can alter the inferred causal structure (Figure 3Biii), although recent techniques can help with these issues (Gong et al., 2015; Hyttinen et al., 2016; Gong et al., 2017). Lastly, when measurements are noisy (Figure 3Biv), Granger causality can assign false interactions and also fail to detect true causality (Newbold, 1978), although some progress has been made on this front (Nalatore et al., 2007).

One might still attempt to infer Granger causality despite the above caveats, especially in situations where caveats can be largely avoided. Linear Granger causality has standard parametric tests: if any of the ${\alpha }_k$ terms in Equation 1 is nonzero, then $X$ linearly Granger-causes $Y$. Parametric tests are computationally inexpensive and available in multiple free and well-documented software packages (Seabold and Perktold, 2010; Barnett and Seth, 2014). These tests assume that time series are ‘covariance-stationary’, which means that certain statistical properties of the series are time-independent (Barnett and Seth, 2014; see also Appendix 3), and can fail when this assumption is violated (Toda and Phillips, 1993; Ohanian, 1988; He and Maekawa, 2001). Additionally, applying linear Granger causality to nonlinear systems can lead to incorrect causal conclusions (Li et al., 2018). One can assess whether the linear model (Equation 1) is a reasonable approximation, for instance by checking whether the model residuals ${\epsilon }_t$ are uncorrelated across time (Feige and Pearce, 1979) as is assumed by Equation 1.

Tests for general Granger causality often use a statistic known as transfer entropy (Papana et al., 2012). Roughly, the transfer entropy from $X$ to $Y$ is the extent to which the entropy (a measurement of uncertainty) of $Y$’s future is reduced when we account for (specifically, condition on) the past of $X$ (Schreiber, 2000; Cover and Thomas, 2006; Montalto et al., 2014; Papana et al., 2017). A significant transfer entropy thus indicates the presence of general Granger causality. Surrogate data are typically used to evaluate significance (Montalto et al., 2014; Papana et al., 2017; Shorten et al., 2021). However, the previously discussed surrogate data procedures are designed to test the null hypothesis of independence, which is different from the null hypothesis of general Granger non-causality (i.e. Equation 2, but replace ‘≠’ with ‘=’). More recent surrogate procedures have been proposed to address this issue (Runge, 2018a; Shorten et al., 2021). Several software implementations of Granger causality tests based on transfer entropy statistics are available (e.g. Montalto et al., 2014; Behrendt et al., 2019; Wollstadt et al., 2019).

Granger causality methods face challenges when datasets have a large number of variables (e.g. in microbial ecology). In this case, the summation in Equation 1 will contain a large number of terms, and so a regression procedure may fail to detect many true interactions (Runge et al., 2019a; Runge et al., 2019b). To handle systems with many variables, one can impose the assumption that only a small number of causal links exist (Gibbons et al., 2017; Mainali et al., 2019). This is sometimes called sparse regression or regularization. Additionally, under certain technical assumptions, it is possible to use a series of logical rules to remove unnecessary terms in a purely data-driven way (Runge et al., 2019b; Runge et al., 2019a). As an example, suppose that we wish to test whether pH is a Granger-cause of chlorophyll concentration in some aquatic environment and we infer based on a prior analysis that chlorophyll concentration is always independent of fluctuations in salinity. Then, most likely, salinity is irrelevant to the pH-chlorophyll relationship and can be safely omitted from our Granger causality analysis. As an aside, this reasoning could theoretically fail in pathological cases where, for instance, the ‘faithfulness’ condition (Appendix 2) is violated (see Example 7 of Runge, 2018b for a worked counterexample). These rules and their associated assumptions are formalized in ‘constraint-based’ causal discovery algorithms (Appendix 2; Peters et al., 2017; Glymour et al., 2019). The development of new causal discovery algorithms, and their application to time series, is a very active area of research (Hyvärinen et al., 2010; Runge et al., 2019b; Runge et al., 2019a; Sanchez-Romero et al., 2019).

Consider the deterministic dynamical system in Figure 4. Here, $Z$ is causally driven by ${\displaystyle X}$ and ${\displaystyle Y}$, but not by $W$ or $V$. We can make a vector out of the current value $Z(t)$ and two past values $Z(t-\tau )$ and $Z(t-2\tau )$, where τ is the time delay and $[Z(t),Z(t-\tau ),Z(t-2\tau )]$ is called a ‘delay vector’ (Figure 4B, red dots). The delay vector can be represented as a single point in the three-dimensional $Z$ ‘delay space’ (Figure 4C, red dot). We then shade each point of the trajectory in $Z$ delay space according to the contemporaneous value of ${\displaystyle Y}$, which causally influences ${\displaystyle Z}$. Since in this example each point of the trajectory in $Z$ delay space corresponds to one and only one $Y(t)$ value, we call this a ‘delay map’ from $Z$ to $Y$. Notice that the $Y(t)$ gradient in this plot looks gradual in the sense that if two points are nearby in the delay space of $Z$, then their corresponding $Y(t)$ shades are also similar. This property is called ‘continuity’ (Appendix 4—figure 1). Overall, there is a continuous map from the $Z$ delay space to $Y$, or more concisely, a ‘continuous delay map’ from $Z$ to $Y$. A similar continuous delay map also exists from $Z$ to its other causer $X$. On the other hand, if we shade the delay space of $Z$ by $W$ or $V$ (neither of which causes $Z$), we do not get a continuous delay map (Figure 4D–E).

Figure 4

Download asset Open asset

In this example, there is a continuous delay map from a causee to a causer, but not the other way around, and also no continuous delay map between causally unrelated variables. If this behavior reflects a broader principle, then perhaps continuous delay maps can be used to infer the presence and direction of causation. Is there in fact a broader principle?

In fact, there is a sort of broader principle, but it may not be fully satisfying for causality testing. The principle stems from a classic theorem due to Takens, 1980. A rough translation of Takens’ theorem is the following: If a particle follows a deterministic trajectory which forms a surface (e.g. an ant crawling all over a doughnut), and if we take one-dimensional measurements of that particle’s position over time (e.g. the distance from the ant’s starting position), then we are almost guaranteed to find a continuous delay map from our measurements (of current distance) to the original surface (the donut), as long as we use enough delays. (We walk through visual examples of these ideas in detail in Appendix 4.) A key result that follows from this theorem is that we can typically (‘generically’) expect to find continuous delay maps from ‘dynamically driven’ variables to ‘dynamically driving’ variables in a coupled deterministic dynamical system, as long as certain technical requirements are met (Cummins et al., 2015). Although the notion of ‘dynamic driving’ (Cummins et al., 2015) differs from our definition of causation, the two are related and we will still use the standard notion of causation when evaluating the performance of SSR methods. In theory, Takens’ theorem says that almost any choice of delay vector should work as long as it contains enough delays. However in practice, with finite noisy data, the behavior of SSR methods can depend on the delay vector selection procedure (Cobey and Baskerville, 2016; see also Appendix 4). Overall, Takens’ theorem and later results (Sauer et al., 1991; Cummins et al., 2015) form the theoretical basis of SSR techniques.

SSR techniques attempt to detect a continuous delay map (or a related feature) between two variables and use this to infer the presence and direction of causation (Sugihara et al., 2012; Ma et al., 2014; Harnack et al., 2017): A continuous delay map from $Y$ to $X$ is taken as an indication that $X$ causes $Y$. The fact that the map points in the opposite direction as the expected causation is potentially counterintuitive. One informal explanation is that the delay vectors of the causee can contain a record of past influence from the causer (Sugihara et al., 2012). As a word of warning, while causation is one possible explanation for a continuous delay map, it is not the only possible explanation. Indeed, we now illustrate scenarios where a causal relationship and a continuous delay map do not coincide.

Figure 5 illustrates four failure modes of SSR. In the first failure mode, which we refer to as ‘nonreverting continuous dynamics’ (top row of Figure 5; see also Appendix 4), a continuous map arises from the delay space of $X$ to $Z$ because a continuous map can be found from the delay space of $X$ to time (‘nonreverting $X$’) and from time to $Z$ (‘continuous $Z$’). This pathology leads to false causal conclusions and may explain apparently causal results in some early works where SSR methods were applied to data with a clear temporal trend. We are not aware of statistical tests for this problem, but Clark et al., 2015 recommend shading points in the delay space with their corresponding time to visually check for a time trend. In the second failure mode (Figure 5, second row; see also Sugihara et al., 2012), one variable drives another variable in such a way that the dynamics of the two variables are synchronized. Consequently, although the true causal relationship is unidirectional, bidirectional causality is inferred. Although the ‘prediction lag test’ (Figure 6B right panel) can sometimes alleviate this problem (Ye et al., 2015; Cobey and Baskerville, 2016), it is not foolproof as we demonstrate in Appendix 4. In the third failure mode (Figure 5 third row), $X$ and $Z$ both oscillate and $X$’s period is an integer multiple of $Z$’s period. In this case, $Z$ is inferred to cause $X$ even though they are causally unrelated (see also Cobey and Baskerville, 2016). In the fourth failure mode (Figure 5, bottom row), SSR gives a false negative error due to ‘pathological symmetry’, although this may be rare in practice.

Figure 5

Download asset Open asset

Figure 6

Download asset Open asset

SSR causal discovery methods require testing for the existence of continuous delay maps between variables. However, testing for continuity in real data is complicated by noise and discrete sampling (Figure 4, compare panels C and F; see also Appendix 4—figure 1).

Several methods have been used to detect SSR causal signals by detecting approximate continuity (Cummins et al., 2015) or related properties (Sugihara et al., 2012; Ma et al., 2014; Harnack et al., 2017). The most popular is convergent cross mapping (CCM), which has been applied to nonlinear (Sugihara et al., 2012) or linear deterministic systems (Barraquand et al., 2020). CCM is based on a statistic called ‘cross map skill’ that quantifies how well a causer can be predicted from delay vectors of its causee (Figure 6A), conceptually similar to checking for gradual transitions when shading the causee delay space by causer values (Figure 4). Four criteria have been proposed to infer causality (Sugihara et al., 2012; Ye et al., 2015; Cobey and Baskerville, 2016; Figure 6B): First, the cross map skill must be positive. Second, the cross map skill must be significant according to some surrogate data test. Third, the cross map skill must increase with an increasing amount of training data. Lastly, the cross map skill must be greater when predicting past values of the causer than when predicting future values of the causer (the prediction lag test [Ye et al., 2015; Cobey and Baskerville, 2016] in the right panel of Figure 6B, but see Appendix 4 for caveats of this test). In practice, many if not most CCM analyses use only a subset of these four criteria (Sugihara et al., 2012; Brookshire and Weaver, 2015; Cramer et al., 2017; Wang et al., 2018). Other approaches to detect various aspects of continuous delay maps have also been proposed (Ma et al., 2014; Cummins et al., 2015; Harnack et al., 2017; Leng et al., 2020). We do not know of a systematic comparison of these alternatives.

The concepts of dependence and independence between random variables are central to many statistical methods, including those that concern causality. A random variable is a variable whose values or experimental measurements depend on outcomes of a random phenomenon and follow a particular probability distribution. Reichenbach’s common cause principle states that if $X$ and $Y$ are random variables with a statistical dependence (such as a nonzero covariance), then one or more of three statements is true: $X$ causes $Y$, $Y$ causes $X$, or a third variable $Z$ causes both $X$ and $Y$. The common cause principle cannot be proven from the axioms of probability; rather, the principle is itself a fundamental assumption that supports much of the modern statistical theory of causality (Section 1.4.2 of Pearl, 2000).

As an example, consider the size and length of a bacterial cell. If a larger cell tends to be longer, then cell volume and cell length covary and are thus dependent. A mathematical definition of dependence (and its opposite, independence) is presented in Appendix 1—figure 1B.

Appendix 1—figure 1

Download asset Open asset

Dependence can be readily generalized from the definition in Appendix 1—figure 1 to become a property between two vectors of random variables. (Note that a time series can be viewed as a vector of random variables.) For example, suppose that we measure two variables $X$ and $Y$ over two days. Our (very short) time series are then $[X_1,X_2]$ and $[Y_1,Y_2]$ where the subscript index denotes the day of measurement. Similar to Appendix 1—figure 1B, we would say that our two time series are independent if

for all choices of $x_1,x_2,y_1,y_2$.

Many statistical techniques require repeated measurements that can be modeled as independent and identically distributed (IID) random variables, and passing non-IID data (such as time series) into such techniques can lead to spurious results (e.g. Figure 2; see also Koplenig and Müller-Spitzer, 2016). Random variables are IID if they have the same probability distribution and are independent (Appendix 1—figure 1). In Appendix 1—figure 2 we give examples of pairs of random variables that are (or are not) identically distributed, and that are (or are not) independent. Note that two dependent random variables can be linearly correlated (Appendix 1—figure 2, 3rd column), or not (Appendix 1—figure 2, 4th column).

Random sampling from a population with replacement is one way to produce “IID data” (which we use as a shorthand for “data which can be modeled as IID random variables”). For example, repeatedly rolling a standard die can be thought of as randomly sampling from the set ${\displaystyle \{1,2,3,4,5,6\}}$ with replacement: if the first trial registers 1, then the second trial can register one as well. Otherwise, if sampling was done without replacement, then the second trial must not register 1, which means that the outcome of the second trial would depend on the outcome of the first trial.

Appendix 1—figure 2

Download asset Open asset

Since the IID concept is so central to statistical analysis, we wish to further clarify one conceptual difficulty that may arise. To set the stage, suppose that a scientist measures the levels of voluntary physical activity in a collection of mice that includes both males and females. Also suppose that female mice tend to be more physically active than male mice (Rosenfeld, 2017). Since this dataset now contains measurements from both the less active males and the more active females, we might naively think that these data cannot be IID.

In fact, such a dataset still might be IID, but this depends on how the scientist chooses which mice to measure. To illustrate this fact, consider the highly simplified scenario in which only two mice are assayed for physical activity. Let ${\displaystyle X_1}$ and ${\displaystyle X_2}$ be random variables that describe the activity levels of these two mice. We consider three different ways that the scientist might select which mice to assay. Only one of these ways will result in an IID dataset.

First, suppose that the scientist chooses to measure X₁ from a male mouse and X₂ from a female mouse. In this case, to see whether X₁ and X₂ are IID, we can use the same visualization strategy as in Appendix 1—figure 1. That is, we imagine many possible ‘parallel universes’, each with a different possible two-mouse dataset (left panel of Appendix 1—figure 3). This allows us to visualize the joint distribution of X₁ and X₂. We can then see that X₁ and X₂ are independent, but not identically distributed.

Appendix 1—figure 3

Download asset Open asset

Second, suppose that the scientist again selects exactly one mouse of each sex, but randomizes the order so that both X₁ and X₂ have an equal chance of being measured from a male or female mouse (middle panel of Appendix 1—figure 3). We can now see that X₁ and X₂ are identically distributed, but not independent.

Lastly, suppose that the scientist selects mice randomly, and without any information about whether a mouse is male or female. In this case, the two-mouse sample might be all male, all female, or have one of each. Once again we plot the joint distribution of X₁ and X₂ by imagining their values across many different parallel universes (right panel of Appendix 1—figure 3). We then see that that X₁ and X₂ are finally independent identically distributed. Overall, a set of measurements can be IID even if they are taken from a mixed population, as long as they are sampled randomly from among different subpopulations.

Here, we first restate the concept of independence in terms of statistical conditioning, and then introduce the related concept of conditional independence.

It is intuitive that two variables are independent if knowledge of one variable tells us nothing about the other. The statistical notion of independence captures this intuition: Random variables $X$ and $Y$ are independent if the conditional distribution of $X$ given $Y$ is always equal to the marginal distribution of $X$. For discrete random variables, this condition can be written

or equivalently written ${\displaystyle P(X=x,Y=y)=P(X=x)P(Y=y)}$ for all $x$ and $y$. For continuous random variables, independence can be written in terms of probability density functions as $f_X(x|y)=f_X(x)$ or equivalently, $f_{X,Y}(x,y)=f_X(x)f_Y(y)$ where $f_X(x|y)$ is the conditional density of $X$ given $Y$, $f_{X,Y}(x,y)$ is the joint density of $X$ and $Y$, and $f_X(x)$ and $f_Y(y)$ are the marginal densities of $X$ and $Y$, respectively.

The statement “$X$ and $Y$ are conditionally independent given $Z$” intuitively means that $X$ and $Y$ are independent when we only analyze outcomes where $Z$ has a certain value. For discrete random variables, this condition is written $P(X=x|Y=y,Z=z)=P(X=x|Z=z)$, or equivalently, $P(X=x,Y=y|Z=z)=P(X=x|Z=z)P(Y=y|Z=z)$, for all $x$, $y$, and $z$. For continuous random variables, we have a similar formulation except that the probability $P$ is replaced by the probability density $f$ (i.e. $f_{X,Y}(x,y|z)=f_X(x|z)f_Y(y|z)$ for all $x,y,z$). If $X$ and $Y$ are not conditionally independent given $Z$, then $X$ and $Y$ are conditionally dependent given $Z$.

One could be forgiven for worrying about the feasibility of testing for dependence between long time series. This is because as a time series grows longer, the amount of data needed to get a sense of its probability distribution would seem to grow extremely rapidly. Thus, when $X$ and $Y$ are vectors that represent long time series, estimating the distributions in Equation 3 seems unrealistic. However, establishing that two time series are dependent only requires that we show that the distributions on the left and right sides of Equation 3 differ. Showing that two distributions differ can be much easier than actually estimating those distributions. For instance, if we know that the averages of two univariate distributions are different, then we immediately know that the two distributions are not the same, even if we know nothing about their shapes. Indeed, Appendix 1—figure 4 demonstrates a way to test for dependence between time series with only a moderate number of replicates, and without any assumptions about the shapes of the distributions. Additionally, surrogate data methods can be used to test for dependence with only one replicate of each time series, as discussed in the main text.

Many theoretical results and data-driven methods for causal analysis begin by representing causal relationships as a directed acyclic graph (DAG). That is, one makes a graph by represeting random variables as nodes and by drawing a directed edge from each direct cause (or parent) to its causee (or child), as in Figure 1B; additionally, the graph is acyclic, meaning that that it does not contain any directed paths from any variable back to itself. The acyclicity condition is often required for nice theoretical properties and ease of analysis (Spirtes and Zhang, 2016). Additionally, when data are temporal, a particular node in the graph commonly refers to a particular variable measured at a particular time (e.g. chapter 10 of Peters et al., 2017). If we follow this convention and note that causation cannot flow backward in time, and if we additionally exclude instantaneous causation, then our causal graph will be acyclic, even for systems with feedback (Appendix 2—figure 4).

DAGs are useful visual tools in their own right, but for many purposes we need to be more mathematically precise about what we mean when we draw an edge from one variable to another. Thus, often one interprets a causal DAG as corresponding to a set of equations with the following two conditions: First, each variable can be written as a function of (only) the variable’s direct causers and a random process noise term unique to the variable. Models that satisfy this condition are called structural equation models (SEMs) (Hitchcock, 2020b). Second, all process noise terms are (jointly) independent of one another. SEMs that satisfy this second condition are called Markovian and have a useful property called the ‘causal Markov condition’ (Pearl, 2000). (Some authors [Peters et al., 2017], but not all [Hitchcock, 2020b], require that all SEMs be Markovian by definition.) The causal Markov condition, along with the related ‘causal faithfulness condition’ are key assumptions that allow one to connect statistical structure to causal structure and infer aspects of causal structure from data, even in observational settings.

The causal Markov condition states that if there is no path from ${\displaystyle X}$ to ${\displaystyle Y}$ in a DAG (i.e. we cannot go from ${\displaystyle X}$ to ${\displaystyle Y}$ by following a sequence of edges in the forward direction), then ${\displaystyle X}$ and ${\displaystyle Y}$ are conditionally independent given ${\displaystyle X}$'s parents (Pearl, 2000; Zhang and Spirtes, 2008). In this context ${\displaystyle Y}$ can be either a variable or a set of variables. As an example, consider the boxed DAG in Appendix 2—figure 2B. Here, ${\displaystyle X}$ and ${\displaystyle Y}$ share the common cause ${\displaystyle Z}$. Each variable depends on its parents, and on its own process noise term. Although ${\displaystyle X}$ and ${\displaystyle Y}$ are dependent, the causal Markov condition expresses the intuitive idea that if we were to control for ${\displaystyle Z}$, then ${\displaystyle X}$ and ${\displaystyle Y}$ would become independent. Note that if ${\displaystyle X}$ does not have any parents, then the statement “${\displaystyle X}$ and ${\displaystyle Y}$ are conditionally independent given ${\displaystyle X}$'s parents” reduces to “${\displaystyle X}$ and ${\displaystyle Y}$ are independent”.

The causal faithfulness condition is, like the Markov condition, very useful in causal discovery and often quite reasonable. However, faithfulness is more difficult to state precisely and concisely without first introducing technical notation such as ‘$d$-separation’ (as in definition 6.33 of Peters et al., 2017). We attempt to give the gist of the idea here and direct readers to other sources (Peters et al., 2017; Zhang and Spirtes, 2008) for more precise definitions. The causal faithfulness condition is a kind of converse to the causal Markov condition. Recall that the causal Markov condition requires certain conditional (or unconditional) independence relationships based on the causal graph structure. Let us call any other independence relationships (i.e. those not required directly or indirectly by the causal Markov condition) ‘extra’ independence relationships. The (joint) probability distribution of random variables is causally faithful to the DAG if no ‘extra’ independence relationships exist (Hitchcock, 2020b). An imprecise shorthand for the faithfulness condition is ‘independence relationships indicate the absence of certain causal relationships’. The faithfulness condition can be violated when two effects precisely cancel each other (Appendix 2—figure 1).

Existing observational causal discovery methods for the IID (e.g. non-temporal) setting are diverse. Such methods can differ greatly in the assumptions they make (e.g. whether there are hidden variables or ‘unknown shift interventions’), the reasoning they employ, and the resolution of causal detail they provide (e.g. a unique causal graph versus a set of several plausible graphs) (Heinze-Deml et al., 2018). We will briefly introduce two classes of causal methods: (1) constraint-based search and (2) structural equation models (SEMs) with assumptions about the functional forms of equations (Spirtes and Zhang, 2016). However, these two classes, while illustrative of different modes of causal discovery, are far from an exhaustive list (Heinze-Deml et al., 2018).

Constraint-based search uses independence and dependence relationships (and their conditional counterparts) to narrow down the scope of possible causal graphs without exhaustively checking all possibilities (which can be enormous in number even for a handful of variables). The PC algorithm (named after its inventors Peter Spirtes and Clark Glymour) and the fast causal inference algorithm are examples of constraint-based search methods (Glymour et al., 2019). However, constraint-based methods often find multiple graphs that are consistent with the same set of data (e.g. Appendix 2—figure 2Biii, see legend; see also Spirtes and Zhang, 2016).

Functional form-based (or SEM-based) approaches to causal discovery begin by assuming a particular functional form for causal relationships, and then assess a given causal hypothesis by inspecting the joint distribution between a potential causer and its potential causee (Spirtes and Zhang, 2016). These methods rely on the fact that in a Markovian SEM, each variable has a noise term that is independent of the noise terms of all other variables (Peters et al., 2017). Given two dependent variables with no hidden common causes, one can use an appropriate regression to estimate values of a proposed causee based on the proposed causer (Spirtes and Zhang, 2016). If the residuals of this regression are independent from the proposed causer, then the proposed causal direction is consistent with the data (Hoyer et al., 2008). Crucially, theoretical results indicate that for a fairly wide variety of scenarios (e.g. linear non-Gaussian and post-nonlinear models), we can expect the data to be consistent with only one causal direction, thus enabling unambiguous identification of the causal direction (Spirtes and Zhang, 2016). An illustrative graphic example is given in Figure 3 of Spirtes and Zhang, 2016 and also in Figure 3 of Glymour et al., 2019. Similar ideas can be applied to multivariate systems (Hoyer et al., 2008; Peters et al., 2012).

Appendix 2—figure 1

Download asset Open asset

Appendix 2—figure 2

Download asset Open asset

Appendix 2—figure 3

Download asset Open asset

Appendix 2—figure 4

Download asset Open asset

Many methods for time series analysis require that data satisfy some stationarity condition, meaning that certain statistical properties of the data must remain constant across time. Two important stationarity conditions are strong stationarity and covariance-stationarity. A stochastic process ${\displaystyle X_t}$ is called strongly stationary if the joint distribution of any ${\displaystyle k}$ consecutive values ${\displaystyle (X_t,X_{t+1},\cdots ,X_{t+k-1})}$ is independent of time ${\displaystyle t}$(Definition 20.1 of Greene, 2012). A stochastic process ${\displaystyle X_t}$ is covariance-stationary (or weakly stationary) if: (1) the ensemble mean ${\displaystyle \mathbb{E}[X_t]}$ is finite and does not depend on $t$; (2) the variance $\mathrm{Var}[X_t]$ is finite and does not depend on $t$; (3) for all choices of $h$, the covariance $\mathrm{Cov}(X_t,X_{t+h})$ is finite and does not depend on $t$ (Definition 20.2 of Greene, 2012).

As an illustrated example of a covariance-stationary process, consider a population whose dynamics are governed by death and stochastic migration (similar to Figure 2C):

Here, X_t is the population size at time ${\displaystyle t}$, $a$ is the proportion of the population lost due to death during one time step, $c$ is the average number of individuals migrating into the population during one time step, and ${ϵ}_t$ is a random variable with a mean of zero which represents temporal fluctuations in the number of migrants. Suppose that we observed the dynamics of 10 populations governed by Equation 4 such that the populations all have the same parameters, but are independent (Appendix 3—figure 2A). Then, at each time point $t$, we will have some distribution of values of X_t. In fact, if we have not just 10, but 1,200 replicates, we can see that the distribution of values of X_t does not appear to depend on time (Appendix 3—figure 2B, top). Furthermore, the covariance between X_t and $X_{t+1}$ does not appear to depend on time either (Appendix 3—figure 2B, bottom).

Appendix 3—figure 2

Download asset Open asset

Although it is common to talk about a time series being stationary or nonstationary, this is technically a slight abuse of language. Just as the mean and variance are properties of a random variable (and not of any single data point obtained from that random variable), stationarity is a property of a stochastic process (and not of any single time series produced by that process). This fact is illustrated by comparing the middle and bottom rows of Appendix 3—figure 3. If we examine any one time series from the middle or bottom rows (e.g. the black curves in each), we see that they have essentially the same dynamics (i.e. they are sine waves with the same frequency). However, the process shown in the middle row is covariance-stationary (as shown below), whereas the process shown in the bottom row is not since its mean changes over time.

Appendix 3—figure 3

Download asset Open asset

To see that the middle row of Appendix 3—figure 3 shows a covariance-stationary process, we can show that the mean, variance, and covariance of the process are independent of time:

A deterministic time series from a system with many variables can be approximated as stochastic. This is illustrated below in Appendix 3—figure 4. When we track the trajectory of a particle in a box with 99 other particles (Appendix 3—figure 4 bottom row), the observed trajectory appears random, even though the governing equations of motion are deterministic. In particular, the motion of our particle over each time step can be approximated as having a random component. Note that this flavor of randomness is in general different from the phenomenon called chaos. In chaotic dynamics, each time step needs not be random, but small changes in initial conditions lead to large changes at later times.

Appendix 3—figure 4

Download asset Open asset

For panel B, we simulated the random walk system

where the $ϵ(t)$ terms were drawn independently from a normal distribution with mean of 0 and standard deviation of 1. We simulated this system from the initial condition of $X(1)=0$ through 999 subsequent steps. For panel C, we simulated the autoregressive system

where the $ϵ(t)$ terms were again drawn independently from a normal distribution with mean of 0 and standard deviation of 1. We simulated this system from the initial condition of $X(1)=40$ for 1999 subsequent steps. We only used the final 1000 steps for computing the correlation between two time series.

To compute the significance of the Pearson correlation between two time series, we used surrogate data generated by either permutation or the random phase procedure. Permutation surrogate time series were generated by randomly shuffling data. Random phase surrogate time series were generated by Ebisuzaki’s random phase method (Ebisuzaki, 1997) as implemented in the rEDM (version 1.5) function make_surrogate_data. For a pair of time series $[X_1(1),X_1(2),\mathrm{\dots },X_1(1000)],[X_2(1),X_2(2),\mathrm{\dots },X_2(1000)]$, we first computed the Pearson correlation $\widehat{\rho }$ between the two time series. We then replaced the X₂ values with surrogate time series and recomputed the Pearson correlation as $\tilde{\rho }$. We computed this shuffled correlation 9,999 times (permutation) or 499 times (random phase) to get a null distribution $[{\tilde{\rho }}_1,{\tilde{\rho }}_2,\mathrm{\dots },{\tilde{\rho }}_n]$. Following Schreiber and Schmitz, 2000, we computed the $p$ value as

where $N_{surr}$ is the number of surrogates, $N_{stronger}$ is the number of surrogate correlations $\tilde{\rho }$ whose magnitude was greater than or equal to the magnitude of the original correlation $\widehat{\rho }$, and the “$+1$” terms account for the original correlation $\widehat{\rho }$.

The system of equations was numerically integrated using the ode45 method in Matlab from $t=0$ to $t=200$ in time steps of 0.03, and plotted in the delay space $Z$ with $\tau =3.6$. The initial condition for all state variables ($V$, $W$, $X$, $Y$, $Z$, $\frac{dX}{dt}$, $\frac{dY}{dt}$, and $\frac{dV}{dt}$) was 1. For panel F, measurement noise was added to $Y(t)$. Specifically, noisy data were generated as:

where $\mathrm{Unif}(a,b)$ is a uniform random variable bounded by $a$ and $b$, and ${\mathrm{\Delta }}_Y$ is the difference between the maximum and minimum values of $Y(t)$ between $t=0$ and $t=200$. These noise parameters are chosen so that $Y^{obs}(t)$ is centered at $Y(t)$ and has a standard deviation of $0.15{\mathrm{\Delta }}_Y$ .

The dynamics in the top row of Figure 5 were generated from the equations:

This continuous-time system was discretized from $t=1$ to $t=20$ on an evenly spaced grid of 400 data points for visualizing delay spaces where the time delay is 50 time points (i.e. $\tau =50(20-1)/(400-1)$).

The dynamics in the second row of Figure 5 were generated from the equations:

with initial conditions of $X(1)=0.4$ and $Z(1)=0.7$. For this system, $\tau =1$ and $t=1,2,\mathrm{\dots },2000$ were used to make the plots of delay spaces.

The dynamics in the third row of Figure 5 were generated from the equations:

with initial conditions of $X(1)=X^{\prime }(1)=Z(1)=Z^{\prime }(1)=1$. For this system, $\tau =0.9$ was used for delay spaces. This continuous-time system was numerically integrated using the ode45 method in Matlab from $t=0$ to $t=13.998$ on a grid of 4,667 evenly-spaced time points for plotting dynamics, and time points $t=0.003$ through $t=7.698$ were used for visualizing delay spaces.

The dynamics in the bottom row of Figure 5 were generated from the classic Lorenz attractor equations:

with initial conditions of $X(0)=Y(0)=Z(0)=1$. A delay of $\tau =0.14$ was used to make delay spaces. This continuous-time system was numerically integrated using the ode45 method in Matlab from $t=0$ to $t=399.98$ on an evenly spaced grid of 5715 data points for visualizing delay spaces.

For each combination of ${\sigma }_m$ and ${\sigma }_p$ (the standard deviation of measurement noise and process noise, respectively), we generated 1000 time series for S₁ and S₂ as described above. For each replicate pair of time series, we used Granger causality and CCM to infer whether S₁ causes S₂ (it does) and whether S₂ causes S₁ (it does not).

We used the multivariate Granger causality Matlab package (MVGC, Barnett and Seth, 2014). We used the following settings:

• regmode = ’OLS’ (We fit the autoregressive model by the ordinary least squares method).

• icregmode = ’LWR’ (We determined the information criterion using the LWR algorithm. This is the default setting).

• morder = ’AIC’ (We used Akaike information criterion to determine the number of lags in the autoregressive model).

• momax = 50 (We used a maximum of 50 lags in the autoregressive model).

• tstat = ” (We used Granger’s F-test for statistical significance. This is the default setting).

We inferred the presence of a causal link if the p-value was less than or equal to 0.05. We inferred no causal link otherwise. When ${\sigma }_m$ and ${\sigma }_p$ were both 0, the MVGC package (correctly) exited with an error on most trials. We reported this as ‘unsuitable data’ in Figure 7D & E.

When ${\sigma }_m$ and ${\sigma }_p$ are both 0, the inferred spectral radius of the stochastic process is close to 1, and the MVGC routines can be prohibitively slow (i.e. when running 1,000 trials, the program would hang at an early stage for hours). In this case, the authors note that switching from the package’s default single-regression mode to an alternative dual-regression mode may improve runtime (Barnett and Seth, 2014). We thus switched to the dual-regression mode when the spectral radius was between 0.9999 and 1 (a spectral radius of 1 or more causes an error). This fix had no effect on benchmark results as long as at least one of ${\sigma }_m$ and ${\sigma }_p$ was not 0.

Convergent cross mapping looks for a delay map from $X$ to $Y$. That is, CCM looks for a map from $[X(t),X(t-\tau ),X(t-2\tau ),\mathrm{\dots },X(t-(E-1)\tau )]$ to $Y(t)$. Thus in order to apply CCM one needs to choose the delay τ and the vector length (dimension of the delay space) E. The parameters ${\displaystyle E}$ and τ should ideally be ‘generic’ in the sense of Takens’ theorem: we want to avoid line-crossing (such as the symbol ‘$\mathrm{\infty }$’) in the delay space, because otherwise, ${\mathrm{\Phi }}^{-1}$ in Appendix 4—figure 2 does not exist. There are different ways to do this, but no method is obviously the best (Harnack et al., 2017; Cobey and Baskerville, 2016).

Following Cobey and Baskerville, 2016 and Sugihara et al., 2012 we chose τ and $E$ to maximize univariate one-step-ahead forecast of the putative causee $X.$ That is, for $X(n)$, we try to predict $X(n+1)$ using the simplex projection method by finding delay vectors in the training data of $X$ that are most similar to $[X(n),X(n-\tau ),X(n-2\tau ),\mathrm{\dots },X(n-(E-1)\tau )]$, and take weighed average of their $X$ values one step in the future (i.e. Figure 6A where $X=Y$ and the prediction lag is 1). If the delay space has a line crossing, then at the cross-point, a one-step-ahead forecast may have more than one possible outcome and thus perform poorly. In more detail, we made one-step-ahead forecasts within the time range 201–400 (we did not use time range 1–200 to avoid transient dynamics). As per the field standard, we used leave-one-out cross-validation to do simplex projection. That is, when making a forecast for a time $t$, we used all times within 201–400 other than $t$ as training data (200 time points). We performed a grid search, varying τ from 1 to 6 and varying $E$ from 1 to 6. We then used the combination of τ and $E$ that maximized the forecast skill (the Pearson correlation between forecasts and true values) for subsequent CCM analysis. Additionally, following (Sugihara et al., 2012), if the optimal combination of τ and $E$ failed to give a significantly positive forecast skill, we did not report CCM results for that trial and reported the trial as “unsuitable data”. To test whether forecast skill is “significantly positive”, we ask whether it is robust to small changes in the training dataset. To do so, we used a naive bootstrap approach to create different versions of training libraries composed of randomly chosen delay vectors (sampling with replacement: some vectors may not be sampled and others may be sampled more than once) from the original training data using the ’random_libs’ setting in the rEDM (version 1.5) ccm method. The training library size (the number of delay vectors in the library) was chosen to be 200. We then calculated forecast skills with 300 such libraries and considered the forecast skill “significant” if at least least 95% gave a forecast skill greater than 0.

Having chosen τ and $E$, we checked three CCM criteria to infer causality (criteria 1–3 in Figure 6) using rEDM version 1.5. We did not use the fourth criterion (the prediction lag test) since its interpretation is unclear for periodic systems (Appendix 4—figure 5). For all three criteria, we used the same cross-validation setting that we used to choose τ and $E$. The first CCM criterion is that cross map skill is greater than 0. Thus, we computed cross map skill using the maximum possible number of distinct delay vectors ($200-(E-1)\tau$) and compared this value to 0.

The second CCM criterion is that the cross map skill from causee to causer with real data must be greater than the cross map skill when the putative causer is replaced with surrogate data. To test this criterion, we first computed cross map skill using the same training and testing time points as before to obtain a single cross map skill value. We then repeatedly (1000 times) computed cross map skill in the same way, but now with the putative causer time series replaced with random phase surrogate data. Random phase surrogate data were generated by Ebisuzaki’s method as implemented in the rEDM function make_surrogate_data. We then computed the p-value according to Equation 5. A putative causal link would pass this criterion if the p-value was less than or equal to 0.05.

The third CCM criterion is that cross map skill increases with more training data. Following Cobey and Baskerville, 2016, we again used a naive bootstrap approach to test for this criterion. Specifically, we computed the cross map skill with a training library composed of randomly chosen delay vectors sampled with replacement from the original training data time points. We used either a large library with $200-(E-1)\tau$ available training vectors as used previously, or a small library with 15 training vectors. For each of 1000 bootstrap trials, we compared the cross map skill from a randomly chosen small library to the cross map skill from a randomly chosen large library. We said that the cross map skill increased with training data if the cross map skill of the large library was greater than that of the small library in at least 95% of the 1000 bootstrap trials.

For ‘alternative’ CCM testing, we only changed how the third CCM criterion (cross map skill increases with more training data) were tested. Here, instead of using the bootstrap test of Cobey and Baskerville, 2016, we tested the third CCM criterion using Kendall’s τ test as suggested in Chang et al., 2017. To do this, we varied the library size from a minimum of 15 vectors to the the maximum library size ($200-(E-1)\tau$), in increments of 3 vectors. For each library size, we computed cross map skill using 50 libraries randomly sampled without replacement (e.g. the 50 libraries would be identical at the maximal library size). We then computed the median cross map skill for each library size. Finally we ran a 1-tailed Kendall’s τ test for a positive association between library size and median cross map skill. We used the function stats.kendalltau from the Python package SciPy to compute a 2-tailed p-value, and then divided this p-value by two to get a 1-tailed p-value. We said that cross map skill increased with training data if the τ statistic was positive and the 1-tailed p-value was $\le 0.05$.

The original subpopulation distributions are normal distributions with standard deviation of 10 and mean of 100 (male) or 130 (female). Each sampling plot shows 300 random samples.

To generate data for panels C-J, the system of panel A was numerically integrated using the ode45 method in Matlab with a time step of 0.005 and with the initial condition that $X$, $Y$, $Z$, $\frac{dX}{dt}$, $\frac{dY}{dt}$ were all set to one at $t=0$. Panels C, E, F, G, and I show data from a single period. For panel H the system was integrated for about five periods to more clearly visualize the lack of a continuous delay map. For panel J, the system was integrated for 1 period for the main panel and about 12 periods (to increase the sampling density) for the inset. This allows us to better see the separated legs of the curve upon zooming in. Panels C, D, F, H, and J were colored $\mathrm{mod}(t,T)$. That is, they were colored by the remainder of $t$ (time) after dividing by $T$ (here $T=2\pi$). $\tau =3.6$ was used for all delay spaces.

All systems were discretized from $t=1$ to $t=20$ on an evenly spaced grid of 200 points for visualizing delay spaces.

The dynamics in the top row were generated from the equations:

A delay time of 12 time indices (i.e. $\tau =12(20-1)/(200-1)$) was used for constructing delay spaces.

The dynamics in the second row were generated from the equations:

with $\delta =(20-1)/(200-1)$ and the initial condition $Y(1)=0.3$. A delay time of 25 time indices (i.e. $\tau =25(20-1)/(200-1)$) was used for constructing delay spaces.

In the third row the dynamics are identical to the first row, except that $X$ and $Y$ are switched, and $\tau =25$ time indices was used for constructing delay spaces.

Top row: For this system, we used the initial conditions $W(0)=\dot{W}(0)=X(0)=Y(0)=\dot{Y}(0)=1$. We numerically integrated this system using ode45 in Matlab with a time step of 0.03. We composed delay vectors of length $E=3$ with a delay of $\tau =3.6$. We visualized the delay space using data from $t=0$ through $t=29.97$ (time indices 1–1000). For CCM with temporally separate training and testing sets, we used data from $t=0$ through $t=14.97$ (time indices 1–500) for training data and data from $t=15$ through $t=29.97$ (time indices 501–1000) for testing. Specifically, in the rEDM (version 0.7.2) ccm method we set the “lib” argument to “c(1, 500)” and set the “pred” argument to “c(501, 1000)”. We used rEDM version 0.7.2 for this analysis because we found that it more easily produced distinct training and test sets than later versions (on a computer running MacOS 11.6 and R version 4.0.2). For CCM with temporally interspersed training and testing sets, we set both the lib and pred arguments to “c(1,1000)”. This setting instructs rEDM to use leave-one-out cross-validation.

Second row: Ground truth data generation and analysis were the same as in the top row, except that the roles of $X$ and $Y$ were swapped.

Third row: For this system, we used the initial conditions $W(0)=0,\dot{W}(0)=1/1.6,X(0)=1.3,Y(0)=1.5$. We numerically integrated this system using ode45 in Matlab with a time step of 0.1. We visualized the delay space using data from $t=0$ through $t=40$ (time indices 1–401). We used the delay vector parameters $(E=3,\tau =3.0)$. For CCM with temporally separate training and testing sets, we used data from $t=0$ through $t=19.9$ (time indices 1–200) for training data and data from $t=20$ through $t=39.9$ (time indices 201–400) for testing. For CCM with temporally interspersed training and testing sets, we used cross-validation over the entire range $t=0$ through $t=39.9$.

Bottom row: We discretized this system with a time step of 0.05. We visualized the delay space using data from $t=0$ through $t=20$ (time indices 1–401). We used the delay vector parameters $(E=2,\tau =2.5)$. For CCM with temporally separate training and testing sets, we used data from $t=0$ through $t=9.95$ (time indices 1–200) for training data and data from $t=10$ through $t=19.95$ (time indices 201–400) for testing. For CCM with temporally interspersed training and testing sets, we used cross-validation over the entire range $t=0$ through $t=19.95$.

For convergent cross mapping, we used the same τ and $E$ as for visualizing delay spaces (see above). “Training data size” on the horizontal axis is the number of delay vectors in the training library. Each dot in these CCM plots represents the average forecast skill over 300 randomly chosen libraries. Error bars represent the 95% confidence interval as calculated by the bias-corrected and accelerated bootstrap (1,000 bootstraps) as implemented in Matlab’s bootci function. Error bars are the same color as the dots and so are not visible when they fit inside the dots.

Top row: For this system, we used the initial conditions $X(1)=0.2,Y(1)=0.4$ and composed delay vectors of length $E=2$ with a delay of $\tau =2$. We visualized the delay space using data from time points 501–2000. We used points 801–1000 for training data and points 1001–2000 for testing cross map predictions.

Second row: For this system, we used the initial conditions $X(1)=0.4,Y(1)=0.2$ and the delay vector parameters $(E=2,\tau =1)$. We visualized the delay space using data from time points 501–2000. We used points 801–1000 for training data and points 1001–2000 for testing cross map predictions.

Third row: For this system, we used the initial conditions $X(1)=0.2,Y(1)=0.4$ and the delay vector parameters $(E=3,\tau =2)$. We visualized the delay space using data from time points 501–2000 (time points 1-$6\times {10}^5$ for the zoomed-in inset). We used points 801–1000 for training data and points 1001–2000 for testing cross map predictions.

Fourth row: For this system, we used the initial conditions $W(0)=Y(0)=0$ and $X(0)=\dot{W}(0)=\dot{Y}(0)=1$. We numerically integrated this system using ode45 in Matlab with a time step of 0.1. We visualized the delay space using data from $t=50.1$ through $t=200$ (time indices 501–2000). We used the delay vector parameters $(E=3,\tau =7.2)$. We used data from $t=70.1$ through $t=100$ (time indices 701–1000) for training data and data from $t=100.1$ through $t=200$ (time indices 1001–2000) for testing cross map predictions.

Fifth row: For this system, we used the initial conditions $X(1)=0.2,Y(1)=0$ and composed delay vectors of length $E=2$ with a delay of $\tau =2$. We visualized the delay space using data from time points 501–2000. We used points 801–1000 for training data and points 1001–2000 for testing cross map predictions.

Sixth row: For this system, the “initial” conditions specified the first three time points since we included a lag of 3. Thus, for $k=1,2,3$, $W(k)=0.2$, $X(k)=0.4$, and $Y(k)=0.3$. We composed delay vectors of length $E=3$ with a delay of $\tau =1$. We visualized the delay space using data from time points 501–2000. We used points 801–1000 for training data and points 1001–2000 for testing cross map predictions.

For convergent cross mapping (in rEDM version 0.7.2), we used the same τ and $E$ as for visualizing delay spaces. The training data size is the number of delay vectors in the training library. For the plots in the fourth column, we chose 300 random libraries of training delay vectors with variable training data size, and used the standard prediction lag of 0. Delay vectors were chosen without replacement. Note that at large training data size, some or all of the 300 random libraries can be identical. Each dot in these CCM plots represents the average forecast skill over all 300 randomly-chosen libraries. Error bars represent the 95% confidence interval as calculated by the bias-corrected and accelerated bootstrap (1,000 bootstraps) as implemented in Matlab’s bootci function. Error bars are the same color as the dots and so are not visible when they fit inside the dots.

In all rows, the cross map skill for the putative causer $Y$ was greater than for at least 95% of random phase surrogate time series (purple dot). The 5% cutoff value was computed for the maximum library size (156 for row four and ∼200 for all other rows) by running the CCM procedure after replacing the putative causer $Y$ with 500 random phase surrogate time series generated using the rEDM function make_surrogate_data.

For the plots in the fifth column we used the full library contained within the training data window (156 delay vectors for row four and ∼200 for all other rows) and varied the prediction lag. We did not use random libraries for these plots.

To generate randomized parameter sets, we randomly selected R_X, R_Y, $A_{XY}$ and $A_{YX}$ from uniform distributions. We also randomly selected the initial conditions $X(1)$ and $Y(1)$ from uniform distributions. To make systems in the ‘friendly’ parameter regime, we drew R_X and R_Y independently from the range 3.7–3.9, we drew $A_{XY}$ and $A_{YX}$ independently from the range 0.05–0.1, and we drew $X(1)$ and $Y(1)$ independently from the range 0.01–0.99. These are the same parameters used in the randomized numerical simulations of Ye et al., 2015. Next, to make systems in the ‘pathological’ parameter regime, we drew R_X from the range 3.7–3.9, we drew R_Y from the range 3.1–3.3, we drew $A_{XY}$ and $A_{YX}$ independently from the range 0.15–0.2, and we drew $X(1)$ and $Y(1)$ independently from the range 0.01–0.99. For both parameter regimes we randomly chose 250 sets of parameters and ran the system for 3,000 time points. Occasionally a randomly chosen system would leave the basin of attraction and reach large values, represented on the computer as positive or negative infinity, or ‘not a number’. When this occurred, we discarded the data and resampled parameters.

To apply CCM (in rEDM version 0.7.2) on each system, we generated a training library of delay vectors of $X$ by randomly selecting 200 vectors from among time points 100–2000. We then evaluated cross map skill from delay vectors of $X$ to values of $Y$ at points 2001–3000. Following Ye et al., 2015, we used delay vectors of length $E=2$ and a delay duration of $\tau =1$. We evaluated cross map skill with a prediction horizon of -8 through 8.