# Causal Discovery from Subsampled Time Series
with Proxy Variables

Mingzhou Liu1,2
&Xinwei Sun3
&Lingjing Hu4
&Yizhou Wang1,2,5
Correspondence to<EMAIL_ADDRESS>
(

1 School of Computer Science, Peking University
2 Center on Frontiers of Computing Studies (CFCS), Peking University
3 School of Data Science, Fudan University
4 Yanjing Medical College, Capital Medical University
5 Institute for Artificial Intelligence, Peking University
)

###### Abstract

Inferring causal structures from time series data is the central interest of
many scientific inquiries. A major barrier to such inference is the problem of
subsampling, _i.e._ , the frequency of measurement is much lower than that of
causal influence. To overcome this problem, numerous methods have been
proposed, yet either was limited to the linear case or failed to achieve
identifiability. In this paper, we propose a constraint-based algorithm that
can identify the entire causal structure from subsampled time series, without
any parametric constraint. Our observation is that the challenge of
subsampling arises mainly from hidden variables at the unobserved time steps.
Meanwhile, every hidden variable has an observed proxy, which is essentially
itself at some observable time in the future, benefiting from the temporal
structure. Based on these, we can leverage the proxies to remove the bias
induced by the hidden variables and hence achieve identifiability. Following
this intuition, we propose a proxy-based causal discovery algorithm. Our
algorithm is nonparametric and can achieve full causal identification.
Theoretical advantages are reflected in synthetic and real-world experiments.
Our code is available at https://github.com/lmz123321/proxy_causal_discovery.

##  1 Introduction

Temporal systems are the primary subject of causal modeling in many sciences,
such as pathology, neuroscience, and economics. For these systems, a common
issue is the difficulty of collecting sufficiently refined timescale data.
That is, we can only observe a _subsampled_ version of the true causal
interaction. This can cause serious problems for causal identification, since
full observation (_i.e._ , causal sufficiency [1]) is believed to be an
important condition for causal identification, and the violation of which
(_e.g._ , the existence of hidden variables) can induce strong bias [2].

######  _Example_ 1.1.

Take the pathology study of Alzheimer’s disease (AD) as an example. In AD,
patients suffer from memory loss due to the atrophy of memory-related brain
regions such as the Hippocampus [3]. To monitor such atrophy, the standard
protocol is to perform a MRI examination on the brain every six months [4].
However, many studies have shown that the disease can progress much more
rapidly than this [5, 6]. For this reason, when recovering the causal
interactions between brain regions, those at the unobserved time steps
constitute hidden variables and induce spurious edges in the causal graph, as
shown in Figure 1 (b).

Figure 1: Recovered causal pathways in Alzheimer’s disease, with (a) our
method that explicitly models the subsampling and (b) Dynotears [7] that
neglects subsampling. Spurious edges that contradict clinical studies [8, 9]
are marked red. Please refer to Section 5.2 and Figure 7 for details.

To resolve this problem, many attempts have been made to recover the correct
causal relations from subsampled data. However, none of them achieved general
identifiability results. For example, in [10, 11, 12], identifiability was
achieved only for linear data. As for nonlinear data, only a small part of the
causal information (_i.e._ , an equivalence class) could be identified [13,
14, 15, 16]. Indeed, encountering such difficulties is no coincidence,
considering the fact that there are many hidden variables in the system.
Particularly, the hidden mediators prevent us from distinguishing the direct
causation from possible mediation and hence only allow us to recover some
ancestral information.

In this paper, we propose a constraint-based algorithm that can identify the
entire causal structure from subsampled data, without any parametric
assumption. Our departure here is that in time series, every hidden variable
has an observed descendent, which is essentially itself at some observable
time in the future111For example, the hidden variable $X(t_{u})$ at the
unobserved time step $t_{u}$ with $X(t_{u})\to X(t_{u}+1)\to
X(t_{u}+2)\small\to\\!\cdots\\!\to X(t_{o})$ has an observable descendent
$X(t_{o})$ at some observable time $t_{o}$ in the future.. Therefore, we can
use the observed descendent as its _proxy variable_ [17] and adjust for the
bias induced by it. For this purpose, we first represent the ancestral
information, which is composed of both direct causation and hidden mediation,
with the Maximal Ancestral Graph (MAG) [18]. Then, to distinguish causation
from mediation, we use the observed descendants of the hidden mediators as
their proxy variables and adjust for their induced mediation bias. We show
that our strategy can achieve a complete identification of the causal
relations entailed in the time series. It is essentially nonparametric and
only relies on common continuity and differentiability of the structural
equation to conduct valid proxy-based adjustment [19].

Back to the AD example, as shown in Figure 1, our method can recover causal
pathways that align well with existing clinical studies [8, 9], while compared
baselines fail to do so.

Our contributions are summarized below:

  1. 1.

We establish nonparametric causal identification in subsampled time series.

  2. 2.

We propose an algorithm to practically recover the causal graph.

  3. 3.

We achieve more accurate causal identification than others on synthetic and
real data.

The rest of the paper is organized as follows. First, in Section 2, we
introduce the problem setting, basic assumptions, and related literature
needed to understand our work. Second, in Section 3, we establish theories on
structural identifiability. Then, based on these theories, in Section 4, we
introduce the proposed causal discovery algorithm and verify its effectiveness
in Section 5. Finally, we conclude the paper and discuss future works in
Section 6.

##  2 Preliminary

We start by introducing the causal framework and related literature needed to
understand our work.

Time series model. Let $\mathbf{X}(t):=\left[X_{1}(t),...,X_{d}(t)\right]$ be
a multivariate time series with $d$ variables defined at discrete time steps
$t=1,...,T$. We assume the data is generated by a first-order structural
vector autoregression (SVAR) process [16]:

$X_{i}(t)=f_{i}(\mathbf{PA}_{i}(t-1),N_{i}),$ (1)

where $f_{i}$ is the structural function, $\mathbf{PA}_{i}$ is the parent of
$X_{i}$, and $N_{i}$ is the exogenous noise.

Implicit in (1) is the assumptions that cause precedes effect, and that
causation is invariant to time, _i.e._ , the structural function $f_{i}$ and
the causal parents $\mathbf{PA}_{i}$ keeps unchanged across time steps. These
assumptions carry the fundamental beliefs of temporal precedence and stability
on causality [1], therefore are widely adopted by existing works [10, 11, 12,
14, 16].

Subsampling. The subsampling problem means that model (1) can be only observed
every $k$ steps [10, 12, 16]. That is, we can only observe
$\mathbf{X}(1),\mathbf{X}(k+1),...,\mathbf{X}(\lfloor\frac{T}{k}\rfloor k+1)$.

Figure 2: Graph terminologies: (a) Summary graph. (b) Full time DAG with
subsampling factor $k=2$. Variables at observed time steps $t_{1},t_{3},t_{5}$
are marked gray, while those at unobserved time steps $t_{2},t_{4}$ are
dashed. (c) MAG over observed variables. The edges $A(t_{1})\to M(t_{3})$ and
$A(t_{3})\leftrightarrow M(t_{3})$ are (_resp._) induced by the inducing paths
$A(t_{1})\to A(t_{2})\to M(t_{3})$ and $A(t_{3})\leftarrow A(t_{2})\to
M(t_{3})$ in the full time DAG. Please refer to the appendix for detailed
explanations.

Graph terminologies. Causal relations entailed in model (1) can be represented
by causal graphs (Figure 2). We first introduce the _full time Directed
Acyclic Graph (DAG)_ [20], which provides a complete description of the
dynamics in the system.

######  Definition 2.1 (Full time DAG).

Let $G:=(\mathbf{V},\mathbf{E})$ be the associated full time DAG of model (1).
The vertex set $\mathbf{V}:=\\{\mathbf{X}(t)\\}_{t=1}^{T}$, the edge set
$\mathbf{E}$ contains $X_{i}(t-1)\to X_{j}(t)$ iff $X_{i}\in\mathbf{PA}_{j}$.

We assume the full time DAG is Markovian and faithful to the joint
distribution $\mathbb{P}(\mathbf{X}(1),...,\mathbf{X}(T))$.

######  Assumption 2.2 (Markovian and faithfulness).

For disjoint vertex sets
$\mathbf{A},\mathbf{B},\mathbf{Z}\subseteq\mathbf{V}$,
$\mathbf{A}\perp\\!\\!\\!\perp\mathbf{B}|\mathbf{Z}\Leftrightarrow\mathbf{A}\perp\\!\\!\\!\perp_{G}\mathbf{B}|\mathbf{Z}$,
where$\perp\\!\\!\\!\perp_{G}$ denotes _d_ -separation in $G$.

In practice, it is often sufficient to know the causal relations between time
series as a whole, without knowing precisely the relations between time
instants [21]. In this regard, we can summarize the causal relations with the
_summary graph_ [20].

######  Definition 2.3 (Summary graph).

Let $G:=(\mathbf{V},\mathbf{E})$ be the associated summary graph of model (1).
The vertex set $\mathbf{V}:=\mathbf{X}$, the edge set $\mathbf{E}$ contains
$X_{i}\to X_{j}(i\neq j)$ iff $X_{i}\in\mathbf{PA}_{j}$.

