導(dǎo)讀

意義：功能性近紅外光譜（fNIRS）數(shù)據(jù)的靜息態(tài)功能連接（RSFC）分析揭示了整個大腦的皮層連接和網(wǎng)絡(luò)。fNIRS信號中的運動偽影和系統(tǒng)生理學(xué)提出了獨特的分析挑戰(zhàn)，目前已有研究探索了控制系統(tǒng)生理噪聲的方法。但是這些相同的方法在應(yīng)用于靜息態(tài)fNIRS（RS-fNIRS）數(shù)據(jù)時是否需要修改仍不清楚。

目的：本研究系統(tǒng)地考察了幾種RSFC分析管道的敏感性和特異性，以確定校正RS-fNIRS數(shù)據(jù)中全局系統(tǒng)生理信號的最佳方法。

方法：使用數(shù)值模擬的RS-fNIRS數(shù)據(jù)，比較了幾種連通性分析管道的真陽性率和假陽性率。使用受試者工作特征分析對其性能進(jìn)行評分。管道包括偏相關(guān)和多變量格蘭杰因果關(guān)系，有/沒有短距測量，以及包括非傳統(tǒng)零滯后交叉項的修正多變量因果關(guān)系模型。研究者還考察了預(yù)白化和穩(wěn)健統(tǒng)計估計對性能的影響。

結(jié)果：與以往關(guān)于雙變量相關(guān)模型的研究一致，本研究結(jié)果表明穩(wěn)健統(tǒng)計和預(yù)白化是校正fNIRS時間序列中運動偽影和自相關(guān)的有效方法。此外，當(dāng)兩個感興趣的信號同步波動時，使用從短距fNIRS通道提取的主成分作為偏相關(guān)模型的一部分進(jìn)行預(yù)濾波，在減少由于共有的系統(tǒng)生理學(xué)引起的偽相關(guān)方面最有效。然而，當(dāng)信號之間存在時間滯后時，結(jié)合短距通道的多變量格蘭杰因果檢驗效果更好。由于不知道實驗數(shù)據(jù)中是否存在這樣的滯后，研究者提出了一種修正版本的格蘭杰因果關(guān)系，其中包括非傳統(tǒng)的零滯后項作為折中解決方案。

結(jié)論：建議結(jié)合預(yù)白化、穩(wěn)健統(tǒng)計方法和處理管道中的偏相關(guān)性，以減少自相關(guān)、運動偽影和全局生理噪聲，從而獲得具有統(tǒng)計有效性的RS-fNIRS連通性指標(biāo)。未來的研究應(yīng)使用人類數(shù)據(jù)來驗證這些方法的有效性。

前言

自發(fā)或靜息態(tài)功能連接（RSFC）已成為表征人腦網(wǎng)絡(luò)結(jié)構(gòu)的重要工具。Biswal等人首先指出，功能性磁共振成像（fMRI）血氧水平依賴性（BOLD）信號中的同步自發(fā)波動被組織成不同的網(wǎng)絡(luò)，這些網(wǎng)絡(luò)在被試之間保持不變。25年后，這些結(jié)果在數(shù)百項研究中被復(fù)制，并擴展到人類大腦功能、發(fā)育和疾病的各個方面。靜息態(tài)網(wǎng)絡(luò)也使用其他神經(jīng)成像方法進(jìn)行了觀察，包括EEG和MEG，以及fNIRS。靜息態(tài)范式的一個主要優(yōu)點是不需要執(zhí)行任務(wù)，這在嬰兒、老年人、患者等人群的研究中特別有用。

目前已經(jīng)開發(fā)了幾種分析方法來量化靜息態(tài)連通性。對于fMRI和fNIRS，這些方法通常關(guān)注大腦區(qū)域之間的零滯后關(guān)系。零滯后關(guān)系是指兩個信號同時發(fā)生信號波動（在測量采樣時間內(nèi)）的關(guān)系。估計功能連通性的最常用方法可能是計算兩個腦區(qū)時間序列之間的Pearson相關(guān)系數(shù)。有幾種方法可以在時域中表征RSFC，包括（i）基于種子的連接，（ii）成對連接，以及（iii）獨立成分分析（ICA）。在基于種子的連接中，選擇一個種子區(qū)域，并評估其與所有其他腦區(qū)域的連通性。在成對連接中，計算所有可能通道對之間的相關(guān)性。ICA將fNIRS通道數(shù)據(jù)的協(xié)方差分解為幾個統(tǒng)計上獨立的時空成分，其中一些成分表示每個被試的大腦功能網(wǎng)絡(luò)，或者表示跨被試連接的功能網(wǎng)絡(luò)。

格蘭杰因果關(guān)系是“有效連接”方法的一個例子，該方法依賴于fNIRS通道在時域上的滯后關(guān)系。另一種常用于模擬fNIRS中功能連通性的方法是小波變換相干（WTC），該方法在時頻域中模擬通道之間的關(guān)系。WTC依賴于連續(xù)小波變換（CWT），CWT將時間序列數(shù)據(jù)轉(zhuǎn)換為時頻域。WTC估計為時間序列在時頻域中的互相關(guān)。

功能性近紅外光譜

