基于 FieldTrip 的置换检验与统计绘图

一、脑电统计的核心矛盾：一次实验，上万次检验

模块 3 你画过两个条件的 ERP 差值图，肉眼看到 P3 更负；模块 4 你看过两块时频图，肉眼觉得 alpha 能量不同。但"肉眼显著"到"统计显著"之间隔着一道鸿沟，而且这道鸿沟比多数新手想象的深得多。

算一笔账：64 通道、500 Hz 采样、分析窗口 -0.2 到 1.0 秒（601 个时间点）。 如果你逐个"通道 × 时间点"做检验，一共要做 64 × 601 = 38464 次。 就算两个条件的数据完全来自同一个分布（真实差异为零），按每次 α = 0.05 算， 期望也会有约 1923 个"显著"点——它们全是噪声，却会画满你的显著图。 加上时频的频率维度，检验次数轻松冲上十几万。

四种常见应对方式放在一起看，就知道为什么脑电界最后集体选了最后一种：

Bonferroni 的问题在于它假设四万次检验彼此独立——但脑电相邻时间点的相关高达 0.9 以上， 相邻通道也高度相关。真正"独立"的检验远没有四万次。cluster-based 方法干脆承认 这种相关结构，把它变成优势：显著的效应从来不是孤立的一个点，而是连成片的时空区域。

二、置换检验的直觉：洗牌发牌

先把 cluster 放一边，讲它的地基——置换检验。想象你怀疑牌局有鬼： 某位玩家这几局赢得太多。怎么判断是实力还是运气？把整副牌重新洗 10000 次、重新发 10000 次， 数一数"纯运气"下能出现多大的领先——这就得到了零分布。 如果真实牌局的领先幅度超过了其中 99% 的洗牌结果，那基本不是运气能解释的。

换到脑电上完全一样：零假设是"条件标签与数据无关"——同一个被试的两组数据， 标签换一下同样说得通。于是：

1. 按真实标签算一次统计量（如配对 t）；
2. 随机打乱标签，重新算一次；
3. 重复 1000 次，得到零分布；
4. p 值 = 零分布中比真实值更极端的比例。数据自己告诉你什么是显著，不需要任何正态分布假设。

% 手工置换检验：20 个被试的条件差值（想象成 B-A 的 alpha 功率）
rng(1);                                 % 固定随机种子，结果可复现
diff = 0.5 + randn(20, 1);              % 真实均值 0.5
t_obs = mean(diff) / (std(diff) / sqrt(numel(diff)));   % 单样本 t

% 零假设：条件无差别 → 每个差值的正负号可以随便换
nPerm = 1000;
t_null = zeros(nPerm, 1);
for p = 1:nPerm
    flip = randi([0 1], numel(diff), 1) * 2 - 1;    % 随机生成 +1/-1
    d = diff .* flip;
    t_null(p) = mean(d) / (std(d) / sqrt(numel(d)));
end

% p 值 = 零分布中与观察值同样极端的比例
pval = sum(abs(t_null) >= abs(t_obs)) / nPerm;
fprintf('观察 t = %.2f，置换 p = %.3f\n', t_obs, pval);

% 画出零分布，标出真实值——显著 = 真实值落在分布的尾巴上
figure; histogram(t_null, 40); hold on;
xline(t_obs, 'r', 'LineWidth', 2);       % 需要 R2018b+

这个一维小例子就是 5.2 全部内容的缩影：脑电的真实场景只是把"一个数"换成 "通道 × 时间（× 频率）的矩阵"，把"打乱标签"交给 FieldTrip，其余逻辑分毫未变。 注意 p 值的分辨率由置换次数决定：1000 次置换时最小 p 值约为 1/1001， 想报告 p < 0.001 就得加码（见第八节）。

三、cluster-based 置换检验：四步拆解

现在把"多重比较"和"置换检验"组装起来，就是论文里最常见的 cluster-based permutation test：

① 逐点检验：每个 (通道, 时间点) 上算一次 t 值
② 阈值聚类：t 值超过 cfg.clusteralpha 的点，按"时间相邻 + 空间相邻"连成 cluster
③ 取统计量：对每个 cluster 求 t 值之和，只保留所有 cluster 中的最大值
④ 置换比较：打乱标签重复 ①-③，用"最大 cluster 统计量"的零分布给真实值算 p

可以把聚类想成把散落的珠子串成项链：过阈值的数据点是珠子，邻接关系是线， 一整条项链的"重量"（t 值之和）才是被检验的对象。一条又长又重的项链， 远比一颗孤零零的大珠子更难用噪声解释。

两个容易误解的参数，先划重点：

• cfg.clusteralpha（默认 0.05）只是珠子入串的门槛：阈值松则 cluster 大而少、紧则小而碎，它不影响最终的假阳性率——假阳性由置换次数和第④步的比较控制；
• cfg.minnbchan（常用 2）要求 cluster 至少含 2 个邻接通道，防止单个噪声通道自成一簇。

这套方法的完整理论出处是 Maris & Oostenveld（2007, Journal of Neuroscience Methods），也是 FieldTrip 官方教程反复引用的文献——写论文时值得放进参考文献列表。