In this paper, our goal is to recover the summary graph, given data observed
at $t=1,k+1,...,\lfloor\frac{T}{k}\rfloor k+1$.

For this purpose, we need a structure to represent the (marginal) causal
relations between observed variables and hence link the observation
distribution to the summary graph. Here, we use the _Maximal Ancestral Graph
(MAG)_ [22], since it can represent casual relations when unobserved variables
exist.

Specifically, for the variable set
$\mathbf{V}:=\\{\mathbf{X}(t)\\}_{t=1}^{T}$, let
$\mathbf{O}:=\\{\mathbf{X}(1),...,\mathbf{X}(\lfloor\frac{T}{k}\rfloor
k+1)\\}$ be the observed subset and
$\mathbf{L}:=\mathbf{V}\backslash\mathbf{O}$ be the unobserved subset. Let
$\mathbf{An}(A)$ be the ancestor set of $A$. Then, given any full time DAG $G$
over $\mathbf{V}$, the corresponding MAG $M_{G}$ over $\mathbf{O}$ is defined
as follows:

######  Definition 2.4 (MAG).

In the MAG $M_{G}$, two vertices $A,B\in\mathbf{O}$ are adjacent iff there is
an inducing path. 222A path $p$ between $A,B\in\mathbf{O}$ is an inducing path
relative to $\mathbf{L}$ if every non-endpoint vertex on $p$ is either in
$\mathbf{L}$ or a collider, and every collider on $p$ is an ancestor of either
$A$ or $B$. relative to $\mathbf{L}$ between them in $G$. Edge orientation is:

  1. 1.

$A\to B$ in $M_{G}$ if $A\in\mathbf{An}(B)$ in $G$;

  2. 2.

$A\leftarrow B$ in $M_{G}$ if $B\in\mathbf{An}(A)$ in $G$;

  3. 3.

$A\leftrightarrow B$ in $M_{G}$ if $A\not\in\mathbf{An}(B)$ and
$B\not\in\mathbf{An}(A)$ in $G$.

######  Remark 2.5.

In the MAG, the directed edge $A\to B$ means $A$ is the ancestor of $B$. The
bidirected edge $A\leftrightarrow B$ means there is an unobserved confounder
$U$ between $A$ and $B$.

Proximal causal discovery. Recent works [23, 19] found that one could use the
descendant (_i.e._ , proxy) of an unobserved variable $M$ to differentiate
direct causation $A\to B$ from hidden mediation $A\to M\to B$. Specifically,
suppose that:

######  Assumption 2.6 (TV Lipschitzness).

For each causal pair $A\to B$, the map $a\mapsto\mathbb{P}(B|A=a)$ is
Lipschitz with respect to the Total Variation (TV) distance, that is, there
exists a constant $L$ such that:

$\mathrm{TV}(\mathbb{P}(B|A=a),\mathbb{P}(B|A=a^{\prime}))\leq
L|a-a^{\prime}|.$

######  Assumption 2.7 (Completeness).

For each causal pair $A\to B$, the conditional distribution $\mathbb{P}(B|A)$
is complete, that is, for any function $g$, we have:

$E\\{g(b)|a\\}=0\,\,\text{almost surely}\,\,\text{if and only
if}\,\,g(b)=0\,\,\text{almost surely}.$

######  Remark 2.8.

Roughly, Assumption 2.6 requires the continuity of the structural function
$f_{i}$ and Assumption 2.7 requires that the characterization function of the
exogenous noise $N_{i}$ is non-zero. Please refer to [19] for a detailed
discussion.

Under the above assumptions, [19] showed that:

######  Theorem 2.9 ([19]).

Suppose that $A$ and $B$ are mediated by an unobserved variable $M$, and that
$M$ has a proxy variable $M^{\prime}$ satisfying
$A\perp\\!\\!\\!\perp_{G}M^{\prime}|M$. Then, we can identify whether $A\to B$
by testing whether $pr(b|a)\sim pr(m^{\prime}|a)$ has a linear relation. If
the linear relation exists, then $A\not\to B$, otherwise, we have$A\to B$.

######  Remark 2.10.

The linear relation can be tested by checking the least square residual of
regressing $pr(b|a)$ on $pr(m^{\prime}|a)$. Please refer to [19] for details.

To ensure the existence of the proxy variable, we require the following self
causation assumption, which states that each variable $X_{i}(t)$ is influenced
by its past $X_{i}(t-1)$.

######  Assumption 2.11 (Self causation).

In model (1), $X_{i}\in\mathbf{PA}_{i}$ for each $i=1,...,d$.

######  Remark 2.12.

In time series statistics, this is also known as the autocorrelation
assumption [24, 25].

##  3 Structural identifiability

In this section, we establish the identifiability of the summary graph. For
this purpose, we connect the observational distribution to the summary graph
with the bridge of the MAG.

Specifically, our analysis consists of three progressive results: Proposition
3.1, Proposition 3.3, and Theorem 3.5. First, in Proposition 3.1, we prove
that the MAG can be identified from observational data. Then, in Proposition
3.3, we show that the identified MAG _almost uniquely_ reflects the structure
of the summary graph, in the sense that only edges connecting vertices and
their ancestors cannot be determined due to hidden mediators. Finally, in
Theorem 3.5, we identify these underdetermined edges with the proxies of the
hidden mediators and therefore identify the whole summary graph.

Next, we first introduce Proposition 3.1, which shows that the MAG is
identifiable.

######  Proposition 3.1 (Identifiability of the MAG).

Assuming model (1) and Assumption 2.2, then the MAG over the observed variable
set $\mathbf{O}$ is identifiable, i.e., its skeleton and edge orientations can
be uniquely derived from the joint distribution $\mathbb{P}(\mathbf{O})$.

Since the identifiability of the MAG’s skeleton under the faithfulness
assumption is a well-known result [18], we focus on explaining the
identification of edge orientations. Specifically, we will divide the edges
into two classes: the instantaneous edges and the lagged edges, and discuss
their orientations respectively.

The _instantaneous edge_ , _e.g._ , $A(t_{3})\leftrightarrow M(t_{3})$, is an
edge that connects two vertices at the same time. Since we assume cause must
precede effect, the instantaneous edge does not represent causation but latent
confounding, therefore is bidirected according to Definition 2.4. For example,
in Figure 2 (c), the instantaneous edge $A(t_{3})\leftrightarrow M(t_{3})$
represents the latent confounding $A(t_{3})\leftarrow A(t_{2})\to M(t_{3})$
between $A(t_{3})$ and $M(t_{3})$. On the other hand, the _lagged edge_ ,
_e.g._ , $A(t_{1})\to M(t_{3})$, is an edge that connects two vertices at
different time. The lagged edge represents ancestral information and hence is
directed ($\to$) from the past to the future.

To connect the identified MAG to the summary graph, we first define the
following graph structures, examples of which are shown in Figure 3:

######  Definition 3.2 (Graph structures).

In the summary graph,

  1. 1.

A _directed path_ $p_{AB}$ from $A$ to $B$ with length $l$ is a sequence of
distinct vertices $A,V_{1},...,V_{l},B$ where each vertex points to its
successor. Two directed paths $p_{AB_{1}},p_{AB_{2}}$ are called _disjoint_ if
they do not share any non-startpoint vertex.

  2. 2.

A _confounding structure_ $c_{AB}$ between $A$ and $B$ with lengths $(r,q)$
consists of a vertex $U$, a directed path $p_{UA}$ from $U$ to $A$ with length
$r$, and a directed path $p_{UB}$ from $U$ to $B$ with length $q$, where
$p_{UA}$ and $p_{UB}$ are disjoint.

Figure 3: Illustration of Definition 3.2. (a) A directed path from $A$ to $B$
with length $l$. (b) A confounding structure between $A$ and $B$ with lengths
$(r,q)$.

In the following, we explain the relationship between the identified MAG and
the summary graph. In particular, we will discuss what the directed and
bidirected edges in the MAG (_resp._) imply about the summary graph.

According to Definition 2.4, the _directed edge_ $A(t)\to B(t+k)$ in the MAG
means $A$ is the ancestor of $B$. Therefore, in the summary graph, there is
either $A\to B$ or a directed path from $A$ to $B$. On the other hand, the
_bidirected edge_ $A(t+k)\leftrightarrow B(t+k)$ in the MAG means there is a
latent confounder between them, which can be either $A(t+1)$, _i.e._ ,
$A(t+k)\small\leftarrow\\!\cdots\\!\leftarrow A(t+1)\small\to\\!\cdots\\!\to
B(t+k)$, or a third variable $U(t+1)$, _i.e._ ,
$A(t+k)\small\leftarrow\\!\cdots\\!\leftarrow U(t+1)\small\to\\!\cdots\\!\to
B(t+k)$. Hence, in the summary graph, there is either a directed path from $A$
to $B$ or a confounding structure between them.

Summarizing the above observations, we have the following proposition:

######  Proposition 3.3 (MAG to summary graph).

If there are $A(t)\to B(t+k)$ and $A(t+k)\leftrightarrow B(t+k)$ in the MAG,
then, in the summary graph, there is either

1. $A\to B$,  2. a directed path from $A$ to $B$ with length $l\leq k-2$, or

  1. 3.

a directed path from $A$ to $B$ with length $l=k-1$ and a confounding
structure between them with lengths $(r\leq k-2,q\leq k-2)$.

