High-recall causal discovery for autocorrelated time series with latent confounders

Introduction

Observational causal discovery [Spirtes et al., 2000, Peters et al., 2017] from time series is a challenge of high relevance to many fields of science and engineering if experimental interventions are infeasible, expensive, or unethical. Causal knowledge of direct and indirect effects, interaction pathways, and time lags can help to understand and model physical systems and to predict the effect of interventions [Pearl, 2000]. Causal graphs can also guide interpretable variable selection for prediction and classification tasks. Causal discovery from time series faces major challenges [Runge et al., 2019a] such as unobserved confounders, high-dimensionality, and nonlinear dependencies, to name a few. Few frameworks can deal with these challenges and we here focus on constraint-based methods pioneered in the seminal works of Spirtes, Glymour, and Zhang [Spirtes et al., 2000, Zhang, 2008]. We demonstrate that existing latent causal discovery methods strongly suffer from low recall in the time series case where identifying lagged and contemporaneous causal links is the goal and autocorrelation is an added, ubiquitous challenge. Our main theoretical contributions lie in identifying low effect size as a major reason why current methods fail and in introducing a novel sound, complete, and order-independent causal discovery algorithm that yields strong gains in recall for autocorrelated continuous data. Our practical contributions lie in extensive numerical experiments that can serve as a future benchmark and in open-source Python implementations of our and major previous time series causal discovery algorithms. The paper is structured as follows: After briefly introducing the problem and existing methods in Sec. 2, we describe our method and theoretical results in Sec. 3. Section 4 provides numerical experiments followed by a discussion of strengths and weaknesses as well as an outlook in Sec. 6. The paper is accompanied by Supplementary Material (SM).

Time series causal discovery in the presence of latent confounders

We assume the reader is familiar with the Fast Causal Inference (FCI) algorithm [Spirtes et al., 1995, Spirtes et al., 2000, Zhang, 2008] and related graphical terminology, see Secs. S1 and S2 of the SM for a brief overview. Importantly, the MAGs (maximal ancestral graphs) considered in this paper can contain directed (${\rightarrow}$) and bidirected (${\leftrightarrow}$) edges (interchangeably also called links). The associated PAGs (partial ancestral graphs) may additionally have edges of the type ${\circ\!{\rightarrow}}$ and ${\circ\!{\--}\!\circ}$.

2 Existing methods

The tsFCI algorithm [Entner and Hoyer, 2010] adapts the constraint-based FCI algorithm to time series. It uses time order and stationarity to restrict conditioning sets and to apply additional edge orientations. SVAR-FCI [Malinsky and Spirtes, 2018] uses stationarity to also infer additional edge removals. There are no assumptions on the functional relationships or on the structure of confounding. Granger causality [Granger, 1969] is another common framework for inferring the causal structure of time series. It cannot deal with contemporaneous links (known as instantaneous effects in this context) and may draw wrong conclusions in the presence of latent confounders, see e.g. [Peters et al., 2017] for an overview. The ANLTSM method [Chu and Glymour, 2008] restricts contemporaneous interactions to be linear, and latent confounders to be linear and contemporaneous. TS-LiNGAM [Hyvärinen et al., 2008] is based on LiNGAM [Shimizu et al., 2006] that is rooted in the structural causal model framework [Peters et al., 2017, Spirtes and Zhang, 2016]. It allows for contemporaneous effects, assumes linear interactions with additive non-Gaussian noise, and might fail in the presence of confounding. The TiMINo [Peters et al., 2013] method restricts interactions to an identifiable function class or requires an acyclic summary graph. Yet another approach are Bayesian score-based or hybrid methods [Chickering, 2002, Tsamardinos et al., 2006]. These often become computationally infeasible in the presence of unobserved variables, see [Jabbari et al., 2017] for a discussion, or make restrictive assumptions about functional dependencies or variable types.

In this paper we follow the constraint-based approach that allows for general functional relationships (both for lagged and contemporaneous interactions), general types of variables (discrete and continuous, univariate and multivariate), and that makes no assumption on the structure of confounding. The price of this generality is that we will not be able to distinguish all members of a Markov equivalence class (although time order and stationarity allow to exclude some members of the equivalence class). Due to its additional use of stationarity we choose SVAR-FCI rather than tsFCI as a baseline and implement the method, restricted to no selection variables, in Python. As a second baseline we implement SVAR-RFCI, which is a time series adaption of RFCI along the lines of SVAR-FCI (also restricted to no selection variables). The RFCI algorithm [Colombo et al., 2012] is a modification of FCI that does not execute FCI’s potentially time consuming second edge removal phase.

3 On maximum time lag, stationarity, soundness, and completeness

In time series causal discovery the assumption of stationarity and the length of the chosen time lag window $t-\tau_{\max}\leq t^{\prime}\leq t$ play an important role. In the causally sufficient case ($\mathbf{X}=\mathbf{V})$ the causal graph stays the same for all $\tau_{\max}\geq p_{ts}$. Not so in the latent case: Let $\mathcal{M}(\mathcal{G})^{\tau_{\max}}$ be the MAG obtained by marginalizing over all unobserved variables and also all generally observed variables at times $t^{\prime}<t-\tau_{\max}$. Then, increasing the considered time lag window by increasing $\tau_{\max}$ may result in the removal of edges that are fully contained in the original window, even in the case of perfect statistical decisions. In other words, $\mathcal{M}(\mathcal{G})^{\tau_{\max,1}}$ with $\tau_{\max,1}<\tau_{\max,2}$ need not be a subgraph of $\mathcal{M}(\mathcal{G})^{\tau_{\max,2}}$. Hence, $\tau_{\max}$ may be regarded more as an analysis choice than as a tunable parameter. For the same reason stationarity also affects the definition of MAGs and PAGs that are being estimated. For example, SVAR-FCI uses stationarity to also remove edges whose separating set extends beyond the chosen time lag window. It does, therefore, in general not determine a PAG of $\mathcal{M}(\mathcal{G})^{\tau_{\max}}$. To formalize this let $i)$ $\mathcal{M}(\mathcal{G})_{statA}^{\tau_{\max}}$ be the MAG obtained from $\mathcal{M}(\mathcal{G})^{\tau_{\max}}$ by enforcing repeating adjacencies, let $ii)$ $\mathcal{P}(\mathcal{G})_{statA}^{\tau_{\max}}$ be the maximally informative PAG for the Markov equivalence class of $\mathcal{M}(\mathcal{G})_{statA}^{\tau_{\max}}$, which can be obtained from running the FCI orientation rules on $\mathcal{M}(\mathcal{G})_{statA}^{\tau_{\max}}$, and let $iii)$ $\mathcal{P}(\mathcal{G})_{statAO}^{\tau_{\max}}$ be the PAG obtained when additionally enforcing time order and repeating orientations at each step of applying the orientation rules. Note that $\mathcal{P}(\mathcal{G})_{statAO}^{\tau_{\max}}$ may have fewer circle marks, i.e., may be more informative than $\mathcal{P}(\mathcal{G})_{statA}^{\tau_{\max}}$. Our aim is to estimate $\mathcal{P}(\mathcal{G})_{statAO}^{\tau_{\max}}$. We say an algorithm is sound if it returns a PAG for $\mathcal{M}(\mathcal{G})_{statA}^{\tau_{\max}}$, and complete if it returns $\mathcal{P}(\mathcal{G})_{statAO}^{\tau_{\max}}$. Below we write $\mathcal{M}(\mathcal{G})=\mathcal{M}(\mathcal{G})_{statA}^{\tau_{\max}}$ and $\mathcal{P}(\mathcal{G})=\mathcal{P}(\mathcal{G})_{statAO}^{\tau_{\max}}$ for simplicity.

4 Motivational example

