王芳娟

(浙江師范大學(xué)數(shù)理與信息工程學(xué)院，浙江金華 321004)

0 引 言

格蘭杰因果關(guān)系在經(jīng)濟(jì)、生物、計算神經(jīng)科學(xué)等領(lǐng)域已有廣泛的應(yīng)用.其中，對于二維的時間序列，可以直接用兩變量格蘭杰因果關(guān)系分析數(shù)據(jù)之間的內(nèi)在聯(lián)系.但對于多維的時間序列，用兩變量格蘭杰因果關(guān)系已不能正確判斷變量之間的關(guān)系是直接的還是間接的，因此提出了條件格蘭杰因果關(guān)系［1］.在實(shí)際應(yīng)用中，若把外部公共輸入和隱變量的影響因素考慮在內(nèi)，則要用到偏相關(guān)格蘭杰因果關(guān)系才能真正刻畫變量之間的內(nèi)在聯(lián)系［1］.

格蘭杰因果關(guān)系是一種很常用的分析方法，在研究變量之間的因果關(guān)系時，通常情況下要求變量的時間序列的長度遠(yuǎn)遠(yuǎn)大于變量個數(shù)［1-2］，此時，格蘭杰因果關(guān)系法效果顯著.但在實(shí)際應(yīng)用中，可以提取到的時間序列長度是有限的，對于數(shù)以萬計的蛋白質(zhì)、基因變量，每個變量可以測量到的時間序列長度相當(dāng)短.此時，變量維數(shù)比時間序列長度要大，若仍然用通常的格蘭杰因果關(guān)系法研究變量之間的內(nèi)在因果關(guān)系，就會得到錯誤的結(jié)論.為了解決這個問題，將提出一種新的方法，在兩變量格蘭杰因果關(guān)系、條件格蘭杰因果關(guān)系和偏相關(guān)格蘭杰因果關(guān)系的基礎(chǔ)上，把簡單的格蘭杰因果關(guān)系應(yīng)用到復(fù)雜網(wǎng)絡(luò)中.從模擬數(shù)據(jù)的結(jié)果可以看出:該方法能比較有效地揭示復(fù)雜網(wǎng)絡(luò)中變量之間真實(shí)的內(nèi)在聯(lián)系.

1 預(yù)備知識

1.1 線性偏相關(guān)格蘭杰因果關(guān)系

兩變量格蘭杰因果關(guān)系和條件格蘭杰因果關(guān)系廣泛應(yīng)用于實(shí)際數(shù)據(jù)處理過程中.然而，這2種方法都沒有考慮到外部輸入和隱變量的影響.如果在實(shí)際分析時不將這些因素考慮在內(nèi)，那么得到的網(wǎng)絡(luò)結(jié)構(gòu)可能是不正確的.為了更好地處理實(shí)際數(shù)據(jù)中出現(xiàn)的外部輸入和隱變量的影響，提出了線性偏相關(guān)格蘭杰因果關(guān)系［1］.

考慮2個時間序列Xt和Zt，它們的聯(lián)合自回歸表示為

