Detecting causal associations in large nonlinear time series datasets

Introduction

• Observational causal discovery is needed because experiments may be expensive, unethical, or impractical, while climate-model results can be limited by unrealistic modeled processes.The introduction notes the growing availability of observational data and climate-model output.Source: Introduction · p. 1
• Testing an observed correlation between two processes, such as ENSO and North American temperatures, requires considering other relevant variables, quickly creating a high-dimensional discovery problem.The example illustrates why pairwise association alone may not resolve competing hypotheses.Source: Introduction · p. 1
• Causal discovery must recover time-lagged links while distinguishing direct effects from autocorrelation, indirect links, and common drivers that can create spurious correlations.The goal is high detection power while controlling false positives; reconstructed causal networks are therefore sparser than correlation networks.Source: Introduction · p. 1, Figure 1 · p. 2
• Standard regression implementations of Granger causality have low detection power in high-dimensional datasets with only a few hundred samples, while bivariate analyses cannot account for indirect links or common drivers [3–10].The cited example is monthly data spanning 30 years of satellite observations.Source: Introduction · pp. 1–2
• Regularized regression and causal-discovery algorithms address some high-dimensional problems, but strong spatio-temporal interdependencies make reliable false-positive control difficult, especially with nonlinear estimators [11–21].The introduction characterizes the resulting network reconstructions as brittle and calls for more reliable methodology.Source: Introduction · p. 2
• The proposed method targets datasets with tens to hundreds of variables and few-hundred-or-more samples, combining linear or nonlinear tests to improve detection power and false-positive control without significantly increasing computational demand.The authors support these claims with analytical results and extensive numerical experiments.Source: Introduction · p. 2

Motivating example from climate science

• In the climate example, lagged correlation suggested influence in both directions, while FullCI found Nino → BCT at lag 2 with partial correlation 0.1 (p = 0.037), but not the reverse.The established causal direction is Nino → BCT, with tropical Pacific temperatures affecting British Columbia temperatures.Source: Motivating example from climate science · p. 3, Figure 2 · p. 3
• At 5% significance, FullCI detected the true Nino at lag 2 → BCT link in 53% of realizations; under higher dimensionality, detection fell to 40%, while effect size remained 0.09.Source: Motivating example from climate science · p. 4
• Because Z contains information about Nino, conditioning on it reduced the partial correlation from 0.1 to 0.09 and detection power to 53%, versus 85% when Z was independent of Nino.The passage attributes the main power loss to the reduced effect size, with increased test dimensionality contributing only partly.Source: Motivating example from climate science · p. 4
• For W pairs with the same coupling coefficient, partial correlations fell from 0.15 to 0.13 to 0.11 as autocorrelation increased, alongside lower power and more spurious correlations.Conditioning on other lags explains away information, reducing effect size and power for more autocorrelated pairs.Source: Motivating example from climate science · p. 4
• PCMCI uses PC1 to select relevant conditions, then MCI to test links; MCI accounts for autocorrelation to control false positives in interdependent time series.PC1 uses a liberal significance level to retain likely parents, while MCI adds conditions on the tested variable’s parents.Source: Motivating example from climate science · pp. 4–5
• PCMCI’s PC1 and MCI stages can use linear partial correlation or nonlinear GPDC and CMI tests, but GPDC and CMI have lower power for linear relationships with small samples.GPDC targets a broad class of nonlinear dependencies with additive noise, while CMI is nonparametric and accommodates almost any dependency type.Source: Motivating example from climate science · pp. 5–6

Experiments with linear relationships

• Across all considered network sizes, PCMCI’s average power exceeded FullCI and PC; at N = 10, 99% of links had detection power >70%.PCMCI maintained high power even when Nτmax exceeded the sample size.Source: Figure 4 · p. 6, Results · p. 8

Theoretical properties of PCMCI

• Under standard causal-discovery assumptions, PCMCI is consistent in the infinite-sample limit, and MCI quantifies the effect of a hypothetical perturbation on a target variable.The MCI statistic can therefore rank causal links by an interpretable causal-strength measure.Source: Results · p. 7

PCMCI on real climate example

• In the climate example, PCMCI detected the Ninot−2 →BCTt link with >80% power even in the high-dimensional case and estimated similar effects across three cases.The link’s partial-correlation effect size was ≈0.10 in all three cases; PCMCI also estimated causal effect strength ≈0.14 for adjacent W pairs despite differing autocorrelations.Source: Results · p. 7