fNIRS是一種光學(xué)成像方法，使用近紅外波長范圍（650至950nm）的光來估計氧合血紅蛋白（HbO）和脫氧血紅蛋白（HbR）濃度的時間變化。fNIRS測量大腦的血流動力學(xué)變化，反映血流量、體積和血氧合的變化，因此，它是一種間接測量神經(jīng)活動的方法，就像fMRI BOLD信號一樣。與fMRI相比，fNIRS具有一些優(yōu)勢，包括成本更低，對頭部運動的彈性更大，以及更高的數(shù)據(jù)采樣率。fNIRS也可以便攜，能夠執(zhí)行更靈活或生態(tài)有效性的任務(wù)。然而，fNIRS的缺點是其空間分辨率較低，深入大腦的深度有限，相對于大腦本身，對頭部表層的敏感性增加。對于fNIRS測量來說，后一個問題尤其令人困惑，因為頭皮的系統(tǒng)生理信號，包括心臟、呼吸和血壓波動，通常會導(dǎo)致空間上的全局噪聲，其振幅通常大于感興趣的潛在大腦信號。因此，在使用fNIRS的RSFC實驗中，將表層生理噪聲和真正的潛在神經(jīng)相關(guān)腦網(wǎng)絡(luò)進(jìn)行分離是具有挑戰(zhàn)性的。

Obrig等人首次證明了fNIRS在捕獲極低頻和低頻振蕩方面的可行性，這些振蕩有助于RSFC的研究。之后，幾項用于考察RSFC可行性的早期研究使用fNIRS和擴散光學(xué)層析成像（DOT）發(fā)現(xiàn)了可靠的和預(yù)期的視覺、感覺運動和語言網(wǎng)絡(luò)。從那時起，研究人員使用靜息態(tài)fNIRS（RS-fNIRS）和DOT研究嬰兒發(fā)育、前額葉皮層的性別差異、早產(chǎn)、語言網(wǎng)絡(luò)、自閉癥譜系障礙、情感障礙和衰老。檢查RSFC的重測信度發(fā)現(xiàn)，在從一小時到一周的多個sessions中，個體和群體水平的RSFC圖之間存在良好到極好的組內(nèi)相關(guān)性。

RS-fNIRS信號的統(tǒng)計特性

鑒于RS-fNIRS的信號和噪聲特性，與RS-fMRI相比，RS-fNIRS具有一些獨特的挑戰(zhàn)。由于fNIRS信號的采集方式，它對頭皮的非神經(jīng)腦外信號變化很敏感。RS-fNIRS信號由三部分組成：（i）低頻神經(jīng)振蕩（感興趣信號），（ii）源自大腦的腦內(nèi)生理噪聲，以及（iii）源自頭皮和非腦組織的腦外生理噪聲。鑒于靜息態(tài)范式的性質(zhì)，所有信號都是非誘發(fā)的。因此，fNIRS信號同時受到腦外信號變化和非神經(jīng)腦信號變化的影響，包括心臟（約1-1.2Hz）、呼吸（0.3-0.6Hz）和血壓/Mayer波（0.1Hz）波動引起的系統(tǒng)生理噪聲。生理噪聲對RS-fNIRS信號有兩種影響：①誘導(dǎo)時間自相關(guān)和②增加大腦通道之間的空間協(xié)方差。

先前的研究表明，高采樣率和系統(tǒng)生理信號(心率、呼吸和血壓)，結(jié)合潛在的血流動力學(xué)反應(yīng)，可以共同導(dǎo)致時序相關(guān)，從而在RS-fNIRS時間序列中會產(chǎn)生自相關(guān)。去除時間序列中的自相關(guān)以獲得統(tǒng)計有效的關(guān)系估計，可以使用預(yù)白化的方法。預(yù)白化是通過去除時間序列中的序列相關(guān)來實現(xiàn)的。自回歸(AR)模型可用于數(shù)據(jù)預(yù)白化和去除序列相關(guān)。具體來說，對于fNIRS，使用預(yù)白化需要使用p階AR模型來模擬當(dāng)前時間點作為前一個時間點的函數(shù)。AR模型中去除時間相關(guān)性后的殘差項被用來計算通道之間的Pearson相關(guān)系數(shù)，而不是使用自相關(guān)通道數(shù)據(jù)。當(dāng)將預(yù)白化納入到預(yù)處理管道中用于功能連接估計時，可以減少I型錯誤。盡管對fNIRS的時間自相關(guān)進(jìn)行了校正，但如果沒有考慮由于系統(tǒng)生理機能導(dǎo)致的空間協(xié)方差增加，這也可能會導(dǎo)致FDR的增加。

由于系統(tǒng)生理信號被認(rèn)為在整個大腦中是相對均勻的，所以由所有通道共有的系統(tǒng)波動引起的信號變化可以增加通道中的空間協(xié)方差。圖1顯示，不校正fNIRS數(shù)據(jù)中增加的空間協(xié)方差可能會導(dǎo)致I型誤差增加和特異性降低，因為基于地面真值的連接路徑遠(yuǎn)比預(yù)期的多。

圖1.連接圖顯示，在模擬的fNIRS數(shù)據(jù)中，加入模擬AR平滑和共享協(xié)方差的生理噪聲，Pearson相關(guān)系數(shù)、AR相關(guān)和AR偏相關(guān)可以有效地恢復(fù)真實的連接。

方法

為了表征每種分析方法的敏感性和特異性，本研究根據(jù)具有已知地面真值的數(shù)值模擬數(shù)據(jù)，探索了模型在幾種模擬情況中的性能，從簡單正態(tài)分布隨機噪聲（其中本研究模型的所有統(tǒng)計假設(shè)都有效）到具有時間自相關(guān)噪聲的生理信號，再到全局空間生理噪聲。最后，加入運動偽影。每個數(shù)據(jù)模擬都代表了一個更具挑戰(zhàn)性的統(tǒng)計問題，通過模擬數(shù)據(jù)中的真實和零相關(guān)性，使用受試者工作特征（ROC）分析來檢驗這個問題。所有方法均在MATLAB（R2019b）的NIRS AnalyzIR工具箱中實現(xiàn)。