We illustrate the challenge posed by unobserved variables with the example of Fig. 1. SVAR-FCI with the partial correlation (ParCorr) CI test correctly identifies the auto-links but misses the true lagged link $Y_{t-1}{\rightarrow}Z_{t}$ and returns a false link $Y_{t-2}{\rightarrow}Z_{t}$ instead. In most realizations the algorithm fails to detect the contemporaneous adjacency $X_{t}{\leftrightarrow}Y_{t}$ and, if detected, fails to orient it as bidirected. The reason are wrong CI tests in its edge removal and orientation phases. When it iterates through conditioning sets of cardinality $p=0$ in the edge removal phase, the correlation $\rho(X_{t};Y_{t})$ is non-significant in many realizations since the high autocorrelation of both $X$ and $Y$ increases their variance and decreases their signal-to-noise ratio (the common signal due to the latent confounder). Further, for $p=1$ also the lagged correlation $\rho(Y_{t-1};Z_{t}|Y_{t-2})$ often is non-significant and the true link $Y_{t-1}{\rightarrow}Z_{t}$ gets removed. Here conditioning away the autocorrelation of $Y_{t-1}$ decreases the signal while the noise level in $Z_{t}$ is still high due to $Z$’s autocorrelation. This false negative has implications for further CI tests since $Y_{t-1}$ won’t be used in subsequent conditioning sets: The path $Y_{t-2}{\rightarrow}Y_{t-1}{\rightarrow}Z_{t}$ can then not be blocked anymore and the false positive $Y_{t-2}{\rightarrow}Z_{t}$ remains even after the next removal phase. In the orientation phase of SVAR-FCI rule $\mathcal{R}1$ yields tails for all auto-links. Even if the link $X_{t}{\circ\!{\--}\!\circ}Y_{t}$ is detected, it is in most cases not oriented correctly. The reason again lies in wrong CI tests: In principle the collider rule $\mathcal{R}0$ should identify $X_{t}{\leftrightarrow}Y_{t}$ since the middle node of the triple $X_{t-1}{\circ\!{\rightarrow}}X_{t}{\circ\!{\--}\!\circ}Y_{t}$ does not lie in the separating set of $X_{t-1}$ and $Y_{t}$ (and similarly for $X$ and $Y$ swapped). In practice $\mathcal{R}0$ is implemented with the majority rule [Colombo and Maathuis, 2014] to avoid order-dependence, which involves further CI test given subsets of the adjacencies of $X_{t-1}$ and $Y_{t}$. SVAR-FCI here finds independence given $Y_{t-1}$ (correct) but also given $X_{t}$ (wrong, due to autocorrelation). Since the middle node $X_{t}$ is in exactly half of the separating sets, the triple is marked as ambiguous and left unoriented. The same applies when $X$ and $Y$ are swapped.

Autocorrelation is only one manifestation of a more general problem we observe here: Low signal-to-noise ratio due to an ‘unfortunate’ choice of conditioning sets that leads to low effect size (here partial correlation) and, hence, low statistical power of CI tests. Wrong CI tests then lead to missing links, and these in turn to false positives and wrong orientations. In the following we analyze effect size more theoretically and suggest a general idea to overcome this issue.

Latent PCMCI

The detection power of a true link $X^{i}_{t-\tau}{\ast\!{\rightarrow}}X^{j}_{t}$, where below we write $A=X^{i}_{t-\tau}$ and $B=X^{j}_{t}$ to emphasize that the discussion also applies to the non-time series case, quantifies the probability of the link not being erroneously removed due to a wrong CI test. It depends on $i)$ the sample size (usually fixed), $ii)$ the CI tests’ significance level $\alpha$ (fixed by the researcher as the desired false positives level), $iii)$ the CI tests’ estimation dimensions (kept at a minimum by SVAR-FCI’s design to preferentially test small conditioning sets), and $iv)$ the effect size. We here define effect size as the minimum of the CI test statistic values $I(A;B|\mathcal{S})$ taken over all conditioning sets $\mathcal{S}$ that are being tested (for fixed $A$ and $B$). As observed in the motivating example, this minimum can become very small and hence lead to low detection power. The central idea of our proposed method Latent PCMCI (LPCMCI) is to increase effect size by $a)$ restricting the conditioning sets $\mathcal{S}$ that need to be tested in order to remove all wrong links, and by $b)$ extending those sets $\mathcal{S}$ that do need to be tested with so called default conditions $\mathcal{S}_{def}$ that increase the CI test statistic values and at the same time do not induce spurious dependencies. Regarding $a)$, Lemma S5 proves that it is sufficient to only consider conditioning sets that consist of ancestors of $A$ or $B$ only. Regarding $b)$, and well-fitting with $a)$, Lemma S4 proves that no spurious dependencies are introduced if $\mathcal{S}_{def}$ consist of ancestors of $A$ or $B$ only. Further, the following theorem shows that taking $\mathcal{S}_{def}$ as the union of the parents of $A$ and $B$ (without $A$ and $B$ themselves) improves the effect size of LPCMCI over that of SVAR-FCI. This generalizes the momentary conditional independence (MCI) idea that underlies the PCMCI and PCMCI+ algorithms [Runge et al., 2019b, Runge, 2020] to causal discovery with latent confounders. We state the theorem in an information theoretic framework, where $I$ denotes (conditional) mutual information and $\mathcal{I}(A;B;C|D)\equiv I(A;B|D)-I(A;B|C\cup D)$ the interaction information.

Let $A{\ast\!{\rightarrow}}B$ (with $A=X^{i}_{t-\tau}$ and $B=X^{j}_{t}$) be a link (${\rightarrow}$ or ${\leftrightarrow}$) in $\mathcal{M}(\mathcal{G})$. Consider the default conditions $\mathcal{S}_{def}=pa(\{A,B\},\mathcal{M}(\mathcal{G}))\setminus\{A,B\}$ and denote $\mathbf{X}^{*}=\mathbf{X}\setminus\mathcal{S}_{def}$. Let $\mathbf{S}=\arg\min_{\mathcal{S}\subseteq\mathbf{X}^{*}\setminus\{A,B\}}I(A;B|\mathcal{S}\cup\mathcal{S}_{def})$ be the set of sets that define LPCMCI’s effect size. If $i)$ there is $\mathcal{S}^{*}\in\mathbf{S}$ with $\mathcal{S}^{*}\subseteq adj(A,\mathcal{M}(\mathcal{G}))\setminus\mathcal{S}_{def}$ or $\mathcal{S}^{*}\subseteq adj(B,\mathcal{M}(\mathcal{G}))\setminus\mathcal{S}_{def}$ and $ii)$ there is a proper subset $\mathcal{Q}\subset\mathcal{S}_{def}$ such that $\mathcal{I}(A;B;\mathcal{S}_{def}\setminus\mathcal{Q}|\mathcal{S}^{*}\cup\mathcal{Q})<0$, then

If the assumptions are not fulfilled, then (trivially) "$\geq$" holds in eq. (2).

The second assumption only requires that any subset $\mathcal{S}_{def}\setminus{Q}$ of the parents contains information that increases the information between $A$ and $B$. A sufficient condition for this is detailed in Corollary S1.

These considerations lead to two design principles behind LPCMCI: First, when testing for conditional independence of $A$ and $B$, discard conditioning sets that contain known non-ancestors of $A$ and $B$. Second, use known parents of $A$ and $B$ as default conditions. Unless the higher effect size is overly counteracted by the increased estimation dimension (due to conditioning sets of higher cardinality), this leads to higher detection power and hence higher recall of true links. While we do not claim that our choice of default conditions as further detailed in Sec. 3.4 is optimal, our numerical experiments in Sec. 4 and the SM indicate strong increases in recall for the case of continuous variables with autocorrelation. In [Runge et al., 2019b, Runge, 2020] it is discussed that, in addition to higher effect size, conditioning on the parents of both $A$ and $B$ also leads to better calibrated tests which in turn avoids inflated false positives. Another benefit is that fewer conditioning sets need to be tested, which is also the motivation for a default conditioning on known parents in [Lee and Honavar, 2020].

The above design principles are only useful if some (non-)ancestorships are known before all CI test have been completed. LPCMCI achieves this by entangling the edge removal and edge orientation phases, i.e., by learning ancestral relations before having removed all wrong links. For this purpose we below develop novel orientation rules. These are not necessary in the causally sufficienct setting considered by PCMCI+ [Runge, 2020] because there the default conditions need not be limited to ancestors of $A$ or $B$ (although PCMCI+ tries to keep the number of default conditions low). While not considered here, background knowledge about (non-)ancestorships can easily be incorporated.

2 Introducing middle marks and LPCMCI-PAGs

To facilitate early orientation of edges we give an unambiguous causal interpretation to the graph at every step of the algorithm. This is achieved by augmenting edges with middle marks. Using generic variable names $A$, $B$, and $C$ indicates that the discussion also applies to the non-time series case.