Experiments with linear relationships

• Correlation produced very high false-positive rates; FullCI power fell from around 80% at N = 5 to 40% at N = 20 and was inapplicable when Nτmax > T = 150.PC did not correctly control false positives for some highly autocorrelated links, while PCMCI maintained similar power across autocorrelation strengths through N = 20. Runtime was of similar order across methods except FullCI; PCMCI was faster than Lasso for smaller networks.Source: Results · pp. 7–8

Experiments with nonlinear relationships

• With GPDC on nonlinear models, PCMCI had the highest power, while false positives were mostly controlled but slightly inflated for larger networks.With CMI, PCMCI controlled false positives and had higher power than PC except for some nonlinear link types; FullCI had almost no power, especially for nonlinear links.Source: Figure 5 · p. 8, Results · p. 9

Causal strength experiments

• For links with the same causal effect, correlation values ranged from near zero to high, whereas MCI partial correlations were larger than FullCI values and estimated causal strength.PCMCI’s higher effect size than FullCI, and often than correlation, helps explain its greater detection power than FullCI.Source: Figure 6 · p. 9, Results · pp. 9–10

Further experiments

• Supplementary experiments found robust results across larger samples, denser networks, and observational noise up to 25% of the dynamical-noise standard deviation.At noise levels on the order of the dynamical noise, performance degraded more strongly and false positives were no longer correctly controlled.Source: Results · p. 10

• Irrelevant variables can reduce power and weaken causal interpretation in high-dimensional analyses; condition selection removes them, while MCI is designed for highly interdependent time series.Source: Discussion and conclusion · p. 10
• MCI controls false-positive rates under high autocorrelation, improves power over classical techniques even with few variables, and supports meaningful causal-strength rankings.These rankings can help identify the strongest inferred links in large datasets.Source: Discussion and conclusion · p. 10
• PCMCI accommodates diverse conditional-independence tests, including for discrete or continuous series, and multivariate graph nodes for studying networks across interdependent layers [38].Source: Discussion and conclusion · p. 10
• The method focuses on time-lagged dependencies, where cause-effect direction is unambiguous; approaches that infer causality without lags may help determine direction for contemporaneous links [39,40].Source: Discussion and conclusion · p. 10
• Causal claims rely on assumptions including causal sufficiency and the Markov condition, so inferred links are relative to the included variables [18].Unobserved variables may still cause a link, making inferred links hypotheses for follow-up analyses or simulations; absence of a link relies on weaker assumptions and is more robust when the data faithfully represent the process and nonlinearities are adequately captured [37].Source: Discussion and conclusion · p. 10
• With data availability creating opportunities for discovery across disciplines, the method enables multivariate analyses that can reconstruct causal networks more credibly than mostly bivariate approaches.Source: Discussion and conclusion · p. 10

• The novel method is described in the main text, with further material available in the Supplement.Source: Materials and Methods · p. 10
• Software to reproduce the examples and numerical experiment results is available online with comprehensive method documentation.Source: Materials and Methods · pp. 10–11

Funding

• J.R. received a James S. McDonnell Foundation postdoctoral award, and the project received ERDF, German federal, and Land Brandenburg support for high-performance computing resources.The resources were provided for the high-performance computer system at the Potsdam Institute for Climate Impact Research; the authors declared no conflicts of interest.Source: Funding · p. 12

S1.1 Causal discovery method (PCMCI)

• The time-series graph represents lag-specific directed links, while contemporaneous dependencies are left undirected; omitting lagged variables can leave associations confounded by past common drivers.The graph is estimated only up to a chosen maximum lag, and the method assumes stationarity.Source: S1.1 Causal discovery method (PCMCI) · p. 14
• PCMCI first selects a superset of each variable’s parents, then uses those parents as conditions in MCI tests of lagged variable pairs.The tests can use different conditional-independence statistics; the method tests links up to a specified maximum time lag.Source: S1.1 Causal discovery method (PCMCI) · p. 14
• PCMCI’s α acts as a regularization parameter: too-low values can omit true parents and increase false positives, while too-high values enlarge conditioning sets, reduce power, and increase runtime.The implementation does not account for reduced degrees of freedom from condition selection, though the authors report no sign of inflated false positives in their numerical experiments.Source: S1.1 Causal discovery method (PCMCI) · p. 15

S1.2.1 FullCI; S1.2.2 Lasso

