1默赂、文獻信息?

復旦大學林偉教授：新型因果網(wǎng)絡(luò)辨識算法可以用于疫情防控

近日逗物，復旦大學數(shù)學科學學院、類腦智能科學與技術(shù)研究院的林偉教授團隊宴猾，與中國科學院陡舅、蘇州大學抵乓、日本東京大學等團隊合作，提出了數(shù)據(jù)驅(qū)動的因果網(wǎng)絡(luò)辨識的新型算法靶衍。該方法可以用于大規(guī)模復雜動力系統(tǒng)內(nèi)蘊因果網(wǎng)絡(luò)的復現(xiàn)灾炭，有助于解析實際系統(tǒng)演化的本質(zhì)機制和規(guī)律。該研究成果不久前以《偏交叉映射排除間接因果影響》為題在線發(fā)表于綜合類學術(shù)期刊《自然-通訊》

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

2颅眶、文獻簡介

由于因果傳遞的影響蜈出，因果檢測可能會將間接因果誤判為直接因果。盡管在傳統(tǒng)框架上提出了一些方法來避免這種誤判涛酗，任然缺少從間接因果中識別出直接因果關(guān)系的方法铡原，特別是當變量隱藏在動力系統(tǒng)中且是相互不分離的偷厦，相互之間有弱關(guān)系這樣情況會具有更大挑戰(zhàn)。這里燕刻，我們給出一種基于數(shù)據(jù)的與模型無關(guān)的偏交叉映射的方法來解決此問題只泼，該方法基于非線性動力學和數(shù)據(jù)的三種工具的整合：相空間重構(gòu)，相互交叉映射卵洗，和偏相關(guān)请唱。我們使用來自不同模型和真實系統(tǒng)的數(shù)據(jù)驗證了我們的方法。因為直接因果是很多復雜動力學基礎(chǔ)的關(guān)鍵过蹂，我們預計我們的方法在解鎖和解密不同學科真實系統(tǒng)的內(nèi)部機制方面必不可少籍滴。

3、研究方法與結(jié)果

3.1 研究模型

3.1.1 不可分離變量概念

我們引入了不可分離變量概念榴啸，它是通過一個連續(xù)的時間動力系統(tǒng)得到： $\dot{x}=F(x)$ 孽惰，其中狀態(tài)變量 $\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}$ 的盒維數(shù)勋功。初始值為 $x_0\in \mathcal{M_x}$ 的動力學由 $x(t)=\varphi_t(x_0)$ 得到，其中 $\varphi_t(\cdot)$ 是流形 $\mathcal{M}_x$ 上的流库说。根據(jù)Takens –Ma?é's的嵌入理論以及它的份行概述狂鞋，可以以概率重構(gòu)出具有一個正延遲 $\tau$ 和平滑觀測函數(shù) $h:\mathcal{M}_x\to \mathbb{R}$ 的系統(tǒng)，?只要 $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}$ 通常是一個嵌入圖骚揍。為了直接說明，我們?nèi)∮^測函數(shù)為一個簡單的坐標函數(shù) $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))$

一方面锰什，系統(tǒng) $\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)$ 式中的關(guān)系暗示了另一個事實，即只要存在嵌入圖 $\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)$ 是可能的谢谦，且這種預測與使用系統(tǒng)的所有變量信息 $x_1(t),x_2(t),...,x_n(t)$ 的預測一樣完美释牺。因此萝衩，Takens–Ma?é的嵌入理論揭示了，在這種非線性的系統(tǒng)中没咙，整個動力系統(tǒng)的信息可以一般的注入到一個變量中猩谊，因此可以通過該變量的感測數(shù)據(jù)進行重構(gòu)。因此祭刚，這引起了不可分離的概念牌捷，也即對動力系統(tǒng)進行任何預測時，通常不能從其它變量中刪除某些變量的信息涡驮。這也表明基于預測框架的方法(如Granger因果暗甥，傳遞熵及其擴展)從數(shù)學上不適合處理非線性系統(tǒng)產(chǎn)生的時間序列數(shù)據(jù)，而非線性系統(tǒng)之間始終存在不可分性捉捅。簡例見附件撤防。