######  Remark 3.4.

In Proposition 3.3, the maximum length of the directed path is $k-1$, this is
because any path longer than this can not be an inducing path and therefore
can not induce $A(t)\to B(t+k)$ in the MAG. Please refer to the appendix for
details.

Proposition 3.3 provides a necessary condition for having $A\to B$ in the
summary graph, _i.e._ , there are both $A(t)\to B(t+k)$ and
$A(t+k)\leftrightarrow B(t+k)$ in the MAG. It also inspires us that, to make
the condition sufficient, we need to distinguish the direct effect ($A\to B$)
from the indirect one ($A\small\to\\!\cdots\\!\to B$).

For this purpose, we propose to use the proximal causal discovery method.
Before preceding any technical detail, we first brief our idea below. The
indirect effect, unlike the direct one, relies on mediation and therefore can
be _d_ -separated by the (unobserved) mediator. To test the _d_ -separation,
we can use the observed proxy of the mediator, thanks to the self causation
assumption in Assumption 2.11. For example, in Figure 4 (b), the indirect
effect $A(t_{1})\to M(t_{2})\to B(t_{3})$ from $A$ to $B$ can be _d_
-separated by the mediator $M(t_{2})$, who has an observed proxy $M(t_{3})$.

For the general case, the separation set is a bit more complicated to ensure
all paths except the one representing direct effect are _d_ -separated.
Specifically, for two vertices $A(t)$ and $B(t+k)$, the separation set is the
union of two sets $\mathbf{M}(t+1)\cup\mathbf{S}(t)$. The set
$\mathbf{M}(t+1)$ contains possible mediators, namely $A(t+1)$ and any vertex
$M_{i}(t+1)$ ($M_{i}\neq B$) such that $A(t)\to M_{i}(t+k)$ in the MAG and
$M_{i}$ is not $B$’s descendant333This can be justified from the MAG.
Specifically, if $M_{i}$ is the descendant of $B$, then there is $B(t)\to
V_{1}(t+k),V_{1}(t)\to V_{2}(t+k),...,V_{l}(t)\to M_{i}(t+k)$ in the MAG.. The
set $\mathbf{S}(t)$ is used to _d_ -separate possible back-door paths between
$A(t)$ and $B(t+k)$, it contains any vertex $S_{i}(t)$ ($S_{i}\neq A$) such
that $S_{i}(t)\to B(t+k)$ or $S_{i}(t)\to M_{j}(t+k)$ for some
$M_{j}\in\mathbf{M}$ in the MAG.

Equipped with the separation set, we then have the following identifiability
result:

Figure 4: Distinguishing the direct effect ($A\to B$) from the indirect one
($A\small\to\\!\cdots\\!\to B$) with proxy variables. Note that though (a) and
(b) have the same MAG, only the indirect effect $A(t_{1})\to M(t_{2})\to
B(t_{3})$ in (b) can be _d_ -separated by $M(t_{2})$. To test the _d_
-separation, we can use $M(t_{3})$ as the proxy variable.

######  Theorem 3.5 (Identifiability of the summary graph).

Assuming model (1), Assumption 2.2, and Assumptions 2.6, 2.7, 2.11, then the
summary graph is identifiable. Specifically,

  1. 1.

There is $A\to B$ in the summary graph iff there are $A(t)\to
B(t+k),A(t+k)\leftrightarrow B(t+k)$ in the MAG, and the set
$\mathbf{M}(t+1)\cup\mathbf{S}(t)$ is not sufficient to d-separate
$A(t),B(t+k)$ in the full time DAG.

  2. 2.

The condition “the set $\mathbf{M}(t+1)\cup\mathbf{S}(t)$ is not sufficient to
d-separate $A(t),B(t+k)$ in the full time DAG” can be tested by the proxy
variable $\mathbf{M}(t+k)$ of the unobserved set $\mathbf{M}(t+1)$.

Theorem 3.5 provides a necessary and sufficient condition for having $A\to B$
in the summary graph. This condition is testable with the proxy
$\mathbf{M}(t+k)$ of the unobserved set $\mathbf{M}(t+1)$ ($\mathbf{S}(t)$ is
observable). The validity of such proxying, _i.e._ ,
$A(t)\perp\\!\\!\\!\perp_{G}\mathbf{M}(t+k)|\mathbf{M}(t+1),\mathbf{S}(t)$ is
proved in the appendix.

##  4 Discovery algorithm

Equipped with the identifiability results in Proposition 3.1 and Theorem 3.5,
we can practically discover the summary graph with Algorithm 1. Specifically,
given the observed data points, Algorithm 1 first recovers the MAG according
to Proposition 3.1. Then, for each pair of vertices in the summary graph, it
first uses the necessary condition of “having $A(t)\to B(t+k)$ and
$A(t+k)\leftrightarrow B(t+k)$ in the MAG” to construct underdetermined edges.
Finally, among these underdetermined edges, it uses the proxy variable
approach to identify direct causation and remove the indirect one, according
to Theorem 3.5.

Below, we first introduce the _Partially Determined Directed Acyclic Graph
(PD-DAG)_ , which will be used to represent the intermediate result of the
algorithm.

######  Definition 4.1 (PD-DAG).

A PD-DAG is a DAG with two kinds of directed edges: solid ($\to$) and dashed
($\dashrightarrow$). A solid edge represents a determined causal relation,
while a dashed one means the causal relation is still underdetermined. We use
a meta-symbol, asterisk ($*\\!\\!\to$), to denote any of the two edges. In a
PD-DAG $G=(\mathbf{V},\mathbf{E})$, the vertex $A$ points to the vertex $B$ if
the edge $A*\\!\\!\to B\in\mathbf{E}$. The definition of the directed path and
confounding structure are the same as in the summary graph.

Data: Data observed at $t=1,k+1,...,\lfloor\frac{T}{k}\rfloor k+1$

Result: The summary graph

1 Construct the skeleton of the MAG with [26], orient edges according to
Proposition 3.1;

For every vertex pair $A,B$, set $A*\\!\\!\to B$ if there are $A(t)\to B(t+k)$
and $A(t+k)\leftrightarrow B(t+k)$ in the MAG. Called the resulted PD-DAG
$G_{0}$ ;

/* Proposition 3.3 */

2

Set $G=G_{0}$ and iteratively execute:

  1. (a)

For every $A*\\!\\!\to B$, check

    1. (i)

a directed path from $A$ to $B$ with length $l\leq k-2$;

    2. (ii)

a directed path from $A$ to $B$ with length $l\leq k-1$ and a confounding
structure between $A$ and $B$ with lengths $(r\leq k-2,q\leq k-2)$;

If $\exists\,$ (i) or (ii), pass. Otherwise, set $A\to B$ ;

/* Proposition 3.3 */

* (b)

Randomly pick a dashed edge $A*\\!\\!\to B$, check whether
$A(t)\perp\\!\\!\\!\perp_{G}B(t+k)|\mathbf{M}(t+1),\mathbf{S}(t)$;

If the _d_ -separation holds, remove $A*\\!\\!\to B$. Otherwise, set $A\to B$
;

3 /* Theorem 3.5 */

until all edges in $G$ are solid ones; return $G$

Algorithm 1 Discover the summary graph

######  Remark 4.2.

Algorithm 1 can be generalized to cases where the subsampling factor $k$ is
unknown, by skipping the step (a) below line 3.

##  5 Experiment

In this section, we evaluate our method on synthetic data and a real-world
application, _i.e._ , discovering causal pathways in Alzheimer’s disease.

Compared baselines. _Methods that account for subsampling_ :  1. SVAR-FCI [16]
that extended the FCI algorithm to time series and recovered a MAG over
observed variables;  2. NG-EM [10] that achieved identifiability with linear
non-Gaussianity and used Expectation Maximization (EM) for estimation.
_Methods that neglect subsampling_ :  1. Dynotears [7] that extended the
score-based Notears [27] algorithm to the time series;  2. PC-GCE [21] that
modified the PC algorithm for temporal data with a new information-theoretic
conditional independence (CI) test.

Metrics. We use the $\mathrm{F}_{1}$-score, precision, and recall, where
precision and recall (_resp._) measure the accuracy and completeness of
identified causal edges, and
$\mathrm{F}_{1}:=2\cdot\frac{\mathrm{precision}\cdot\mathrm{recall}}{\mathrm{precision}+\mathrm{recall}}$.

Implementation details. The significance level is set to $0.05$. For the MAG
recovery, we use the FCI algorithm implemented in the $\mathrm{causallearn}$
package444https://github.com/py-why/causal-learn. Temporal constraints, _e.g._
, causal precedence, time invariance, and self causation are added as the
background knowledge.

###  5.1 Synthetic study

Data generation. We generate radnom summary graphs with the Erdos-Renyi model
[28], where the vertex number is set to $5$, the probability of each edge is
set to $0.3$. For each graph, we generate temporal data with the structural
equation $X_{i}(t)=\sum_{j\in\mathbf{PA}_{i}}f_{ij}(X_{j}(t-1))+N_{i}$, where
the function $f_{ij}$ is randomly chosen from $\\{linear,sin,tanh,sqrt\\}$,
the exogenous noise $N_{i}$ is randomly sampled from
$\\{uniform,gauss.,exp.,gamma\\}$. We consider different subsampling factors
$k=\\{2,3,4,5\\}$ and sample sizes $n=\\{600,800,1000,1200\\}$. For each
setting555We also consider different graph scales $d=\\{5,15,25,35,45\\}$.
Please refer to the appendix for details., we replicate over $100$ random
seeds.