• For FullCI, linear Granger tests use a VAR model, while adaptive Lasso addresses high-dimensional regression through iterative feature reweighting and OLS-based p-values for selected variables.The adaptive Lasso’s regularization parameter was chosen by time-series cross-validation; some AIC-selected tests did not perform well.Source: S1.2.1 FullCI · p. 16, S1.2.2 Lasso · p. 16

S1.2.3 PC algorithm

• The PC algorithm’s aggregated p-values tended to be over-conservative, with outliers under strong autocorrelation, and false-positive rates were not reliably controlled just below a desired threshold.The implementation uses qmax = 10 and aggregates a link’s p-value as the maximum across tested condition sets.Source: S1.2.3 PC algorithm · pp. 16–17

S1.2.4 PCMCI variants

• PC1+MCI0 and its pre-whitened variant both failed to control false positives in the reported numerical experiments.PC1+MCI0 retests links removed during condition selection; the pre-whitened variant removes estimated lag-1 autocorrelation before testing.Source: S1.2.4 PCMCI variants · p. 17

S1.2.5 BivCI; S1.2.6 Unconditional pairwise measures

• The bivariate test is a pairwise lag-specific Granger-causality or transfer-entropy measure that only partly reduces autocorrelation and does not remove common drivers or indirect dependencies; unconditional measures include correlation, distance correlation, and mutual information.BivCI did not correctly control false positives in the analyzed experiments.Source: S1.2.5 BivCI · p. 17, S1.2.6 Unconditional pairwise measures · p. 17

S2.1 ParCorr

• ParCorr regresses X and Y on Z, then tests correlation between residuals; it is fast but limited to linear dependencies in the multivariate Gaussian case.Its null distribution under conditional independence is known analytically.Source: S2.1 ParCorr · p. 18

S2.2 GPDC

• For a quadratic dependence, GPDC identifies X ⊥⊥ Y | Z, whereas ParCorr does not.GPDC uses Gaussian-process regression followed by a distance-correlation test on uniformized residuals; it assumes additive functional dependence.Source: S2.2 GPDC · p. 18

S2.3 CMI

• For continuous data, CMI is estimated with a nearest-neighbor method and tested using local permutations, enabling data-adaptive estimation of conditional dependence.The implementation uses kCMI = 50, kperm = 5, and B = 500 surrogates; finite-sample convergence rates for the estimator are unknown.Source: S2.3 CMI · pp. 18–19

S2.4 Discrete data

• For discrete data, CMIsymb estimates CMI from symbol frequencies, and a permutation test is recommended because its null distribution has no analytical results.Source: S2.4 Discrete data · p. 19

S2.5 Multivariate extensions

• Although the focus is on univariate X and Y, the method can extend to multivariate time series and structured or non-Euclidean variables with suitable tests.CMI extends directly to multivariate X and Y; vector-valued regression with kernel conditional-independence tests is another option.Source: S2.5 Multivariate extensions · p. 19

S3.1 Computational complexity

• In the worst case, PCMCI’s complexity is polynomial, O(N^3τmax^2 + N^2τmax) without α-optimization, while sparser dependencies speed PC1 convergence.Runtime also depends on how the conditional-independence test scales with dimensionality and time-series length T.Source: S3.1 Computational complexity · p. 19

S3.2 Consistency

• Under causal sufficiency, faithfulness, the causal Markov condition, and consistent tests, population PCMCI recovers the true graph as sample size tends to infinity.The consistency proof also assumes no contemporaneous causal effects and stationarity; with the Markov condition, PC1 estimates the true parents.Source: S3.2 Consistency · p. 20
• Population consistency is weaker than uniform consistency: convergence can be arbitrarily slow for nearly unfaithful distributions, and uniform guarantees require further assumptions.The population-version results for MCI also do not by themselves establish properties of empirical estimators.Source: S3.2 Consistency · p. 21, S3.3 False positive control, effect… · p. 23

S3.3 False positive control, effect size, and causal strength

