1、文献信息

复旦大学林伟教授：新型因果网络辨识算法可以用于疫情防控

近日，复旦大学数学科学学院、类脑智能科学与技术研究院的林伟教授团队，与中国科学院、苏州大学、日本东京大学等团队合作，提出了数据驱动的因果网络辨识的新型算法。该方法可以用于大规模复杂动力系统内蕴因果网络的复现，有助于解析实际系统演化的本质机制和规律。该研究成果不久前以《偏交叉映射排除间接因果影响》为题在线发表于综合类学术期刊《自然-通讯》

Leng S , Ma H F , Kurths J , et al. Partial cross mapping eliminates indirect causal influences[J]. Nature Communications, 2020, 11:2632.

2、文献简介

由于因果传递的影响，因果检测可能会将间接因果误判为直接因果。尽管在传统框架上提出了一些方法来避免这种误判，任然缺少从间接因果中识别出直接因果关系的方法，特别是当变量隐藏在动力系统中且是相互不分离的，相互之间有弱关系这样情况会具有更大挑战。这里，我们给出一种基于数据的与模型无关的偏交叉映射的方法来解决此问题，该方法基于非线性动力学和数据的三种工具的整合：相空间重构，相互交叉映射，和偏相关。我们使用来自不同模型和真实系统的数据验证了我们的方法。因为直接因果是很多复杂动力学基础的关键，我们预计我们的方法在解锁和解密不同学科真实系统的内部机制方面必不可少。

3、研究方法与结果

3.1 研究模型

3.1.1 不可分离变量概念

我们引入了不可分离变量概念，它是通过一个连续的时间动力系统得到： $\dot{x}=F(x)$ ，其中状态变量 $\mathbf{x}(t)=\left[x_{1}(t), x_{2}(t), \ldots, x_{n}(t)\right]^{\top}$ 在紧流形中演化，形成维度为 $d_\mathcal{A}$ 的吸引子 $\mathcal{A}$ 。这里可以将 $d_\mathcal{A}$ 计算为 $\mathcal{A}$ 的盒维数。初始值为 $x_0\in \mathcal{M_x}$ 的动力学由 $x(t)=\varphi_t(x_0)$ 得到，其中 $\varphi_t(\cdot)$ 是流形 $\mathcal{M}_x$ 上的流。根据Takens –Mañé's的嵌入理论以及它的份行概述，可以以概率重构出具有一个正延迟 $\tau$ 和平滑观测函数 $h:\mathcal{M}_x\to \mathbb{R}$ 的系统，只要 $L>2d_{\mathcal{A}}$ ，延迟坐标图 $\Gamma_{h, \varphi, \tau}(\mathbf{x})=\left[h(\mathbf{x}), h\left(\varphi_{-\tau}(\mathbf{x})\right), h\left(\varphi_{-2 \tau}(\mathbf{x})\right), \ldots, h\left(\varphi_{-(L-1) \tau}(\mathbf{x})\right)\right]^{\top}$ 通常是一个嵌入图。为了直接说明，我们取观测函数为一个简单的坐标函数 $h(\mathbf{x})=x_i$ ，其中 $x_i$ 为 $\mathbf{x}$ 的第i个分量。因此，我们有 $\mathbf{y}(t)=[x_i(t),x_i(t-\tau),...,x_i(t-(L-1)\tau)]^\top$ 并且通过嵌入图 $\Gamma$ 将流形 $\mathcal{M}_x$ 映射到阴影流形 $\mathcal{M}_y$ 。因为嵌入图是一对一的，因此阴影流形 $\mathcal{M}_y$ 上的动力学 $\psi_\tau$ 和 $\mathcal{M}_x$ 上的动力学 $\varphi_\tau$ 在拓扑上是共轭的，也即：

$\mathbf{y}(t+\tau)=\psi_{\tau}(\mathbf{y}(t))=\Gamma \circ \varphi_{\tau} \circ \Gamma^{-1}(\mathbf{y}(t))$

一方面，系统 $\dot{x}=F(x)$ 隐含着这样一个事实，即某个特定组件（例如 $xj,j =(or≠)i$ ）的未来动力学受以下因素支配：

$\left[\varphi_{\tau}\right]_{j}: \mathbf{x}(t)=\left[x_{1}(t), x_{2}(t), \ldots, x_{n}(t)\right]^{\top} \rightarrow x_{j}(t+\tau)$