Middle marks are denoted above the link symbol and can be ‘?’, ‘L’, ‘R’, ‘!’, or ‘’ (empty). The ‘L’ (‘R’) on $A{\mathbin{\ensurestackMath{\stackon[-1pt]{{\ast\!{\--}\!\ast}}{{\scriptscriptstyle L}}}}}B$ ($A{\mathbin{\ensurestackMath{\stackon[-1pt]{{\ast\!{\--}\!\ast}}{{\scriptscriptstyle R}}}}}B$) asserts that if $A<B$ ($B<A$) then $B\notin an(A,\mathcal{G})$ or there is no $\mathcal{S}\subseteq pa(A,\mathcal{M}(\mathcal{G}))$ that m-separates $A$ and $B$ in $\mathcal{M}(\mathcal{G})$. Here $<$ is any total order on the set of variables. Its choice is arbitrary and does not influence the causal information content, the sole purpose being to disambiguate $A{\mathbin{\ensurestackMath{\stackon[-1pt]{{\ast\!{\--}\!\ast}}{{\scriptscriptstyle L}}}}}B$ from $A{\mathbin{\ensurestackMath{\stackon[-1pt]{{\ast\!{\--}\!\ast}}{{\scriptscriptstyle R}}}}}B$. Moreover, ‘$\ast$’ is a wildcard that may stand for all three edge marks (tail, head, circle) that appear in PAGs. Further, the ‘!’ on $A{\mathbin{\ensurestackMath{\stackon[-2pt]{{\ast\!{\--}\!\ast}}{{\scriptscriptstyle!}}}}}B$ asserts that both $A{\mathbin{\ensurestackMath{\stackon[-1pt]{{\ast\!{\--}\!\ast}}{{\scriptscriptstyle L}}}}}B$ and $A{\mathbin{\ensurestackMath{\stackon[-1pt]{{\ast\!{\--}\!\ast}}{{\scriptscriptstyle R}}}}}B$ are true, and the empty middle mark on $A{\ast\!{\--}\!\ast}B$ says that $A\in adj(B,\mathcal{M}(\mathcal{G}))$. Lastly, the ‘?’ on $A{\mathbin{\ensurestackMath{\stackon[-2pt]{{\ast\!{\--}\!\ast}}{{\scriptscriptstyle?}}}}}B$ doesn’t promise anything. Non-circle edge marks (here potentially hidden by the ‘$\ast$’ symbol) still convey their standard meaning of ancestorship and non-ancestorship, and the absence of an edge between $A$ and $B$ still asserts that $A\notin adj(B,\mathcal{M}(\mathcal{G}))$. We call a PAG $\mathcal{C}(\mathcal{G})$ whose edges are extended with middle marks a LPCMCI-PAG for $\mathcal{M}(\mathcal{G})$, see Sec. S3 in the SM for a more formal definition. The ‘$\ast$’ symbol is also used as a wildcard for the five middle marks.

Note that we are not changing the quantity we are trying to estimate, this is still the PAG $\mathcal{P}(\mathcal{G})$ as explained in Sec. 2.3. The notion of LPCMCI-PAGs is used in intermediate steps of LPCMCI and has two advantages. First, $A{\ast\!{\--}\!\ast}B$ is reserved for $A\in adj(B,\mathcal{M}(\mathcal{G}))$ and thus has an unambiguous meaning at every point of the algorithm, unlike for (SVAR-)FCI and (SVAR-)RFCI. In fact, even if LPCMCI is interrupted at any arbitrary point it still yields a graph with unambiguous and sound causal interpretation. Second, middle marks carry fine-grained causal information that allows to determine definite adjacencies early on:

In LPCMCI-PAG $\mathcal{C}(\mathcal{G})$ one may replace 1.) $A{\mathbin{\ensurestackMath{\stackon[-0pt]{{\rightarrow}}{{\scriptscriptstyle!}}}}}B$ by $A{\rightarrow}B$, 2.) $A{\mathbin{\ensurestackMath{\stackon[-0pt]{{\rightarrow}}{{\scriptscriptstyle L}}}}}B$ for $A>B$ by $A{\rightarrow}B$, and 3.) $A{\mathbin{\ensurestackMath{\stackon[-0pt]{{\rightarrow}}{{\scriptscriptstyle R}}}}}B$ for $A<B$ by $A{\rightarrow}B$.

When LPCMCI has converged all middle marks are empty and hence $\mathcal{C}(\mathcal{G})$ is a PAG. We choose a total order consistent with time order, namely $X^{i}_{t-\tau}<X^{j}_{t}$ iff $\tau>0$ or $\tau=0$ and $i<j$. Lagged links can then be initialized with edges ${\mathbin{\ensurestackMath{\stackon[-1pt]{{\circ\!{\rightarrow}}}{{\scriptscriptstyle L}}}}}$ (contemporaneous links as ${\mathbin{\ensurestackMath{\stackon[-2pt]{{\circ\!{\--}\!\circ}}{{\scriptscriptstyle?}}}}}$).

3 Orientations rules for LPCMCI-PAGs

We now discuss rules for edge orientation in LPCMCI-PAGs. For this we need a definition:

In MAG $\mathcal{M}(\mathcal{G})$ let $A$ and $B$ be m-separated by $\mathcal{S}$. The set $\mathcal{S}$ is a weakly minimal separating set of $A$ and $B$ if $i)$ it decomposes as $\mathcal{S}=\mathcal{S}_{1}\dot{\cup}\,\mathcal{S}_{2}$ with $\mathcal{S}_{1}\subseteq an(\{A,B\},\mathcal{M}(\mathcal{G}))$ such that $ii)$ if $\mathcal{S}^{\prime}=\mathcal{S}_{1}\dot{\cup}\,\mathcal{S}_{2}^{\prime}$ with $\mathcal{S}_{2}^{\prime}\subseteq\mathcal{S}_{2}$ m-separates $A$ and $B$ then $\mathcal{S}_{2}^{\prime}=\mathcal{S}_{2}$. The pair $(\mathcal{S}_{1},\mathcal{S}_{2})$ is called a weakly minimal decomposition of $\mathcal{S}$.

This generalizes the notion of minimal separating sets, for which additionally $\mathcal{S}_{1}=\emptyset$. Since LPCMCI is designed to extend conditioning sets by known ancestors, the separating sets it finds are in general not minimal. However, they are still weakly minimal. The following Lemma, a generalization of the unshielded triple rule [Colombo et al., 2012], is central to orientations in LPCMCI-PAGs:

Let $A{\mathbin{\ensurestackMath{\stackon[-2pt]{{\ast\!{\--}\!\ast}}{{\scriptstyle\ast}}}}}B{\mathbin{\ensurestackMath{\stackon[-2pt]{{\ast\!{\--}\!\ast}}{{\scriptstyle\ast}}}}}C$ be an unshielded triple in LPCMCI-PAG $\mathcal{C}(\mathcal{G})$ and $\mathcal{S}_{AC}$ the separating set of $A$ and $C$. 1.) If $i)$ $B\in\mathcal{S}_{AC}$ and $ii)$ $\mathcal{S}_{AC}$ is weakly minimal, then $B\in an(\{A,C\},\mathcal{G})$. 2.) Let $\mathcal{T}_{AB}\subseteq an(\{A,B\},\mathcal{M}(\mathcal{G}))$ and $\mathcal{T}_{CB}\subseteq an(\{C,B\},\mathcal{M}(\mathcal{G}))$ be arbitrary. If $i)$ $B\notin\mathcal{S}_{AC}$, $ii)$ $A$ and $B$ are not m-separated by $\mathcal{S}_{AC}\cup\mathcal{T}_{AB}\setminus\{A,B\}$, $iii)$ $C$ and $B$ are not m-separated by $\mathcal{S}_{AC}\cup\mathcal{T}_{CB}\setminus\{C,B\}$, then $B\notin an(\{A,C\},\mathcal{G})$. The conditioning sets in $ii)$ and $iii)$ may be intersected with the past and present of the later variable.

