letswave教程：腦電資料的時頻分析組平均與統計分析

Hello，

這裡是

行上行下

，我是

喵君姐姐

~

上一期，我們學習了

如何用Letswave進行資料的預處理和ERP分析

，包括letswave7的安裝、資料集匯入；預處理中常用的幾種降噪方法；ICA分解；ERP分析。

這一期，我們就來教大家學習如何用Letswave對單個主題進行

時頻分析

以及對多個主題進行

平均

和

統計分析

。

下一期，將進行

圖形生成

和

指令碼處理

，還請持續關注我們，敬請期待喲~

溫馨提示

：以下閱讀大概需要8分鐘左右，若是需要相關軟體與其他相關資料，也可直接在後臺回覆“

letswave

“獲取喲~

首先，我們接著上一期內容繼續對單主題進行分析。時頻分析描述了非平穩EEG訊號隨時間增加在頻域上的變化。Letswave7提供了兩種時頻分析方法，分別是

STFT

（短時傅立葉變換）和

CWT

（連續小波變換）。

STFT的原理是

：把整個時域過程分解成無數個等長的小過程，每個小過程近似平穩，再進行傅立葉變換，便可得到某個時間點出現的頻率。

CWT的原理是

：將傅立葉變換中無限長的三角函式基換成有限長的會衰減的小波基。這樣不僅能夠獲取頻率，還可以定位到時間。

具體連續小波變換和短時傅立葉變換的區別，請見：

傅立葉分析和小波分析的通俗含義

。

通常，我們可以在時域平均之前或之後進行時頻分析。對於在時域平均之前的時頻分解，可以觀察到非鎖相活動，但訊號會產生更多的噪聲。時域平均後的時頻分解通常具有較好的訊號噪聲比，但會丟失非鎖相活動。

1.時頻分析

1.1

連續小波變換

在這裡，我們將展示在時域平均之前，使用CWT進行時頻分析。

首先，

選擇資料集

“bl reref ep_S 9 sp_filter ica chan_interp butt sel_chan sub093”

和

“bl reref ep_S 10 sp_filter ica chan_interp butt sel_chan sub093”

，然後在選單中單擊”

Plugins-> my_tfa-> Averaged CWT"

。

其次，

將dafault引數設定保留在批處理模組中，然後單擊批處理模組底部的

“Run”

按鈕以進行時頻分析。

最後，

名稱為

“avg cwt bl reref ep_S 9 sp_filter ica chan_interp butt sel_chan sub093”

和

“avg cwt bl reref ep_S 10 sp_filter ica chan_interp butt sel_chan sub093”

的兩個新資料集將出現在管理器模組的資料列表中。

由於時頻分析非常耗時，並且還需要較大的儲存空間，因此記憶體較小的計算機很容易出現“記憶體不足”的錯誤。因此，上述步驟實際上是將時頻分析和求平均的步驟結合在一起，以節省儲存空間。

1.2 基準線校正

在進行時頻分析之後，我們還需要對時頻域進行基線校正。

首先，

選擇資料集

"avg cwt bl reref ep_S 9 sp_filter ica chan_interp butt sel_chan sub093 ”

和

“ avg cwt bl reref ep_S 10 sp_filter ica chan_interp butt sel_chan sub093”

，然後在選單中單擊

“Process- >Baseline corrections->Baseline correction ”

。

其次，

建議在時頻域中縮小基線校正的間隔。因此，我們將基準設定為

-0.75

到

-0.25s

。在時頻域中，可以使用幾種方法（例如減法，ER百分比，除法和 z分數）作為基線校正。在這裡，我們使用預設方法

減法

。單擊批處理模組底部的

“Run”

按鈕以進行基線校正。

最後，

名稱為

“bl avg cwt bl reref ep_S 9 sp_filter ica chan_interp butt sel_chan sub093”

和

“bl avg cwt bl reref ep_S 10 sp_filter ica chan_interp butt sel_chan sub093”

的兩個新資料集將出現在管理器模組的資料列表中。

1.3 檢視結果