數(shù)據(jù)模擬

采用具有指定空間協(xié)方差矩陣的多變量正態(tài)分布方法，對由16個長距通道（27mm）和9個短距通道（7.1mm）組成的探測器進(jìn)行了fNIRS數(shù)據(jù)模擬。模擬了近紅外波長為690和830nm的數(shù)據(jù)。這種已知的空間協(xié)方差在數(shù)據(jù)中產(chǎn)生潛在的“真實”相關(guān)性，在此基礎(chǔ)上評估各種算法。這些真實相關(guān)性是通過通道之間的零滯后（無時間偏移）或時間滯后（數(shù)學(xué)因果模型）產(chǎn)生的。僅針對長距通道在通道空間中模擬目標(biāo)真實相關(guān)性。為了給模型添加額外的噪聲，使用半無限均勻平板模型來計算光學(xué)正向模型。使用光學(xué)正向模型將這些體素投影到通道空間，并添加到基本真值相關(guān)信號。

對于每個模擬，以4Hz采樣率生成300 s的靜息態(tài)數(shù)據(jù)。為了生成ROC曲線，每個測試總共生成了200個數(shù)據(jù)集。在ROC分析中，在探測器中任何兩個通道之間的所有可能連接中，模擬了10%左右的可能連接路徑的真實連接，確保全局系統(tǒng)噪聲是主要的空間噪聲特征。因此，用于生成模擬數(shù)據(jù)的協(xié)方差矩陣中只有10%的非對角元素為非零。由于時間自相關(guān)、全局系統(tǒng)生理信號或運動偽影，其余90%路徑之間的任何顯著連接都是假陽性。為了生成真和假陽性估計，從總鄰接矩陣中隨機抽取相等數(shù)量的真和空連接。

fNIRS實驗數(shù)據(jù)

采集

研究者還測試了本研究的功能連接方法對實驗靜息態(tài)數(shù)據(jù)的有效性和適用性。使用商業(yè)NIRScout-2（NIRx Medizintechnik GmbH，柏林，德國）連續(xù)波fNIRS系統(tǒng)記錄NIRS數(shù)據(jù)。共有50個通道（42個長距通道，8個短距通道）分布在雙側(cè)額葉和感覺運動皮層（見圖2）。如圖2所示，長距通道由放置在頭皮上的16個源光電二極管(橙色圓圈)和13個探測器(藍(lán)色矩形)組成。8個探測器（綠色菱形）放置在頭皮的8個短距通道上。藍(lán)色實線表示長距通道，綠色虛線表示短距通道。以7.8125Hz的采樣率記錄兩個波長（760nm和850nm）的數(shù)據(jù)。定位頭帽后，使用NIRx Aurora軟件優(yōu)化信號質(zhì)量。

圖2.fNIRS數(shù)據(jù)采集示意圖。

參與者和任務(wù)

約24名健康被試參與了實驗（9名男性和15名女性；平均年齡27.4歲，SD=8.7歲；23名右利手）。告知被試實驗情況，并獲得書面同意。這項研究得到了匹茲堡大學(xué)機構(gòu)審查委員會的批準(zhǔn)。每個被試進(jìn)行一次靜息態(tài)掃描和其他任務(wù)掃描。指示被試在靜息態(tài)掃描期間盡量減少身體運動，保持坐姿放松5分鐘，而不使用任何腦力勞動。

數(shù)據(jù)預(yù)處理

將模擬的fNIRS信號強度轉(zhuǎn)換為光密度的變化，然后使用修正的Beer-Lambert定律（MBLL）將光密度單位轉(zhuǎn)換為HbO和HbR濃度，兩種波長的微分路徑長度因子為6，部分體積校正均為60。除了對實驗fNIRS時間序列進(jìn)行降采樣外，以與模擬數(shù)據(jù)類似的方式對數(shù)據(jù)進(jìn)行預(yù)處理。使用MBLL將原始信號強度轉(zhuǎn)換為HbO和HbR濃度后，將實驗靜息態(tài)fNIRS數(shù)據(jù)從7.8125Hz降采樣至4Hz，以匹配模擬數(shù)據(jù)的采樣率。在本文的其余部分，分析的重點在于HbO濃度的變化，結(jié)果也應(yīng)適用于HbR。

圖3顯示了三種情況下模擬數(shù)據(jù)的fNIRS時間序列示例，（i）僅時間自相關(guān)，（ii）具有時間自相關(guān)和全局系統(tǒng)生理信號，以及（iii）具有時間自相關(guān)、全局系統(tǒng)生理和運動偽影。圖3(左)顯示了在使用MBLL將原始模擬fNIRS信號轉(zhuǎn)換為HbO信號變化后，三種情況下的兩個通道的時間軌跡。圖3(右)顯示了通道的相關(guān)結(jié)構(gòu)。正如預(yù)期的那樣，隨著全局系統(tǒng)生理信號的引入，這些通道中的協(xié)方差增加。隨著在模擬數(shù)據(jù)中引入運動偽影，隨著通道之間的相關(guān)值接近1，協(xié)方差進(jìn)一步急劇增加。

圖3.顯示了模擬數(shù)據(jù)中兩個代表性通道HbO變化的時間軌跡以及通道的歸一化協(xié)方差矩陣。

預(yù)白化