其中:εit為白噪聲;為外部輸入;Bi(L為隱變量;i=1，2.

令ui(t)=εit++Bi(L)，i=1，2，則誤差項協(xié)方差矩陣為

若考慮3個變量Xt，Yt和Zt，它們之間的聯(lián)合自回歸模型為

對S，若不考慮u2(t)=的影響，則u1t的方差可以表示為

對Σ，類似地可以考慮去除u5(t)=ε5t+εE5t+B5(L)εL5t的影響，則u3t的方差可以表示為

定義1 Z條件下，Y對X的偏相關(guān)格蘭杰因果關(guān)系為

1.2 Bootstrap 方法

Bootstrap法類似重復(fù)抽樣法，先對原始的數(shù)據(jù)擬合向量自回歸模型［3］，利用最小二乘估計法［3］計算其系數(shù)矩陣A和誤差項方差E;再利用已知條件，構(gòu)造系數(shù)矩陣為A、誤差項均值為0、方差為E的多元正態(tài)分布白噪聲模型，共模擬M(M=500)次，每次計算相應(yīng)的因果關(guān)系或偏相關(guān)因果值.其中，3倍的σ(均方差)作為置信區(qū)間，若置信區(qū)間下限小于0，則表示因果關(guān)系不存在，否則，存在因果關(guān)系.

2 格蘭杰因果關(guān)系在復(fù)雜網(wǎng)絡(luò)中的應(yīng)用

一般情況下，格蘭杰因果關(guān)系主要用于研究簡單的網(wǎng)絡(luò)結(jié)構(gòu).在計算過程中，用最小二乘法估計多變量向量自回歸模型(MVAR)的參數(shù)時，在理論上有一個約束條件，即每個時間序列長度要大于模型中估計的參數(shù)個數(shù)，對于一個一階的模型，時間序列長度要大于變量個數(shù)的平方［2］.但在實(shí)際應(yīng)用中，對于數(shù)以萬計的蛋白質(zhì)、基因變量，每個變量可以測量到的時間序列長度相當(dāng)短.當(dāng)變量個數(shù)超過時間序列長度時，就會出現(xiàn)“維數(shù)災(zāi)難”，通常的格蘭杰因果關(guān)系就會失效.為了解決這個問題，提出了一種利用格蘭杰因果關(guān)系法進(jìn)行迭代計算的復(fù)雜網(wǎng)絡(luò)研究方法.

2.1 有關(guān)概念

在有向網(wǎng)絡(luò)圖中，將對目標(biāo)節(jié)點(diǎn)T有直接或間接作用的所有節(jié)點(diǎn)統(tǒng)稱為節(jié)點(diǎn)T的源節(jié)點(diǎn);并將所有對目標(biāo)節(jié)點(diǎn)T有直接作用的節(jié)點(diǎn)統(tǒng)稱為節(jié)點(diǎn)T的父節(jié)點(diǎn).對任何一個給定的目標(biāo)節(jié)點(diǎn)T，若可以找到所有對它有直接作用的節(jié)點(diǎn)(即父節(jié)點(diǎn))，則可構(gòu)造出整個網(wǎng)絡(luò).因此，對于給定一個目標(biāo)節(jié)點(diǎn)T，希望找到它所有的父節(jié)點(diǎn).圖1(a)表示一個包含目標(biāo)節(jié)點(diǎn)T和它所有源節(jié)點(diǎn)的網(wǎng)絡(luò).假設(shè){X1，X2，…，Xn}中的每個點(diǎn)到目標(biāo)節(jié)點(diǎn)T只有1條路徑，{Y1，Y2，…，Yn}中的每個點(diǎn)到目標(biāo)節(jié)點(diǎn)T有2條不同的路徑.圖1(b)中通過1條路徑連接的間接點(diǎn)都被刪除了，剩下的間接點(diǎn)Y1，Y2，…，Yn通過2條不同路徑連接.圖1(c)中通過1條路徑和2條不同路徑連接的間接點(diǎn)都被刪除了，只剩下目標(biāo)節(jié)點(diǎn)T的父節(jié)點(diǎn)T1，T2，T3.

圖1 目標(biāo)節(jié)點(diǎn)T的父節(jié)點(diǎn)的篩選步驟圖

2.2 具體算法

首先，應(yīng)用兩變量格蘭杰因果關(guān)系，尋找目標(biāo)節(jié)點(diǎn)T的所有源節(jié)點(diǎn)，用點(diǎn)集A0(T)表示.在圖1(a)中，A0(T)={T1，T2，T3，X1，X2，…，Xn，Y1，Y2，…，Yn}.

其次，進(jìn)一步檢測第1步中找到的連接是直接的還是間接的.執(zhí)行如下迭代過程［4］:

1)對 A0(T)中的每個節(jié)點(diǎn) a∈{T1，T2，T3，X1，X2，…，Xn，Y1，Y2，…，Yn}，計算所有 a→T|b(b∈A0(T)，b≠a)的偏相關(guān)格蘭杰因果關(guān)系值.若經(jīng)計算存在某個b∈A0(T)，b≠a，使得a→T|b的偏相關(guān)格蘭杰因果關(guān)系值用Bootstrap法(M=500)計算得到的3δ置信區(qū)間下限小于0，則表示a→T的連接是間接的，它們之間不存在直接連接，把a(bǔ)從A0(T)中刪除.重復(fù)上述步驟，刪除所有通過1條路徑連接的間接點(diǎn).圖1(a)中，{X1，X2，…，Xn}是通過單一路徑連接的間接點(diǎn)集，將其從A0(T)中刪除，剩下的點(diǎn)構(gòu)成一個新的點(diǎn)集，記為 A1(T)={T1，T2，T3，Y1，Y2，…，Yn}.