3.1.2 傳遞引起間接因果關(guān)系

為了說明傳遞引起間接因果關(guān)系，我們考慮一個3物種的啟發(fā)logistic模型有如下聯(lián)系：

$\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}$ 是非零的。現(xiàn)在无牵，我們將上述模型的第二個方程移動一個時間步長漾肮，然后代入模型的最后一個方程，得到:

$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(jié)的單向因果關(guān)系克懊。但是，這種因果關(guān)系是間接的七蜘，是由傳遞性引起的谭溉，且這個印象對離散動力學系統(tǒng)具有時間延遲的影響。

3.1.3 一階和高階PCM方法

現(xiàn)在我們正式指定PCM框架（見附件圖1）.第一步是將時間序列 $Y=\{y_t\}$ 以時間步長 $\tau_i(i=1,2,...,m)$ 轉(zhuǎn)換生成m個變量 $Y_{\tau_i}=\{y_{t+\tau_i}\}$ 崔梗。對于時間序列對 $Y_{\tau_i}$ 和 $Z$ 夜只，我們應用MCM方法從 $Y_{\tau_i}$ 獲得 $\hat{Z}^{Y_{\tau_i}}$ 垒在，并計算其相關(guān)系數(shù) $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)$ ∩■辏現(xiàn)在得到的 $\hat{X}^{\hat{Z}^Y}$ 表示間接信息流。通過直接應用MCM模型到 $Y_{\tau_i}$ 和 $X$ 签舞，我們可以用 $\hat{X}^Y$ 表示從X到Y(jié)的全部信息秕脓，這是 $\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)$ 。我們引入相關(guān)性指標 $\mathcal{Q}_{\mathrm{D}}=\left|\operatorname{Pcc}\left(X, \hat{X}^{Y} \mid \hat{X}^{\hat{Z}^{Y}}\right)\right|$ 吠架，其中 $Pcc(\cdot,\cdot|\cdot)$ 是一個描述剔除第三個變量信息的前兩個變量聯(lián)系程度偏相關(guān)系數(shù)芙贫。我們來回顧下偏相關(guān)系數(shù)的定義，對于時間序列X,Y以及 $Z^1,...,Z^s$ 傍药，那么在 $Z^1$ 的條件下X和Y之間的相關(guān)系數(shù)為：

$\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的偏相關(guān)系數(shù)為：

$\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)}}$

可以遞歸定義在更多條件下的偏相關(guān)系數(shù)磺平。為了提供我們方法的詳細說明，我們總結(jié)了實際步驟：

過程A：MCM檢測從 $U=\{u_t\}^L_{t=1}$ 到 $V=\{v_t\}^L_{t=1}$ 的因果關(guān)系拐辽。（1）通過使用U和V時間序列的延遲坐標嵌入來重構(gòu)相空間拣挪，可以通過使用FNN算法和DMI方法選擇重構(gòu)參數(shù)(嵌入維度 $E_u$ , $E_v$ ,時滯 $\tau_u$ , $\tau_v$ ；見附件Note5)俱诸；（2）對于每個時間索引t菠劝，找到 $\mathbf{v}_t$ 的鄰居節(jié)點集合 $\mathcal{N}(\mathbf{v}_t)$ ( $E_v+1$ 個最近鄰居被使用，因為它是 $E_v$ 維空間中有界單純性所需的最小數(shù)目)睁搭；（3）在 $M_U$ 中找出哪些與 $\mathcal{N(\mathbf{v}_t)}$ 具有相同時間索引的相應點闸英，并計算它們的加權(quán)平均得到估計值 $\hat{\mathbf{u}}^\mathbf{v}_t$ （這里的權(quán)重取決于 $\mathcal{N}(\mathbf{v}_t)$ 中的節(jié)點與 $\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(jié)的直接因果辙喂。（1）選擇不同的時間延遲 $\tau_i(i=1,2,...,m)$ 將時間序列Y 轉(zhuǎn)化生成 $Y_{\tau_i}=\{y_{t+\tau_i}\}$ ；（2）對于每對 $Z$ 到 $Y_{\tau_i}$ 鸠珠，執(zhí)行過程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}$ ，執(zhí)行過程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}$ 可缚，執(zhí)行過程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(jié)的直接因果關(guān)系知给。

