采用基于频率簇(Cluster)的置换检验(Permutation)方法选取感兴趣频段

1 min read

作者:北京师范大学 龙宇航,[email protected]
代码来源(见本页底部):周思远

在使用wtc计算脑间神经同步后,我们需要在多个频率段、多个通道组合上对神经同步值进行统计检验,因此当进行频段选择时,面临多重比较的问题。为了解决多重比较的问题,可以采取基于参数或非参数检验的多重比较矫正的方法。由于基于非参数检验的多重比较矫正对数据的分布形态没有严格要求,因此具有更广泛的应用场景 (Maris and Oostenveld, 2007)。本文即介绍基于随机置换的非参数检验的方法 (Zheng et al., 2020; Long et al., 2021)。

在寻找感兴趣的效应时,我们采取了基于频率簇(Cluster)的方法,即在频率方向寻找连续显著的Cluster,该方法比基于最强效应点的方法具有更为优秀的抗噪音能力。值得注意的是,我们并没有沿着通道的方向去寻找连续显著的通道簇,这是因为沿着通道方向寻找Cluster容易受到生理噪音的影响。

下面进入具体的实操部分。假设本例招募了22对组1被试及22对组2被试,每对被试分别进行3种条件的任务,因此本例是2(组别,被试间因素)*3(条件,被试内因素)的实验设计。本例对神经同步值进行2*3的混合方差分析,并关注交互作用。

具体来讲,进行置换检验需要进行以下几个步骤:1. 重采样;2. 对随机样本进行计算及统计;3. 计算真实样本的统计量;4. 真实样本与随机样本的对比。下面依次进行介绍。

1. 重采样

首先,在每种实验条件下都进行重采样。置换检验采取的是不放回抽样,一般来讲至少需要1000个随机样本(不放回随机采样)。以组1在条件1下的重采样为例,本例中的每对被试由一男一女组成,表1中Original Sample (OS) 指组1在条件1时的配对情况,即原始配对。而每个Permutation Sample (PS) 中的22对被试并不是真实被试,如在PS1的第一对被试中,第7号男性在条件1下的血氧信号与第20号女性在条件1下的血氧信号进行了匹配,因此并不是真实互动条件下的神经信号。同样地,对每个水平都进行重采样。

2. 对随机样本进行计算及统计

2.1 对每个随机配对进行WTC的计算得到神经同步值。

2.2 对每个随机样本(PS)的神经同步值进行所关心的统计分析。在本例中我们得到交互作用的P值矩阵,其中一个维度为频率,一个维度为通道。假设第1个随机样本我们得到了以下统计P值矩阵(表2),表中显示有400个通道,以及109个频率点。

表格 2 某个随机样本的P值矩阵,其中行代表频率,列代表通道。

2.3 按频率方向挑选最长连续的P < 0.05的Cluster,得到该Cluster所对应的频段以及通道。具体来讲,在表3中标红的为P < 0.05的值,我们可以看到被框起的红色部分为表中按频率方向(纵向)最长的Cluster,该Cluster所对应的频率段为104 – 109,通道组合为6。至此,我们已经找到了该随机样本最长的Cluster及其所对应的统计值、对应的通道、以及对应的频率段。寻找最长连续P < 0.05 Cluster的代码见本页底部。

表格 3 某个随机样本的P值矩阵,绿色背景部分为最长Cluster所对应的频段及通道。

2.4 得到最终用于画图的统计值。有3种方法可以得到该值,如下图,第一种为将该Cluster所对应的F值平均。第二种为将该Cluster所对应的F值求和。第三种为将该Cluster所对应的频段以及通道的WTC值进行平均,对平均后的WTC值再次进行统计分析。

图 1 三种用于得到最终统计值的方法。

2.5 对每个随机样本(PS)都进行上述2.1-2.4的操作,最终得到1000个随机样本(PS)统计值的95百分位数,如下图(图中为示意图,统计值与上表并不对应),图中横轴代表最终用于画图的F值,纵轴代表样本数量。在计算完第一个随机样本后,我们发现其F值为8.97(图2中左上角子图);依此类推,当计算完1000个随机样本(PS)后,我们能够画出1000个随机样本的分布情况(图2中右下角子图),并根据该分布得到大于95百分位数区域(灰色区域)。

图 2 不同随机样本数目的统计值分布情况。其中右下子图表示1000个随机样本(PS)的分布情况,其中的阴影表示位于1000个随机样本(PS)统计值分布的95百分位数以上。

3. 计算真实样本的统计量

3.1 对真实样本进行WTC计算。

3.2 对真实样本的神经同步值进行统计分析得到二维的统计值矩阵(一个维度为频率一个维度为通道)。