2)對 A1(T)中的每個節(jié)點(diǎn) a∈{T1，T2，T3，Y1，Y2，…，Yn}，按上述方法計算所有 a→T|b，c(b，c∈A1(T)，c≠a;b≠c)的偏相關(guān)格蘭杰因果關(guān)系值，其中，b和c的組合共C2n+2種.從而刪除通過2條路徑連接的間接點(diǎn).圖1(b)中，{Y1，Y2，…，Yn}是通過2條不同路徑連接的間接點(diǎn)集，將其從A1(T)中刪除，剩下的新點(diǎn)集為 A2(T)={T1，T2，T3}.

3)繼續(xù)上述步驟，直到不能再從點(diǎn)集Ak(T)中移去任何點(diǎn)為止［4］.

2.3 基本原理

若計算得到的兩變量格蘭杰因果關(guān)系值FY→X很大，但在第3個變量Z條件下，偏相關(guān)格蘭杰因果關(guān)系值FY→X|Z顯著降低，則表示Y→X的連接是間接的，應(yīng)刪除.利用這個原理，可以尋找每個節(jié)點(diǎn)的父節(jié)點(diǎn).

第1階段，先用兩變量格蘭杰因果關(guān)系法，找出使得Y→X的因果值比較大的所有節(jié)點(diǎn)Y，構(gòu)成候選集A0.

第2階段，過濾點(diǎn)集A0，保留Y∈A0，滿足在任意一個剩余變量Z∈A0，Z≠Y條件下的偏相關(guān)格蘭杰因果關(guān)系值FY→X|Z都很大，得到剩余點(diǎn)集A1.

第3階段，利用新點(diǎn)集A1，再次過濾點(diǎn)集A1，保留Y∈A1，滿足在任意2個剩余變量Z1，Z2∈A0(Z1，Z2≠Y;Z1≠Z2)條件下的偏相關(guān)格蘭杰因果關(guān)系值FY→X|Z1，Z2都很大，得到剩余點(diǎn)集A2.

第4階段，再利用點(diǎn)集A2，以其中任意3個不同剩余變量為條件，繼續(xù)上述步驟，直到這種迭代進(jìn)行不下去為止.

通過上述迭代步驟，將格蘭杰因果關(guān)系法應(yīng)用到復(fù)雜網(wǎng)絡(luò)中，使得每個節(jié)點(diǎn)對應(yīng)一列父節(jié)點(diǎn)，就可以得到復(fù)雜網(wǎng)絡(luò)的因果關(guān)系圖.

3 數(shù)值模擬

下面驗證本文方法在復(fù)雜網(wǎng)絡(luò)中應(yīng)用的有效性(能否用格蘭杰因果關(guān)系法正確刻畫長度有限的多變量之間的因果關(guān)系).首先，利用一個模擬的例子(由12個時間序列組成)來進(jìn)一步說明.其次，利用例子擬合的多維數(shù)據(jù)，將本文的新方法和簡單網(wǎng)絡(luò)中的偏相關(guān)格蘭杰因果關(guān)系法的計算結(jié)果進(jìn)行比較.

例1 假設(shè)12個變量的時間序列由如下的動力方程確定:

其中，εi(i=1，2，…，12)是相互獨(dú)立、均值為0、方差為1的白噪聲.由這個模型直接可得12個變量之間的因果關(guān)系，如圖2(b)所示，共有18條連接邊.取模型中每個變量的200個時間點(diǎn)作為擬合數(shù)據(jù)，所得數(shù)據(jù)軌跡如圖2(a)所示.為了增加可視性，將X1，X2，…，X12進(jìn)行了平移，依次從下到上排列.

圖2 數(shù)據(jù)的軌跡及相應(yīng)的網(wǎng)絡(luò)結(jié)構(gòu)

利用上述方法研究該擬合數(shù)據(jù)的網(wǎng)絡(luò)結(jié)構(gòu).先對數(shù)據(jù)計算兩變量格蘭杰因果關(guān)系值，用Bootstrap法得到的連接圖如圖3(a)所示，共有95條連接邊，比真實(shí)的網(wǎng)絡(luò)多出很多邊，原因是兩變量格蘭杰因果關(guān)系不能區(qū)分連接是直接的還是間接的，在其中混雜了許多間接連接的邊.

