基于格蘭杰因果的效應(yīng)性連接分析方法綜述

尚志剛，沈曉陽，李蒙蒙，萬 紅

(1.鄭州大學(xué) 電氣工程學(xué)院，河南 鄭州 450001； 2.鄭州大學(xué) 河南省腦科學(xué)與腦機接口技術(shù)重點實驗室，河南 鄭州 450001)

0 引言

生物的感知與認(rèn)知普遍依賴不同腦區(qū)間的相互作用，深入了解多個腦區(qū)間的作用關(guān)系有助于闡明正常或病理狀態(tài)下特定腦功能的神經(jīng)機制[1]。目前對大腦不同腦區(qū)間交互性連接關(guān)系研究主要有3類分析視角：結(jié)構(gòu)連接、功能性連接和效應(yīng)性連接。結(jié)構(gòu)連接一般指解剖連接的網(wǎng)絡(luò)，研究不同的大腦結(jié)構(gòu)之間如何通過直接或間接的神經(jīng)纖維通路產(chǎn)生相互影響，常采用體內(nèi)侵入性同位素標(biāo)記示蹤技術(shù)或無創(chuàng)的磁共振擴散加權(quán)成像方法；功能性連接是指大腦信息處理過程中各區(qū)域間的神經(jīng)活動存在時間上的相關(guān)性，通過相關(guān)性度量來探究腦區(qū)間的功能聯(lián)系，常用的方法有時域相關(guān)、譜相干、鎖相值等[2]；效應(yīng)性連接被定義為一個神經(jīng)系統(tǒng)對另一個神經(jīng)系統(tǒng)施加的直接或間接影響，側(cè)重描述大腦各區(qū)域間的動態(tài)定向相互作用與因果關(guān)系[3]。考慮到不同腦區(qū)間的交互作用是定向的，效應(yīng)性連接分析越來越受到研究者的關(guān)注，出現(xiàn)了多種典型算法，在使用時應(yīng)對其原理與特點有所了解，才能合理選擇和有效使用。

效應(yīng)性連接可直接基于指定因果關(guān)系的模型估計，也可從數(shù)據(jù)驅(qū)動的角度直接從信號本身估計，其中基于格蘭杰因果關(guān)系(Granger causality，GC)派生出的一類效應(yīng)性連接分析的方法得到了廣泛應(yīng)用。代表性算法包括：定向傳遞函數(shù)(directed transfer function，DTF)、直接定向傳遞函數(shù)(direct DTF，dDTF)、偏定向相干(partial directed coherence，PDC)、廣義偏定向相干(generalized PDC，GPDC)。在對包括腦電圖(electroencephalogram，EEG)、腦磁圖(magnetoencephalography，MEG)、局部場電位(local field potential，LFP)、鋒電位(spike)、功能磁共振成像(functional magnetic resonance imaging，fMRI)等不同模態(tài)與尺度的神經(jīng)數(shù)據(jù)分析中取得了良好效果[4-8]。首先對此類常用算法的基本原理及其特點進行了系統(tǒng)介紹，接著對實際應(yīng)用中的模型擬合、非平穩(wěn)數(shù)據(jù)處理、閾值估計等注意事項進行總結(jié)，最后以較具代表性的廣義正交化的偏定向相干(generalized orthogonalized PDC，GOPDC)算法為例，展示了此類方法的應(yīng)用效果。

1 方法原理及功能特性

1.1 二元因果連接分析

格蘭杰因果的基本思想可追溯到維納提出的概念[9]：若對時間序列x(t)的預(yù)測可以通過納入時間序列y(t)的信息而得到改進，則表明y(t)會對x(t)產(chǎn)生因果影響。然而維納的概念缺乏能夠?qū)嶋H執(zhí)行的機制，格蘭杰后來在線性回歸模型的背景下將該預(yù)測思想進行了形式化表示：

(1)

(2)