需要注意的是上述MCM過程中，我們搜索的是不同候選時間延遲下的強因果關(guān)系。為了一致性涩赢，在整個研究中戈次，所有MCM結(jié)果都基于該策略。此外筒扒，可以在時間延遲的分布（即因果譜）上表征變量之間的因果關(guān)系朝扼。這種完整的因果關(guān)系將包含在我們未來的工作中。

正如上述描述霎肯，一階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|$

在一個復雜的動力學網(wǎng)絡(luò)中，直接因果也能夠通過不止一個變量傳遞（如 $X\to Z^1\to Z^2\to Y$ ）观游。高階PCM方法被用于這種特殊情況搂捧。特別的，我們通過刪除從s個變量（ $Z^1,...,Z^s$ ）中任意兩個變量的交叉映射變量的信息懂缕，可以計算 $X$ 和 $\hat{X}^Y$ 相關(guān)系數(shù)和它們間的偏相關(guān)系數(shù)：

$\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.$

該式是區(qū)分直接因果的一種有效方法允跑，表示從X到Y(jié)通過兩個變量傳遞的間接因果。類似的搪柑，可以通過定義 $\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|$

結(jié)合 $\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$ ，用其反映這些系數(shù)的接近程度工碾，我們得到了高階PCM方法弱睦，來檢測大型網(wǎng)絡(luò)的因果聯(lián)系。然而渊额，對于較大的n階况木，n個中間變量的可能組合非常大，我們將在未來的工作中研究高階方法的計算和應用旬迹，本文只考慮一階問題火惊。

在實踐中，如果網(wǎng)絡(luò)規(guī)模相對較大奔垦，則部分相關(guān)過程將遇到計算問題屹耐，因此應考慮使用較大的條件集。在這種情況下椿猎，我們要選擇一些節(jié)點 $Z^i$ 來最大化 $\mathcal{Q}_\mathrm{C}^{X\to Z^i}+\mathcal{Q}_\mathrm{C}^{ Z^i\to Y}$ 惶岭，這就意味著很大概率上存在通過 $Z^i$ 的間接因果，且以這些點為條件鸵贬。此外俗他，如果我們的先驗知識表明網(wǎng)絡(luò)是稀疏的，也即阔逼，間接因果是很少的，我們也可以一個個以 $Z^1,...,Z^s$ 為條件地沮，取 $\mathcal{Q}_\mathrm{D}^{X\to Y|Z^i}$ 的最小值作為最終結(jié)果嗜浮。

此外羡亩，通過將偏相關(guān)性替換為表征條件依賴性的其他可能的度量，可以進一步發(fā)展或改變PCM思想危融。例如畏铆，確定相關(guān)系數(shù)（表示為 $r^2$ ）是一種可能的選擇，用作從交叉圖鄰居通過分配每個影響因素的效應大小直接估計一個指標吉殃。另一種啟發(fā)性的思路是對于間接因果影響 $X\to Z\to Y$ 辞居，切斷 $X\to Y$ 或 $Z\to Y$ 都可以削除整個間接信息流，這也提供了PCM框架的變化蛋勺。這些進一步的變化將包含在未來工作中瓦灶。

3.2 研究結(jié)果