因此，取决于所有组件 $x_1,x_2,…,x_n$ 的历史值。另一方面， $\mathbf{y}(t+\tau)$ 式中的关系暗示了另一个事实，即只要存在嵌入图 $\Gamma$ ， $x_j$ 的未来动力学也受下式影响：

$\left[\Gamma^{-1} \circ \psi_{\tau}\right]_{j}: \mathbf{y}(t)=\left[x_{i}(t), x_{i}(t-\tau), \ldots, x_{i}(t-(L-1) \tau)\right]^{\top} \rightarrow x_{j}(t+\tau)$

因此仅取决于变量 $x_i$ 的历史和嵌入图 $\Gamma$ 。

一般的，仅通过观察一个变量来预测 $x_j(t+\tau)$ 是可能的，且这种预测与使用系统的所有变量信息 $x_1(t),x_2(t),...,x_n(t)$ 的预测一样完美。因此，Takens–Mañé的嵌入理论揭示了，在这种非线性的系统中，整个动力系统的信息可以一般的注入到一个变量中，因此可以通过该变量的感测数据进行重构。因此，这引起了不可分离的概念，也即对动力系统进行任何预测时，通常不能从其它变量中删除某些变量的信息。这也表明基于预测框架的方法(如Granger因果，传递熵及其扩展)从数学上不适合处理非线性系统产生的时间序列数据，而非线性系统之间始终存在不可分性。简例见附件。

3.1.2 传递引起间接因果关系

为了说明传递引起间接因果关系，我们考虑一个3物种的启发logistic模型有如下联系：

$\begin{array}{l}x_{t}=x_{t-1}\left(\alpha_{x}-\alpha_{x} x_{t-1}\right) \\z_{t}=z_{t-1}\left(\alpha_{z}-\alpha_{z} z_{t-1}-\beta_{z x} x_{t-1}\right) \\y_{t}=y_{t-1}\left(\alpha_{y}-\alpha_{y} y_{t-1}-\beta_{y z} z_{t-1}\right)\end{array}$

其中这三个物种在因果链中相互作用，表示为 $X\to Z\to Y$ ，且耦合强度 $\beta_{zx}$ 和 $\beta_{yz}$ 是非零的。现在，我们将上述模型的第二个方程移动一个时间步长，然后代入模型的最后一个方程，得到:

$y_{t}=y_{t-1}\left[\alpha_{y}-\alpha_{y} y_{t-1}-\beta_{y z} z_{t-2}\left(\alpha_{z}-\alpha_{z} z_{t-2}-\beta_{z x} x_{t-2}\right)\right]$

同样最后一个方程可变化为：

$\begin{array}{l}z_{t-1}=\frac{1}{\beta_{y z}}\left(\alpha_{y}-\alpha_{y} y_{t-1}-y_{t} / y_{t-1}\right) \\z_{t-2}=\frac{1}{\beta_{y z}}\left(\alpha_{y}-\alpha_{y} y_{t-2}-y_{t-1} / y_{t-2}\right)\end{array}$

因此有：

$\begin{aligned}y_{t}=& y_{t-1}\left\{\alpha_{y}-\alpha_{y} y_{t-1}-\beta_{y z} \frac{1}{\beta_{y z}}\left(\alpha_{y}-\alpha_{y} y_{t-2}-y_{t-1} / y_{t-2}\right)\right.\\&\left.\times\left[\alpha_{z}-\alpha_{z} \frac{1}{\beta_{y z}}\left(\alpha_{y}-\alpha_{y} y_{t-2}-y_{t-1} / y_{t-2}\right)-\beta_{z x} x_{t-2}\right]\right\}\end{aligned}$

因此这个方程和模型的第一个方程，形成了从X到Y的单向因果关系。但是，这种因果关系是间接的，是由传递性引起的，且这个印象对离散动力学系统具有时间延迟的影响。

3.1.3 一阶和高阶PCM方法