與自回歸模型(2)中x(t)的誤差項相比，若回歸模型(1)中x(t)的預(yù)測誤差項ε1的方差顯著減小，則認(rèn)為時間序列y(t)對x(t)有因果影響。結(jié)合式(1)和(2)中誤差項的方差可對因果影響進行量化，量化后的值稱為格蘭杰因果指數(shù)GCI(Granger causality index)，可表示為[10]：

(3)

同理也可判斷是否存在反方向因果影響的問題。

GCI直接在時域度量因果關(guān)系，Geweke首次提出二元情況下的頻域格蘭杰因果關(guān)系。對式(1)兩端進行傅里葉變換(Fourier transform，F(xiàn)T)：

(4)

依據(jù)式(4)可將兩通道間的譜矩陣分解為：

S(w)=H(w)ΣH*(w)，

(5)

式中：Σ為式(1)中二者殘差項的協(xié)方差矩陣；*表示該矩陣的復(fù)共軛轉(zhuǎn)置。

結(jié)合式(4)與式(5)將不同頻率w下序列y(t)對x(t)產(chǎn)生的因果影響度量定義為[11]：

GCy→x(w)=

(6)

式中：Sxx為序列x(t)自身的功率譜；Σxx=var(ε1(t))，Σyy=var(ε2(t))，二者之間的協(xié)方差為Σxy= cov(ε1(t),ε2(t))。

在特定頻率w下，若序列y(t)對x(t)無因果影響，則對應(yīng)序列x(t)自譜成分分解中不包含序列y(t)的影響，GCy→x(w)為0。反之，GCy→x(w)值越大，表明在此頻率成分序列y(t)對x(t)的影響越大。同理，可定義序列x(t)對y(t)在不同頻率w下產(chǎn)生的因果影響[12]。

1.2 多元因果連接分析

格蘭杰因果關(guān)系最初只涉及二元概念，格蘭杰本人在后續(xù)研究中也指出[13]：只有在沒有其他變量影響這一二元判定過程的情況下，因果關(guān)系原則才成立。然而多通道神經(jīng)信號之間通常存在復(fù)雜的相互作用，對信號進行兩兩分析無法有效區(qū)分直接和間接效應(yīng)而可能導(dǎo)致錯誤的結(jié)論[14]。為分析多通道神經(jīng)信號中的整體多元結(jié)構(gòu)，主要采用多變量自回歸模型(multivariate autore-gressive，MVAR)解決這一問題。

用p階MVAR模型表示一組n通道同步記錄的信號X(t)=[x1(t),…,xn(t)]T：

(7)

式中：C(r)為滯后r時的系數(shù)矩陣,大小為n×n；E(t)∈Rn為一個均值為零的多變量高斯白噪聲序列。

系數(shù)矩陣C包含的因果耦合信息可以采用時域或頻域分析提取，鑒于大腦中普遍存在的節(jié)律與振蕩活動，以下重點介紹所涉及的頻域分析方法。對MVAR模型兩端同時進行FT可得到對應(yīng)關(guān)系：

(8)

E(w)=C(w)X(w)，

(9)

X(w)=C-1(w)E(w)=H(w)E(w)，

(10)

式中：I為單位矩陣；H(w)為系統(tǒng)轉(zhuǎn)移矩陣。

基于C(w)或H(w)，國外學(xué)者提出了多種具有不同功能特點的多元因果連接分析方法。

1.2.1 DTF

基于MVAR模型中轉(zhuǎn)移矩陣H(w)的性質(zhì)，Kaminski等[15]在1991年提出了DTF的分析方法：

(11)

式(11)為廣泛應(yīng)用的歸一化形式的DTF(w)，取值范圍為[0,1]。它描述了在頻率為w時通道j對通道i產(chǎn)生的影響，從信息傳遞角度可理解為通道j流入通道i的信息量占所有通道流入通道i信息量的比例。DTF無法區(qū)分兩通道信號間所有直接與間接因果連接的關(guān)系。

1.2.2 dDTF