3.3 按频率方向挑选所有连续的P < 0.05的Cluster所对应的频段以及通道。

3.4 将每一个符合要求的Cluster所对应的F值进行平均或求和或对应的WTC值进行平均后再次进行统计分析。

4. 真实样本与随机样本的对比

看真实样本中哪个cluster能够落在灰色区域。在本例中,红色线条为真实样本的某一个Cluster所对应的F值,落在了灰色区域内,表明在0.05水平显著。至此,我们已基于频率簇(Cluster)的置换检验(Permutation)方法选取出了感兴趣的频段。

图 3 真实样本中的Cluster高于随机样本95百分位数。红色线表明真实样本所代表的统计值,落在了阴影区域内。

参考文献:
1. Long Y, Zheng L, Zhao H, Zhou S, Zhai Y, Lu C (2021) Interpersonal Neural Synchronization during Interpersonal Touch Underlies Affiliative Pair Bonding between Romantic Couples. Cerebral cortex 31:1647-1659.
2. Maris E, Oostenveld R (2007) Nonparametric statistical testing of EEG- and MEG-data. J Neurosci Methods 164:177-190.
3. Zheng L, Liu W, Long Y, Zhai Y, Zhao H, Bai X, Zhou S, Li K, Zhang H, Liu L, Guo T, Ding G, Lu C (2020) Affiliative bonding between teachers and students through interpersonal synchronisation in brain activity. Social cognitive and affective neuroscience 15:97-109.

MatLab 代码:

function T_selected = bandSelect(pmap, threshold_peak, threshold_bound)
% Dependency: findContineElement
% Input: pmap: the p-value map (fs x ch)
%        threshold_peak: the threshold of the peak point
%        threshold_bound: the threshold to determine adjacent points around
%        peak point
% Output: T_selected: first column is channel pair ID, second column is fsband
% By Siyuan Zhou,2020/10
% Bug Fix: 2021/3

[t2_pos1,t2_pos2]=find(pmap<=threshold_peak);
source_x = t2_pos2;
source_y = t2_pos1;
selected = {};
value = {};

for ii = 1:size(source_x)
    
    pos0 = source_y(ii);
    raw = pmap(:,source_x(ii));
    select_chwise = findContinueElement(raw, pos0, threshold_bound);
    selected{ii}=select_chwise;
    value{ii} = raw(select_chwise);
end

T = table(source_x,selected',value','VariableNames',{'ch' 'fs' 'value'});

% convert fs and value to string to remove duplicated data
tempT = T(:,{'ch','fs'});
tempT.fs = cellfun(@num2str,T.fs,'UniformOutput',false);
[B, ci]=unique(tempT(:,[1 2]));

T_selected= T(ci,:);

    function continual_set = findContinueElement(raw, pos0, threshold)
        % find continual elements in a vector based on a certain threshold
        %   Input: raw: the vector to find continual elements
        %          pos0: the start position
        %          threshold: to determine whether an element belongs to the continual vector
        %   Output: continual_set: the positions of continual elements
        % By Siyuan Zhou, 2020/10
        
        pos = pos0;
        idx = 1;
        pos_idx = [];
        
        while(pos >0 & pos<=numel(raw) & raw(pos)<threshold)
            pos_idx(idx) = pos;
            pos = pos +1;
            idx = idx+1;
        end
        pos = pos0 - 1;
        
        while(pos >0 & pos<=numel(raw) & raw(pos)<threshold)
            pos_idx(idx) = pos;
            pos = pos -1;
            idx = idx+1;
        end
        
        continual_set = sort(pos_idx);
    end
end



写作助手,把中式英语变成专业英文


Want to receive new post notification? 有新文章通知我

第五十八期fNIRS Journal Club视频 王硕教授团队

Youtube: https://youtu.be/Qx9zMaIzAwo 优酷:https://v.youku.com/v_show/id_XNjQ0MzM3MDE1Mg==.html 理解噪音中的
Wanling Zhu
14 sec read

第五十八期fNIRS Journal Club通知2024/12/07, 10am 王硕教授团队

理解噪音中的言语对老年听力损失患者来说是一个重大挑战。来自首都医科大学附属北京同仁医院耳鼻咽喉科研究所王硕教授团队的助理研究员王松建将为大家介绍他们采用同步EEG-fNIRS技术,从神经与血流动力学两
Wanling Zhu
10 sec read

第五十七期fNIRS Journal Club视频 王欣悦博士

Youtube: https://youtu.be/vyo-kECC2Ps 优酷:https://v.youku.com/v_show/id_XNjQzNTA0ODIwMA==.html 肢体语言——
Wanling Zhu
20 sec read

Leave a Reply

Your email address will not be published. Required fields are marked *