现在我们正式指定PCM框架（见附件图1）.第一步是将时间序列 $Y=\{y_t\}$ 以时间步长 $\tau_i(i=1,2,...,m)$ 转换生成m个变量 $Y_{\tau_i}=\{y_{t+\tau_i}\}$ 。对于时间序列对 $Y_{\tau_i}$ 和 $Z$ ，我们应用MCM方法从 $Y_{\tau_i}$ 获得 $\hat{Z}^{Y_{\tau_i}}$ ，并计算其相关系数 $Corr(Z,\hat{Z}^{Y_{\tau_i}})$ 。为了简化，我们用 $\hat{Z}^Y$ 表示 $\hat{Z}^{Y_{\tau_{{i_1}}}}$ ，其中 $i_{1}=\operatorname{argmax}_{1 \leq i \leq m} \operatorname{Corr}\left(Z, \hat{Z}^{Y_{\tau_{i}}}\right)$ 。下一步重复上述过程得到从 $\hat{Z}^{Y}_{\tau_i}$ 到 $\hat{X}^{\hat{Z}_{\tau_i}^Y}$ 的映射，并用 $\hat{X}^{\hat{Z}^Y}$ 表示 $\hat{X}^{\hat{Z}_{\tau_2}^Y}$ ，其中 $i_{2}=\operatorname{argmax}_{1 \leq i \leq m} \operatorname{Corr}\left(X, \hat{X}^{\hat{Z}_{i_{i}}^{\mathrm{Y}}}\right)$ 。现在得到的 $\hat{X}^{\hat{Z}^Y}$ 表示间接信息流。通过直接应用MCM模型到 $Y_{\tau_i}$ 和 $X$ ，我们可以用 $\hat{X}^Y$ 表示从X到Y的全部信息，这是 $\hat{X}^{Y_{\tau_{i_3}}}$ 的简化，其中 $i_{3}=\operatorname{argmax}_{1 \leq i \leq m} \operatorname{Corr}\left(X, \hat{X}^{Y_{\tau_{i}}}\right)$ 。我们引入相关性指标 $\mathcal{Q}_{\mathrm{D}}=\left|\operatorname{Pcc}\left(X, \hat{X}^{Y} \mid \hat{X}^{\hat{Z}^{Y}}\right)\right|$ ，其中 $Pcc(\cdot,\cdot|\cdot)$ 是一个描述剔除第三个变量信息的前两个变量联系程度偏相关系数。我们来回顾下偏相关系数的定义，对于时间序列X,Y以及 $Z^1,...,Z^s$ ，那么在 $Z^1$ 的条件下X和Y之间的相关系数为：

$\operatorname{Pcc}\left(X, Y \mid Z^{1}\right)=\frac{\operatorname{Corr}(X, Y)-\operatorname{Corr}\left(X, Z^{1}\right) \operatorname{Corr}\left(Y, Z^{1}\right)}{\sqrt{\left(1-\operatorname{Corr}\left(X, Z^{1}\right)^{2}\right)\left(1-\operatorname{Corr}\left(Y, Z^{1}\right)^{2}\right)}}$

在 $Z^1$ 和 $Z^2$ 的条件下X和Y的偏相关系数为：

$\operatorname{Pcc}\left(X, Y \mid Z^{1}, Z^{2}\right)=\frac{\operatorname{Pcc}\left(X, Y \mid Z^{1}\right)-\operatorname{Pcc}\left(X, Z^{2} \mid Z^{1}\right) \operatorname{Pcc}\left(Y, Z^{2} \mid Z^{1}\right)}{\sqrt{\left(1-\operatorname{Pcc}\left(X, Z^{2} \mid Z^{1}\right)^{2}\right)\left(1-\operatorname{Pcc}\left(Y, Z^{2} \mid Z^{1}\right)^{2}\right)}}$

可以递归定义在更多条件下的偏相关系数。为了提供我们方法的详细说明，我们总结了实际步骤：

过程A：MCM检测从 $U=\{u_t\}^L_{t=1}$ 到 $V=\{v_t\}^L_{t=1}$ 的因果关系。（1）通过使用U和V时间序列的延迟坐标嵌入来重构相空间，可以通过使用FNN算法和DMI方法选择重构参数(嵌入维度 $E_u$ , $E_v$ ,时滞 $\tau_u$ , $\tau_v$ ；见附件Note5)；（2）对于每个时间索引t，找到 $\mathbf{v}_t$ 的邻居节点集合 $\mathcal{N}(\mathbf{v}_t)$ ( $E_v+1$ 个最近邻居被使用，因为它是 $E_v$ 维空间中有界单纯性所需的最小数目)；（3）在 $M_U$ 中找出哪些与 $\mathcal{N(\mathbf{v}_t)}$ 具有相同时间索引的相应点，并计算它们的加权平均得到估计值 $\hat{\mathbf{u}}^\mathbf{v}_t$ （这里的权重取决于 $\mathcal{N}(\mathbf{v}_t)$ 中的节点与 $\mathbf{v}_t$ 间的距离，定义为操作 $\mathbb{E}[\cdot]$ ）；（4）使用一个合适的指标（如 $\mathcal{Q}_\mathrm{C}=|Corr(\mathbf{u}_t,\hat{\mathbf{u}}^\mathbf{v}_t)|$ ）来表征时间序列 $\hat{U}^{V_0}$ 的估计（下表0表示这里没有V的变换，下文含义相同）且初始时间序列 $U$ ，它衡量从 $U$ 到 $V$ 的因果。