預(yù)白化通常使用AR模型來去除時間序列中的自相關(guān)，并對信號的頻率內(nèi)容進(jìn)行白化。預(yù)白化大大降低了連接模型中的FDR。這里，使用AR模型實現(xiàn)預(yù)白化，如下所示：

fNIRS時間進(jìn)程可表示為：

其中εt∈N(0，σ2)，yt是fNIRS的時間序列，ai是AR參數(shù)，p是模型階數(shù)，由貝葉斯信息準(zhǔn)則（BIC）確定，εt是新息項，用于計算連通性指標(biāo)，而不是時間序列本身。圖4顯示了三組模擬數(shù)據(jù)預(yù)白化前后代表性通道的（a）功率譜密度和（b）自相關(guān)函數(shù)圖，包括：（i）時間自相關(guān)，（ii）時間自相關(guān)和全局系統(tǒng)生理信號，（iii）時間自相關(guān)、全局系統(tǒng)生理和運動偽影。正如預(yù)期的那樣，預(yù)白化步驟使功率譜密度圖變得平坦，并降低了fNIRS時間序列中的時間自相關(guān)。

圖4.顯示了三組模擬數(shù)據(jù)預(yù)白化前后代表性通道的(a)功率譜密度和(b)自相關(guān)函數(shù)圖。

穩(wěn)健方法

在對時間序列進(jìn)行預(yù)白化后，新息項包含在該時間點添加的信息。由于運動偽影被視為異常值，它們可能會影響正常性假設(shè)。通常，這些高影響點會對連通性指標(biāo)的估計產(chǎn)生重大影響。如前所述，運動偽影對連接性指標(biāo)的影響不同，這取決于運動是跨多個通道共享還是僅限于幾個通道。跨多個通道共享的運動損壞時間點可以通過增加它們之間的協(xié)方差來影響功能連接。因此，為了識別統(tǒng)計異常值，研究者使用了之前在Santosa等人研究中引入的聯(lián)合加權(quán)函數(shù)的多變量版本。聯(lián)合加權(quán)函數(shù)可以通過計算時間序列的幾何長度來識別異常值，并根據(jù)它們偏離平均值的方式對它們進(jìn)行適當(dāng)?shù)慕禉?quán)，將運動損壞的異常值分配較低的權(quán)重，使它們對相關(guān)和回歸方法的影響較低。

根據(jù)所有通道（包括短距通道）的新息項計算幾何長度函數(shù)rt，如下所示：

其中εi，t表示時間點t的第i個通道的新息項，c是通道總數(shù)。

根據(jù)離群值偏離平均值的程度對離群值進(jìn)行加權(quán)的加權(quán)函數(shù)(S)，可使用Tukey雙平方函數(shù)的平方根，對所有時間點進(jìn)行如下計算：