首先，

選擇資料集

“bl avg cwt bl reref ep_S 9 sp_filter ica chan_interp butt sel_chan sub093”

和

“bl avg cwt bl reref ep_S 10 sp_filter ica chan_interp butt sel_chan sub093”

。單擊右鍵選單中的

View

，我們可以在多檢視器中檢視波形的時頻域結果。

其次，

同時選擇兩個資料集，選擇通道Pz，並將顏色範圍設定為

-10

至

10

；開啟

topograph

並將

cursor

設定為

x = 0.35，y = 3

，我們可以觀察到topograph上的鎖相P300活動。

另外，

將

cursor

設定為

x = 1，y = 4

，可以觀察到事件相關去同步（ERD）活動，這是非鎖相響應。

到目前為止，我們已經完成了分析單個被試時頻分析的教程。

1. 總體平均

1.1

降低取樣率

考慮到要上傳和下載的示例資料集的大小，我們將資料集的取樣範圍從1000Hz降至250Hz。

首先，

選擇資料集

“ avg bl reref ep_S 9 sp_filter ica chan_interp butt sel_chan sub093 ”

和

“ avg bl reref ep_S 10 sp_filter ica chan_interp butt sel_chan sub093 ”

，然後選擇

“ Edit->Resample Signals->Downsample signals（integer radio)”

。

其次，

在批處理模組中，將向下取樣率設定為

4

。單擊按鈕“

Run"

以進行向下取樣。

最後，

名稱為

“ds avg bl reref ep_S 9 sp_filter ica chan_interp butt sel_chan sub093”

和

“ds avg bl reref ep_S 10 sp_filter ica chan_interp butt sel_chan sub093”

的兩個新資料集將出現在管理器模組的資料列表中。

1.2 重新命名

在單個被試的分析中，對於每個步驟，將使用字首生成新的資料集。雖然包含了基本資訊，但名稱過於冗長。因此，在進行組分析之前，我們需要重新命名資料集。

首先，

選擇資料集

“ ds avg bl reref ep_S 9 sp_filter ica chan_interp butt sel_chan sub093 ”

，在右鍵單擊的選單中按

Rename

。將資料集

“ ds avg bl reref ep_S 9 sp_filter ica chan_interp butt sel_chan sub093 ”

重新命名為

“ Sub093 P300 target ”

。

同樣，

將資料集

“ avg bl reref ep_S 10 sp_filter ica chan_interp butt sel_chan sub093 ”

重新命名為

“ Sub093 P300 nontarget ”

。

1.3

資料集複製和合並

對於組分析，我們需要將時域分析結果複製到新資料夾中。例如，資料夾rawdata2中的資料集是所有93個被試的時域分析結果。解壓縮

rawdata2.zip

檔案，在93 * 2檔案中有93個數據集。

在單個被試的分析中，結果是針對目標和非目標時期的所有試次的平均值。對於組分析，我們首先需要將93個被試的資料集合併為一個數據集，其中新資料集中的每個被試都被視為一個epoch。

首先，

在管理器模組左側的Include列表框中選擇標籤

“ nontarget ”

，然後在管理器模組右側的資料集列表框中選擇所有資料集。選擇

“Edit->Arrange signals->Merge dataset epochs,channels,indexes”

。

其次，

將預設設定保留在批處理模組中，然後單擊

“ Run ”

按鈕以在非目標條件下合併所有選定的資料集。

以同樣的方式，我們可以在目標條件下合併所有資料集，從管理器模組左側的Include列表框中選擇標籤

“ target ”

。

最後，

產生兩個新的資料集，

“ merge_epoch Sub001 P300 nontarget ”

和

“ merge_epoch Sub001 P300 target ”

，每個都有93個epoch，其中包含93個被試在“ target”和“ nontarget”條件下的時域分析結果。

1.4 平均

首先，

在管理器模組左側的Include列表框中選擇標籤

“ merge_epoch ”

，然後在管理器模組右側的資料集列表框中選擇所有資料集。與單個被試的平均操作類似，選擇