过程B：PCM方法检测在条件Z下从X到Y的直接因果。（1）选择不同的时间延迟 $\tau_i(i=1,2,...,m)$ 将时间序列Y 转化生成 $Y_{\tau_i}=\{y_{t+\tau_i}\}$ ；（2）对于每对 $Z$ 到 $Y_{\tau_i}$ ，执行过程A获得 $Corr(Z,\hat{Z}^{Y_{\tau_i}})$ ，且用 $\hat{Z}^Y$ 表示 $\hat{Z}^{Y_{\tau_{i_1}}}$ ，其中 $\tau_{i_1}$ 可使 $Corr(Z,\hat{Z}^{Y_{\tau_i}})$ 达到最大；（3）改变时间序列 $\hat{Z}^Y$ 的时间延迟 $\tau_i(i=1,2,...,m)$ 得到 $\hat{Z}^Y_{\tau_i}$ ；（4）对于每对 $X$ 到 $\hat{Z}^Y_{\tau_i}$ ，执行过程A得到 $Corr(X,\hat{X}^{\hat{Z}_{\tau_i}^Y})$ ，且用 $\hat{X}^{\hat{Z}^Y}$ 表示 $\hat{X}^{\hat{Z}_{\tau_{i_2}}^Y}$ ，其中时间延迟 $\tau_{i_2}$ 使 $Corr(X,\hat{X}^{\hat{Z}^Y_{\tau_i}})$ 达到最大；（5）对于每对 $X$ 到 $Y_{\tau_i}$ ，执行过程A的到 $Corr(X,\hat{X}^{Y_{\tau_i}})$ ，且用 $\hat{X}^Y$ 表示 $\hat{X}^Y_{\tau_{i_3}}$ ，其中时间延迟 $\tau_{i_3}$ 使得 $Corr(X,\hat{X}^{Y_{\tau_i}})$ 达到最大；（6）使用 $\mathcal{Q}_\mathrm{D}=\left|\operatorname{Pcc}\left(X, \hat{X}^{Y} \mid \hat{X}^{\hat{Z}^{Y}}\right)\right|$ 来衡量在条件Z下X到Y的直接因果关系。

需要注意的是上述MCM过程中，我们搜索的是不同候选时间延迟下的强因果关系。为了一致性，在整个研究中，所有MCM结果都基于该策略。此外，可以在时间延迟的分布（即因果谱）上表征变量之间的因果关系。这种完整的因果关系将包含在我们未来的工作中。

正如上述描述，一阶PCM模型可以按下述定义估计超过3个相互作用变量，高阶模型可写为：

$\mathcal{Q}_{\mathrm{D}_{1}}=\left|\operatorname{Pcc}\left(X, \hat{X}^{Y}\left|\left\{\hat{X}^{\hat{Z}^{i^{Y}}} \mid i=1, \ldots, s\right\}\right)\right.\right|$

在一个复杂的动力学网络中，直接因果也能够通过不止一个变量传递（如 $X\to Z^1\to Z^2\to Y$ ）。高阶PCM方法被用于这种特殊情况。特别的，我们通过删除从s个变量（ $Z^1,...,Z^s$ ）中任意两个变量的交叉映射变量的信息，可以计算 $X$ 和 $\hat{X}^Y$ 相关系数和它们间的偏相关系数：

$\mathcal{Q}_{\mathrm{D}_{2}}=\left|\operatorname{Pcc}\left(X, \hat{X}^{Y}\left|\left\{\hat{X}^{\hat{Z}^{i^{\hat{Z}^{jY}}}} \left| i \neq j, i, j \in\{1, \ldots, s\}\right\}\right) \right| \right.\right.$

该式是区分直接因果的一种有效方法，表示从X到Y通过两个变量传递的间接因果。类似的，可以通过定义 $\mathcal{Q}_{\mathrm{D}_n}$ 来表示 $Z^1,...,Z^s$ 中任意n个变量的组合传递因果的情况。

$mathcal{Q}_{\mathrm{D}_{n}}=\left|\operatorname{Pcc}\left(X, \hat{X}^{Y}\left|\left\{\hat{X}^{\hat{Z}^{i_{1}} \cdot 2^{i_{n} Y}} \mid\left(i_{1}, \ldots, i_{n}\right) \text { is an } n-\text { combination from }\{1, \ldots, s\}\right\}\right) \right.\right|$