Part 2.) of this Lemma generalizes the FCI collider rule $\mathcal{R}0$ to rule $\mathcal{R}0^{\prime}$ (of which there are several variations when restricting to particular middle marks), and part 1.) generalizes $\mathcal{R}1$ to $\mathcal{R}1^{\prime}$. Rules $\mathcal{R}2$ and $\mathcal{R}8$ generalize trivially to triangles in $\mathcal{C}(\mathcal{G})$ with arbitrary middle marks, giving rise to $\mathcal{R}2^{\prime}$ and $\mathcal{R}8^{\prime}$. Rules $\mathcal{R}3$, $\mathcal{R}9$ and $\mathcal{R}10$ are generalized to $\mathcal{R}3^{\prime}$, $\mathcal{R}9^{\prime}$ and $\mathcal{R}10^{\prime}$ by adding the requirement that the middle variables of certain unshielded colliders are in the separating set of the two outer variables, and that these separating sets are weakly minimal. Since there are no selection variables, rules $\mathcal{R}5$, $\mathcal{R}6$ and $\mathcal{R}7$ are not applicable. Rule $\mathcal{R}4^{\prime}$ generalizes the discriminating path rule [Colombo et al., 2012] of RFCI. These rules are complemented by the replacements specified in Lemma 1 and a rule for updating middle marks. Precise formulations of all rules are given in Sec. S4 of the SM.

We stress that these rules are applicable at every point of the algorithm and that they may be executed in any order. This is different from the (SVAR-)FCI orientation phase which requires that prior to orientation a PAG has been found. Also (SVAR-)RFCI orients links only once an RFCI-PAG has been determined, and both (SVAR-)FCI and (SVAR-)RFCI require that all colliders are oriented before applying their other orientation rules.

The relevance of the novel orientation rules bears on them allowing to determine ancestorships and non-ancestorships already after only few CI tests have been performed. This is utilized in LPCMCI by entangling the edge removal and edge orientation phases, which then allows to implement the idea of the PCMCI and PCMCI+ algorithms [Runge et al., 2019b, Runge, 2020] to increase the effect sizes of CI tests (and hence the recall of the algorithm, see Theorem 1 and the subsequent discussion) by conditioning on known parents also in the causally insufficient case considered here (where latent confounders are allowed). The aspect of determining ancestorships with only few CI tests is similar in spirit to an approach taken in the recent work [Mastakouri et al., 2020], which considers the narrower but important task of causal feature selection in time series with latent confounders: The SyPI algorithm introduced there does not aim at finding the full PAG $\mathcal{P}(\mathcal{G})$ but rather at finding ancestors of given target variables. Under several assumptions on the connectivity pattern of the time series graph $\mathcal{G}$ the work [Mastakouri et al., 2020] presents conditions that are sufficient for a given variable (the potential cause) to be an ancestor of another given variable (the potential effect). For certain types of ancestors, namely for parents that belong to a time series which in the summary graph is not confounded with the time series of the target variable, these conditions are even necessary. These findings allow SyPI to determine (some) ancestorships with only two CI tests per pair of potential cause and potential effect. It would be interesting to investigate whether the problem of causal feature selection as framed in [Mastakouri et al., 2020] can benefit from some of the ideas presented here, for example from the novel orientation rules (which do not require restrictions on the connectivity pattern of $\mathcal{G}$) or the idea to increase the effect sizes of CI tests by conditioning on known parents. Similarly, it would be interesting to explore whether the ideas behind SyPI can be utilized to further improve the statistical performance of algorithms that approach the more general task of finding the full PAG.

4 The LPCMCI algorithm

LPCMCI is a constraint-based causal discovery algorithm that utilizes the findings of Sec. 3.1 to increase the effect size of CI tests. High-level pseudocode is given in Algorithm 1. After initializing $\mathcal{C}(\mathcal{G})$ as a complete graph, the algorithm enters its preliminary phase in lines 2 to 4. This involves calls to Algorithm S2 (pseudocode in Sec. S5 of the SM), which removes many (but in general not all) false links and, while doing so, repeatedly applies the orientation rules introduced in the previous section. These rules identify a subset of the (non-)ancestorships in $\mathcal{G}$ and accordingly mark them by heads or tails on edges in $\mathcal{C}(\mathcal{G})$. This information is then used as prescribed by the two design principles of LPCMCI that were explained in Sec. 3.1: The non-ancestorships further constrain the conditioning sets $\mathcal{S}$ of subsequent CI tests, the ancestorships are used to extend these sets to $\mathcal{S}\cup\mathcal{S}_{def}$ where $\mathcal{S}_{def}=pa(\{X^{i}_{t-\tau},X^{j}_{t}\},\mathcal{C}(\mathcal{G}))$ are the by then known parents of those variables whose independence is being tested. All parentships marked in $\mathcal{C}(\mathcal{G})$ after line 3 are remembered and carried over to an elsewise re-initialized $\mathcal{C}(\mathcal{G})$ before the next application of Alg. S2. Conditioning sets can then be extended with known parents already from the beginning. The purpose of this iterative process is to determine an accurate subset of the parentships in $\mathcal{G}$. These are then passed on to the final phase in lines 5 - 6, which starts with one final application of Alg. S2. At this point there may still be false links because Alg. S2 may fail to remove a false link between variables $X^{i}_{t-\tau}$ and $X^{j}_{t}$ if neither of the two is an ancestor of the other. This is the purpose of Algorithm S3 (pseudocode in Sec. S5 of the SM) that is called in line 6, which thus plays a similar role as the second removal phase in (SVAR-)FCI. Algorithm S3 repeatedly applies orientation rules and uses identified (non-)ancestorships in the same way as Alg. S2. As stated in the following theorems, LPCMCI will then have found the PAG $\mathcal{P}(\mathcal{G})$. Moreover, its output does not depend on the order of the $N$ time series variables $X^{j}$. The number $k$ of iterations in the preliminary phase is a hyperparameter and we write LPCMCI($k=k_{0}$) when specifying $k=k_{0}$. Stationarity is enforced at every step of the algorithm, i.e., whenever an edge is removed or oriented all equivalent time shifted edges (called ‘homologous’ in [Entner and Hoyer, 2010]) are removed too and oriented in the same way.

Assume that there is a process as in eq. (1) without causal cycles, which generates a distribution $P$ that is faithful to its time series graph $\mathcal{G}$. Further assume that there are no selection variables, and that we are given perfect statistical decisions about CI of observed variables in $P$. Then LPCMCI is sound and complete, i.e., it returns the PAG $\mathcal{P}(\mathcal{G})$.

The output of LPCMCI does not depend on the order of the $N$ time series variables $X^{j}$ (the $j$-indices may be permuted).

5 Back to the motivational example in Fig. 1

The first iteration ($l=0$) of LPCMCI also misses the links $Y_{t-1}{\rightarrow}Z_{t}$ and finds $X_{t}{\ast\!{\--}\!\ast}Y_{t}$ in only few realizations (we here suppress middle marks for simpler notation), but orientations are already improved as compared to SVAR-FCI. Rule $\mathcal{R}1^{\prime}$ applied after $p=1$ orients the auto-links $X_{t-1}{\rightarrow}X_{t}$ and $Y_{t-1}{\rightarrow}Y_{t}$. This leads to the parents sets $pa(X_{t},\mathcal{C}(\mathcal{G}))=\{X_{t-1}\}$ and $pa(Y_{t},\mathcal{C}(\mathcal{G}))=\{Y_{t-1}\}$, which are then used as default conditions in subsequent CI tests. This is relevant for orientation rule $\mathcal{R}0^{\prime}$ that tests whether the middle node of the unshielded triple $X_{t-1}{\circ\!{\rightarrow}}X_{t}{\circ\!{\--}\!\circ}Y_{t}$ does not lie in the separating set of $X_{t-1}$ and $Y_{t}$. Due to the extra conditions the relevant partial correlation $\rho(X_{t-1};Y_{t}|X_{t},X_{t-2},Y_{t-1})$ now correctly turns out significant. This identifies $X_{t}$ as collider and (since the same applies with $X$ and $Y$ swapped) the bidirected edge $X_{t}{\leftrightarrow}Y_{t}$ is correctly found. The next iteration ($l=1$) then uses the parents obtained in the $l=0$ iteration, here the autodependencies plus the (false) link $Y_{t-2}{\rightarrow}Z_{t}$, as default conditions already from the beginning for $p=0$. While the correlation $\rho(X_{t};Y_{t})$ used by SVAR-FCI is often non-significant, the partial correlation $\rho(X_{t};Y_{t}|X_{t-1},Y_{t-1})$ is significant since the autocorrelation noise was removed and effect size increased (indicated as link color in Fig. 1) in accord with Theorem 1. Also the lagged link is correctly detected because $\rho(Y_{t-1};Z_{t}|Y_{t-2},Z_{t-1})$ is larger than $\rho(Y_{t-1};Z_{t}|Y_{t-2})$. The false link $Y_{t-2}{\rightarrow}Z_{t}$ is now removed since the separating node $Y_{t-1}$ was retained. This wrong parentship is then also not used for default conditioning anymore. Orientations of bidirected links are facilitated as before and $Y_{t-1}{\rightarrow}Z_{t}$ is oriented by rule $\mathcal{R}1^{\prime}$.