圖1：直接因果聯(lián)系和間接因果聯(lián)系

3.2.1 Partial cross mapping

為了方便描述偏交叉映射的方法，我們考慮三個變量在單向鏈中因果關(guān)系的簡單情況抱完。令 $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\}$ 。這些嵌入大小和時差參數(shù)可以分別通過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$ 的最鄰近節(jié)點的固定值（通常取 $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)$ 變?yōu)?img class="math-inline" src="https://math.jianshu.com/math?formula=%7B%5Cmathcal%7BN%7D%7D%5Cleft(%5Cboldsymbol%7B%5Ceta%7D_%7Bt%7D%5Cright)" alt="{\mathcal{N}}\left(\boldsymbol{\eta}_{t}\right)" mathimg="1">鄰居的交叉映射逛薇。從 ${\mathcal{N}}\left(\boldsymbol{\eta}_{t}\right)$ 到 $\hat{\mathcal{N}}^{\xi}\left(\boldsymbol{\eta}_{t}\right)$ 的獨立性表征了從變量產(chǎn)生 $\eta_t$ 到變量產(chǎn)生 $\xi_t$ 的因果影響捺疼。先前用于量化這種獨立性和因果影響的啟發(fā)式方法構(gòu)成了MCM框架。我們利用 $\eta_t$ 和 $\hat{\eta_t}^\xi=\mathbb{E}[\hat{\mathcal{N}}^\xi(\eta_t)]$ 之間的相關(guān)系數(shù)永罚，其中 $\hat{\eta_t}^\xi$ 是 $\xi_t$ 的一個映射啤呼， $\mathbb{E}[\cdot]$ 是一個對給定集合中所有的點取近似加權(quán)平均的操作。特別的呢袱，如果相關(guān)系數(shù) $\mathcal{Q}_C=|Corr(x_t,\hat{x}_t^y)|$ 大于一個經(jīng)驗閾值T官扣，這個MCM方法將認為從X到Y(jié)存在因果關(guān)系。MCM對成對不可分離系統(tǒng)中的因果關(guān)系領(lǐng)域做了補充羞福。然而惕蹄，由于因果關(guān)系的傳遞性，MCM得到的因果連接可以是直接的，也可以是間接的（如圖2a所示）卖陵。另外遭顶，因為因果關(guān)系會在一個特定的時間延遲后顯示其影響，我們尋找一個最優(yōu)的時間延遲來極大化X和Y的因果關(guān)系（如泪蔫，相關(guān)系數(shù) $\mathcal{Q}_C$ ）棒旗。

圖2：PCM框架的基本原理