Comparison with baselines. Figure 5 shows the performance of our method and
baselines under different subsampling factors (upper row) and sample sizes
(lower row). As we can see, our method significantly outperforms the others in
all settings. Specifically, compared with methods that account for subsampling
(SVAR-FCI and NG-EM), our method achieves both higher precision and recall,
indicating less detection error and missing edges. This advantage can be
attributed to the fact that our method enjoys an identifiability guarantee
(_v.s._ SVAR-FCI) and meanwhile requires no parametric assumption on the
structural model (_v.s._ NG-EM). Compared with methods that ignore subsampling
(Dynotears and PC-GCE), our method achieves higher $\mathrm{F}_{1}$-score,
precision, and is comparable in recall. This result shows that our method can
effectively reduce the spurious detection induced by unobserved time steps.
Besides, we can observe that Dynotears and PC-GCE slightly outperform SVAR-FCI
and NG-EM, which again demonstrates the necessity of establishing
nonparametric identification in subsampling problems.

Figure 5: Performance of our method and baselines under different subsampling
factors (upper row) and sample sizes (lower row). Figure 6: Evaluation of
intermediate results. We report $F_{1}$-score, precision, and recall of the
recovered MAG, edges identified by proxies, and the recovered summary graph.

Intermediate results. We evaluate the intermediate results of Algorithm 1
(MAG, edges identified by proxies, summary graph) and report the results in
Figure 6. As shown, both the recovered MAG and edges identified by proxies
have high accuracy under moderate subsampling factors and sample sizes, hence
explaining the effectiveness of our algorithm in recovering the summary graph.
Meanwhile, we can also observe that the performance slightly decreases when
the subsampling factor $k$ is large. This is explained as follows. When $k$ is
large, there are many unobserved steps between two observations, which weakens
the correlation pattern in data and therefore breaks our faithfulness
assumption.

###  5.2 Discovering causal pathways in Alzheimer’s disease

Background. Alzheimer’s disease (AD) is one of the most common neuro-
degenerative diseases. In AD, patients suffer from memory loss due to atrophy
of memory-related brain regions, such as the Hippocampus [3] and Temporal lobe
[29]. A widely accepted explanation for such atrophy is: the disease releases
toxic proteins, _e.g._ , $A\beta$ [30] and ${tau}$ [31]; along anatomical
pathways, these proteins spread from one brain region to another, eventually
leading to atrophy of the whole brain [32]. Recovering these anatomical
pathways, _i.e._ , underlying causal mechanisms, will benefit the
understanding of AD pathology and inspire potential treatment methods.

Dataset and preprocessing. We consider the Alzheimer’s Disease Neuroimaging
Initiative (ADNI) dataset [4], in which the imaging data is acquired from
structural Magnetic Resonance Imaging (sMRI) scans. We apply the Dartel VBM
[33] for preprocessing and the Statistical Parametric Mapping (SPM) [34] for
segmenting brain regions. Then, we implement the Automatic Anatomical Labeling
(AAL) atlas [35] to partition the whole brain into $90$ regions. In total, we
use $n=558$ subjects with baseline and month-6 follow-up visits enrolled in
ADNI-GO/1/2/3 periods.

Figure 7: Recovered summary graph in Alzheimer’s disease.

Results. Figure 7 shows the recovered summary graph over $90$ brain regions.
For illustration, in Figure 1 (a), we further group the brain regions into $8$
meta-regions according to their anatomical structures. We can observe that i)
most causal pathways between meta-regions are unidirectional; ii) the
identified sources of atrophy are the Hippocampus, Amygdala, and Temporal
lobe, which are all early-degenerate regions [3, 29]; iii) the topology order
induced from our result, _i.e._ ,
$\\{Hip.,Amy.,Tem.,Cin.,Ins.,Fro.,Par.,Occ.\\}$ coincides with the temporal
degeneration order found in [8, 9]. These results constitute an important
finding that can be supported by existing studies: the atrophy (releasing of
toxic proteins) is sourced from the Hippocampus and gradually propagates to
other brain structures along certain anatomical pathways [32, 36].

In contrast, in Figure 1 (b), the causal relations recovered by the Dynotears
baseline666Please refer to the appendix for results of the other baselines.
are less clear. For example, most of the identified interactions are
bidirectional, which may be due to the spurious edges induced by subsampling.
Besides, the identified atrophy source is the Frontal lobe, with AD-related
regions identified as outcomes, which is inconsistent with clinical studies.
These results, from another perspective, show the importance of modeling
subsampling in time series causal discovery.

##  6 Conclusion

In this paper, we propose a causal discovery algorithm for subsampled time
series. Our method leverages the recent progress of proximal causal discovery
and can achieve complete identifiability without any parametric assumption.
The proposed algorithm can outperform baselines on synthetic data and recover
reasonable causal pathways in Alzheimer’s disease.

Limitation and future work. Our method relies on an accurate test of the
conditional independence (CI) and therefore may suffer from low recall when
the CI patterns in data are weak. To solve this problem, we will investigate
theories on high-efficiency CI tests and pursue a dedicated solution.

## Acknowledgments and Disclosure of Funding

This work was supported by National Key R&D Program of China (2022ZD0114900).