Numerical experiments

We here compare LPCMCI to the SVAR-FCI and SVAR-RFCI baselines with CI tests based on linear partial correlation (ParCorr), for an overview of further experiments presented in the SM see the end of this section. To limit runtime we constrain the cardinality of conditioning sets to $3$ in the second removal phase of SVAR-FCI and in Alg. S3 of LPCMCI (excluding the default conditions $\mathcal{S}_{def}$, i.e., $|\mathcal{S}|\leq 3$ but $|\mathcal{S}\cup\mathcal{S}_{def}|>3$ is allowed). We generate datasets with this variant of the SCM in eq. (1):

In Fig. 2A we show LPCMCI for $k=0,\ldots,4$ against increasing autocorrelation $a$. Note that $a=0$ implies a different true PAG than $a>0$. The largest gain, both in recall and precision, comes already from $k=0$ to $k=1$. For higher $k$ LPCMCI maintains false positive control and orientation precision, and improves recall before converging at $k=4$. The gain in recall is largely attributable to improved effect size. On the downside, larger $k$ increase cardinality (estimation dimension) and runtime. However, the runtime increase is only marginal because later $l$-steps converge faster and the implementation caches CI test results. Fig. 2B shows a comparison of LPCMCI with SVAR-FCI and SVAR-RFCI against autocorrelation, depicting LPCMCI for $k=0$ and $k=4$. Already LPCMCI($k=0$) has higher adjacency and orientation recall than SVAR-FCI and SVAR-RFCI for increasing autocorrelation while they are on par for $a=0$. This comes at the price of precision, especially lagged orientation precision. LPCMCI($k=4$) has more than 0.4 higher contemporaneous orientation recall and still 0.1 higher lagged orientation recall than SVAR-FCI and SVAR-RFCI. Lagged precision is higher for high autocorrelation and contemporaneous precision is slightly lower. LPCMCI($k=4$) maintains high recall for increasing autocorrelation $a\geq 0.5$ while SVAR-FCI and SVAR-RFCI’s recall sharply drops. These results can be explained by improved effect size while the increased cardinality ($\approx 5$) of separating sets is still moderate compared to the sample size $T=500$. LPCMCI($k=0$) has similar low runtime as SVAR-RFCI, for LPCMCI($k=4$) it is comparable to that of SVAR-FCI. In Fig. 2C we show results for different numbers of variables $N$. As expected, all methods have decreasing adjacency and orientation recall for higher $N$, but LPCMCI starts at a much higher level. For $N=3$ both SVAR-FCI and SVAR-RFCI cannot control false positives for lagged links while for larger $N$ false positives become controlled. The reason is the interplay of ill-calibrated CI tests for smaller $N$ due to autocorrelation (inflating false positives) with sequential testing for larger $N$ (reducing false positives), as has been discussed in [Runge et al., 2019b, Runge, 2020] for the similar PC algorithm [Spirtes and Glymour, 1991]. LPCMCI better controls false positives here, its decreasing recall can be explained by decreasing effect size and increasing cardinality. Runtime becomes slightly larger than that of SVAR-FCI for larger $N$. Fig. 2D shows results for different maximum time lags $\tau_{\max}$. Note that these imply different true PAGs, especially since further lagged links appear for larger $\tau_{\max}$. All methods show a decrease in lagged recall and precision, whereas contemporaneous recall and precision stay almost constant. For SVAR-FCI there is an explosion of runtime for higher $\tau_{\max}$ due to excessive searches of separating sets in its second removal phase. In LPCMCI this is partially overcome since the sets that need to be searched through are more restricted.

In Sec. S9 in the SM we present further numerical experiments. This includes more combinations of model parameters $N,\,a,\,\lambda,\,T$, nonlinear models together with the nonparametric GPDC CI test [Runge et al., 2019b], and a comparison to a residualization approach. In these cases the results are largely comparable to those above regarding relative performances. For non-time series models we find that, although all findings of Secs. 3.1 through 3.4 still apply, LPCMCI($k$) is on par with the baselines for $k=0$ while it shows inflated false positives for $k=4$. Similarly, for models of discrete variables together with a $G$-test of conditional independence LPCMCI($k$) performs comparable to the baselines for $k=0$ and gets worse with increasing $k$. A more detailed analysis of LPCMCI’s performance in these two cases, non-time series and discrete models, is subject to future research.

Application to real data

We here discuss an application of LPCMCI to average daily discharges of rivers in the upper Danube basin, measurements of which are made available by the Bavarian Environmental Agency at https://www.gkd.bayern.de. We consider measurements from the Iller at Kempten ($X$), the Danube at Dillingen ($Y$), and the Isar at Lenggries ($Z$). While the Iller discharges into the Danube upstream of Dillingen with the water from Kempten reaching Dillingen within about a day, the Isar reaches the Danube downstream of Dillingen. We thus expect a contemporaneous link $X_{t}{\rightarrow}Y_{t}$ and no direct causal relationships between the pairs $X,Z$ and $Y,Z$. Since all variables may be confounded by rainfall or other weather conditions, this choice allows to test the ability of detecting and distinguishing directed and bidirected links. To keep the sample size comparable with those in the simulation studies we restrict to the records of the past three years (2017-2019). We set $\tau_{\max}=2$ and apply LPCMCI($k$) for $k=0,\ldots,4$ and $\alpha=0.01$ with ParCorr CI tests. Restricting the discussion to contemporaneous links, LPCMCI correctly finds $X_{t}{\rightarrow}Y_{t}$ for $k=1,\ldots,4$ and for $k=0$ wrongly finds $X_{t}{\leftrightarrow}Y_{t}$. For all $k$ it infers the bidirected link $X_{t}{\leftrightarrow}Z_{t}$, which is plausible due to confounding by weather. For $k=3,4$ LPCMCI wrongly finds the directed link $Z_{t}{\rightarrow}Y_{t}$, which should either be absent or bidirected. The results are similar for $\alpha=0.05$, with the difference that LPCMCI then always correctly finds $X_{t}{\rightarrow}Y_{t}$ but wrongly infers $Z_{t}{\rightarrow}Y_{t}$ also for $k=1,2$. In comparison, SVAR-FCI with ParCorr CI tests finds the contemporaneous adjacencies $Y_{t}{\circ\!{\--}\!\circ}X_{t}{\circ\!{\--}\!\circ}Z_{t}$ for $\alpha=0.01,0.03,0.05,0.08,0.1,0.3,0.5$ and $Y_{t}{{\leftarrow}\!\circ}X_{t}{\circ\!{\--}\!\circ}Z_{t}$ for $\alpha=0.8$. The estimated PAGs are shown in Sec. S10 of the SM.

We note that since the discharge values show extreme events caused by heavy rainfall, the assumption of stationarity is expected to be violated. For other analyses of the dataset of average daily discharges see [Asadi et al., 2015, Engelke and Hitz, 2020, Mhalla et al., 2020, Gnecco et al., 2020]. More detailed applications to and analyses of LPCMCI on real data are subject to future research.

Discussion and future work