$\mathcal{Q}_C$ 按照上述定義，表示整個空間中 $X$ 和 $\hat{X}^Y$ 夾角的余弦值撩荣，如圖2b所示铣揉。為了區(qū)分因果傳遞的存在，我們考慮將 $\mathcal{Q}_C$ 投影到與因果傳遞性引起的間接信息正交的信息空間上餐曹。為此逛拱，我們給出我們的PCM框架。首先凸主，對于一個時間序列 $Z$ 和平移量 $Y_{\tau_i}=\{y_{t+\tau_i}\}$ (候選時間延遲為 $\tau_i(i=1,2,...,m)$ )，我們應用傳統(tǒng)的MCM方法去決定最優(yōu)的時間延遲 $\tau_i=\tau_{i_1}$ 旁舰，就是要最大化相關(guān)系數(shù) $Corr(Z,\hat{Z}^{Y_{\tau_i}})$ 嗡官。相應的，從 $Y_{\tau_{i_1}}$ 得到的映射 $\hat{Z}^{Y_{\tau_{i_1}}}$ 可用 $\hat{Z}^Y$ 簡化表示磺樱。接下來重復這個過程婆咸，產(chǎn)生時間序列對 $X$ 和translated $\hat{Z}^Y_{\tau_i}$ 以獲得最優(yōu)的延遲 $\tau_{i_2}$ ，同樣從 $\hat{Z}^Y_{\tau_{i_2}}$ 得到的映射 $\hat{X}^{\hat{Z}^Y_{\tau_{i_2}}}$ 块差，這樣最大化了相關(guān)系數(shù) $Corr(X,\hat{X}^{\hat{Z}^Y_{\tau_i}})$ 。定義獲得的映射為 $\hat{X}^{\hat{Z}^\hat{Y}}$ 倔丈，它是從一個連續(xù)的MCM過程獲得的憨闰，且表征了一個流經(jīng)Z的間接信息流，然后就得到了 $\hat{X}^Y$ 需五，它表征了所有從X到Y(jié)的因果信息鹉动，通過對時間序列對X和translated $Y_{\tau_i}$ 重復上述過程，我們就引入了相關(guān)指標： $\mathcal{Q}_{\mathrm{D}}=\left|\operatorname{Pcc}\left(X, \hat{X}^{Y} \mid \hat{X}^{\hat{Z}^{Y}}\right)\right|$ 來測量從X到Y(jié)的直接因果關(guān)系(通過Z的間接因果為條件)宏邮，其中 $Pcc(\cdot,\cdot|\cdot)$ 是一個偏相關(guān)系數(shù)用來描述刪除第三個變量信息后前兩個變量相關(guān)性的程度泽示，與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}$ 变姨。設(shè)定一個經(jīng)驗閾值 $1>T\gg0$ ，這些相關(guān)指標有三種大小情況： $\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}$ ，相應的怒竿，這三種因果關(guān)系分別為：從X到Y(jié)的直接因果砍鸠，從X到Y(jié)的間接因果，從X到Y(jié)沒有因果關(guān)系耕驰。指標 $\mathcal{Q}_\mathrm{D}$ 表征削除間接聯(lián)系后直接因果聯(lián)系的程度爷辱。對于圖2a中的例子，這里X和Y的因果關(guān)系屬于第二種情況媒抠，可以從其相關(guān)指標推斷得到阀趴。在真實的應用中刘急，可能會發(fā)生因果信號不夠強，導致 $\mathcal{Q}_\mathrm{C}\gtrsim T$ , $\mathcal{Q}_\mathrm{D}$ 趨于T的情況攻柠。在這種情況下瑰钮，直接因果的檢測對T的值更加敏感开睡。為了克服這個問題篇恒，我們引入 $\gamma =\mathcal{Q}_\mathrm{D}/\mathcal{Q}_\mathrm{C}$ 來衡量兩個指標的接近度胁艰。越接近于1腾么，則越是一種直接因果聯(lián)系解虱。進行了多次測試保證了統(tǒng)計上的可靠性。