接著，按上述迭代步驟，依次計算1個、2個變量條件下的偏相關(guān)格蘭杰因果關(guān)系值，分別用于除去通過1條、2條路徑連接的間接點(diǎn)，用Bootstrap法后所得的結(jié)果如圖3(b)和(c)所示，分別有19條和18條邊.它們和圖3(a)相比已經(jīng)去掉了很大一部分間接連接邊，比較接近真實(shí)結(jié)構(gòu).說明該模擬數(shù)據(jù)利用上述迭代步驟可以比較正確地判別出網(wǎng)絡(luò)結(jié)構(gòu).

最后，作為比較，對12個變量兩兩之間直接計算偏相關(guān)格蘭杰因果關(guān)系值(對研究的2個變量，以其他剩余的所有變量為條件)，用Bootstrap法后所得的網(wǎng)絡(luò)結(jié)構(gòu)如圖3(d)所示，與真實(shí)的網(wǎng)絡(luò)結(jié)構(gòu)相比缺少一些連接邊.該方法所得結(jié)果和真實(shí)結(jié)構(gòu)相比誤差比較大.這是由于模型的階為4，在模型中要估計的參數(shù)總數(shù)是12×12×4=576，而已有的時間序列長度只有200，不足以用它來估計所有參數(shù)，導(dǎo)致計算得到錯誤的因果連接.

圖3(A)，(B)，(C)，(D)表示網(wǎng)絡(luò)圖 3(a)，(b)，(c)，(d)的關(guān)系圖.若某種關(guān)系(行指標(biāo)表示原因，列指標(biāo)表示結(jié)果)經(jīng)Bootstrap法得到的置信區(qū)間下限小于0，則為黑色，表示變量之間沒有因果關(guān)系;若下限大于0，則為白色，表示變量之間有因果關(guān)系存在.

圖3 不同方法所得的網(wǎng)絡(luò)結(jié)構(gòu)圖

4 結(jié) 語

在實(shí)際應(yīng)用中，研究的網(wǎng)絡(luò)通常是多變量復(fù)雜網(wǎng)絡(luò)，此時直接應(yīng)用偏相關(guān)格蘭杰因果關(guān)系得到的結(jié)果往往是不穩(wěn)定的，與真實(shí)的網(wǎng)絡(luò)結(jié)構(gòu)存在很大偏差.為了解決“維數(shù)災(zāi)難”這個難題，把格蘭杰因果關(guān)系法應(yīng)用到復(fù)雜網(wǎng)絡(luò)中，基于兩變量格蘭杰因果關(guān)系、條件格蘭杰因果關(guān)系、偏相關(guān)格蘭杰因果關(guān)系，用迭代法一步步除去通過1條單一路徑、2條不同路徑等的間接連接，直至除去所有間接連接.最后剩下每個節(jié)點(diǎn)的父節(jié)點(diǎn)，從而構(gòu)造出相應(yīng)復(fù)雜網(wǎng)絡(luò)的結(jié)構(gòu).格蘭杰因果關(guān)系在復(fù)雜網(wǎng)絡(luò)中的迭代應(yīng)用是對通常格蘭杰因果關(guān)系的優(yōu)化，彌補(bǔ)了其對維數(shù)限制的不足.該方法對今后利用格蘭杰因果關(guān)系法研究其他各種復(fù)雜網(wǎng)絡(luò)結(jié)構(gòu)有很大的幫助.

［1］Guo Shuixia，Seth A K，Kendriek K M，et al.Partial granger causality-eliminating exogenous inputs and latent variables［J］.Journal of Neuroscience Methods，2008，172(1):79-93.

［2］Gopikrishna D，Stephan L C，George A J，et al.Multivariate granger causality analysis of fMRI data［J］.Human Brain Mapping，2009，30(4):1361-1373.

［3］何書元.應(yīng)用時間序列分析［M］.北京:北京大學(xué)出版社，2003:300-305.

［4］Zou Cunlu，Christophe L，Guo Shuixia，et al.Identifying interactions in the time and frequency domains in local and global networks:A granger causality approach［J］.BMC Bioinformatics，2010，11:337-341.