Major strengths of LPCMCI lie in its significantly improved recall as compared to the SVAR-FCI and SVAR-RFCI baselines for autocorrelated continuous variables, which grows with autocorrelation and is particularly strong for contemporaneous links. At the same time LPCMCI (for $k>0$) has better calibrated CI test leading to better false positive control than the baselines. We cannot prove false positive control, but are not aware of any such proof for other constraint-based algorithms in the challenging latent, nonlinear, autocorrelated setting considered here. A general weakness, which also applies to (SVAR-)FCI and (SVAR-)RFCI, is the faithfulness assumption. If violated in practice this may lead to wrong conclusions. We did not attempt to only assume the weaker form of adjacency-faithfulness [Ramsey et al., 2006], which to our knowledge is however generally an open problem in the causally insufficient case. Moreover, like all constraint-based methods, our method cannot distinguish all members of Markov equivalence classes like methods based on the SCM framework such as e.g. TS-LiNGAM [Hyvärinen et al., 2008] and TiMINo [Peters et al., 2013] do. These, however, restrict the type of dependencies. Concluding, this paper shows how causal discovery in autocorrelated time series benefits from increasing the effect size of CI tests by including causal parents in conditioning sets. The LPCMCI algorithm introduced here implements this idea by entangling the removal and orientation of edges. As demonstrated in extensive simulation studies, LPCMCI achieves much higher recall than the SVAR-FCI and SVAR-RFCI baselines for autocorrelated continuous variables. We further presented novel orientation rules and an extension of graphical terminology by the notions of middle marks and weakly minimal separating sets. Code for all studied methods is provided as part of the tigramite Python package at https://github.com/jakobrunge/tigramite. In future work one may relax assumptions of LPCMCI to allow for selection bias and non-stationarity. Background knowledge about (non-)ancestorships may be included without any conceptual modification. Since the presented orientation rules are applicable at any point and thus able to determine (non-)ancestorships already after having performed only few CI tests, the rules may also be useful for causal feature selection in the presence of hidden confounders, a task that for time series has recently been considered in [Mastakouri et al., 2020]. Lastly, it would be interesting to combine the ideas presented here with ideas from the structural causal model framework.

Broader Impact

Observational causal discovery is especially important for the analysis of systems where experimental manipulation is impossible due to ethical reasons, e.g., in climate research or neuroscience. Our work focuses on the challenging time series case that is of particular relevance in these fields. Understanding causal climate mechanisms from large observational satellite datasets helps climate researchers in understanding and modeling climate change as a main challenge of humanity. Since all code will be published open-source, our methods can be used by anyone. Causal discovery is a rather fundamental topic and we deem the potential for misuse as low.

Acknowledgments and Disclosure of Funding

We thank the anonymous referees for considered and helpful comments that helped to improve the paper. Thanks also goes to Christoph Käding for proof-reading.

DKRZ (Deutsches Klimarechenzentrum) provided computational resources (grant no. 1083).

References

Supplementary material

In this supplementary material we present a brief overview of the FCI algorithm and related graphical terminology as well as details, proofs, further simulation studies, and figures for illustrating the application to the real data example that have been omitted from the main text for reasons of space.

We always assume that there are no selection variables. When saying that $\mathcal{S}$ is a separating set of $A$ and $B$ the exclusions $A\notin\mathcal{S}$ and $B\notin\mathcal{S}$ are implicit. The term subset without the attribute proper refers to both proper subsets and the original set itself, although in formulas we make this explicit by using the symbol $\subseteq$ instead of $\subset$. We switch between using variable names such as $X^{i}_{t-\tau}$ and $X^{j}_{t}$ that make the time structure explicit, and generic names such as $A$ and $B$ that do not make this explicit (using generic names does, however, not imply that there is no time structure). The precise configurations of numerical experiments are given in the respective panel label and figure caption.

S1 Relevant graphical terminology and notation

The structural causal model (SCM) in eq. (1) can be graphically represented by its time series graph (also known as full time graph) $\mathcal{G}$ [Spirtes et al., 2000, Pearl, 2000, Peters et al., 2017]. This graph contains a node for each variable in the SCM (we use the words node and variable interchangeably in this context) and an edge (link, words again used interchangeably) $X^{i}_{t-\tau}{\rightarrow}X^{j}_{t}$ if and only if $X^{i}_{t-\tau}\in pa(X^{j}_{t})$. It can be understood as a directed acyclic graph (DAG) with infinite extension and repeating structure along the time axis. The parents $pa(X^{j}_{t},\mathcal{G})=pa(X^{j}_{t})$ of $X^{j}_{t}$ are the set of nodes $X^{i}_{t-\tau}$ with $X^{i}_{t-\tau}{\rightarrow}X^{j}_{t}$ in $\mathcal{G}$, the ancestors $an(X^{j}_{t},\mathcal{G})$ are the set of nodes connected to $X^{j}_{t}$ by a directed path in $\mathcal{G}$ together with $X^{j}_{t}$ itself (so every node is an ancestor of itself), and the adjacencies $adj(X^{j}_{t},\mathcal{G})$ the set of nodes connected to $X^{j}_{t}$ by any edge in $\mathcal{G}$. Parents are a special case of ancestors. We call $X^{j}_{t}$ a descendant of $X^{i}_{t-\tau}$ if $X^{i}_{t-\tau}$ is an ancestor of $X^{j}_{t}$ (this implies that every node is a descendant of itself). A link between $X^{i}_{t-\tau}$ and $X^{j}_{t}$ is lagged if $\tau>0$, contemporaneous if $\tau=0$, for $i=j$ we speak of an autodependency link, and for $i\neq j$ of a cross link.

In the presence of unobserved variables so called maximal ancestral graphs (MAGs) [Richardson and Spirtes, 2002] provide an appropriate graphical language for representing causal relationships. Since in this paper we assume the absence of selection variables, the relevant MAGs $\mathcal{M}$ contain two types of edges: directed ‘${\rightarrow}$’ and bidirected ‘${\leftrightarrow}$’. These edges are interpreted as composite objects constituted by the symbols at their ends (edge marks), which can be an (arrow-)head (‘>’ or ‘<’) or a tail (‘-’). These edge marks carry a causal meaning: Tails convey ancestorships in $\mathcal{G}$, i.e., $X^{i}_{t-\tau}{\rightarrow}X^{j}_{t}$ in $\mathcal{M}$ asserts that $X^{i}_{t-\tau}\in an(X^{j}_{t},\mathcal{G})$; heads convey non-ancestorships in $\mathcal{G}$, i.e., $X^{i}_{t-\tau}{\rightarrow}X^{j}_{t}$ and $X^{i}_{t-\tau}{\leftrightarrow}X^{j}_{t}$ in $\mathcal{M}$ say that $X^{j}_{t}\notin an(X^{i}_{t-\tau},\mathcal{G})$. As an immediate consequence of time order there cannot be a link $X^{i}_{t-\tau}{\leftarrow}X^{j}_{t}$ for $\tau>0$ (an effect cannot precede its cause). Parents, ancestors and adjacencies are defined in the same way as for DAGs, and the spouses $sp(X^{j}_{t},\mathcal{M})$ of $X^{j}_{t}$ are the set of nodes $X^{i}_{t-\tau}$ with $X^{i}_{t-\tau}{\leftrightarrow}X^{j}_{t}$ in $\mathcal{M}$. Two variables are connected by an edge in $\mathcal{M}$ if and only if they cannot be d-separated by a subset of observed variables in $\mathcal{G}$, and d-separation in $\mathcal{G}$ restricted to observed variables is equivalent to m-separation in $\mathcal{M}$ [Pearl, 1988, Verma and Pearl, 1990, Richardson and Spirtes, 2002]. The parents (ancestors, adjacencies, spouses) of a set of variables are defined as the union of parents (ancestors, adjacencies, spouses) of the individual variables. Example: $pa(\{A,B\},\cdot)=pa(A,\cdot)\cup pa(B,\cdot)$.

The Markov equivalence class of a MAG is the set of all MAGs that yield the exact same set of m-separations [Zhang, 2008]. These are graphically represented by partial ancestral graphs (PAGs), in which the set of allowed edge marks is extended by the circle mark ‘$\circ$’ [Zhang, 2008]. Such a graph is said to be a PAG for MAG $\mathcal{M}$ if $i)$ it has the same nodes and adjacencies as $\mathcal{M}$ and if $ii)$ all its non-circle edge marks are shared by all members in the Markov equivalence class of $\mathcal{M}$. It is further said to be maximally informative if for all its circle marks there is some member of the equivalence class in which there is a tail instead and some other member in which there is a head instead. The wildcard symbol ‘$\ast$’ may stand for all three possible edge marks (head, tail, circle). This is a notational device only, there are no ‘$\ast$’ marks in PAGs.