為區(qū)分直接和間接的交互影響，Korzeniewska等[16]在2003年將重新定義的全頻定向傳遞函數(shù)(full-frequency DTF，ffDTF)與偏相干(partial coherence，pCOH)相結(jié)合提出了dDTF，用于刻畫多通道情形下通道間的直接因果連接關(guān)系：

dDTFj→i(w)=ffDTFj→i(w)pCOHij(w)。

(12)

ffDTF側(cè)重在全頻段分析確定通道j對i在不同頻帶上施加的因果作用，而pCOH在多通道情況下分析兩通道間的無向耦合，排除了其他通道的影響。

1.2.3 PDC

對以上基于DTF及其改進方法的計算過程都需對矩陣求逆，帶來的計算精度問題會對結(jié)果造成一定影響。2001年Baccala等基于系數(shù)矩陣C(w)的角度對pCOH進行因式分解，將pCOH的思想進行擴展提出了PDC分析方法[17]：

(13)

PDC度量了通道情形下兩兩通道之間的直接因果連接關(guān)系，從信息傳遞角度可理解為從通道j流出到通道i的信息量占通道j所有流出信息量的比例。

1.2.4 GPDC

(14)

GPDC通過引入目標(biāo)通道對應(yīng)的誤差方差項進行量級的規(guī)范，大大改善了PDC的性能，成為多元因果連接分析方法中的代表。

1.2.5 GOPDC

當(dāng)采用以上方法對神經(jīng)信號分析時，不可避免會受到容積傳導(dǎo)效應(yīng)的影響，多個潛在的信號源因距離或傳導(dǎo)介質(zhì)差異會導(dǎo)致各通道不同程度的混疊，MVAR模型參數(shù)對這種潛在的涂抹效應(yīng)較敏感，可能導(dǎo)致電極間出現(xiàn)虛假連接[19]。研究表明:這種虛假連接能通過正交化信號功率進行消除，而且譜相干方法中虛部不受容積傳導(dǎo)偽跡的影響[20]。Omidvarnia等[21]結(jié)合這兩種處理方式在2014年對PDC作出了進一步改進，提出正交化的偏定向相干(orthogonalized PDC，OPDC)，并借鑒GPDC推導(dǎo)，進一步提出了GOPDC。

(15)

該方法是在MVAR模型參數(shù)水平上進行正交化，Omidvarnia通過模型仿真進行方法對比，結(jié)果表明，GOPDC在減輕容積傳導(dǎo)效應(yīng)偽跡影響的同時，仍能抑制由通道誤差項幅度量級不同而產(chǎn)生的影響。

2 實際應(yīng)用需考慮的問題

實際應(yīng)用中，模型擬合、非平穩(wěn)數(shù)據(jù)的分析和非零值對應(yīng)閾值的估計是需要注意的問題。

2.1 模型擬合

2.1.1 數(shù)據(jù)平穩(wěn)性判斷

以上因果連接分析方法均要求數(shù)據(jù)為“弱平穩(wěn)”或協(xié)方差平穩(wěn)[22]，對數(shù)據(jù)平穩(wěn)性一般采取單位根檢驗。但由于最終目的是用數(shù)據(jù)擬合對應(yīng)的模型，因此數(shù)據(jù)擬合的角度更適用，首先要求MVAR模型中系數(shù)矩陣C平方可和[23]：

(16)

若MVAR模型系數(shù)矩陣C進一步滿足[24]：

(17)

將此條件進一步簡化可表示為：

(18)

式(18)中矩陣M為系數(shù)矩陣C的增廣矩陣。

若矩陣M所有特征值的模均小于1，則認(rèn)為模型系數(shù)表征了協(xié)方差平穩(wěn)的過程。在M的構(gòu)造過程中涉及模型系數(shù)矩陣C的求取，目前對MVAR模型系數(shù)矩陣C及誤差協(xié)方差Σ的估計算法較多，其中，Nuttall-Strand方法在不同數(shù)據(jù)集上均表現(xiàn)更好[25]。

2.1.2 模型階數(shù)選擇