## References

  * [1] Judea Pearl.  Causality.  Cambridge University Press, 2009.
  * [2] Patrick Suppes.  Causal analysis of hidden variables.  In PSA: Proceedings of the Biennial Meeting of the Philosophy of Science Association, volume 1980, pages 563–571. Cambridge University Press, 1980.
  * [3] Mark J West, Paul D Coleman, Dorothy G Flood, and Juan C Troncoso.  Differences in the pattern of hippocampal neuronal loss in normal ageing and alzheimer’s disease.  The Lancet, 344(8925):769–772, 1994.
  * [4] Ronald Carl Petersen, Paul S Aisen, Laurel A Beckett, Michael C Donohue, Anthony Collins Gamst, Danielle J Harvey, Clifford R Jack, William J Jagust, Leslie M Shaw, Arthur W Toga, et al.  Alzheimer’s disease neuroimaging initiative (adni): clinical characterization.  Neurology, 74(3):201–209, 2010.
  * [5] Rachelle Smith Doody, Paul Massman, and J Kay Dunn.  A method for estimating progression rates in alzheimer disease.  Archives of neurology, 58(3):449–454, 2001.
  * [6] Craig J Thalhauser and Natalia L Komarova.  Alzheimer’s disease: rapid and slow progression.  Journal of the Royal Society Interface, 9(66):119–126, 2012.
  * [7] Roxana Pamfil, Nisara Sriwattanaworachai, Shaan Desai, Philip Pilgerstorfer, Konstantinos Georgatzis, Paul Beaumont, and Bryon Aragam.  Dynotears: Structure learning from time-series data.  In International Conference on Artificial Intelligence and Statistics, pages 1595–1605. PMLR, 2020.
  * [8] Alexandra L Young, Razvan V Marinescu, Neil P Oxtoby, Martina Bocchetta, Keir Yong, Nicholas C Firth, David M Cash, David L Thomas, Katrina M Dick, Jorge Cardoso, et al.  Uncovering the heterogeneity and temporal complexity of neurodegenerative diseases with subtype and stage inference.  Nature communications, 9(1):4273, 2018.
  * [9] Jacob W Vogel, Alexandra L Young, Neil P Oxtoby, Ruben Smith, Rik Ossenkoppele, Olof T Strandberg, Renaud La Joie, Leon M Aksman, Michel J Grothe, Yasser Iturria-Medina, et al.  Four distinct trajectories of tau deposition identified in alzheimer’s disease.  Nature medicine, 27(5):871–881, 2021.
  * [10] Mingming Gong, Kun Zhang, Bernhard Schoelkopf, Dacheng Tao, and Philipp Geiger.  Discovering temporal causal relations from subsampled data.  In International Conference on Machine Learning, pages 1898–1906. PMLR, 2015.
  * [11] Mingming Gong, Kun Zhang, Bernhard Schölkopf, Clark Glymour, and Dacheng Tao.  Causal discovery from temporally aggregated time series.  In Uncertainty in artificial intelligence: proceedings of the… conference. Conference on Uncertainty in Artificial Intelligence, volume 2017\. NIH Public Access, 2017.
  * [12] Alex Tank, Emily B Fox, and Ali Shojaie.  Identifiability and estimation of structural vector autoregressive models for subsampled and mixed-frequency time series.  Biometrika, 106(2):433–452, 2019.
  * [13] David Danks and Sergey Plis.  Learning causal structure from undersampled time series.  2013\.
  * [14] Sergey Plis, David Danks, Cynthia Freeman, and Vince Calhoun.  Rate-agnostic (causal) structure learning.  Advances in neural information processing systems, 28, 2015.
  * [15] Antti Hyttinen, Sergey Plis, Matti Järvisalo, Frederick Eberhardt, and David Danks.  Causal discovery from subsampled time series data by constraint optimization.  In Conference on Probabilistic Graphical Models, pages 216–227. PMLR, 2016.
  * [16] Daniel Malinsky and Peter Spirtes.  Causal structure learning from multivariate time series in settings with unmeasured confounding.  In Proceedings of 2018 ACM SIGKDD workshop on causal discovery, pages 23–47. PMLR, 2018.
  * [17] Manabu Kuroki and Judea Pearl.  Measurement bias and effect restoration in causal inference.  Biometrika, 101(2):423–437, 2014.
  * [18] P Spirtes, C Meek, and T Richardson.  Computation, causation and discovery, 1999.
  * [19] Mingzhou Liu, Xinwei Sun, Yu Qiao, and Yizhou Wang.  Causal discovery with unobserved variables: a proxy variable approach.  arXiv preprint arXiv:2305.05281, 2023.
  * [20] Chang Gong, Di Yao, Chuzhe Zhang, Wenbin Li, and Jingping Bi.  Causal discovery from temporal data: An overview and new perspectives.  arXiv preprint arXiv:2303.10112, 2023.
  * [21] Charles K Assaad, Emilie Devijver, and Eric Gaussier.  Discovery of extended summary graphs in time series.  In Uncertainty in Artificial Intelligence, pages 96–106. PMLR, 2022\.
  * [22] Jiji Zhang.  On the completeness of orientation rules for causal discovery in the presence of latent confounders and selection bias.  Artificial Intelligence, 172(16-17):1873–1896, 2008.
  * [23] Wang Miao, Zhi Geng, and Eric J Tchetgen Tchetgen.  Identifying causal effects with proxy variables of an unmeasured confounder.  Biometrika, 105(4):987–993, 2018.
  * [24] Maurice S Bartlett.  On the theoretical specification and sampling properties of autocorrelated time-series.  Supplement to the Journal of the Royal Statistical Society, 8(1):27–41, 1946.
  * [25] ZA Lomnicki and SK Zaremba.  On the estimation of autocorrelation in time series.  The Annals of Mathematical Statistics, 28(1):140–158, 1957.
  * [26] Peter Spirtes, Clark N Glymour, Richard Scheines, and David Heckerman.  Causation, prediction, and search.  MIT press, 2000.
  * [27] Xun Zheng, Bryon Aragam, Pradeep K Ravikumar, and Eric P Xing.  Dags with no tears: Continuous optimization for structure learning.  Advances in neural information processing systems, 31, 2018.
  * [28] Paul Erdős, Alfréd Rényi, et al.  On the evolution of random graphs.  Publ. Math. Inst. Hung. Acad. Sci, 5(1):17–60, 1960.
  * [29] Clifford R Jack, Ronald C Petersen, Yuecheng Xu, Peter C O’Brien, Glenn E Smith, Robert J Ivnik, Eric G Tangalos, and Emre Kokmen.  Rate of medial temporal lobe atrophy in typical aging and alzheimer’s disease.  Neurology, 51(4):993–999, 1998.
  * [30] Frank M LaFerla, Kim N Green, and Salvatore Oddo.  Intracellular amyloid-$\beta$ in alzheimer’s disease.  Nature Reviews Neuroscience, 8(7):499–509, 2007.
  * [31] Michala Kolarova, Francisco García-Sierra, Ales Bartos, Jan Ricny, and Daniela Ripova.  Structure and pathology of tau protein in alzheimer disease.  International journal of Alzheimer’s disease, 2012, 2012.
  * [32] George S Bloom.  Amyloid-$\beta$ and tau: the trigger and bullet in alzheimer disease pathogenesis.  JAMA neurology, 71(4):505–508, 2014.
  * [33] John Ashburner.  Vbm tutorial.  Wellcome Trust Centre for Neuroimaging, London, UK, 2010.
  * [34] Karl J Friston.  Statistical parametric mapping.  Neuroscience databases: a practical guide, pages 237–250, 2003\.
  * [35] Nathalie Tzourio-Mazoyer, Brigitte Landeau, Dimitri Papathanassiou, Fabrice Crivello, Octave Etard, Nicolas Delcroix, Bernard Mazoyer, and Marc Joliot.  Automated anatomical labeling of activations in spm using a macroscopic anatomical parcellation of the mni mri single-subject brain.  Neuroimage, 15(1):273–289, 2002.
  * [36] Mingzhou Liu, Xiangyu Zheng, Xinwei Sun, Fang Fang, and Yizhou Wang.  Which invariance should we transfer? A causal minimax learning approach.  In Proceedings of the 40th International Conference on Machine Learning, volume 202, pages 22488–22527. PMLR, 23–29 Jul 2023.

###### Appendix

  1. 1 Introduction
  2. 2 Preliminary
  3. 3 Structural identifiability
  4. 4 Discovery algorithm
  5. 5 Experiment
    1. 5.1 Synthetic study
    2. 5.2 Discovering causal pathways in Alzheimer’s disease
  6. 6 Conclusion
  7. A Preliminary
    1. A.1 Explanation of Figure 2(c): MAG
    2. A.2 Proof of Remark 2.5: Edge orientations in the MAG
  8. B Structural identifiability
    1. B.1 Proof of Proposition 3.1: Identifiability of the MAG
    2. B.2 Proof of Proposition 3.3: MAG to summary DAG
    3. B.3 Proof of Theorem 3.5: Identifiability of the summary DAG
  9. C Experiment
    1. C.1 Extra results of Section 5.1: Synthetic study
    2. C.2 Extra results of Section 5.2: Discovering causal pathways in Alzheimer’s disease

##  Appendix A Preliminary

###  A.1 Explanation of Figure 2(c): MAG

Table 1: Detailed explanation of edges in the MAG Edge in the MAG | Inducing path in the full time DAG
---|---
$A(t_{1})\to A(t_{3})$ | $A(t_{1})\to A(t_{2})\to A(t_{3})$
$A(t_{1})\to M(t_{3})$ | $A(t_{1})\to A(t_{2})\to M(t_{3})$
$A(t_{1})\to M(t_{2})\to M(t_{3})$
$A(t_{1})\to B(t_{3})$ | $A(t_{1})\to M(t_{2})\to B(t_{3})$
$M(t_{1})\to M(t_{3})$ | $M(t_{1})\to M(t_{2})\to M(t_{3})$
$M(t_{1})\to B(t_{3})$ | $M(t_{1})\to M(t_{2})\to B(t_{3})$
$M(t_{1})\to B(t_{2})\to B(t_{3})$
$B(t_{1})\to B(t_{3})$ | $B(t_{1})\to B(t_{2})\to B(t_{3})$
$A(t_{3})\leftrightarrow M(t_{3})$ | $A(t_{3})\leftarrow A(t_{2})\to M(t_{3})$
$M(t_{3})\leftrightarrow B(t_{3})$ | $M(t_{3})\leftarrow M(t_{2})\to B(t_{3})$

###  A.2 Proof of Remark 2.5: Edge orientations in the MAG

Remark 2.5._In the MAG, the directed edge $A\to B$ means $A$ is the ancestor
of $B$. The bidirected edge $A\leftrightarrow B$ means there is an unobserved
confounder $U$ between $A$ and $B$._

###### Proof.

Since the meaning of the directed edge is directly induced from the
definition, we focus on explaining the bidirected edge.

Specifically, $A\leftrightarrow B$ means there is an inducing path $p$
relative to $\mathbf{L}$ between $A$ and $B$. Since $A\not\in\mathbf{An}(B)$
and $B\not\in\mathbf{An}(A)$, the inducing path $p$ must contains non-
mediators, _i.e._ , colliders or confounders. Suppose that there are $r$
confounders on $p$, since between every two confounders there is a collider,
the number of colliders on $p$ is $r-1$. Denote the confounders on $p$ as
$U_{1},...,U_{r}$, the colliders as $C_{1},...,C_{r-1}$, we have
$p=A\small\leftarrow\\!\cdots\\!\leftarrow U_{1}\small\to\\!\cdots\\!\to
C_{1}\small\leftarrow\\!\cdots\\!\leftarrow\cdots\small\to\\!\cdots\\!\to
C_{h-1}\small\leftarrow\\!\cdots\\!\leftarrow U_{h}\small\to\\!\cdots\\!\to
B$. According to the definition of inducing path, each collider $C_{i}$ on $p$
is the ancestor of either $A$ or $B$. Based on this, the following algorithm
is assured to find a latent confounder between $A$ and $B$.

1 for  _$i=1,...,r-1$_ do

2 if  _$C_{i}\in\mathbf{An}(B)$_ then

return $U_{i}$ ;

/* $A\small\leftarrow\\!\cdots\\!\leftarrow
C_{i-1}\small\leftarrow\\!\cdots\\!\leftarrow U_{i}\small\to\\!\cdots\\!\to
C_{i}\small\to\\!\cdots\\!\to B$ */

3

4 end if

5

6 end for

return $U_{r}$ ;

/* $A\small\leftarrow\\!\cdots\\!\leftarrow
C_{r-1}\small\leftarrow\\!\cdots\\!\leftarrow U_{r}\small\to\\!\cdots\\!\to B$
*/

Algorithm 2 Find the latent confounder

∎

##  Appendix B Structural identifiability

###  B.1 Proof of Proposition 3.1: Identifiability of the MAG