S2 Some background on FCI

The Fast Causal Inference (FCI) algorithm is an algorithm for constraint-based causal discovery in the presence of unobserved variables [Spirtes et al., 1995, Spirtes et al., 2000, Zhang, 2008]. It allows for both latent confounders and selection variables, although in this paper we assume the absence of selection variables. Under the assumptions of faithfulness [Spirtes et al., 2000], acyclicity, and the existence of an underlying SCM the algorithm determines the maximally informative PAG from perfect statistical decisions of conditional independencies in the distribution $P$ generated by the SCM. The algorithm is based on the following fact:

Let $A$ and $B$ be two nodes such that $A\notin adj(B,\mathcal{M})$ and $B\notin an(A,\mathcal{M})$, then they are m-separated by some subset of $\text{D-Sep}(B,A,\mathcal{M})$. Here:

Node $V\in\mathcal{M}$ is in $\text{D-Sep}(B,A,\mathcal{M})$ if and only if $i)$ it is not $B$ and $ii)$ there is a path $p_{V}$ between $B$ and $V$ such that $iia)$ all nodes on $p_{V}$ are in $an(\{A,B\},\mathcal{M})$ and $iib)$ all non end-point nodes on $p_{V}$ are colliders on $p_{V}$.

A node $B$ is a collider on a path $p$ if the two edges on $p$ involving $B$ both have a head at $B$, as e.g. in $A{\ast\!{\rightarrow}}B{{\leftarrow}\!\ast}C$, otherwise it is a non-collider. Together with acyclicity Proposition S1 guarantees that non-adjacent variables $A$ and $B$ are m-separated by a subset of $\text{D-Sep}(B,A,\mathcal{M})$ or a subset of $\text{D-Sep}(A,B,\mathcal{M})$. However, $\mathcal{M}$ is initially unknown and the D-Sep sets cannot be determined without prior work. Therefore, starting from the complete graph over the set of variables, FCI first performs tests of CI given subset of $pa(B,\mathcal{M}^{\prime})$ and $pa(A,\mathcal{M}^{\prime})$ where $\mathcal{M}^{\prime}$ is the (changing) graph that the algorithm operates on. Whenever two variables are found to be conditionally independent given some subset of variables, the edge between them is removed and their separating set is remembered. This removes some, but in general not all false links. Second, the algorithm orients all resulting unshielded triples $A{\ast\!{\--}\!\ast}B{\ast\!{\--}\!\ast}C$ in $\mathcal{M}^{\prime}$ (these are triples $A{\ast\!{\--}\!\ast}B{\ast\!{\--}\!\ast}C$ such that $A$ and $C$ are not adjacent) as colliders $A{\ast\!{\rightarrow}}B{{\leftarrow}\!\ast}C$ if $B$ is not in the separating set of $A$ and $C$ (rule $\mathcal{R}0$). We note that at this point head marks are not guaranteed to convey non-ancestorships, but those unshielded triples in $\mathcal{M}^{\prime}$ that are part of $\mathcal{M}$ are oriented correctly. This is enough to determine the Possible-D-Sep sets, see [Spirtes et al., 2000], which are supersets of the D-Sep sets define above. Third, FCI performs tests of CI given subsets of $\text{Possible-D-Sep}(B,A,\mathcal{M}^{\prime})$ and $\text{Possible-D-Sep}(A,B,\mathcal{M}^{\prime})$. This removes all false links. Fourth, all previous orientations are undone, $\mathcal{R}0$ is applied once more and then followed by exhaustive application of the ten rules $\mathcal{R}1$ through $\mathcal{R}10$. Tests of CI are preferentially made given smaller conditioning sets $\mathcal{S}$, i.e., FCI first tests sets with $|\mathcal{S}|=p=0$, then those with $|\mathcal{S}|=p=1$ and so on.

S3 LPCMCI-PAGs

Section 3.2 introduced middle marks and LPCMCI-PAGs. We here give a more formal definition of these notions. Recall that we assume the absence of selection variables.

Consider a simple graph $\mathcal{C}(\mathcal{G})$ over the same set of variables as $\mathcal{M}(\mathcal{G})$ with edges of the type ${\mathbin{\ensurestackMath{\stackon[-0pt]{{\rightarrow}}{{\scriptstyle\ast}}}}}$, ${\mathbin{\ensurestackMath{\stackon[-0pt]{{\leftrightarrow}}{{\scriptstyle\ast}}}}}$, ${\mathbin{\ensurestackMath{\stackon[-2pt]{{\circ\!{\rightarrow}}}{{\scriptstyle\ast}}}}}$, and ${\mathbin{\ensurestackMath{\stackon[-2pt]{{\circ\!{\--}\!\circ}}{{\scriptstyle\ast}}}}}$ where the wildcard ‘$\ast$’ can stand for the five possible middle marks ‘?’, ‘L’, ‘R’, ‘!’, or ‘’ (empty). Such $\mathcal{C}(\mathcal{G})$ is a LPCMCI-PAG for $\mathcal{G}$ with respect to total order $<$ if for any probability distribution $P$ that is Markov relative and faithful to $\mathcal{G}$ the following seven conditions hold:

If $A{\mathbin{\ensurestackMath{\stackon[-2pt]{{\ast\!{\rightarrow}}}{{\scriptstyle\ast}}}}}B$ in $\mathcal{C}(\mathcal{G})$, then $B\notin an(A,\mathcal{G})$.

If $A{\mathbin{\ensurestackMath{\stackon[-0pt]{{\rightarrow}}{{\scriptstyle\ast}}}}}B$ in $\mathcal{C}(\mathcal{G})$, then $A\in an(B,\mathcal{G})$.

If $A\notin adj(B,\mathcal{C}(\mathcal{G}))$, then $A\notin adj(B,\mathcal{M}(\mathcal{G}))$.

If $A{\mathbin{\ensurestackMath{\stackon[-1pt]{{\ast\!{\--}\!\ast}}{{\scriptscriptstyle L}}}}}B$ in $\mathcal{C}(\mathcal{G})$ for $A<B$, then $B\notin an(A,\mathcal{G})$ or there is no $\mathcal{S}\subseteq pa(A,\mathcal{M}(\mathcal{G}))$ that m-separates $A$ and $B$ in $\mathcal{M}(\mathcal{G})$.

If $A{\mathbin{\ensurestackMath{\stackon[-1pt]{{\ast\!{\--}\!\ast}}{{\scriptscriptstyle R}}}}}B$ in $\mathcal{C}(\mathcal{G})$ for $A<B$, then $A\notin an(B,\mathcal{G})$ or there is no $\mathcal{S}\subseteq pa(B,\mathcal{M}(\mathcal{G}))$ that m-separates $A$ and $B$ in $\mathcal{M}(\mathcal{G})$.

If $A{\mathbin{\ensurestackMath{\stackon[-2pt]{{\ast\!{\--}\!\ast}}{{\scriptscriptstyle!}}}}}B$ in $\mathcal{C}(\mathcal{G})$, then both $A{\mathbin{\ensurestackMath{\stackon[-1pt]{{\ast\!{\--}\!\ast}}{{\scriptscriptstyle L}}}}}B$ and $A{\mathbin{\ensurestackMath{\stackon[-1pt]{{\ast\!{\--}\!\ast}}{{\scriptscriptstyle R}}}}}B$ would be correct.

If $A{\ast\!{\--}\!\ast}B$ in $\mathcal{C}(\mathcal{G})$, then $B\in adj(A,\mathcal{M}(\mathcal{G}))$.

The first two points give the same causal meaning to head and tail edge marks as they have in MAGs and PAGs. We repeat that while this definition involves a fixed total order $<$ , its choice is arbitrary and without influence on the conveyed causal information. Moreover, the definition does not depend on time order. Also note that if all middle marks in $\mathcal{C}(\mathcal{G})$ are empty, then $\mathcal{C}(\mathcal{G})$ is a PAG for $\mathcal{M}(\mathcal{G})$ (guaranteed by the first, second, third, and seventh point). Parents, ancestors, descendants, spouses, and adjacencies in $\mathcal{C}(\mathcal{G})$ are defined (and denoted) in the same way as for MAGs and PAGs, i.e., without being influenced by middle marks.

