之前用神經網絡煉丹的時候,受制于硬體條件,訓練起來極為痛苦,尤其是一開始我頭鐵一波直接用全腦的4D圖像作為輸入資料,網絡複雜度根本上不去,隻要結構稍微複雜一點就必定炸顯存。畢竟單獨一個時間點上的3D資料大小就能達到幾十萬(如果體素大小3毫米的話normalise之後就是61×73×61)。那就要考慮一些資料降維的方法了,先行提取一些特征出來(答辯的時候還有個老師問我為啥不直接用全腦資料來着,我說窮啊,要是硬體條件允許我也想試試)。
功能連接配接(Functional Connectivity)矩陣
先說說
功能連接配接矩陣吧,這應該是用得最多而且也比較簡單的一種方法了。思路也比較簡單:
fMRI是四維的時間序列,也就是在一段時間内不同時間點上的三維圖像。那麼,隻要提取出某一個體素在全部時間上的資料也就得到了
這個體素的時間序列。這個
時間序列也就代表着這個體素之内
血氧水準随時間的變化,也就反映了功能的表現。 是以,從
時間序列之間的相關性,就可以得到功能的相關性。(當然這裡存在問題,這也是功能連接配接以及下面有效連接配接這兩種方法的缺點,因為這兩種方法兩兩計算腦區連接配接度,那麼就考慮不到是否存在互相作用的第三點的問題。)
計算功能連接配接一般使用
皮爾森相關系數。也有研究使用其他方法,以期得到更稀疏的功能連接配接(但我忘了具體咋算的了)。兩組時間序列的皮爾森相關系數由下式得,也就是這兩段時間序列的
協方差與除以兩段時間序列的标準差之積:
如果用Python寫的話可以直接使用numpy提供的函數:
corr = np.corrcoef(X, Y)
不過很多工具箱也提供了直接計算功能連接配接矩陣的功能(如
DynamicBC和
GRETNA,這些同樣也可以用來計算有效連接配接等其他腦網絡),就不需要自己手動寫腳本提取了。
皮爾森相關系數的結果
取值範圍為
。當
時,表示兩組時間序列正相關,反之則為負相關。當等于零時,代表兩組時間序列互相獨立,無相關關系。由此可以得出某兩組腦區在某一段時間内的功能是互相協同的還是互相拮抗的。
通過
計算各個腦區的平均時間序列,
兩兩計算相關系數,就可以得到這段時間内全腦的相關性矩陣,即功能連接配接矩陣。
功能連接配接矩陣,使用AAL90模闆,共90個腦區,是以連接配接矩陣為90x90
另外,還可以使用動态視窗法在不同的時間窗内計算功能連接配接矩陣,得到随時間變化的功能連接配接矩陣。
有效連接配接(Effective Connectivity)矩陣
有效連接配接與功能連接配接相似,但是計算的是
腦區之間的因果關系,也就是腦區之間的互相影響。與功能連接配接不同,
功能連接配接是
無向有權圖,得到的矩陣是
對稱矩陣;而
有效連接配接使用因果關系,是有方向的,也就是說會得到
有向有權圖,是以有效連接配接矩陣是
不對稱的。
常用的方法有
動态因果模型(Dynamic Causal Modelling, DCM)和
格蘭傑因果關系(Granger Causality)。我常用的方法是格蘭傑因果。其定義為:對于兩個平穩的時間序列X和Y,如果通過X和Y過去的值來預測X的目前值要比僅使用X自身的過去值預測X的目前值要更準确,則認為Y是X的格蘭傑原因。
打個比方就是,
半個小時前天氣很好很暖和,但是突然開始下雪,大雪下了半個小時,現在天氣很冷。假設
氣溫是X,
天氣狀況是Y,那麼隻使用
X的過去值(暖和)來預測
現在(冷)明顯是不準确的。那麼
同時考慮X和Y,因為
下雪是以
使暖和的溫度變冷了,這就要更為準确,是以我們可以認為
天氣狀況Y是氣溫X的原因。
那麼建構因果關系首先要建構X和Y的自回歸模型:
其中a和b是模型的自回歸系數,是常數;p為模型階數,即使用目前時刻前1到p個時間點來預測目前值;ε是預測誤差。因為現在的自回歸模型是
單變量的自回歸模型,即隻使用X來預測X(或隻使用Y來預測Y),是以
預測誤差隻取決于X(或Y)的過去值。我們可以使用預測誤差的方差來評估模型精度。
接下來咱們考慮
X和Y互相影響的情況,這時候在上式中X加入一項Y,在Y中加入一項X,構造雙變量自回歸模型:
c和d同樣是自回歸系數。這時
預測誤差ε會同時受到X和Y的影響。同樣使用預測誤差的方差來評估模型精度。那麼,
當雙變量自回歸模型的精度更高時,例如:
。
這個關系的度量通過下式計算:
反之,考慮X對Y的因果關系也同理。
計算的方法與功能連接配接相似,提取出腦區的平均時間序列兩兩計算。得到的矩陣也類似。
獨立成分分析(Independent Component Analysis,ICA)
這個方法就不太一樣了,
獨立成分分析可以把全腦信号分解成
若幹個互相獨立的信号源産生的
信号的疊加。例如,某個腦區同時受到了多個腦區産生的信号的影響,此時這個區域的信号是這些原始信号的疊加,那麼我們可以通過這種方法分離出原始信号源。
全腦的 fMRI 資料可以認為是由若幹組互相獨立的信号的線性疊加。是以利 用 ICA 方法可以提取出空間上獨立的信号源。(也可以進行時間獨立的分解。)
獨立成分分析的問題可以描述為如下:
簡單來說,
就是找到一個矩陣W,能夠将原始信号分解為若幹獨立的信号。
那麼,用什麼方法可以找到W呢?在這裡介紹一種
Infomax方法,按照資訊極大化原則進行。它的思想是,當
系統的輸出熵極大的時候(熵極大也就代表
輸出中包含的資訊量極大),這個時候輸出的
各分量之間的互資訊最小(也就是各分量中包含的
重複資訊得到了去除),這個時候可以認為各個輸出成分是互相獨立的。
Infomax算法原理圖
Infomax的實作和神經網絡非常相似,也就是
用W将初始信号分解為若幹分量後,把這若幹個輸出經過一個
非線性函數得到結果,
計算結果的總熵。然後在不斷的疊代中
調節W的值,使得總熵極大。
當然,我們還是可以使用MATLAB工具箱(例如Group ICA of fMRI Toolbox,GIFT)來幫助我們進行分析。
可以使用工具箱先對資料進行成分數量的估計,然後按照估計的結果設定要分離出多少個獨立成分。那麼最後我們可以得到一個尺寸為
獨立成分數x時間序列長度的二維矩陣(也就是分離出的
各個獨立信号源腦區的平均時間序列)。還可以通過GIFT工具箱檢視各個獨立成分的詳細資訊,例如位置、成分内各個體素位置的信号等等。
總結
以上大概介紹了三種特征提取的方式(也基本就是我畢設裡用的玩意了),其實适用于fMRI資料的降維手段還有很多,例如張量的方法,像Tucker分解或者CP分解(張量方法我還不熟悉,可能以後會再了解了解),當然也有一些研究選擇了直接用全腦的資料去煉丹……這樣會帶來巨大的時空開銷,我也沒實際仔細調試過,不好評價效果如何。
總之就是根據實際情況,我個人覺得ICA的方法更靈活,但也要更加麻煩一些。三種方法應用都很廣泛,根據實際情況選擇吧~