ERP 统计与报告写作

%% 配对 t 检验：同一批被试、两个条件
% ampA / ampB：每个被试一个数的平均幅值（如顶叶 ROI 250-350ms）
[h, p, ci, stats] = ttest(ampA, ampB);
fprintf('t(%d) = %.2f, p = %.4f\n', stats.df, stats.tstat, p);
fprintf('均值差 95%% CI: [%.2f, %.2f] uV\n', ci(1), ci(2));

% 效应量：配对版 Cohen d = 差值均值 / 差值标准差
d = mean(ampA - ampB) / std(ampA - ampB);
fprintf("Cohen's d = %.2f\n", d);

%% 重复测量 ANOVA：三个及以上条件（条件为被试内因素）
% ampA / ampB / ampC：每个条件一个 nSub x 1 的向量
nSub = numel(ampA);
y    = [ampA; ampB; ampC];                     % 所有观测摞成一列
cond = [repmat({'A'}, nSub, 1); repmat({'B'}, nSub, 1); repmat({'C'}, nSub, 1)];
sub  = repmat((1:nSub)', 3, 1);                % 每个被试出现三次

% 被试设为随机因素：条件 x 被试交互项自动成为正确的误差项
[pVal, anovaTbl] = anovan(y, {cond, sub}, 'model', 2, ...
                          'random', 2, 'varnames', {'Cond', 'Sub'});
% 看 anovaTbl 中 Cond 一行的 p 值；显著后做事后两两配对比较（记得再校正）

%% 通道 x 时间的大规模检验 + FDR 校正
% ERP_A / ERP_B：被试 x 通道 x 时间 的矩阵（每人每条件一条 ERP）
[nSub, nChan, nTime] = size(ERP_A);
pmap = ones(nChan, nTime);                     % 每个通道-时间点的 p 值
for ch = 1:nChan
    for tt = 1:nTime
        [~, pmap(ch, tt)] = ttest(squeeze(ERP_A(:, ch, tt)), ...
                                   squeeze(ERP_B(:, ch, tt)));
    end
end

%% Benjamini-Hochberg FDR：手工实现，一共几行
pv = sort(pmap(:));  m = numel(pv);  q = 0.05; % 目标 FDR 水平
thresh = q * (1:m)' / m;                       % 第 i 小 p 值对应的阈值
ok = pv <= thresh;
if any(ok)
    pCrit  = pv(find(ok, 1, 'last'));          % 能通过 FDR 的最大 p 值
    sigMap = pmap <= pCrit;                    % 校正后仍显著的通道-时间点
else
    pCrit = 0;  sigMap = false(size(pmap));    % 一个显著点都没有
end
fprintf('FDR(q=0.05) 临界 p 值: %.4f，显著点占比: %.1f%%\n', ...
        pCrit, 100 * mean(sigMap(:)));

%% ERP 指标与行为指标的相关（被试级）
% amp：每个被试的 P300 平均幅值；rt：每个被试的平均反应时 (ms)
[r, p]    = corr(amp, rt);                     % 默认皮尔逊（线性关系）
[rho, p2] = corr(amp, rt, 'Type', 'Spearman'); % 秩相关：对离群值和非线性更稳健
fprintf('Pearson  r   = %.2f, p = %.3f\n', r, p);
fprintf('Spearman rho = %.2f, p = %.3f\n', rho, p2);

% 散点图必画：小样本里一个离群点就能"造出"一个显著相关
figure; scatter(rt, amp, 'filled'); hold on;
b  = polyfit(rt, amp, 1);                      % 一元线性拟合
xv = linspace(min(rt), max(rt), 100);
plot(xv, polyval(b, xv), 'r-');
xlabel('平均反应时 (ms)'); ylabel('P300 平均幅值 (\muV)');
title(sprintf('被试级相关（n = %d）', numel(amp)));