MVAR模型階數(shù)p的選取非常重要，若p過小，模型不足以包含原數(shù)據(jù)的信息，反之會增加計算量，同時也易造成過擬合。最優(yōu)階數(shù)的選取常常依據(jù)一些常用的準(zhǔn)則，包括Akaike信息準(zhǔn)則(Akaike information criterion，AIC)與貝葉斯信息準(zhǔn)則(Bayesian information criterion，BIC)。已有研究表明，BIC準(zhǔn)則比AIC準(zhǔn)則更適合時間序列分析[22]。此外，估計理論要求已知樣本數(shù)量大于待估計未知參數(shù)數(shù)量[26]，因此p也存在上限，假設(shè)所分析數(shù)據(jù)的長度為L，則要求nL≥n2p；否則MVAR模型參數(shù)估計誤差將增大，導(dǎo)致模型合理性檢驗不能通過。

2.1.3 模型合理性檢驗

模型合理性檢驗是為確定選定階數(shù)后的模型是否對原數(shù)據(jù)進行充分地擬合，一般從一致性、穩(wěn)定性和殘差白噪聲性三方面進行檢驗。

模型一致性檢驗描述了原數(shù)據(jù)中所存在的相關(guān)結(jié)構(gòu)能被模型表征的程度[27]，對原始數(shù)據(jù)和模型產(chǎn)生的數(shù)據(jù)的所有通道分別做自相關(guān)與互相關(guān)處理，用Rr和Rs分別表示真實數(shù)據(jù)和模型產(chǎn)生數(shù)據(jù)得出的相關(guān)性向量，則一致性百分比可表示為：

(19)

當(dāng)PC值大于80%時，表示模型充分體現(xiàn)了原始數(shù)據(jù)的相關(guān)結(jié)構(gòu)[23]。

模型穩(wěn)定性與數(shù)據(jù)平穩(wěn)性判斷中所要求的條件相同，若系數(shù)矩陣C構(gòu)成的增廣矩陣M所有特征值的絕對值均小于1，則認(rèn)為擬合模型穩(wěn)定。

模型殘差白噪聲檢驗的目的是為了判斷模型殘差中是否仍存在一些未被描述的相關(guān)結(jié)構(gòu)，若模型充分?jǐn)M合了原始數(shù)據(jù)，則殘差項間近似于白噪聲。一般采用Durbin-Watson統(tǒng)計檢驗[28]對模型的殘差E(t)進行相關(guān)性判斷。

以上三方面的檢驗均滿足的情況下，由數(shù)據(jù)所擬合的MVAR模型才可用于下一步的因果分析。

2.2 非平穩(wěn)數(shù)據(jù)處理

實際應(yīng)用中EEG、LFP等神經(jīng)數(shù)據(jù)多是高度非平穩(wěn)的，且與腦區(qū)間因果連接關(guān)系的變化有關(guān)[22]，而對數(shù)據(jù)進行“平穩(wěn)化”操作的方式往往會破壞其中所蘊含的結(jié)構(gòu)，得出的結(jié)果往往與所采取的處理方式有關(guān)，具有一定的虛假性。為研究非平穩(wěn)數(shù)據(jù)中所蘊含的時變因果連接變化，目前時變MVAR模型已成為主流，其中應(yīng)用較廣泛的兩種方式包括[29]：數(shù)據(jù)截段滑窗法和狀態(tài)空間法。

數(shù)據(jù)截段滑窗法最早由Ding等[27]在2000年提出，主要策略是認(rèn)為原數(shù)據(jù)在足夠小的時間窗口內(nèi)近似平穩(wěn)，從而對原數(shù)據(jù)進行高度重疊的窗口截段，可從這些高度重疊的數(shù)據(jù)段內(nèi)得到隨時間變化的系數(shù)矩陣C等其他參數(shù)，進而得到隨時間變化的因果連接關(guān)系的度量。該方法需要考慮數(shù)據(jù)窗口大小的選擇，一般依據(jù)以下4個原則進行折中[30]：