S4 Orientation rules for LPCMCI-PAGs

The following is a list of rules for orienting edges in LPCMCI-PAGs. These are extensions of the standard FCI rules [Zhang, 2008] as well as the unshielded triple rule and discriminating path rule of RFCI [Colombo et al., 2012]. If a rule proposes to orient the same edge mark as both tail and head, this is resolved by putting a conflict mark ‘x’ instead. The edge mark wildcard ‘$\ast$’ is redefined to stand for the circle, head, tail or conflict mark; the second wildcard symbol ‘$\star$’ excludes the conflict mark. For two reasons we explicitly present and prove also those rules that generalize without much modification: To demonstrate their validity for LPCMCI-PAGs, and to show in which cases the rules also apply to structures with conflict marks.

If $X{\mathbin{\ensurestackMath{\stackon[-2pt]{{\ast\!{\--}\!\ast}}{{\scriptstyle\ast}}}}}Y{\mathbin{\ensurestackMath{\stackon[-2pt]{{\ast\!{\--}\!\ast}}{{\scriptstyle\ast}}}}}Z$ is an unshielded triple we write $\mathcal{S}_{XZ}$ for the separating set of $X$ and $Z$. Many rules require that $\mathcal{S}_{XZ}$ be weakly minimal and $Y\in\mathcal{S}_{XZ}$. In all these case the requirement of weak minimality can be dropped if $X{\ast\!{\--}\!\ast}Y{\ast\!{\--}\!\ast}Z$, i.e., if both middle marks on $X{\mathbin{\ensurestackMath{\stackon[-2pt]{{\ast\!{\--}\!\ast}}{{\scriptstyle\ast}}}}}Y{\mathbin{\ensurestackMath{\stackon[-2pt]{{\ast\!{\--}\!\ast}}{{\scriptstyle\ast}}}}}Z$ are empty. For this reason the standard FCI orientation rules are implied as special cases.

$\mathbf{\mathcal{R}0^{\prime}a}$: For all unshielded triples $A{\mathbin{\ensurestackMath{\stackon[-2pt]{{\ast\!{\--}\!\ast}}{{\scriptstyle\ast}}}}}B{\mathbin{\ensurestackMath{\stackon[-2pt]{{\ast\!{\--}\!\ast}}{{\scriptstyle\ast}}}}}C$: If $ia)$ $A{\ast\!{\--}\!\ast}B$ or $ib)$ $A$ and $B$ are conditionally dependent given [\mathcal{S}_{AC}\cup pa(\{A,B\},\mathcal{C}(\mathcal{G}))]\setminus\{A,B,\text{nodes in the future of bothA$and$B}\}, $iia)$ $C{\ast\!{\--}\!\ast}B$ or $iib)$ $C$ and $B$ are conditionally dependent given [\mathcal{S}_{AC}\cup pa(\{C,B\},\mathcal{C}(\mathcal{G}))]\setminus\{C,B,\text{nodes in the future of bothC$and$B}\}, $iii)$ none of the edge mark‘$\ast$’s at $B$ on $A{\mathbin{\ensurestackMath{\stackon[-2pt]{{\ast\!{\--}\!\ast}}{{\scriptstyle\ast}}}}}B{\mathbin{\ensurestackMath{\stackon[-2pt]{{\ast\!{\--}\!\ast}}{{\scriptstyle\ast}}}}}C$ is ‘-’ or ‘x’, and $iv$) $B\notin\mathcal{S}_{AC}$, then mark the unshielded triple for orientation as collider $A{\mathbin{\ensurestackMath{\stackon[-2pt]{{\ast\!{\rightarrow}}}{{\scriptstyle\ast}}}}}B{\mathbin{\ensurestackMath{\stackon[-2pt]{{{\leftarrow}\!\ast}}{{\scriptstyle\ast}}}}}C$. Condition $ib)$ need only be checked if not $ia)$, $iib)$ need only be checked if not $iia)$, and $iv)$ need only be checked if all previous conditions are true. If $ib)$ or $iib)$ find a conditional independence, mark the corresponding edge(s) for removal.

$\mathbf{\mathcal{R}0^{\prime}b}$: For all unshielded triples $A{\mathbin{\ensurestackMath{\stackon[-2pt]{{\ast\!{\rightarrow}}}{{\scriptstyle\ast}}}}}B{\mathbin{\ensurestackMath{\stackon[-2pt]{{\circ\!{\--}\!\star}}{{\scriptscriptstyle!}}}}}C$ and for all unshielded triples $A{\mathbin{\ensurestackMath{\stackon[-2pt]{{\ast\!{\rightarrow}}}{{\scriptstyle\ast}}}}}B{\mathbin{\ensurestackMath{\stackon[-1pt]{{\circ\!{\--}\!\star}}{{\scriptscriptstyle R}}}}}C$ with $B<C$ and for all unshielded triples $A{\mathbin{\ensurestackMath{\stackon[-2pt]{{\ast\!{\rightarrow}}}{{\scriptstyle\ast}}}}}B{\mathbin{\ensurestackMath{\stackon[-1pt]{{\circ\!{\--}\!\star}}{{\scriptscriptstyle L}}}}}C$ with $B>C$: If $ia)$ $A{\ast\!{\rightarrow}}B$ or $ib)$ $A$ and $B$ are conditionally dependent given [\mathcal{S}_{AC}\cup pa(\{A,B\},\mathcal{C}(\mathcal{G}))]\setminus\{A,B,\text{nodes in the future of bothA$and$B}\}, and $ii)$ $B\notin\mathcal{S}_{AC}$, then mark the edge between $B$ and $C$ for orientation as $B{\mathbin{\ensurestackMath{\stackon[-2pt]{{{\leftarrow}\!\star}}{{\scriptstyle\ast}}}}}C$ (the middle mark remains as it was before). Condition $ib)$ need only be checked if not $ia)$. If $ib)$ finds a conditional independence, mark the corresponding edge for removal.

$\mathbf{\mathcal{R}0^{\prime}c}$: For all unshielded triples $A{\ast\!{\--}\!\ast}B{\mathbin{\ensurestackMath{\stackon[-2pt]{{\circ\!{\--}\!\star}}{{\scriptscriptstyle!}}}}}C$ and for all unshielded triples $A{\ast\!{\--}\!\ast}B{\mathbin{\ensurestackMath{\stackon[-1pt]{{\circ\!{\--}\!\star}}{{\scriptscriptstyle R}}}}}C$ with $B<C$ and for all unshielded triples $A{\ast\!{\--}\!\ast}B{\mathbin{\ensurestackMath{\stackon[-1pt]{{\circ\!{\--}\!\star}}{{\scriptscriptstyle L}}}}}C$ with $B>C$: If $B\notin\mathcal{S}_{AC}$, then mark the edge between $B$ and $C$ for orientation as $B{\mathbin{\ensurestackMath{\stackon[-2pt]{{{\leftarrow}\!\star}}{{\scriptstyle\ast}}}}}C$ (the middle mark remains as it was before).

$\mathbf{\mathcal{R}0^{\prime}d}$: For all unshielded triples $A{\ast\!{\--}}\!\!\circ\!B{\circ\!{\--}}\!\!\!\ast C$ and for all unshielded triples $A{\ast\!{\rightarrow}}B{\circ\!{\--}}\!\!\!\ast C$: If $B\notin\mathcal{S}_{AC}$, then mark the unshielded triple for orientation as collider $A{\ast\!{\rightarrow}}B{{\leftarrow}\!\ast}C$.

$\mathbf{\mathcal{R}1^{\prime}}$: For all unshielded triples $A{\mathbin{\ensurestackMath{\stackon[-2pt]{{\ast\!{\rightarrow}}}{{\scriptstyle\ast}}}}}B{\mathbin{\ensurestackMath{\stackon[-2pt]{{\circ\!{\--}\!\star}}{{\scriptstyle\ast}}}}}C$: If $\mathcal{S}_{AC}$ is weakly minimal and $B\in\mathcal{S}_{AC}$, then mark the edge between $B$ and $C$ for orientation as $B{\mathbin{\ensurestackMath{\stackon[-0pt]{{\rightarrow}}{{\scriptstyle\ast}}}}}C$.