结合 $\mathcal{Q}_\mathrm{C}$ , $\mathcal{Q}_{\mathrm{D}_n}$ （n=1,...,s）和PCM方法得到 $\gamma=(\prod^s_{n=1}\mathcal{Q}_{\mathrm{D}_n})/\mathcal{Q}_{\mathrm{C}}^s$ ，用其反映这些系数的接近程度，我们得到了高阶PCM方法，来检测大型网络的因果联系。然而，对于较大的n阶，n个中间变量的可能组合非常大，我们将在未来的工作中研究高阶方法的计算和应用，本文只考虑一阶问题。

在实践中，如果网络规模相对较大，则部分相关过程将遇到计算问题，因此应考虑使用较大的条件集。在这种情况下，我们要选择一些节点 $Z^i$ 来最大化 $\mathcal{Q}_\mathrm{C}^{X\to Z^i}+\mathcal{Q}_\mathrm{C}^{ Z^i\to Y}$ ，这就意味着很大概率上存在通过 $Z^i$ 的间接因果，且以这些点为条件。此外，如果我们的先验知识表明网络是稀疏的，也即，间接因果是很少的，我们也可以一个个以 $Z^1,...,Z^s$ 为条件，取 $\mathcal{Q}_\mathrm{D}^{X\to Y|Z^i}$ 的最小值作为最终结果。

此外，通过将偏相关性替换为表征条件依赖性的其他可能的度量，可以进一步发展或改变PCM思想。例如，确定相关系数（表示为 $r^2$ ）是一种可能的选择，用作从交叉图邻居通过分配每个影响因素的效应大小直接估计一个指标。另一种启发性的思路是对于间接因果影响 $X\to Z\to Y$ ，切断 $X\to Y$ 或 $Z\to Y$ 都可以削除整个间接信息流，这也提供了PCM框架的变化。这些进一步的变化将包含在未来工作中。

3.2 研究结果

图1：直接因果联系和间接因果联系

3.2.1 Partial cross mapping

为了方便描述偏交叉映射的方法，我们考虑三个变量在单向链中因果关系的简单情况。令 $X=\left\{x_{t}\right\}_{t=1}^{L}, Y=\left\{y_{t}\right\}_{t=1}^{L}$ $, Z=\left\{z_{t}\right\}_{t=1}^{L}$ 是长度为L的时间序列。使用Takens-Mañé的delay-cooordinate embedding，我们可以得到三个shadow manifolds： $X=\left\{x_{t}\right\}_{t=r}^{L}, Y=\left\{y_{t}\right\}_{t=r}^{L}$ $, Z=\left\{z_{t}\right\}_{t=r}^{L}$ 向量表示为：

$\begin{aligned}\mathbf{x}_{t} &=\left(x_{t}, x_{t-\tau_{x}}, \ldots, x_{t-\left(E_{x}-1\right) \tau_{x}}\right) \\\mathbf{y}_{t} &=\left(y_{t}, y_{t-\tau_{y}}, \cdots, y_{t-\left(E_{y}-1\right) \tau_{y}}\right) \\\mathbf{z}_{t} &=\left(z_{t}, z_{t-\tau_{z}}, \cdots, z_{t-\left(E_{z}-1\right) \tau_{z}}\right)\end{aligned}$