“Process->Average->Compute average,std,median across epochs”

。

其次，

將預設設定保留在批處理模組中，然後單擊

“ Run ”

按鈕以對兩個資料集中的所有epoch進行平均，這實際上是在目標和非目標條件下對所有主題進行的總和。

最後，

兩個名為

“ avg merge_epoch Sub001 P300 nontarget ”

和

“ avg merge_epoch Sub001 P300 target ”

的新資料集”將出現在管理器模組的資料列表中。

1.5 檢視結果

首先，

要檢視總體平均值的結果，請選擇兩個資料集

“ avg merge_epoch Sub001 P300 nontarget ”

和

“ avg merge_epoch Sub001 P300 target ”

，然後在右鍵選單中單擊

View

。

其次，

在波形的多畫面檢視器中，選擇通道Pz，以與單個主題相同的順序檢視結果，我們需要更改資料集的順序。選擇資料集

“ avg merge_epoch Sub001 P300 target ”

，然後單擊

“dataset”

按鈕。

接著，

選擇兩個資料集，然後啟用工具欄中的

enable interval selection

按鈕。選擇

0.2

至

0.7s

的間隔引數。

然後，

單擊

“ Table ”

按鈕，彈出用於統計最大值、最小值和平均值的表格，其中我們發現P300的峰值為11。22 V，發生時間在332ms。

最後，

在工具欄中將

cursor

位置設定為

0.332

。開啟topograph。在目標和非目標條件下，P300的峰值時間點都可以看到P300實驗的總體平均分佈。

2.統計分析

有兩種進行統計推斷的方法，分別是

假設驅動

和

資料驅動

。

通常在P300分析中使用假設驅動，我們根據對資料的觀察或先前的經驗提出一個假設，例如假設從預定間隔0。2到0。7s的ERP的最大值或平均值在目標和非目標條件之間有所不同。

因此，我們將進行特定數量的t檢驗或ANOVA來檢驗該假設。假設驅動是基於EEG分析探索神經心理學機制的有效方法。但是，在假設驅動的方法中，結果在很大程度上取決於研究人員的主觀經驗。

資料驅動透過比較詳盡的時間、頻率和空間點上的ERP，為EEG分析提供了更為客觀的方法。此外，逐點分析可以準確說明發生效果的時間和地點。但是，隨著資料驅動中假設檢驗數量的急劇增加，我們必須面對多重比較問題中的

Family-Wise錯誤率（FWER）

。

2.1假設驅動

在本文中，我們透過以下假設顯示了統計推斷：

在目標條件和非目標條件下，ERP在通道Pz處從預定間隔0.2到0.7s的最大值或平均值會有所不同。

（

PS:

常常根據波形圖選時間視窗，地形圖選點，以及前人文獻共同參考）

對93名受試者進行配對樣本t檢驗。因此，對於目標和非目標條件，我們首先需要在通道Pz的0。2s至0。7s區間內獲取最大值和平均值。

首先，

選擇資料集

“merge_epoch Sub001 P300 nontarget ”和“ merge_epoch Sub001 P300 target ”

，然後單擊右鍵選單中的

View

。

其次，

在多畫面的波形中，選擇資料集

“merge_epoch Sub001 P300 target”

和通道的

Pz

，設定疊加波為

epochs

，開啟

enable interval selection

，設定Explore interval為

0.2

至

0.7s

。

接著，

按下

“Table”

按鈕，我們可以在表格中獲得最大值和平均值的統計資訊，並將這些資料複製到外部軟體（例如Excel或SPSS），以進行以下統計分析。

同樣，

選擇資料集

“ merge_epoch Sub001 P300 nontarget ”

，然後按

“Table”

按鈕，在另一個數據集上拾取相同的值，並將其複製到外部軟體。

最後，

我們可以進行統計推斷。以軟體Excel為例。將

目標條件

下的

最大值

和

平均值

複製到A和B列，將

非目標條件

下的

最大值

和

平均值

複製到C和D列。

在E1項中，輸入

“ = T.TEST（A：