這個PCM框架可以拓展到具有任意多交互變量的網(wǎng)絡(luò)系統(tǒng)中悍汛， $X,Y,Z^1,...,Z^s(s\ge2)$ （例如圖1d）员凝。利用 $X$ 和 $\hat{X}^Y$ 所有相關(guān)性，我們可以計算它們的偏相關(guān)系數(shù)為 $\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}}$ 是區(qū)分X到Y(jié)直接和間接因果連接聯(lián)系的一階度量。我們這里強調(diào)非線性系統(tǒng)中的強耦合（同步）變量不再PCM框架內(nèi)馆铁，因為在這種環(huán)境下埠巨，完成的系統(tǒng)將崩潰為一些因果子系統(tǒng)的sub-manifold望侈，并且這種影響變量將稱為一個在因果系統(tǒng)上的可觀測函數(shù)脱衙，其中這種雙向的因果將會從計算中發(fā)現(xiàn)岂丘。另外，我們的PCM框架是基于Takens–Ma?é理論的仪召，它僅能應用到自治系統(tǒng)。從非自治系統(tǒng)種得到的數(shù)據(jù)不能直接用到我們的框架種召娜，但我們的方法可以應用到一些非自治系統(tǒng)中玖瘸。特別的雅倒，它可以用交換系統(tǒng)的數(shù)據(jù)從數(shù)值上被用來檢測分段因果關(guān)系，其中這個轉(zhuǎn)換點可以被定位裁良，且每個連續(xù)轉(zhuǎn)換點的持續(xù)時間足夠長价脾。同樣灭忠，我們的框架也適用于一些強制系統(tǒng)以及一些噪聲較弱或中等的系統(tǒng)弛作，因為一些廣義的嵌入定理可以支持我們框架的健全性映琳。對于一種重要的非自治系統(tǒng)，即具有隨時間變化的耦合函數(shù)或/和各種噪聲的動態(tài)振蕩器谎脯，通過貝葉斯動態(tài)推理和一組精細的函數(shù)庫可以提供非常實用的解決方案源梭。對于未來的調(diào)查主題，可能的調(diào)查包括結(jié)合上述互補方法結(jié)合起來烛愧，用于更一般的動力學系統(tǒng)的因果檢測怜姿，而無需了解明確的模型方程式社牲，但具有高度復雜的交互結(jié)構(gòu)。

3.2.2 確定基準系統(tǒng)中的直接因果關(guān)系

為了能夠驗證我們的方法熟空，我們使用下述三個相互作用的基本系統(tǒng)：

$\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的白噪聲。不同的耦合參數(shù)選擇 $\beta_{xy},\beta_{yx},\beta_{yz},\beta_{zx}$ 導致了不同的相互作用模型迈喉，見圖3a孩革。從這個時間序列中膝蜈，我們可以分別計算出MCM和PCM的指標饱搏， $\mathcal{Q}_\mathrm{C}$ 和 $\mathcal{Q}_\mathrm{D}$ 推沸，以檢測從X到Y(jié) 的直接連邊(如圖3b和3c)坤学。在某些情況下压怠，兩種方法都可以有效的檢測直接因果邊蜗顽，但對于閾值T=0.5時的因果鏈和因果環(huán)結(jié)構(gòu)雇盖，這個PCM方法可以成功的區(qū)分出間接因果關(guān)系崔挖，然而MCM方法由于不能削除因果傳遞的影響而不能區(qū)分。隨著T的變化脓鹃，PCM方法表現(xiàn)的比MCM方法更加魯棒瘸右，使得PCM方法可以在沒有足夠的T先驗信息的情況下更好的應用到真實系統(tǒng)中逞频。這個結(jié)果在圖3b和c中也通過多測實驗矯正得到了驗證苗胀。此外基协，對于其它所有可能的三種相互作用結(jié)構(gòu)，包括代表性的網(wǎng)絡(luò)結(jié)構(gòu)：fan-in,fan-out,和cascading structures杂穷，我們的系統(tǒng)研究表明耐量，PCM方法可以完全準確的實現(xiàn)因果檢測。更重要的是山叮，我們的系統(tǒng)還對Granger因果關(guān)系屁倔，傳遞熵及其所有條件的拓展進行比較锐借，以檢測這三個系統(tǒng)的因果關(guān)系以及測試不同噪聲水平和時間序列長度的魯棒性。正如補充材料Note3所示嗅战，這個PCM方法表現(xiàn)超出了現(xiàn)有方法驮捍，這些方法原則上只適用于變量分離的條件启具。在補充材料Note3中鲁冯，我們還給出了PCM框架和動態(tài)貝葉斯推斷的方法薯演。兩種方法都有它們各自的有點且可以互補的適用跨扮。所有這些結(jié)果系統(tǒng)的證明了我們的方法相較于經(jīng)典方法的普遍性和特殊性衡创，這里動力系統(tǒng)的變量是不可分離的。

圖3：檢測基準模型中從X到Y(jié)的因果連接