四、邻接结构：ft_prepare_neighbours

第②步"空间相邻"需要一个邻接表：谁和谁算邻居。FieldTrip 用 ft_prepare_neighbours 生成， 常用两种方法：

% 方法一：模板法——用官方内置的标准邻接表（推荐，最省心）
cfg = [];
cfg.method   = 'template';
cfg.template = 'elec1020_neighb.mat';   % 10-20 系统；10-10 帽子用 elec1010_neighb.mat
cfg.channel  = 'EEG';
neighbours = ft_prepare_neighbours(cfg);

% 方法二：距离法——按电极物理距离定义（非标准电极帽时用）
cfg = [];
cfg.method        = 'distance';
cfg.neighbourdist = 0.04;               % 40 mm 以内算邻居（默认值）
neighbours = ft_prepare_neighbours(cfg, data);  % data 自带电极坐标

% 一定要画出来检查：孤岛通道进不了任何 cluster，统计力白白流失
cfg = [];
cfg.neighbours = neighbours;
cfg.layout = 'elec1005';
ft_neighbourplot(cfg, data);            % enableedit 设为 yes 还能交互增删邻接边

模板法的前提是电极名规范（又是模块 2.2 强调过的那件事：Fp1、Cz、O2……）； 距离法需要数据里带电极坐标（data.elec），40 mm 这个默认值对多数脑电帽都合理。 无论用哪种，都请用 ft_neighbourplot 画出来看一眼——邻接结构错了，cluster 就错了。

五、ft_timelockstatistics 完整实战

万事俱备。输入是 5.1 批处理保存的每被试、每条件 timelock 结构——先把它们合并成 保留个体值的组结构（统计需要每个被试的 ERP，不能只要组平均），再配统计：

%% 第 0 步：准备输入——保留每个被试个体值的组结构
cfg = [];
cfg.keepindividual = 'yes';
gaA = ft_timelockgrandaverage(cfg, tlkA{:});   % tlkA 来自 5.1 的批处理
gaB = ft_timelockgrandaverage(cfg, tlkB{:});

%% 第 1 步：设计矩阵——统计的"户口本"
subj = numel(tlkA);
design = zeros(2, 2*subj);
design(1,:) = [ones(1,subj) 2*ones(1,subj)];  % 第 1 行：条件标签（1=A，2=B）
design(2,:) = [1:subj 1:subj];                % 第 2 行：被试编号

%% 第 2 步：完整统计 cfg
cfg = [];
cfg.channel = 'EEG';
cfg.latency = [0 1];                 % 检验空间限定在刺激后 0-1 s（按先验假设圈，别圈全段）
cfg.method  = 'montecarlo';          % 蒙特卡洛置换检验
cfg.statistic = 'ft_statfun_depsamplesT';   % 被试内配对 t
cfg.correctm = 'cluster';            % cluster 校正解决多重比较
cfg.clusteralpha = 0.05;             % 聚类阈值（只影响 cluster 形状）
cfg.clusterstatistic = 'maxsum';     % cluster 统计量：簇内 t 值之和的最大值
cfg.minnbchan = 2;                   % cluster 至少含 2 个邻接通道
cfg.neighbours = neighbours;         % 空间邻接结构（上一节准备的）
cfg.numrandomization = 1000;         % 置换次数：论文级 1000 起步
cfg.tail = 0;                        % 双侧检验：不预设方向
cfg.design = design;
cfg.ivar = 1;                         % 条件在第 1 行
cfg.uvar = 2;                         % 被试在第 2 行

stat = ft_timelockstatistics(cfg, gaA, gaB);

%% 第 3 步：看看输出里有什么
stat.stat    % (通道 x 时间) 的 t 值
stat.prob    % 每个数据点的置换 p 值
stat.mask    % 显著性蒙版（0/1 逻辑矩阵），绘图直接用
stat.posclusters   % 正向 cluster 列表，每个元素的 .prob 是它的 p 值
stat.negclusters   % 负向 cluster 列表

cfg 字段一多容易晕，按"四件事"归类记忆：

设计矩阵是最常见的翻车点：ivar 指条件所在行，uvar 指被试所在行； 条件标签必须是 1 和 2，被试编号是 1 到 N 的整数。报错"design 与数据不匹配"时， 先数 design 的列数是否等于"被试数 × 条件数"。

六、时频统计：ft_freqstatistics 同套路

时频数据的统计零新知识：同一个 cfg 换个函数名，输入换成频域的组结构， 再多一个频率窗 cfg.frequency：

% 输入：5.1 批处理算好的每被试 TFR，先合并成保留个体值的组结构
cfg = [];
cfg.keepindividual = 'yes';
gaTFR_A = ft_freqgrandaverage(cfg, tfrA{:});
gaTFR_B = ft_freqgrandaverage(cfg, tfrB{:});