Proposition 3.1. _Assuming model (1) and Assumption 2.2, then the MAG over the
observed variable set $\mathbf{O}$ is identifiable, i.e., its skeleton and
edge orientations can be uniquely derived from the joint distribution
$\mathbb{P}(\mathbf{O})$._

To prove Proposition 3.1, we first introduce the following lemma, which
studies the property of the path between a vertex at $t$ and a vertex at
$t+\tau$ ($\tau\geq 0$).

######  Lemma B.1.

Suppose that there is a path between $A(t),B(t+\tau)$ ($\tau\geq 0$) denoted
as $p=A(t),V_{1},...,V_{r},B(t+\tau)$, if $\exists\,V_{i}$ on $p$ such that
$time(V_{i})\geq t+\tau$, then there is a collider on $p$ at $time\geq
t+\tau$.

###### Proof.

We discuss two possible cases and show that $p$ always contains a collider at
$time\geq t+\tau$.

  1. 1.

$\forall\,V_{i},time(V_{i})\leq t+\tau$. In this case, there is a vertex
$V_{j}$ on $p$ such that $time(V_{j})=t+\tau$. Consider its adjacent vertices
$V_{j-1},V_{j+1}$ on the path. According to our assumption, we have
$time(V_{j-1})=time(V_{j+1})=t+\tau-1$. Hence, according to the causal
precedence assumption, the edges among $V_{j-1},V_{j},V_{j+1}$ are
$V_{j-1}(t+\tau-1)\to V_{j}(t+\tau)\leftarrow V_{j+1}(t+\tau+1)$, which means
$V_{j}(t+\tau)$ is a collider on $p$ at $time=t+\tau$.

  2. 2.

$\exists\,V_{i},time(V_{i})>t+\tau$. Let $j:=\arg\max_{i}time(V_{i})$ and
$t^{*}:=time(V_{j})$. We have $t^{*}>t+\tau$. Consider the vertices
$V_{j-1},V_{j},V_{j+1}$, the edges among them are $V_{j-1}(t^{*}-1)\to
V_{j}(t^{*})\leftarrow V_{j+1}(t^{*}-1)$. Hence, $V_{j}(t^{*})$ is a collider
on $p$ at $time>t+\tau$.

∎

Equipped with Lemma B.1, we now introduce the proof of Proposition 3.1below:

###### Proof of Proposition 3.1.

According to [22], the skeleton of the MAG is identifiable under Assumption
2.2.

The orientation of edges in the MAG is either directed ($\to$) or bi-directed
($\leftrightarrow$). Since we assume that the cause must precedes the effect,
instantaneous edges in the MAG must be bi-directed. In the following, we show
that lagged edges in the MAG must be directed and cross $k$ time steps, _e.g._
, $A(t)\to B(t+k)$.

Suppose that $A(t_{1})$ and $B(t_{2})$ ($t_{2}>t_{1}$) are adjacent in the
MAG. We first prove that $t_{2}=t_{1}+k$. Prove by contradiction. Since
$A(t_{1})$ and $B(t_{2})$ are adjacent in the MAG, there is an inducing path
$p=A(t_{1}),V_{1},...,V_{r},B(t_{2})$ between them in the full time DAG.
Suppose that $t_{2}\neq t_{1}+k$, in other words, $t_{2}\geq t_{1}+2k$. Since
the order of the SVAR process is $1$, there exist $V_{i}$ on $p$ such that
$time(V_{i})=t_{1}+k,time(V_{i+1})=t_{1}+k+1$. Hence, $V_{i}(t_{1}+k)$ is
neither latent nor a collider777$V_{i}(t_{1}+k)\to V_{i+1}(t_{1}+k+1)$ means
$V_{i}(t_{1}+k)$ can not be a collider.. As a result, the path $p$ is not an
inducing path, which is a contradiction.

Next, we prove that the edge between $A(t),B(t+k)$ is directed. We prove this
by showing that $A(t)\leftrightarrow B(t+k)$ can not be true. Prove by
contradiction. Suppose that we have $A(t)\leftrightarrow B(t+k)$ in the MAG,
then there is an inducing path $p=A(t),V_{1},...,V_{r},B(t+k)$ between
$A(t),B(t+k)$, and $A(t)\not\in\mathbf{An}(B(t+k))$888$B(t+k)$ can not be the
ancestor of $A(t)$ according to the causal precedence assumption. in the full
time DAG. There are three kinds of possible inducing paths:

  1. 1.

$\exists\,V_{i}$ such that $time(V_{i})\geq t+k$;

  2. 2.

$\forall\,V_{i}$, $time(V_{i})<t+k$ and $time(V_{i})>t$;

  3. 3.

$\forall\,V_{i}$, $time(V_{i})<t+k$ and $\exists\,V_{j}$ such that
$time(V_{j})\leq t$.

In the following, we will show neither of these paths exist, which is a
contradiction and shows that the edge between $A(t),B(t+k)$ is directed.

  1. 1.

$\exists\,V_{i}$ such that $time(V_{i})\geq t+k$. In this case, according to
Lemma B.1, there is a collider on $p$ at $time\geq t+k$. Since the assumed
causal precedence, the collider can not be the ancestor of $A(t)$ or
$B(t+k+1)$. Hence, the path $p$ is not an inducing path.

  2. 2.

$\forall\,V_{i}$, $time(V_{i})<t+k$ and $time(V_{i})>t$. In other words, we
have $time(V_{i})=t+1,..,t+k-1$. We show that in this case, we have
$A(t)\in\mathbf{An}(B(t+k))$ in the full time DAG, which contradicts with
$A(t)\leftrightarrow B(t+k)$ in the MAG.

Specifically, when $k=2$, the path is $p=A(t)\to V_{1}(t+1)\to B(t+2)$. Hence,
$A(t)\in\mathbf{An}(B(t+k))$. When $k\geq 3$, consider the first
black$\leftarrow$ edge along the path999If all edges among the path are $\to$,
then $A(t)\in\mathbf{An}(B(t+k))$., we have $A(t)\small\to\\!\cdots\\!\to
V_{i-1}(t_{i}+1)$ black$\leftarrow$ $V_{i}(t_{i})$ and $i\geq 3$. Therefore,
$V_{i-1}(t_{i}+1)$ is the descendent of $A(t)$ and a collider on $p$.
According to the definition of the inducing path, $V_{i-1}(t_{i}+1)$ is the
ancestor of $B(t+k)$, which means $A(t)\in\mathbf{An}(B(t+k))$.

  3. 3.

$\forall\,V_{i}$, $time(V_{i})<t+k$ and $\exists\,V_{j}$ such that
$time(V_{j})\leq t$. In this case, we further consider two possibilities,

    1. a.

$\forall\,V_{i},t\leq time(V_{i})<t+k$. In this case, there is a vertex
$V_{j}$ on $p$ such that $time(V_{j})=t$ and we have $V_{j-1}(t+1)\leftarrow
V_{j}(t)\to V_{j+1}(t+1)$. Therefore, $V_{j}(t)$ is neither a latent variable
nor a collider, which means $p$ is not an inducing path.

    2. b.