其中 $E_x,E_y,E_z$ 分别是嵌入大小， $\tau _x,\tau_y,\tau_z$ 是时差， $r=\max _{\xi=x, y, z}\left\{1+\left(E_{\xi}-1\right) \tau_{\xi}\right\}$ 。这些嵌入大小和时差参数可以分别通过false nearest neighbor(FNN)和delayed mutual information(DMI)计算得到。一般的，对于任意一对变量 $\xi$ 和 $\eta \in \{x,y,z\}$ ，我们令 $\hat{\mathcal{N}}^{\xi}\left(\boldsymbol{\eta}_{t}\right)=\left\{\boldsymbol{\eta}_{t^{\prime}} \mid \xi_{t^{\prime}} \in \mathcal{N}\left(\xi_{t}\right)\right\}$ ，其中 $\mathcal{N}(\xi_{t})$ 是一个包含在corresponding shadow manifold中 $\xi_t$ 的最邻近节点的固定值（通常取 $E_\xi +1$ ，这是在 $E_\xi$ 维空间中有界单纯性所需点的最小值）。对于 $\xi=\eta$ ， $\hat{\mathcal{N}}^{\xi}\left(\boldsymbol{\eta}_{t}\right)={\mathcal{N}}\left(\boldsymbol{\eta}_{t}\right)$ 。当 $\xi\ne\eta$ ， $\hat{\mathcal{N}}^{\xi}\left(\boldsymbol{\eta}_{t}\right)$ 变为 ${\mathcal{N}}\left(\boldsymbol{\eta}_{t}\right)$ 邻居的交叉映射。从 ${\mathcal{N}}\left(\boldsymbol{\eta}_{t}\right)$ 到 $\hat{\mathcal{N}}^{\xi}\left(\boldsymbol{\eta}_{t}\right)$ 的独立性表征了从变量产生 $\eta_t$ 到变量产生 $\xi_t$ 的因果影响。先前用于量化这种独立性和因果影响的启发式方法构成了MCM框架。我们利用 $\eta_t$ 和 $\hat{\eta_t}^\xi=\mathbb{E}[\hat{\mathcal{N}}^\xi(\eta_t)]$ 之间的相关系数，其中 $\hat{\eta_t}^\xi$ 是 $\xi_t$ 的一个映射， $\mathbb{E}[\cdot]$ 是一个对给定集合中所有的点取近似加权平均的操作。特别的，如果相关系数 $\mathcal{Q}_C=|Corr(x_t,\hat{x}_t^y)|$ 大于一个经验阈值T，这个MCM方法将认为从X到Y存在因果关系。MCM对成对不可分离系统中的因果关系领域做了补充。然而，由于因果关系的传递性，MCM得到的因果连接可以是直接的，也可以是间接的（如图2a所示）。另外，因为因果关系会在一个特定的时间延迟后显示其影响，我们寻找一个最优的时间延迟来极大化X和Y的因果关系（如，相关系数 $\mathcal{Q}_C$ ）。

图2：PCM框架的基本原理

$\mathcal{Q}_C$ 按照上述定义，表示整个空间中 $X$ 和 $\hat{X}^Y$ 夹角的余弦值，如图2b所示。为了区分因果传递的存在，我们考虑将 $\mathcal{Q}_C$ 投影到与因果传递性引起的间接信息正交的信息空间上。为此，我们给出我们的PCM框架。首先，对于一个时间序列 $Z$ 和平移量 $Y_{\tau_i}=\{y_{t+\tau_i}\}$ (候选时间延迟为 $\tau_i(i=1,2,...,m)$ )，我们应用传统的MCM方法去决定最优的时间延迟 $\tau_i=\tau_{i_1}$ ，就是要最大化相关系数 $Corr(Z,\hat{Z}^{Y_{\tau_i}})$ 。相应的，从 $Y_{\tau_{i_1}}$ 得到的映射 $\hat{Z}^{Y_{\tau_{i_1}}}$ 可用 $\hat{Z}^Y$ 简化表示。接下来重复这个过程，产生时间序列对 $X$ 和translated $\hat{Z}^Y_{\tau_i}$ 以获得最优的延迟 $\tau_{i_2}$ ，同样从 $\hat{Z}^Y_{\tau_{i_2}}$ 得到的映射 $\hat{X}^{\hat{Z}^Y_{\tau_{i_2}}}$ ，这样最大化了相关系数 $Corr(X,\hat{X}^{\hat{Z}^Y_{\tau_i}})$ 。定义获得的映射为 $\hat{X}^{\hat{Z}^\hat{Y}}$ ，它是从一个连续的MCM过程获得的，且表征了一个流经Z的间接信息流，然后就得到了 $\hat{X}^Y$ ，它表征了所有从X到Y的因果信息，通过对时间序列对X和translated $Y_{\tau_i}$ 重复上述过程，我们就引入了相关指标： $\mathcal{Q}_{\mathrm{D}}=\left|\operatorname{Pcc}\left(X, \hat{X}^{Y} \mid \hat{X}^{\hat{Z}^{Y}}\right)\right|$ 来测量从X到Y的直接因果关系(通过Z的间接因果为条件)，其中 $Pcc(\cdot,\cdot|\cdot)$ 是一个偏相关系数用来描述删除第三个变量信息后前两个变量相关性的程度，与MCM指标 $\mathcal{Q}_\mathrm{C}=|Corr(X,\hat{X}^Y)|$ 不同。需要注意的是我们搜索的是在每个MCM过程中的不同候选时间延迟中的最强的因果。因此 $\mathcal{Q}_{\mathrm{D}}$ 可以被直观的认为是 $\mathcal{Q}_\mathrm{C}$ 在与直接信息 $\hat{X}^{\hat{Z}^Y}$ 正交的信息空间上的投影（见图2b），因此消除了间接因果的影响。