這里，σ是分布的標(biāo)準(zhǔn)偏差，根據(jù)中值絕對偏差（MAD）和正態(tài)分布標(biāo)準(zhǔn)偏差（σ）之間的關(guān)系估計(σ≈ 1.4826*MAD）。使用MAD是因為它是更穩(wěn)健的分散統(tǒng)計。常數(shù)κ的值為4.685，這對于沒有任何異常值的標(biāo)準(zhǔn)正態(tài)分布提供了95%的效率。具有對角元素的加權(quán)矩陣S表示所有時間點的加權(quán)函數(shù)可與所有時間序列的新息項相乘，從而降低正態(tài)分布統(tǒng)計異常點的權(quán)值。

其中，Yc是所有通道重加權(quán)后的新息項矩陣，Ec是預(yù)白化后獲得的新息項矩陣。

因此，為了減少運動偽影的影響，對于回歸模型，不使用線性模型和普通最小二乘估計參數(shù)，而是使用穩(wěn)健統(tǒng)計，使用加權(quán)最小二乘（WLS）來降低從學(xué)生化殘差計算出的異常值的權(quán)重。權(quán)重是使用Tukey的雙平方函數(shù)計算的。然后將計算出的權(quán)重應(yīng)用于原始數(shù)據(jù)，并重新估計殘差。重復(fù)該過程直到收斂。穩(wěn)健回歸在內(nèi)置的MATLAB函數(shù)“robustfit”中實現(xiàn)。MATLAB函數(shù)robustfit使用50次迭代或回歸系數(shù)值無可檢測數(shù)值變化的收斂標(biāo)準(zhǔn)來評估系數(shù)的收斂性。由于多變量格蘭杰因果關(guān)系（MVGC）和修正的MVGC在50次迭代中不收斂，研究者將最大迭代限制更改為150次迭代，以確?；貧w參數(shù)收斂。

使用WLS回歸估計的參數(shù)可用如下公式獲得：

這里，是參數(shù)估計，S是加權(quán)矩陣。

用于減少全局信號變化的基于PCA的回歸

如果探針廣泛分布在大腦中，那么跨多個通道共享的大信號方差可能主要反映生理噪聲。因此，PCA的前幾個成分可用于識別由于生理噪聲引起的全局信號變化，并從通道數(shù)據(jù)中回歸其影響。數(shù)據(jù)協(xié)方差矩陣可以使用奇異值分解法進(jìn)行分解

其中C是給定C=YT·Y的數(shù)據(jù)協(xié)方差矩陣，U是包含左奇異向量的矩陣。∑是包含奇異值的對角矩陣，V是包含右奇異向量的矩陣。通常，使用解釋最大方差的前幾個成分。研究者使用最小數(shù)量的主成分來解釋至少90%的通道協(xié)方差。對長距通道、短距通道或兩者進(jìn)行協(xié)方差分解可以提取全局信號。這三種選項都在測試中進(jìn)行了檢查。經(jīng)PCA回歸后的通道數(shù)據(jù)如下：

這里，y是新息項或通道的原始數(shù)據(jù)，取決于分析通道。Z是包含前幾個成分作為列的矩陣。

功能連接指標(biāo)

本研究考察了幾種估計功能連接的方法。在這里，重點關(guān)注時域中的方法，包括相關(guān)性、偏相關(guān)和格蘭杰因果關(guān)系的幾種變式。時頻方法，如小波相干，由于兩個原因沒有在本研究中進(jìn)行檢查。首先，在研究者之前的工作中發(fā)現(xiàn)，在雙變量連接性分析中，小波相干優(yōu)于ROC分析中的標(biāo)準(zhǔn)Pearson相關(guān)，但性能不如預(yù)白化的穩(wěn)健相關(guān)模型。除非標(biāo)準(zhǔn)小波相干模型被修改為包括高階AR項和穩(wěn)?。ó惓Ｖ稻芙^）方法，否則小波相干也具有較高的假陽性率和不受控的I型誤差。經(jīng)過這些修改發(fā)現(xiàn)，與預(yù)白化穩(wěn)健模型相比，小波相干沒有任何優(yōu)勢。第二，小波相干的多變量或偏相關(guān)版本尚未描述，需要開發(fā)一個新的統(tǒng)計模型。

基于相關(guān)性的分析模型

Pearson相關(guān)系數(shù)。兩個時間序列x和y之間的Pearson相關(guān)系數(shù)計算如下：

其中σx，y是時間序列x和y的協(xié)方差，σx和σy分別是x和y標(biāo)準(zhǔn)差。標(biāo)準(zhǔn)Pearson相關(guān)性并不是一個穩(wěn)健的統(tǒng)計估計，Santosa等人在fNIRS的背景下提出并描述了若干穩(wěn)健變量。穩(wěn)健相關(guān)性計算為回歸模型y對x和x對y的標(biāo)準(zhǔn)回歸系數(shù)的幾何平均值，以減少x或y中的異常值，但這仍然會影響相關(guān)估計。步驟如下：

基于上述回歸模型中的學(xué)生化殘差估計加權(quán)矩陣，加權(quán)矩陣Sx→y在上述回歸模型的兩側(cè)相乘，并且更新參數(shù)估計。

同樣，可以在y上擬合x，

加權(quán)矩陣Sy→x根據(jù)學(xué)生化殘差估計，以降低x中的異常值。Sy→x在上述回歸模型的兩側(cè)相乘，并更新參數(shù)估計。

兩個回歸模型交替求解，直到參數(shù)βy→x和βx→y收斂。使用以下公式估計最終穩(wěn)健相關(guān)系數(shù)：

穩(wěn)健相關(guān)系數(shù)估計的符號由回歸系數(shù)的符號確定。有效自由度為

。通過降低數(shù)據(jù)點的權(quán)重，特別是對異常點的權(quán)重，降低了有效自由度。

對于標(biāo)準(zhǔn)和穩(wěn)健相關(guān)估計，在給定零假設(shè)（p值）的情況下，觀察數(shù)據(jù)的概率可以通過使用T統(tǒng)計量估計

其中df為標(biāo)準(zhǔn)Pearson相關(guān)的有效自由度，值為n?2，r為Pearson相關(guān)系數(shù)。這里n是時間點的數(shù)量。

自回歸相關(guān)。AR相關(guān)性是基于預(yù)白化時間序列計算的Pearson相關(guān)系數(shù)。如Santosa等人所述，為每個數(shù)據(jù)通道估計p階AR模型，并用于預(yù)白信號。然后在相關(guān)分析中使用AR模型中的新息項（例如，在每個時間步長添加的新信息）。研究者使用的最大模型階數(shù)為40（10×采樣率），高于模擬fNIRS數(shù)據(jù)時使用的AR濾波器。根據(jù)經(jīng)驗觀察，在實驗靜息態(tài)fNIRS數(shù)據(jù)中，模型階數(shù)通常＜10·Fs（此處Fs為采樣頻率）。

AR偏相關(guān)。偏相關(guān)涉及一個兩步過程，其中首先將其他通道從感興趣的數(shù)據(jù)中投影出來，然后對殘差執(zhí)行相關(guān)。使用基于PCA回歸的投影來部分消除其他通道的影響。PCA回歸用于避免模型的共線不穩(wěn)定性。使用最小數(shù)量的特征向量來解釋至少90%的空間協(xié)方差。對于兩個通道i和j，如果預(yù)白化后的新息項由εi和εj表示，則白化PCA濾波時間序列計算為

現(xiàn)在，當(dāng)從通道中去除全局信號時，AR偏相關(guān)可以計算為兩個殘差之間的Pearson相關(guān)性?；诒Ａ舻闹鞒煞?jǐn)?shù)量進(jìn)一步減少自由度。該模型的輸入數(shù)據(jù)是預(yù)白化新息項時間序列，而不是通道數(shù)據(jù)，因此該方法被稱為AR偏相關(guān)。

基于格蘭杰因果關(guān)系的分析模型

多變量格蘭杰因果關(guān)系。在控制其他時間序列的影響后，MVGC檢驗兩個時間序列之間的滯后關(guān)系。這本質(zhì)上是兩個回歸模型的模型擬合檢驗?？梢允褂枚嘧兞緼R模型（MVAR）控制通道Z后檢驗x“格蘭杰原因”是否為y。對于通道時間序列，可以將當(dāng)前時間點yt建模為限制性模型中y和其他通道Z之前滯后的函數(shù)：