$\forall\,V_{i},time(V_{i})<t+k$ and $\exists\,V_{i}$ such that
$time(V_{i})<t$. In this case, we have $p=A(t),...,V_{i}(t-\tau),...,B(t+k)$
($\tau>0$). As a result,there is a vertex $V_{j}$ on $p$ between
$V_{i}(t-\tau)$ and $B(t+k)$ such that $time(V_{j})=t$ and $time(V_{j+1}=t+1$.
Therefore, $V_{j}(t)$ is neither a latent variable nor a collider, which means
$p$ is not an inducing path.

To conclude, we have shown that if there is an inducing path between
$A(t),B(t+k)$, then $A(t)\in\mathbf{An}(B(t+k))$, which means all lagged edges
in the MAG are directed ones. ∎

###  B.2 Proof of Proposition 3.3: MAG to summary DAG

Proposition 3.3._If there are $A(t)\to B(t+k)$ and $A(t+k)\leftrightarrow
B(t+k)$ in the MAG, then, in the summary DAG, there is either_

  1. 1.

$A\to B$_,_

  2. 2.

_a directed path from_ $A$ to $B$ with length $l\leq k-2$, or

  3. 3.

_a directed path from_ $A$ to $B$ with length $l=k-1$ and a confounding
structure between them with lengths $(r\leq k-2,q\leq k-2)$.

###### Proof.

  1. 1.

We first show that If there is $A\to B$ in the summary DAG, then, there are
$A(t)\to B(t+k)$ and $A(t+k)\leftrightarrow B(t+k)$ in the MAG:

Suppose that there is $A\to B$ in the summary DAG. Then, there is a directed
inducing path $p_{1}=A(t)\to A(t+1)\small\to\\!\cdots\\!\to A(t+k-1)\to
B(t+k)$ from $A(t)$ to $B(t+k)$ in the full time DAG. Thus, we have $A(t)\to
B(t+k)$ in the MAG.

Besides, there is an inducing path $p_{2}=A(t+k)\leftarrow A(t+k-1)\to B(t+k)$
in the full time DAG, which means we have $A(t+k)\leftrightarrow B(t+k)$ in
the MAG.

  2. 2.

We then show that if there is $A(t)\to B(t+k)$ in the MAG, then, there is a
directed path from $A$ to $B$ with length $0\leq l\leq k-1$ (a directed path
with $l=0$ means $A\to B$) in the summary DAG:

Suppose that there is $A(t)\to B(t+k)$ in the MAG. Then, there is a directed
path $A(t)\to V_{1}(t+1)\small\to\\!\cdots\\!\to V_{k-1}(t+k-1)\to B(t+k)$ in
the full time DAG. Denote the number of $V_{i}$ such that $V_{i}=A$ as $r$,
then, we have $0\leq r\leq k-1$. Hence, in the summary DAG, there is a
directed path from $A$ to $B$ with length $l=k-r-1$ and we have $0\leq l\leq
k-1$.

  3. 3.

We finally show that if there is $A(t+k)\leftrightarrow B(t+k)$ in the MAG,
then, in the summary DAG, at least one of the following structures exists:

    1. a.

A directed path101010The path can be from $A$ to $B$ or from $B$ to $A$. $p$
between $A$ and $B$ of length $0\leq l\leq k-2$ (a directed path with $l=0$
means $A\to B$);

    2. b.

A confounding structure $c$ between $A$ and $B$ of length $(r\leq k-2,q\leq
k-2)$.

Suppose that there is $A(t+k)\leftrightarrow B(t+k)$ in the MAG, we will show
that there is a latent confounder $U$ between $A(t+k)$ and $B(t+k)$ such that
$t+1\leq time(U)\leq t+k-1$.

In this regard, if $U=A$ or $U=B$, there is a directed path between $A$ and
$B$ with length $0\leq l\leq k-2$; Otherwise when $U\neq A$ and $U\neq B$,
there is a confounding structure between $A$ and $B$ with length $(r\leq
k-2,q\leq k-2)$.

In the following, we prove that there is a latent confounder $U$ between
$A(t+k)$ and $B(t+k)$ such that $t+1\leq time(U)\leq t+k-1$. Prove by
contradiction. Suppose that $time(U)\leq t$ or $time(U)\geq t+k$. If
$time(U)\leq t$, then on the directed inducing path from $U$ to $A$, there is
a vertex at time $t$, which means the path can not be an inducing path. If
$time(U)\geq t+k$, then there is a collider on the inducing path between $U$
and $A$ according to Lemma B.1, which means the path is not directed and $U$
is not a latent confounder. Hence, the time of the latent confounder $U$ is
between $t+1$ and $t+k-1$.

To conclude, combining the results in 2. and 3., we prove the proposition. ∎

###  B.3 Proof of Theorem 3.5: Identifiability of the summary DAG

For two vertex $A,B$, let the vertex set $\mathbf{M}$ contain $A$ and any
$M_{i}\neq B$ such that $A(t)\to M_{i}(t+k)$ in the MAG and $M_{i}$ is not
$B$’s descendant111111This can be justified from the MAG. Specifically, if
$M_{i}$ is the descendant of $B$, then there is $B(t)\to
V_{1}(t+k),V_{1}(t)\to V_{2}(t+k),...,V_{l}(t)\to M_{i}(t+k)$ in the MAG.. Let
the vertex set $\mathbf{S}$ contains any $S_{i}\neq A$ such that $S_{i}(t)\to
B(t+k)$ or $S_{i}(t)\to M_{j}(t+k)$ for some $M_{j}\in\mathbf{M}$ in the MAG.

Theorem 3.5. _Assuming model (1), Assumption 2.2, and Assumptions
2.6,2.7,2.11, then the summary DAG is identifiable. Specifically,_

  1. 1.

_There is_ $A\to B$ in the summary DAG iff there are $A(t)\to
B(t+k),A(t+k)\leftrightarrow B(t+k)$ in the MAG, and the set
$\mathbf{M}(t+1)\cup\mathbf{S}(t)$ is not sufficient to d-separate
$A(t),B(t+k)$ in the full time DAG.

  2. 2.

_The condition “the set_ $\mathbf{M}(t+1)\cup\mathbf{S}(t)$ is not sufficient
to d-separate $A(t),B(t+k)$ in the full time DAG” can be tested by the proxy
variable $\mathbf{M}(t+k)$ of the unobserved set $\mathbf{M}(t+1)$.

###### Proof.

  1. 1.

$\Rightarrow$. Suppose that there is $A\to B$ in the summary DAG. According to
Proposition 3.3, we have $A(t)\to B(t+k)$ and $A(t+k)\leftrightarrow B(t+k)$
in the MAG. In addition, in the full time DAG, there is a directed path
$p=A(t)\to B(t+1)\small\to\\!\cdots\\!\to B(t+k)$ from $A(t)$ to $B(t+k)$,
which is not _d_ -separated by $\mathbf{M}(t+1)\cup\mathbf{S}(t)$ (because
$B\not\in\mathbf{M}$).

$\Leftarrow$. Suppose that there are $A(t)\to B(t+k)$ and
$A(t+k)\leftrightarrow B(t+k)$ in the MAG. According to Proposition 3.3, at
least one of the following structures exists:

    1. a.

$A\to B$;

    2. b.

A directed path $p_{AB}$ from $A$ to $B$ of length $0<l\leq k-2$;

    3. c.

A directed path $p_{AB}$ from $A$ to $B$ of length $l=k-1$ and a confounding
structure $c_{AB}$ between $A$ and $B$ of length $(r\leq k-2,q\leq k-2)$.

In the following, we show that if the set $\mathbf{M}(t+1)\cup\mathbf{S}(t)$
is not sufficient to _d_ -separate $A(t),B(t+k)$ in the full time DAG, $A\to
B$ must exist in the summary DAG. We prove this by its contrapositive
statement, _i.e._ , if $A\to B$ does not exist, the set
$\mathbf{M}(t+1)\cup\mathbf{S}(t)$ is sufficient to _d_ -separate $A(t)$ and
$B(t+k)$.

Consider the path between $A(t)$ and $B(t+k)$ in the full time DAG
$p=A(t),V_{1},...,V_{l},B(t+k)$. There are three possible cases:

    1. a.

$\exists\,V_{i},time(V_{i})\geq t+k$;

    2. b.

$\forall\,V_{i},time(V_{i})<t+k$ and $\exists\,V_{j}$ such that
$time(V_{j})\leq t$;

    3. c.

$\forall\,V_{i},t<time(V_{i})<t+k$.

Next, we will show that if $A\to B$ does not exist, the set
$\mathbf{M}(t+1)\cup\mathbf{S}(t)$ can _d_ -separate all these three kinds of
paths.

    1. a.

For the path that $\exists\,V_{i},time(V_{i})\geq t+k$. According to Lemma
B.1, there is a collider at $time\geq t+k$. Because the collider and its
descendants are not in $\mathbf{M}(t+1)\cup\mathbf{S}(t)$, the path is _d_
-separated.

    2. b.

For the path that $\forall\,V_{i},time(V_{i})<t+k$ and $\exists\,V_{j}$ such
that $time(V_{j})\leq t$. There is a vertex $V_{j}$ on the path such that
$V_{j}(t)\to V_{j+1}(t+1)\to V_{j+2}(t+2)\to\\!\cdots$121212Strictly, when
$k=2$, $V_{j+2}$ is $B(t+2)$. If the path from $V_{j}(t)$ to $B(t+k)$ is a
directed path, we have $V_{j}\in\mathbf{S}$ and the path is _d_ -separated.
Otherwise, there is a collider on the path at $time\geq t+2$. Neither the
collider nor its descendants are in $\mathbf{M}(t+1)\cup\mathbf{S}(t)$, the
path is also _d_ -separated.

    3. c.

For the path that $\forall\,V_{i},t<time(V_{i})<t+k$. We have $A(t)\to
V_{1}(t+1)\to V_{2}(t+2)\to\\!\cdots$. If the path from $V_{1}(t+1)$ to
$B(t+k)$ is a directed one, $V_{1}\in\mathbf{M}$131313$A(t)\to V_{1}(t+1)$
means we have $A(t)\to V_{1}(t+k)$ in the MAG. $V_{1}$ is the ancestor of $B$,
so $V_{1}\not\in\mathbf{De}(B)$. Hence, $V_{i}\in\mathbf{M}$.. Otherwise,
there is a collider on the path at $time\geq t+2$. Neither the collider nor
its descendants are in $\mathbf{M}(t+1)\cup\mathbf{S}(t)$, the path is also
_d_ -separated.

  2. 2.

In the following, we show that the condition “the set
$\mathbf{M}(t+1)\cup\mathbf{S}(t)$ is not sufficient to _d_ -separate
$A(t),B(t+k)$ in the full time DAG” can be justified by testing
$A(t)\perp\\!\\!\\!\perp B(t+k)|\mathbf{M}(t+1)\cup\mathbf{S}(t)$ with proxy
variables.

    1. a.

We first explain that $A(t),B(t+k),\mathbf{M}(t+1),\mathbf{M}(t+k)$ have the
causal graph shown in Fig. 8.

Figure 8: Causal graph over $A(t),B(t+k),\mathbf{M}(t+1),\mathbf{M}(t+k)$.

       * •

$A(t)\to\mathbf{M}(t+1)$. By definition, $A$ is the ancestor of vertices in
$\mathbf{M}$.

       * •

$\mathbf{M}(t+1)\to\mathbf{M}(t+k)$. By Assumption 2.10the self causation edge
always exists.

       * •

$\mathbf{M}(t+1)\to B(t+k)$. By definition, $B$ is not the ancestor of any
vertex in $\mathbf{M}$.

       * •

$\mathbf{M}(t+1)\leftrightarrow B(t+1)$. According to Proposition 3.1, the
instantaneous edge must be bi-directed.

Indeed, the directed edge $\mathbf{M}(t+1)\to B(t+k)$ and the bi-directed edge
$\mathbf{M}(t+k)\leftrightarrow B(t+k)$ may not exist. However, according to
the theory of single proxy causal discovery, as long as we have
$A(t)\perp\\!\\!\\!\perp\mathbf{M}(t+k)|\mathbf{M}(t+1),\mathbf{S}(t)$ (which
will be proved in the following), we can test whether $A(t)\to B(t+k)$ with
the proxy variable $\mathbf{M}(t+k)$.

    2. b.

We next show that
$A(t)\perp\\!\\!\\!\perp\mathbf{M}(t+k)|\mathbf{M}(t+1),\mathbf{S}(t)$, _i.e._
, $\mathbf{M}(t+k)$ can act as a legal proxy variable of $\mathbf{M}(t+1)$. We
prove this by showing that all path between $A(t)$ to any
$Z_{i}(t+k)\in\mathbf{M}(t+k)$ can be _d_ -separated by
$\mathbf{M}(t+1)\cup\mathbf{S}(t+1)$. Specifically, we consider two different
cases:

       * •

When $A\not\to B$ in the summary DAG. To prove this conclusion, again,
consider three kinds of paths $A(t),V_{1},...,V_{l},Z_{i}(t+k)$ between $A(t)$
and $Z_{i}(t+k)$,

         1. i.

For the path that $\exists\,V_{i},time(V_{i})\geq t+k$. It must contain a
collider at $time\geq t+k$. Hence, neither the collider nor its descendant is
in $\mathbf{M}(t+1)\cup\mathbf{S}(t)$ and the path is _d_ -separated.

         2. ii.

For the path that $\forall\,V_{i},time(V_{i})<t+k$ and $\exists\,V_{j}$ such
that $time(V_{j})\leq t$. Then, there is a vertex $V_{j}$ on the path such
that $V_{j}(t)\to V_{j+1}(t+1)\to V_{j+2}(t+2)\to\\!\cdots$. If the path
between $V_{j}(t)$ and $Z_{i}(t+k)$ is a directed one, then we have
$V_{j}(t)\in\mathbf{S}(t)$ and the path is _d_ -separated. Otherwise, there is
a collider on the path between $V_{j+2}(t+2)$ and $Z_{i}(t+k)$. Thus, neither
the collider nor its descendant is in $\mathbf{M}(t+1)\cup\mathbf{S}(t)$ and
the path is also _d_ -separated.

         3. iii.

For the path that $\forall\,V_{i},t<time(V_{i})<t+k$. In this case, the path
must be $A(t)\to V_{1}(t+1)\to V_{2}(t+2)\to\\!\cdots$. If the path between
$V_{1}(t1)$ and $Z_{i}(t+k)$ is a directed one, we have
$V_{1}(t+1)\in\mathbf{M}(t+1)$, because of the following facts: $V_{1}\neq B$,
$A\to V_{1}$, $V_{1}\not\in\mathbf{De}(B)$141414If $V_{1}\in\mathbf{De}(B)$,
since $V_{1}$ is the ancestor of $Z_{i}$, we have $Z_{i}\in\mathbf{De}(B)$,
which contradicts with the definition of the vertex set $\mathbf{M}$.. As a
result, the path is _d_ -separated. Otherwise, when the path between
$V_{1}(t1)$ and $Z_{i}(t+k)$ is not a directed one, there is a collider on the
path at $time\geq t+2$. Therefore, neither the collider nor its descendant is
in $\mathbf{M}(t+1)\cup\mathbf{S}(t)$ and the path is _d_ -separated.

       * •

When $A\to B$ in the summary DAG. In this case, except for the above analysis,
we need to extra show that every path $p=A(t)\to
B(t+1),V_{1},...,V_{l},Z_{i}(t+k)$ can be _d_ -separated by
$\mathbf{M}(t+1)\cup\mathbf{S}(t)$. Specifically, since we define
$Z_{i}\not\in\mathbf{De}(B)$, $p$ can not be a directed path. Therefore,
starting from $B(t+1)$ along the path, at least one of the edges is
$\leftarrow$. As a result, there is a collider on $p$ which is either $B(t+1)$
itself or at $time\geq t+2$. For both cases, neither the collider nor its
descendant is in $\mathbf{M}(t+1)\cup\mathbf{S}(t)$ and the path is _d_
-separated.

∎

##  Appendix C Experiment

###  C.1 Extra results of Section 5.1: Synthetic study

Figure 9: Performance of our method under different graph scales ($d$ denotes
the variable number) and densities ($\mathrm{prob.}$ denotes the edge
probability in the Erdos-Renyi model).

###  C.2 Extra results of Section 5.2: Discovering causal pathways in
Alzheimer’s disease

Figure 10: Summary DAG over eight meta-regions in AD recovered by the NG-EM
[10] baseline. Figure 11: Summary DAG in AD recovered by the NG-EM [10]
baseline. Figure 12: Summary DAG in AD recovered by the Dynotears [7]
baseline. Table 2: Indices for brain regions partition in Fig. 7.

AAL Index | Abbreviation | Full name | AAL Index | Abbreviation | Full name
---|---|---|---|---|---
1 | PL | Precentral L | 2 | PR | Precentral R
3 | FSL | Frontal Sup L | 4 | FSR | Frontal Sup R
5 | FSOL | Frontal Sup Orb L | 6 | FSOR | Frontal Sup Orb R
7 | FML | Frontal Mid L | 8 | FMR | Frontal Mid R
9 | FMOL | Frontal Mid Orb L1 | 10 | FMOR | Frontal Mid Orb R1
11 | FIOL | Frontal Inf Oper L | 12 | FIOR | Frontal Inf Oper R
13 | FITL | Frontal Inf Tri L | 14 | FITR | Frontal Inf Tri R
15 | FIpL | Frontal Inf Orb L | 16 | FIpR | Frontal Inf Orb R
17 | ROL | Rolandic Oper L | 18 | ROR | Rolandic Oper R
19 | SMAL | Supp Motor Area L | 20 | SMAR | Supp Motor Area R
21 | OL | Olfactory L | 22 | OR | Olfactory R
23 | FSML | Frontal Sup Medial L | 24 | FSMR | Frontal Sup Medial R
25 | FMrL | Frontal Mid Orb L2 | 26 | FMrR | Frontal Mid Orb R2
27 | RL | Rectus L | 28 | RR | Rectus R
29 | IL | Insula L | 30 | IR | Insula R
31 | CAL | Cingulum Ant L | 32 | CAR | Cingulum Ant R
33 | CML | Cingulum Mid L | 34 | CMR | Cingulum Mid R
35 | CPL | Cingulum Post L | 36 | CPR | Cingulum Post R
37 | HL | Hippocampus L | 38 | HR | Hippocampus R
39 | ParL | ParaHippocampal L | 40 | ParR | ParaHippocampal R
41 | AL | Amygdala L | 42 | AR | Amygdala R
43 | CaL | Calcarine L | 44 | CaR | Calcarine R
45 | CuL | Cuneus L | 46 | CuR | Cuneus R
47 | LL | Lingual L | 48 | LR | Lingual R
49 | OSL | Occipital Sup L | 50 | OSR | Occipital Sup R
51 | OML | Occipital Mid L | 52 | OMR | Occipital Mid R
53 | OIL | Occipital Inf L | 54 | OIR | Occipital Inf R
55 | FL | Fusiform L | 56 | FR | Fusiform R
57 | PosL | Postcentral L | 58 | PosR | Postcentral R
59 | PSL | Parietal Sup L | 60 | PSR | Parietal Sup R
61 | PIL | Parietal Inf L | 62 | PIR | Parietal Inf R
63 | SL | SupraMarginal L | 64 | SR | SupraMarginal R
65 | AnL | Angular L | 66 | AnR | Angular R
67 | PreL | Precuneus L | 68 | PreR | Precuneus R
69 | PLL | Paracentral Lobule L | 70 | PLR | Paracentral Lobule R
71 | CL | Caudate L | 72 | CR | Caudate R
73 | PutL | Putamen L | 74 | PutR | Putamen R
75 | PalL | Pallidum L | 76 | PalR | Pallidum R
77 | TL | Thalamus L | 78 | TR | Thalamus R
79 | HesL | Heschl L | 80 | HesR | Heschl R
81 | TSL | Temporal Sup L | 82 | TSR | Temporal Sup R
83 | TPSL | Temporal Pole Sup L | 84 | TPSR | Temporal Pole Sup R
85 | TML | Temporal Mid L | 86 | TMR | Temporal Mid R
87 | TPML | Temporal Pole Mid L | 88 | TPMR | Temporal Pole Mid R
89 | TIL | Temporal Inf L | 90 | TIR | Temporal Inf R