(1)數(shù)據(jù)截段的同時保證窗口內(nèi)數(shù)據(jù)局部平穩(wěn)；

(2)數(shù)據(jù)截段重疊滑窗本質(zhì)上為時域平滑操作，窗口大小選擇應(yīng)在考慮平滑性的同時保留原數(shù)據(jù)中動態(tài)變化的特征；

(3)保證已知數(shù)據(jù)的數(shù)量大于模型中所需估計未知參數(shù)的數(shù)量；

(4)數(shù)據(jù)段對應(yīng)的時間間隔需大于生理意義下該過程動態(tài)變化所預(yù)期的最大時間間隔。

狀態(tài)空間法從參數(shù)自適應(yīng)估計角度出發(fā)，對數(shù)據(jù)是否平穩(wěn)并無要求，其實質(zhì)是將時變MVAR模型參數(shù)用狀態(tài)空間形式進行表示，然后結(jié)合數(shù)據(jù)采用卡爾曼濾波算法進行參數(shù)自適應(yīng)估計。Wan等[31]在2003年提出的對偶擴展卡爾曼濾波(dual extended Kalman filter，DEKF)算法實現(xiàn)了上述過程，該方法采用了兩個相互交錯的狀態(tài)空間方程，分別用于系統(tǒng)狀態(tài)估計和參數(shù)估計，詳細(xì)推導(dǎo)過程可見Omidvarnia等[32]2011年發(fā)表的論文。采用狀態(tài)空間方法可估計出時變MVAR模型中不同時刻的系數(shù)矩陣C等參數(shù)，由此推導(dǎo)出的多種時變因果連接關(guān)系度量更能體現(xiàn)連接關(guān)系的瞬時變化。

2.3 閾值估計問題

實際應(yīng)用中的干擾因素會導(dǎo)致不存在真正相互作用的通道間的因果連接度量在某些頻率上出現(xiàn)非零情況，因此估計不同頻率下因果連接關(guān)系度量閾值具有十分重要的現(xiàn)實意義。

閾值估計一般通過統(tǒng)計檢驗方法來解決，首先，基于無相互因果作用的零假設(shè)，計算通道間不同頻率下的因果連接關(guān)系值的經(jīng)驗分布，然后根據(jù)得出的置信區(qū)間上限確定出所需閾值，其中基于原數(shù)據(jù)生成替代數(shù)據(jù)來確定閾值是國內(nèi)外學(xué)者廣泛采用的方式。由于不同通道間因果連接關(guān)系的存在與其相位有關(guān)，生成替代數(shù)據(jù)的基本思想是將原數(shù)據(jù)對應(yīng)的相位隨機化，即破壞原數(shù)據(jù)間相位的依賴關(guān)系[33]。替代數(shù)據(jù)的生成步驟主要包括：①將原數(shù)據(jù)FT轉(zhuǎn)換到頻域；②對各頻率對應(yīng)的相頻部分施加隨機擾動；③通過傅里葉逆變換得到替代數(shù)據(jù)。替代數(shù)據(jù)產(chǎn)生方式因牽涉到傅里葉變換也被稱為FT surrogate，其產(chǎn)生的替代數(shù)據(jù)僅打亂了原數(shù)據(jù)間的相位信息，但同時保留了時域和幅頻特征，十分適合因果連接閾值估計。

考慮到非任務(wù)態(tài)與任務(wù)態(tài)腦活動中存在因果連接模式的區(qū)別，利用兩種狀態(tài)下通道間因果連接關(guān)系在統(tǒng)計意義上的差異性程度也可達(dá)到去閾值的效果[21]。該方式比FT surrogate所需的計算量少，且可以消除因大腦自發(fā)活動產(chǎn)生的因果連接關(guān)系的影響，并突出差異部分，去閾值效果更好，但該方式的應(yīng)用受限于實驗范式中對照組的設(shè)置，其通用性不如FT surrogate方法。

3 實際數(shù)據(jù)集應(yīng)用