（限制性模型）

在非限制性模型中，添加x的滯后項，以模擬y中的當(dāng)前時間點：

（非限制性模型）

這里，p是模型中包含的滯后項的數(shù)量，使用BIC的擬合測量確定。為了檢驗因果關(guān)系，進(jìn)行了以下統(tǒng)計檢驗，與沒有x（限制性模型）的模型相比，在x（非限制性模型）中添加額外滯后項是否提高了yt的可預(yù)測性。如果x的過去包含關(guān)于當(dāng)前時間點yt的預(yù)測信息，則即使在控制丟失的自由度之后，yt預(yù)測中的誤差也會得到改善。這在統(tǒng)計上可以評估為：

零假設(shè)：x不是y的格蘭杰原因

H0: β1=β2=β3=……βp=0

備擇假設(shè)：x是y的格蘭杰原因。

HA：至少有一個βi為非零

因為有嵌套模型，所以可以使用嵌套F統(tǒng)計量來檢驗格蘭杰因果關(guān)系

這里，SSER和SSEU分別是限制性和非限制性模型的平方誤差之和，n是時間點的數(shù)量。因此，可以在自由度p和[n?(2p+1)]處計算F統(tǒng)計量的p值。格蘭杰因果指標(biāo)G也可以用于評估滯后關(guān)系的強度，如下所示：

在雙變量格蘭杰因果關(guān)系上使用MVGC的目的是控制包括系統(tǒng)生理信號在內(nèi)的第三變量的影響。在MVGC模型中添加多個通道的滯后項，容易使回歸模型無法求解，因為變量的數(shù)量將超過時間點的數(shù)量。因此，與使用通道數(shù)據(jù)本身不同，首先執(zhí)行PCA分析，并像之前使用偏相關(guān)方法一樣，包含解釋至少90%方差的最少成分。然后，在限制性和非限制性模型中包括成分和時滯成分，而不是通道數(shù)據(jù)。

修正的多變量格蘭杰因果關(guān)系。傳統(tǒng)的MVGC僅包括方程右側(cè)的y、x和Z的時間滯后歷史。因此，MVGC僅對非零滯后關(guān)系建模。重要的是要認(rèn)識到零滯后與采樣率有關(guān)，在模擬fNIRS數(shù)據(jù)的4Hz采樣時，一個采樣（250ms）內(nèi)發(fā)生的任何同步信號變化在數(shù)學(xué)上都是零滯后。雖然大腦中的神經(jīng)信號傳遞相當(dāng)快（100ms），但我們不希望假設(shè)所有的相關(guān)性都是零滯后，因而本研究引入了一個改進(jìn)的MVGC模型，以建模通道之間的滯后和零滯后關(guān)系。在修正的MVGC中，將Z的零滯后項包括在限制性和非限制性模型中，以及通道x的零滯后項包括在非限制性模型中。因此，限制性模型可以寫成：

（限制性模型）

類似地，非限制性模型可以寫成：

（非限制性模型）

上述方程還包括Z和x的零滯后項。現(xiàn)在，為了檢驗格蘭杰因果關(guān)系的統(tǒng)計顯著性，檢驗了如下零假設(shè)和備擇假設(shè)：

H0: β1=β2=β3=……βp=0

HA：至少有一個βi為非零

統(tǒng)計顯著性采用嵌套F統(tǒng)計進(jìn)行評估，格蘭杰因果關(guān)系的計算與傳統(tǒng)MVGC類似。該修正模型的一個警告是，假設(shè)不檢驗y和x之間的連接是否是由于零或非零滯后項。我們可以通過修改非限制性模型以包含x的非零滯后項來進(jìn)一步檢驗這一點，使得非限制性模型和限制性模型之間的唯一區(qū)別是在x中包含零滯后項。

穩(wěn)健的多變量格蘭杰因果關(guān)系。最后，研究者開發(fā)了MVGC和修正的MVGC模型的穩(wěn)健統(tǒng)計版本。該多變量模型中離群值的定義具有挑戰(zhàn)性，因為如雙變量相關(guān)模型，離群值可能存在于y、x和/或Z的時間過程中。在雙變量相關(guān)中，使用離群值權(quán)重的并集（例如，y或x中的離群值）。然而結(jié)果發(fā)現(xiàn)，擴展到多個通道的數(shù)據(jù)通常會導(dǎo)致太多的數(shù)據(jù)被向下加權(quán)，因為任何一個通道中的異常點會導(dǎo)致所有通道的時間點被向下加權(quán)。特別是，MVGC中的時滯要求包含該列矩陣的所有滯后列也被向下加權(quán)。

為了降低上述回歸模型中使用多個項的穩(wěn)健回歸方法的復(fù)雜性，研究者使用新息項，而不是單獨使用通道數(shù)據(jù)。穩(wěn)健MVGC的過程如下所述。首先，對通道數(shù)據(jù)（短距或長距）進(jìn)行預(yù)白化，并計算新息項。使用從新息項計算出的加權(quán)矩陣對數(shù)據(jù)進(jìn)行預(yù)加權(quán)。經(jīng)過預(yù)白化和預(yù)加權(quán)后，其他通道（Z）的新息項包括在限制性模型的回歸中，如下所示：

（限制性模型）

這里，εty是在預(yù)白化和預(yù)加權(quán)后的時間點t處通道y的新息項。εjz,t?i是時間點t-i處第j個成分的主成分權(quán)值。nPC表示使用的主成分的數(shù)量。