A，C：

C，2,1）”

表示配對樣本t檢驗的最大值。結果p = 6。80 * 10 ^ -40表示，對於通道Pz上0。2至0。7 s間隔內的EEG訊號最大值，目標和非目標條件之間存在顯著差異。

類似地，在項E2中輸入

“ = T.TEST（B：

B，D：

D，2,1）”

，以對平均值進行配對樣本t檢驗。結果p = 3。2 * 10 ^ -28 表示對於通道Pz上0。2至0。7 s的EEG訊號平均值，目標和非目標條件之間存在顯著差異。

2.2資料驅動

2.2.1逐點T檢驗

對於資料驅動的測試，我們需要對每個時間通道點進行t檢驗，這將導致嚴重的多重比較問題。（

PS:

不過這個也是一種可以快速找到時間窗的方法）

例如，

首先，

選擇資料集

“ merge_epoch Sub001 P300 nontarget ”和“ merge_epoch Sub001 P300 target ”

，然後單擊

“ Process->Compare two datasets （paired sample/ two sample t-test)”

。

然後，

將dafault設定保留在批處理模組中，單擊”

Run"

按鈕以在管理器模組中獲取逐點ttest結果

“ ttest merge_epoch Sub001 P300 target ”

。

最後，

單擊右鍵選單中的

View

在

Pz

通道上觀看結果，我們將

Index

設定為

p-value

。由於我們將有效水平設定為α= 0。05，因此我們可以將Y軸設定為

0

到

0.05

，以觀察p值小於0。05的間隔。

可以發現，除了大約0。2-0。7秒的主要簇以外，在其他時間間隔內甚至在刺激之前還存在幾個簇。因此，顯然FWER很高。由於測試中有750個時間點和63個通道，因此透過Bonferroni方法校正的α值為0。05 / 750/63 = 10 ^ -6。

2.2.2 基於群集的置換測試

基於以上結果，我們需要運用多重比較校正進行基於群集的置換測試。Letswave7提供了兩種基於群集的置換測試的方法，即

單感測器分析

和

多感測器分析

。

單感測器分析

對於基於群集的置換測試，我們可以分別檢測每個通道的群集。

首先，

在逐點t檢驗中，選擇資料集

“ merge_epoch Sub001 P300 nontarget ”

和

“ merge_epoch Sub001 P300 target ”

，然後單擊

“ Process->Compare two datasets （paired sample/ two sample t-test）”

。

其次，

為了獲得更精確的結果，將

置換數

設定為

20000

，這將使計算變得很耗時。單擊

“ Run ”

按鈕，此過程需要10分鐘以上的時間。

然後，

結果將以相同的名稱

“ ttest merge_epoch Sub001 P300 target ”

儲存在資料集中。

最後，

透過單擊右鍵選單中的

View

，開啟資料集

“ ttest merge_epoch Sub001 P300 target ”

。未校正的結果與逐點t檢驗中的結果完全相同。將

Index

設定為

“ cluster p-value ”

，並將

cursor

設定為

0.332

，這是P300的峰值。

可以發現，在基於群集的置換測試之後，僅保留了通道Pz上的主群集，其他群集均被排除。

多感測器分析

首先，

選擇資料集

“ merge_epoch Sub001 P300 nontarget ”

和

“ merge_epoch Sub001 P300 target ”

，然後單擊

“Process->Compare two datasets （paired sample/ two sample t-test)”

。並將

置換數

設定為

20000

。

對於多感測器分析，我們將閾值設定為

0.16

。單擊

“ Run ”

按鈕，計算可能需要1個小時以上。

然後，

具有相同名稱

“ ttest merge_epoch Sub001 P300 target ”

的結果將覆蓋原始資料集。

最後，

用相同的方法在單感測器分析中觀察結果。可以發現，在通道Pz上保留了更多的時間間隔。這是因為在時空聯合域中，它們被視為同一群集。

總結：

到目前為止，我們已經展示了從資料匯入，單個被試和組的ERP以及時頻分析，最後到統計分析的整個過程。