另外拔莱，我們在八個相互作用網(wǎng)絡(luò)中驗證了PCM的有效性。見附件的圖10动看，通過選擇合適的T組菱皆，網(wǎng)路的直接因果邊可被重構(gòu)出來，而成功的削除所有的間接連邊重父。相反方淤，對于相同的T值携茂，MCM方法會產(chǎn)生一個包含直接讳苦，間接甚至錯誤的因果網(wǎng)絡(luò)。我們可以發(fā)現(xiàn)即使閾值參數(shù)T的值較小卿堂，也可以通過適用比值 $\gamma = \mathcal{Q}_\mathrm{D}/\mathcal{Q}_\mathrm{C}$ 來提升檢測準確性(附件Note4)草描。此外，選擇一個實際有效的閾值更容易實現(xiàn)PCM方法的魯棒性(附件圖11和Note5)逛绵。在此模型中术浪，及時數(shù)據(jù)量較小噪聲較強胰苏，PCM對時間序列長度和噪聲標度的魯棒性測試具有良好的效果(附件Note3)硕并。這些結(jié)果證明了我們PCM方法在檢測直接因果以及準確的重建因果網(wǎng)絡(luò)方面具有強大的功能倔毙。

3.2.3 在真實網(wǎng)路中檢測因果關(guān)系

我們在基因表達控制網(wǎng)絡(luò)中測試了我們的結(jié)果。有5個不同合成結(jié)構(gòu)的網(wǎng)絡(luò)毙玻。每個網(wǎng)絡(luò)具有100個基因桑滩。我們使用GeneNetWeaver軟件隨機選擇了20個基因，其中每個基因都有21個基因表達時間序列數(shù)據(jù)的10個實現(xiàn)缭受。圖4a呈現(xiàn)了一種基因規(guī)則網(wǎng)絡(luò)(其它見附件圖12)韭畸。對于每個基因胰丁，我們將所有實現(xiàn)作為一個時間序列來進行相空間重構(gòu)。我們將PCM檢測到的直接因果聯(lián)系與五個網(wǎng)絡(luò)的先驗已知邊緣進行比較甘萧，并計算出各自的ROC（接收器工作特性）曲線（圖4b）。我們發(fā)現(xiàn)ROC曲線下五個區(qū)域的平均值接近?0.75邀泉，這表明了在基因規(guī)則網(wǎng)絡(luò)中（甚至數(shù)據(jù)集很小）具有很高的檢測準確性拔恰。第二個例子考慮食物鏈網(wǎng)絡(luò)，微型藍藻，輪蟲和類胡蘿卜素，其關(guān)系如圖4cd夷恍。第三個例子考慮香港空氣污染和心血管疾病的住院記錄的關(guān)系酿雪，如圖4ef。

圖4：在3個真實網(wǎng)絡(luò)中檢測因果關(guān)系

4、主要結(jié)論

我們的方法主要有兩點優(yōu)勢：基于PCM可以檢測間接因果吓揪；基于Takens–Ma?é的嵌入理論處理變量不可分問題。現(xiàn)有方法要么將間接因果關(guān)系鏈接錯誤識別為直接因果關(guān)系钾腺，要么由于違反可分離性條件而失敗，因此间螟，我們在理論和計算上開發(fā)了一種方法來解決此突出問題厢破，以應對這種情況。現(xiàn)有框架無法有效地解決這些問題见坑。我們的PCM方法通過應用于許多實際系統(tǒng)中得到了驗證嚷掠，從而對這些系統(tǒng)的動力學產(chǎn)生了新的見解。明確消除直接因果聯(lián)系并消除間接因果影響是理解和準確建能衤浚基礎(chǔ)系統(tǒng)的關(guān)鍵不皆，因此我們的框架提供了實現(xiàn)此目標的工具。

5熊楼、后續(xù)討論

霹娄。。孙蒙。

偏交叉映射排除間接因果影響

1默赂、 文獻信息?