• Autocorrelation can inflate false positives by violating iid test assumptions, but MCI conditions on the lagged variable’s parents so tests depend on iid innovations and are well calibrated.Numerical experiments report expected false-positive rates under strong autocorrelation and nonlinear dependencies, while high-dimensional GPDC showed inflated rates; the formal result concerns population MCI and empirical estimators require caution.Source: S3.3 False positive control, effect… · pp. 22–23
• MCI has an effect size at least as large as FullCI and lower dimensionality, properties that explain its empirically higher detection power.The tests are equal when FullCI’s additional conditioning variables are independent of Y given the MCI conditions.Source: S3.3 False positive control, effect… · p. 23
• For linear additive dependencies, MCI reflects the causal coefficient and noise variances rather than dependencies from the parents, whereas correlation can differ substantially from causal effect.The stated theoretical result is for population MCI; nonlinear cases can retain dependencies because the link function mixes the innovation with the lagged variable’s parents.Source: S3.3 False positive control, effect… · p. 24

S4.1 Model setup

• The benchmark uses synthetic networks with known links to model high-dimensionality, nonlinearity, strong autocorrelation, and lagged causal dependencies.The study generated 20 network topologies for each number of variables and coupling strength, unless noted otherwise.Source: S4.1 Model setup · pp. 24–25

S4.2 Performance evaluation

• Performance is assessed on cross-links using true- and false-positive rates across 100 realizations per model, with links stratified by weak versus strong autocorrelation.The reported autocorrelation groups use a mean autocorrelation threshold of 0.7.Source: S4.2 Performance evaluation · p. 25

S5 Tigramite software package

• Tigramite 3.0 provides PCMCI and conditional-independence-test classes, plotting functions for time series, lag functions, and causal graphs, and repository documentation.The package is available at https://github.com/jakobrunge/tigramite.Source: S5 Tigramite software package · p. 25

Algorithm S1

• Algorithm S1 screens candidate parents with conditional-independence tests and is used as a PCMCI pre-selection step with pmax=Nτmax (no maximum-parent restriction) and qmax=1; standalone PC-stable uses qmax=10.Source: S6 Algorithms · p. 26

Algorithm S2

• Algorithm S2 tests MCI links using the sorted parents estimated by S1, optionally adjusts p-values with FDR, and leaves contemporaneous (τ=0) links undirected.It can restrict each variable's conditioning set to its first pX parents.Source: S6 Algorithms · p. 26

Algorithm S3

• Algorithm S3 uses cross-validated LassoCV in an adaptive-lasso procedure, with five time-series splits and kmax=5 in the numerical experiments, then fits OLS on the active set for p-values.The lasso penalty is selected by time-series-based cross-validation.Source: S6 Algorithms · p. 27

S7 Supplementary Tables

• Table S2 surveys conditional-independence tests for continuously valued data and notes that implementations standardize data and fit Gaussian processes using an RBF-plus-white-noise kernel.The sklearn bandwidth is estimated by maximizing marginal likelihood (ML-II), and DZ denotes the cardinality of Z.Source: Table S2 · p. 28
• Table S3 reports model configurations and evaluates false and true positives across 100 time-series realizations, with configurations varying links, functions, sample size, coefficient c, and random networks.The listed coupling functions include f(1)(x) = x, f(2)(x) = (1−4e−x2/2)x, and f(3)(x) = (1−4x3e−x2/2)x.Source: Table S3 · p. 29, S7 Supplementary Tables · p. 29
• Table S4 summarizes the methods compared in numerical experiments, all run with τmax = 5.Source: Table S4 · p. 30

S8 Supplementary Figures

• The supplementary numerical experiments vary linear-model variable counts, link densities, and sample lengths, and also test nonlinear models at different variable counts and sample lengths.Source: Figure S3 · p. 33, Figure S4 · p. 34, Figure S5 · p. 35
• As sample size grows, PCMCI outperforms FullCI, Lasso, and PC, while Lasso still misses some links at the largest tested sample size.For ParCorr, power increases with larger samples; runtime does not rise linearly with sample size because PCMCI and Lasso converge faster.Source: Figure S8 · p. 38
• For nonlinear GPDC experiments, power rises with sample size and FullCI better controls false positives at larger samples, but GPDC runtime scales as ∼T^3.The passage notes that efficient Gaussian-process approximation methods can speed up estimation [31].Source: Figure S10 · p. 40
• In nonlinear CMI experiments, power increases with sample length, while runtime does not increase quadratically like the individual CMI tests.Source: Figure S12 · p. 42
• With observational noise up to 25% of the dynamical-noise standard deviation, all methods show only minor effects; at comparable noise levels, performance degrades and false positives are less controlled.The passage attributes the latter pattern to common drivers no longer being detected well.Source: Figure S13 · p. 43