对于这三个因果变量X,Y和Z，我们有 $\mathcal{Q}_\mathrm{C}\ge \mathcal{Q}_\mathrm{D}$ 。设定一个经验阈值 $1>T\gg0$ ，这些相关指标有三种大小情况： $\mathcal{Q}_\mathrm{C}\ge \mathcal{Q}_\mathrm{D}\ge T$ ， $\mathcal{Q}_\mathrm{C}\ge T \gg \mathcal{Q}_\mathrm{D}$ ， $T>\mathcal{Q}_\mathrm{C}\ge \mathcal{Q}_\mathrm{D}$ ，相应的，这三种因果关系分别为：从X到Y的直接因果，从X到Y的间接因果，从X到Y没有因果关系。指标 $\mathcal{Q}_\mathrm{D}$ 表征削除间接联系后直接因果联系的程度。对于图2a中的例子，这里X和Y的因果关系属于第二种情况，可以从其相关指标推断得到。在真实的应用中，可能会发生因果信号不够强，导致 $\mathcal{Q}_\mathrm{C}\gtrsim T$ , $\mathcal{Q}_\mathrm{D}$ 趋于T的情况。在这种情况下，直接因果的检测对T的值更加敏感。为了克服这个问题，我们引入 $\gamma =\mathcal{Q}_\mathrm{D}/\mathcal{Q}_\mathrm{C}$ 来衡量两个指标的接近度。越接近于1，则越是一种直接因果联系。进行了多次测试保证了统计上的可靠性。

这个PCM框架可以拓展到具有任意多交互变量的网络系统中， $X,Y,Z^1,...,Z^s(s\ge2)$ （例如图1d）。利用 $X$ 和 $\hat{X}^Y$ 所有相关性，我们可以计算它们的偏相关系数为 $\mathcal{Q}_{\mathrm{D}_{1}}=\left| \operatorname{Pcc}\left(X, \hat{X}^{Y} \left\{\hat{X}^{\hat{Z}^{i Y}} \mid i=1, \ldots, s\right\}\right) \right|$ ，通过移除这s个变量 $Z^1,...,Z^s$ 的交叉映射变量信息，其中 $\mathcal{Q}_{\mathrm{D}_{1}}$ 是区分X到Y直接和间接因果连接联系的一阶度量。我们这里强调非线性系统中的强耦合（同步）变量不再PCM框架内，因为在这种环境下，完成的系统将崩溃为一些因果子系统的sub-manifold，并且这种影响变量将称为一个在因果系统上的可观测函数，其中这种双向的因果将会从计算中发现。另外，我们的PCM框架是基于Takens–Mañé理论的，它仅能应用到自治系统。从非自治系统种得到的数据不能直接用到我们的框架种，但我们的方法可以应用到一些非自治系统中。特别的，它可以用交换系统的数据从数值上被用来检测分段因果关系，其中这个转换点可以被定位，且每个连续转换点的持续时间足够长。同样，我们的框架也适用于一些强制系统以及一些噪声较弱或中等的系统，因为一些广义的嵌入定理可以支持我们框架的健全性。对于一种重要的非自治系统，即具有随时间变化的耦合函数或/和各种噪声的动态振荡器，通过贝叶斯动态推理和一组精细的函数库可以提供非常实用的解决方案。对于未来的调查主题，可能的调查包括结合上述互补方法结合起来，用于更一般的动力学系统的因果检测，而无需了解明确的模型方程式，但具有高度复杂的交互结构。

3.2.2 确定基准系统中的直接因果关系

为了能够验证我们的方法，我们使用下述三个相互作用的基本系统：

$\begin{array}{l}x_{t}=x_{t-1}\left(\alpha_{x}-\alpha_{x} x_{t-1}-\beta_{x y} y_{t-1}\right)+\epsilon_{x, t} \\y_{t}=y_{t-1}\left(\alpha_{y}-\alpha_{y} y_{t-1}-\beta_{y x} x_{t-1}-\beta_{y z} z_{t-1}\right)+\epsilon_{y, t} \\z_{t}=z_{t-1}\left(\alpha_{z}-\alpha_{z} z_{t-1}-\beta_{z x} x_{t-1}\right)+\epsilon_{z, t}\end{array}$