算法示例選擇來自于Delorme等[34]2002年研究視覺處理過程中腦電動力學(xué)變化的EEG數(shù)據(jù)，隨機抽選5位受試者看到圖像中包含動物并作出正確反應(yīng)時的腦電數(shù)據(jù)，依據(jù)文中確定的偶極子空間位置及生理視覺通路相關(guān)信息，選取F3、C3、T5、O1、F4、C4、T6、O2共8個通道來研究圖片刺激后腦區(qū)間的信息交互關(guān)系。對數(shù)據(jù)進行去基線漂移、去工頻干擾等操作并剔除眼動干擾試次，選擇受試者各自反應(yīng)時間相對集中的試次(反應(yīng)時間為0.3～0.4 s，每位受試者平均約85個試次)，分別截取刺激前后各0.3 s的數(shù)據(jù)進行分析，其中前0.3 s數(shù)據(jù)作為對照組用于閾值的估計。采用GOPDC方法進行數(shù)據(jù)分析，其中時變MVAR模型中參數(shù)的估計采用DEKF算法，模型階數(shù)依據(jù)BIC準(zhǔn)則進行確定，5位受試者所有試次的平均結(jié)果如圖1所示。

圖1 基于GOPDC的5位受試者的平均結(jié)果Figure 1 Average results for 5 subjects based on GOPDC

從圖1的結(jié)果中可看出，基于DEKF方式的時變GOPDC能夠捕獲非平穩(wěn)腦電信號間的動態(tài)交互變化，并具有良好的時頻分辨率?？梢钥吹剑趫D像呈現(xiàn)后，主要是額葉與頂葉對顳葉和枕區(qū)產(chǎn)生因果影響，這與生理上視覺通路信息的傳遞機制相符。此外容積傳導(dǎo)所產(chǎn)生的影響一般在3 Hz以下[26]，從圖中可以看到3 Hz以下部分均未產(chǎn)生容積傳導(dǎo)中偽跡的影響。結(jié)果表明：采用GOPDC方法能夠有效地削弱這種影響，且所確定的交互頻段集中在θ頻段，與Delorme等采用信號獨立成分建立的等價偶極子模型確定的頻段相符。

4 結(jié)束語

對格蘭杰因果關(guān)系引申出的常用效應(yīng)性連接分析方法進行了概述，分別闡明了這些常用方法的優(yōu)缺點及其在因果連接關(guān)系度量中的意義。

對于實際應(yīng)用中MVAR模型擬合的相關(guān)問題，分別從數(shù)據(jù)平穩(wěn)性的判斷、模型階數(shù)的選擇、擬合后模型合理性驗證3個互相聯(lián)系的方面對需注意的要點進行概述，并建議從模型擬合的角度進行數(shù)據(jù)平穩(wěn)性的判斷。一般選用BIC準(zhǔn)則進行模型階數(shù)的選擇，分別從模型一致性、平穩(wěn)性和殘差白噪聲性等3個方面進行模型合理性驗證。針對非平穩(wěn)數(shù)據(jù)的處理，建議采用時變MVAR模型分析非平穩(wěn)數(shù)據(jù)中所蘊含的因果連接關(guān)系。針對數(shù)據(jù)截段滑窗法和狀態(tài)空間法這兩種常用的時變MVAR模型參數(shù)估計方式，給出數(shù)據(jù)截段滑窗法中窗長的選擇原則，并指出狀態(tài)空間法在數(shù)據(jù)處理和時域解析度上的優(yōu)勢。最后，關(guān)于因果連接分析中必須解決的閾值估計問題，對于無對照組的實驗建議采用國內(nèi)外學(xué)者廣泛使用的FT surrogate方式，而在有對照組的實驗中建議根據(jù)任務(wù)態(tài)與非任務(wù)態(tài)的差異進行閾值估計。

在實際EEG數(shù)據(jù)的應(yīng)用示例中，采用基于DEKF方式的時變GOPDC方法進行分析，結(jié)果證實了該方法能減輕容積傳導(dǎo)效應(yīng)中偽跡的影響。