% 统计 cfg 与 ERP 版完全同一套路
cfg = [];
cfg.method = 'montecarlo';
cfg.statistic = 'ft_statfun_depsamplesT';
cfg.correctm = 'cluster';
cfg.clusteralpha = 0.05;
cfg.clusterstatistic = 'maxsum';
cfg.minnbchan = 2;
cfg.neighbours = neighbours;
cfg.numrandomization = 1000;
cfg.tail = 0;
cfg.design = design;       % 与 ERP 统计用的是同一个 design
cfg.ivar = 1;
cfg.uvar = 2;
cfg.channel   = 'EEG';
cfg.frequency = [4 30];    % 时频统计专属：限定检验的频率窗 (Hz)
cfg.latency   = [0 1];     % 限定时间窗 (s)

stat_tfr = ft_freqstatistics(cfg, gaTFR_A, gaTFR_B);

唯一的新东西是聚类多了一个频率维度：相邻频率bin的点也会被串进同一条项链 （时空谱聚类）。检验空间从"通道 × 时间"变成"通道 × 频率 × 时间"——这也是第一节那笔账 涨到十几万次的原因，也是为什么 cfg.frequency 和 cfg.latency 要按先验假设圈定、而不是全段全频扫。

七、统计绘图：ft_clusterplot 与显著蒙版

统计结果有两类画法。ft_clusterplot 一键出图：按时间点（或频率bin）画一排地形图， 把显著 cluster 高亮标号，适合快速浏览和组会汇报；显著蒙版叠加则把 stat.mask 作为透明度蒙版画在自己指定的图上，是论文配图的主力：

%% 方法一：ft_clusterplot 一键出图
cfg = [];
cfg.layout = 'elec1005';
cfg.alpha = 0.05;                % 只画 p 小于 0.05 的 cluster（上限 0.3）
ft_clusterplot(cfg, stat);       % 时频结果同理：ft_clusterplot(cfg, stat_tfr)

%% 方法二：差值地形图 + 显著蒙版（论文常用）
% 1) 算两条件的组平均差：x1 - x2
cfg = [];
cfg.operation = 'x1-x2';
cfg.parameter = 'avg';
gaDiff = ft_math(cfg, gaA, gaB);

% 2) 把置换检验的显著蒙版挂到差值结构上（0=不显著，1=显著）
gaDiff.mask = stat.mask;

% 3) 画图：不显著的电极自动变透明，一眼看出"效应长在哪里"
cfg = [];
cfg.layout = 'elec1005';
cfg.parameter = 'avg';
cfg.maskparameter = 'mask';     % 引用刚才挂上的蒙版字段
cfg.xlim = [0.3 0.5];            % 显示 300-500 ms 窗口
cfg.zlim = 'maxabs';
cfg.colorbar = 'yes';
ft_topoplotER(cfg, gaDiff);

时频版同理：用 ft_math 算两条件的 TFR 差值，把 stat_tfr.mask 挂上去， ft_singleplotTFR（单通道时频图）和 ft_topoplotTFR（频窗地形图） 都支持 cfg.maskparameter。图注务必写清蒙版的含义（如"透明区域未通过 cluster 校正的置换检验，α = 0.05"）。源空间的统计绘图（ft_sourceplot）属于溯源，模块 7 再讲。

八、结果怎么报告：置换次数、cluster 统计量、校正逻辑

审稿人最常追问的就是这几项，投稿前对照自查：

1. 统计类型与设计：依赖样本 t（被试内）还是独立样本 t（被试间）；检验方向（双侧/单侧）；
2. cluster 参数：cluster 形成阈值（clusteralpha）、cluster 统计量（maxsum）、最小邻接通道数；
3. 置换次数：numrandomization = 多少次。p 值分辨率 = 1/(次数+1)，报告 p < 0.001 至少要 1000 次；
4. 邻接定义：模板法（哪个模板）或距离法（多少 mm）；
5. 检验空间：纳入的时间窗、频率窗、通道集合——它必须由先验假设决定；
6. 结果表述：显著 cluster 的 p 值（保留三位小数）、簇的时空范围用"描述"口吻给出。

九、常见误区：cluster 检验不是免检金牌

1. 先看数据再定检验方向：画完图挑个方向明显的效应再跑单侧检验——这是变相的两次检验，假阳性翻倍。没有预注册的把握，一律双侧；
2. 把 cluster 当定位证据：检验的结论是"存在效应"，不是"恰好这几个电极、这几毫秒显著"。显著 cluster 的边界受噪声影响很大，描述可以，推断不行；
3. 置换次数太少：100 次置换时 p 值分辨率只有 0.01，"p = 0"其实只是小于 0.01。论文级一律 1000 次起步，谨慎者 5000；
4. 用错统计函数或指错行：被试内设计必须用 depsamplesT（被试间用 indepsamplesT），ivar/uvar 指错行会直接报错或悄悄算错；
5. 邻接结构不检查就开跑：改名后的电极进不了模板、孤岛通道进不了 cluster——统计力在你看不见的地方流失。

下一模块（模块 6）我们从"每个通道自己的活动"走向"通道之间的对话"——功能连接与相位同步。 你会发现本课的 cfg 哲学一路通用：连接指标用 ft_connectivityanalysis 算， 统计照旧用置换检验。