同樣，新息項也用于非限制性模型，而不是通道數(shù)據(jù)本身：

（非限制性模型）。

代碼實現(xiàn)

所有方法均在MATLAB中實現(xiàn)。利用MATLAB內(nèi)置函數(shù)和NIRS Brain AnalyzIR工具箱對相關(guān)方法進(jìn)行預(yù)處理和實現(xiàn)。使用自定義的MATLAB腳本實現(xiàn)了偏相關(guān)、格蘭杰因果關(guān)系和穩(wěn)健格蘭杰因果等方法。本文中使用的工具箱和代碼可在線訪問：https://github.com/pradlanka/rsfc-fnirs獲得。

結(jié)果

處理方法的比較

如預(yù)期那樣，隨著模擬的進(jìn)展和噪聲的增加，標(biāo)準(zhǔn)Pearson相關(guān)等基本模型很快就難以為繼了。例如，Pearson相關(guān)假設(shè)獨立的測量，這與加入的時間自相關(guān)噪聲相違背。如圖5(b)所示，與預(yù)白化和MVGC模型相比，不包括AR（預(yù)白化）濾波器的Pearson和穩(wěn)健Pearson相關(guān)模型在ROC曲線上的性能顯著更低。幾乎所有其他方法的AUC接近1。此外，如圖5(a)所示，Pearson和穩(wěn)健Pearson模型具有較大的假陽性率。其他分析模型的I型誤差圖更接近理想（對角線）。AR相關(guān)和穩(wěn)健AR相關(guān)都具有最接近理想的I型誤差控制。AR偏相關(guān)模型略微高估了測試統(tǒng)計的顯著性，而穩(wěn)健的AR偏相關(guān)過校正導(dǎo)致了對測試統(tǒng)計顯著性的低估（低于預(yù)期的假陽性率）。修正MVGC的兩個版本都傾向于略微低估p值，從而導(dǎo)致略高于預(yù)期的假陽性率。在這里，沒有顯示標(biāo)準(zhǔn)MVGC模型的結(jié)果，該模型不包括特定的零項。因為數(shù)據(jù)是在零滯后時模擬的，MVGC沒有建模，所以MVGC沒有顯示任何重要通道。這些結(jié)果支持了校正時間序列中自相關(guān)的必要性。

圖5.包含時間自相關(guān)的模擬數(shù)據(jù)在所有方法上的(a)I型誤差控制圖和(b)ROC曲線。虛線表示方法的非穩(wěn)健版本，而實線表示其穩(wěn)健版本。

第二組模擬（圖6）使用時間自相關(guān)和空間全局系統(tǒng)生理信號生成。正如預(yù)期的那樣，AR相關(guān)性與Pearson相關(guān)系數(shù)相結(jié)合，導(dǎo)致ROC曲線中假陽性增加[圖6(a)]和AUC降低[圖6(b)]。當(dāng)噪聲獨立于通道之間時，使用成對相關(guān)是有效的，但空間噪聲引入了全局假陽性連接和模型顯著性的大量過度報告。AR偏相關(guān)和修正的格蘭杰因果關(guān)系仍然表現(xiàn)良好。在圖5和圖6中，模擬數(shù)據(jù)中沒有添加運動偽影，因此穩(wěn)健和非穩(wěn)健方法在I型誤差控制曲線和ROC曲線中表現(xiàn)相似。盡管使用穩(wěn)健方法可能導(dǎo)致自由度損失，因為數(shù)據(jù)點的權(quán)重可能小于1，但時間點的數(shù)量仍然大得多（＞1000），因此自由度損失幾乎可以忽略不計。

圖6.包含時間自相關(guān)和全局信號的模擬數(shù)據(jù)在所有方法上的(a)I型誤差控制圖和(b)ROC曲線。虛線表示方法的非穩(wěn)健版本，而實線表示其穩(wěn)健版本。

在第三組模擬中，在包含時間自相關(guān)和全局信號的模擬數(shù)據(jù)中加入了以尖峰和位移偽影的頭動。如圖7所示，降低離群點權(quán)重的穩(wěn)健方法比其非穩(wěn)健方法的性能更好，具有穩(wěn)健的修正格蘭杰因果關(guān)系和AR偏相關(guān)，短距通道具有相似的AUC，但AR偏相關(guān)比修正MVGC具有更少的假陽性。值得注意的是，雖然穩(wěn)健的修正格蘭杰因果關(guān)系有效地減少了假陽性，但其仍然有遠(yuǎn)高于預(yù)期的假陽性率（0.05的預(yù)期率下，假陽性率為0.25）。

圖7.包含時間自相關(guān)、全局信號和頭動偽影信號的模擬數(shù)據(jù)在所有方法上的(a) I型誤差控制圖和(b)ROC曲線。

短距通道的有效性比較

本研究深入地探討這些模型，并比較了僅用短距通道、長距通道或同時使用短距通道和長距通道對白化偏相關(guān)方法性能的影響。圖8(a)顯示，在I型誤差控制圖中，僅使用長距通道是最差的，可能會導(dǎo)致檢驗統(tǒng)計量（相關(guān)系數(shù)）顯著性的高估，表明結(jié)果仍然保留了太多的全局連接。僅使用短距通道產(chǎn)生了最佳結(jié)果，盡管它略微低估了穩(wěn)健估計器的顯著性，并導(dǎo)致在估計中略微引入假陰性。最后，短距和長距通道的使用略差于單獨使用短距通道，但優(yōu)于僅使用長距通道。這可能是由于模型中的某些過擬合導(dǎo)致的。