其中 $\alpha_x=3.6$ , $\alpha_y=3.72$ , $\alpha_z=3.68$ , $\epsilon_{i,t}(i\in \{x,y,z\})$ 是零均值标准差为0.005的白噪声。不同的耦合参数选择 $\beta_{xy},\beta_{yx},\beta_{yz},\beta_{zx}$ 导致了不同的相互作用模型，见图3a。从这个时间序列中，我们可以分别计算出MCM和PCM的指标， $\mathcal{Q}_\mathrm{C}$ 和 $\mathcal{Q}_\mathrm{D}$ ，以检测从X到Y 的直接连边(如图3b和3c)。在某些情况下，两种方法都可以有效的检测直接因果边，但对于阈值T=0.5时的因果链和因果环结构，这个PCM方法可以成功的区分出间接因果关系，然而MCM方法由于不能削除因果传递的影响而不能区分。随着T的变化，PCM方法表现的比MCM方法更加鲁棒，使得PCM方法可以在没有足够的T先验信息的情况下更好的应用到真实系统中。这个结果在图3b和c中也通过多测实验矫正得到了验证。此外，对于其它所有可能的三种相互作用结构，包括代表性的网络结构：fan-in,fan-out,和cascading structures，我们的系统研究表明，PCM方法可以完全准确的实现因果检测。更重要的是，我们的系统还对Granger因果关系，传递熵及其所有条件的拓展进行比较，以检测这三个系统的因果关系以及测试不同噪声水平和时间序列长度的鲁棒性。正如补充材料Note3所示，这个PCM方法表现超出了现有方法，这些方法原则上只适用于变量分离的条件。在补充材料Note3中，我们还给出了PCM框架和动态贝叶斯推断的方法。两种方法都有它们各自的有点且可以互补的适用。所有这些结果系统的证明了我们的方法相较于经典方法的普遍性和特殊性，这里动力系统的变量是不可分离的。

图3：检测基准模型中从X到Y的因果连接

另外，我们在八个相互作用网络中验证了PCM的有效性。见附件的图10，通过选择合适的T组，网路的直接因果边可被重构出来，而成功的削除所有的间接连边。相反，对于相同的T值，MCM方法会产生一个包含直接，间接甚至错误的因果网络。我们可以发现即使阈值参数T的值较小，也可以通过适用比值 $\gamma = \mathcal{Q}_\mathrm{D}/\mathcal{Q}_\mathrm{C}$ 来提升检测准确性(附件Note4)。此外，选择一个实际有效的阈值更容易实现PCM方法的鲁棒性(附件图11和Note5)。在此模型中，及时数据量较小噪声较强，PCM对时间序列长度和噪声标度的鲁棒性测试具有良好的效果(附件Note3)。这些结果证明了我们PCM方法在检测直接因果以及准确的重建因果网络方面具有强大的功能。

3.2.3 在真实网路中检测因果关系

我们在基因表达控制网络中测试了我们的结果。有5个不同合成结构的网络。每个网络具有100个基因。我们使用GeneNetWeaver软件随机选择了20个基因，其中每个基因都有21个基因表达时间序列数据的10个实现。图4a呈现了一种基因规则网络(其它见附件图12)。对于每个基因，我们将所有实现作为一个时间序列来进行相空间重构。我们将PCM检测到的直接因果联系与五个网络的先验已知边缘进行比较，并计算出各自的ROC（接收器工作特性）曲线（图4b）。我们发现ROC曲线下五个区域的平均值接近〜0.75，这表明了在基因规则网络中（甚至数据集很小）具有很高的检测准确性。第二个例子考虑食物链网络，微型蓝藻，轮虫和类胡萝卜素，其关系如图4cd。第三个例子考虑香港空气污染和心血管疾病的住院记录的关系，如图4ef。

图4：在3个真实网络中检测因果关系

4、主要结论

我们的方法主要有两点优势：基于PCM可以检测间接因果；基于Takens–Mañé的嵌入理论处理变量不可分问题。现有方法要么将间接因果关系链接错误识别为直接因果关系，要么由于违反可分离性条件而失败，因此，我们在理论和计算上开发了一种方法来解决此突出问题，以应对这种情况。现有框架无法有效地解决这些问题。我们的PCM方法通过应用于许多实际系统中得到了验证，从而对这些系统的动力学产生了新的见解。明确消除直接因果联系并消除间接因果影响是理解和准确建模基础系统的关键，因此我们的框架提供了实现此目标的工具。

5、后续讨论

。。。

偏交叉映射排除间接因果影响

1、 文献信息