圖8(b)顯示了附加運動偽影的模擬數(shù)據(jù)的比較。正如預(yù)期的那樣，穩(wěn)健的統(tǒng)計估計器大大優(yōu)于非穩(wěn)健的統(tǒng)計估計，僅用短距通道模型對這些偽影最敏感。與靜止模擬一致，僅使用穩(wěn)健的短距通道和短距和長距的穩(wěn)健版本給出了最佳結(jié)果。在穩(wěn)健模型中，僅使用長距通道仍然對測試統(tǒng)計的顯著性有大量過度報告。

最后，圖9顯示了，使用修正的MVGC模型代替AR偏相關(guān)性進(jìn)行了相同的分析。圖9(a)顯示了無運動模擬數(shù)據(jù)的結(jié)果。結(jié)果發(fā)現(xiàn)MVGC對模型中使用的通道不太敏感，短距、長距和甚至是兩種通道模型都具有相似的性能。據(jù)觀察，即使在靜止數(shù)據(jù)模擬的情況下，MVGC的穩(wěn)健版本對測試統(tǒng)計數(shù)據(jù)的顯著性也有大量的過度報告。研究者認(rèn)為，這是模型中關(guān)于估計有效自由度的未修正誤差。特別是，由于測量變量y、x和/或Z中的時間點都可能是運動引起的異常值，因此難以在穩(wěn)健估計器中定義統(tǒng)計異常值。該模型中使用的Tukey雙平方權(quán)重假設(shè)離群值是獨立的，但當(dāng)前的模型沒有有效地考慮多個通道的共享離群值時間點，因此自由度稍微高估。在圖9(b)中，針對添加運動偽影的數(shù)據(jù)模擬比較MVGC模型，結(jié)果發(fā)現(xiàn)穩(wěn)健方法優(yōu)于非穩(wěn)健方法。

建模滯后連接

修正后的MVGC同時考慮了零時滯和非零時滯。因此，當(dāng)僅使用零時滯連接模擬數(shù)據(jù)時，AR偏相關(guān)性被確定為最佳。標(biāo)準(zhǔn)MVGC（非零時滯）在前幾節(jié)中沒有出現(xiàn)，是因為它在零非零時滯模擬中沒有發(fā)現(xiàn)任何結(jié)果[圖10(a)]。同樣，當(dāng)數(shù)據(jù)僅以非零時滯項的關(guān)系進(jìn)行模擬時，AR偏相關(guān)模型失敗，如圖10(b)所示。修正后的MVGC是唯一一個在兩種情況下都有效的模型。

在圖11中，在添加了運動偽影的情況下，顯示出了相同的模式，其中ROC曲線對于零時滯模擬的情況下的穩(wěn)健AR偏相關(guān)和穩(wěn)健修正MVGC模型是最佳的，并且對于非零時滯模擬下的穩(wěn)健MVGC和穩(wěn)健修正的MVGC是最佳的。在這兩種情況下，修正后的MVGC都不如AR偏相關(guān)或MVGC好，因此如果知道存在哪種連接，則修正后的MVGC就不是首選。反之，這是兩種情況下唯一可行的模型。

靜息態(tài)fNIRS數(shù)據(jù)的結(jié)果

為了檢驗本研究的功能連接方法對實驗靜息態(tài)數(shù)據(jù)的適用性，研究者比較了Pearson相關(guān)、AR相關(guān)和AR偏相關(guān)的穩(wěn)健版本，以在24名參與者樣本上模擬通道之間的功能連接。圖12顯示了三種方法連接路徑的組水平結(jié)果（q＜0.001，F(xiàn)DR校正）。結(jié)果表明，與AR相關(guān)和Pearson相關(guān)相比，使用AR偏相關(guān)導(dǎo)致更稀疏的網(wǎng)絡(luò)，因為它校正了全局生理信號和自相關(guān)。使用Pearson相關(guān)性揭示了一個非常密集的連接網(wǎng)絡(luò)，但大多數(shù)連接可能是假陽性的。雖然沒有已知的地面真值或黃金標(biāo)準(zhǔn)的RS-fMRI網(wǎng)絡(luò)可用于定量比較這些方法，但AR偏相關(guān)可能是三種方法中假陽性率最低的一種方法。

結(jié)論

復(fù)雜的數(shù)值模擬模型在理論上可以獲取更多關(guān)于fNIRS時間通道之間關(guān)系的信息，但在實踐中很難求解。盡管修正的格蘭杰因果關(guān)系等方法可以對零時滯和時滯關(guān)系進(jìn)行建模，但當(dāng)連通性信息存在零時滯時，預(yù)白化偏相關(guān)性表現(xiàn)最佳，而當(dāng)時間序列之間存在時滯關(guān)系時，格蘭杰因果關(guān)系模型表現(xiàn)最佳。修正的MVGC在兩種情況下都表現(xiàn)良好，但在低信噪比情況下，自由度較低，靈敏度可能更低。同樣，帶通濾波雖然可能會提高信噪比和靈敏度，但也可能導(dǎo)致自由度損失和自相關(guān)增加，從而抵消增益。因此，建議仔細(xì)考慮預(yù)處理流程中的步驟，以最大限度地提高靈敏度，同時將假陽性率降低到期望水平。

原文：Correction of global physiology in resting-state functional near-infrared spectroscopy.

DOI: 10.1117/1.NPh.9.3.035003

小伙伴們點個“在看”，加

(星標(biāo))關(guān)注茗創(chuàng)科技，將第一時間收到精彩內(nèi)容推送哦～