Cluster equivalence check MNE-Python<->FieldTrip

cluster_equiv_ft.m

% FieldTrip reference runs for MNE-Python's mne/stats/tests/test_cluster_equiv.py
% (see mne-tools/mne-python#12663).
%
% The input file cluster_equiv_data.mat is produced by the scenario builders in
% test_cluster_equiv.py (exported with scipy.io.savemat); each scenario holds
% `data` (n_rows x n_chan x n_times), per-row factor labels, and (for
% within-subject designs) 1-based subject indices. This script runs the
% matching FieldTrip cluster-based permutation test for each scenario with
% ft_timelockstatistics (cfg.method = 'montecarlo', cfg.correctm = 'cluster',
% cfg.clusterstatistic = 'maxsum') and saves/prints the observed stat map and
% per-cluster statistics for hard-coding into the Python test.
%
% Tested with fieldtrip-20260812 on MATLAB R2026a.

ft_path = fullfile(getenv('HOME'), 'Desktop', 'fieldtrip-20260812');
data_file = 'cluster_equiv_data.mat';
out_file = 'cluster_equiv_results.mat';

addpath(ft_path);
ft_defaults;
S = load(data_file);
results = struct();

% ---------------------------------------------------------------- one_sample
% pre-subtracted paired diffs: depsamplesT against all-zero partner condition
sc = S.one_sample;
[tls, zero_tls] = deal(make_timelocks(sc), make_zero_timelocks(sc));
n = numel(tls);
cfg = base_cfg(sc, 0);
cfg.statistic = 'depsamplesT';
cfg.design = [1:n, 1:n; ones(1, n), 2 * ones(1, n)];
cfg.uvar = 1;
cfg.ivar = 2;
cfg.numrandomization = 'all';
results.one_sample = run_and_report('one_sample', cfg, [tls, zero_tls]);

% -------------------------------------------------------------------- paired
% 2 conditions x 12 subjects: depsamplesT (a - b)
sc = S.paired;
tls = make_timelocks(sc);
[~, ~, cond] = unique(sc.factor1);  % alphabetical -> a=1, b=2 (matches MNE)
cfg = base_cfg(sc, 0);
cfg.statistic = 'depsamplesT';
cfg.design = [sc.subject; cond'];
cfg.uvar = 1;
cfg.ivar = 2;
cfg.numrandomization = 'all';
results.paired = run_and_report('paired', cfg, tls);

% ----------------------------------------------------------------- between_t
% 2 independent groups (unequal n): indepsamplesT, and indepsamplesF
sc = S.between_t;
tls = make_timelocks(sc);
[~, ~, grp] = unique(sc.factor1);  % ctrl=1, pat=2
cfg = base_cfg(sc, 0);
cfg.statistic = 'indepsamplesT';
cfg.design = grp';
cfg.ivar = 1;
cfg.numrandomization = 20000;
results.between_t_T = run_and_report('between_t (indepsamplesT)', cfg, tls);

cfg = base_cfg(sc, 1);
cfg.statistic = 'indepsamplesF';
cfg.design = grp';
cfg.ivar = 1;
cfg.numrandomization = 20000;
results.between_t_F = run_and_report('between_t (indepsamplesF)', cfg, tls);

% ------------------------------------------------------------- between_anova
% 3 independent groups (unequal n): indepsamplesF
sc = S.between_anova;
tls = make_timelocks(sc);
[~, ~, grp] = unique(sc.factor1);
cfg = base_cfg(sc, 1);
cfg.statistic = 'indepsamplesF';
cfg.design = grp';
cfg.ivar = 1;
cfg.numrandomization = 20000;
results.between_anova = run_and_report('between_anova', cfg, tls);

% ------------------------------------------------------ rm_anova_interaction
% 2x2 within-subject interaction: depsamplesT on the per-subject double
% difference (a1b1 - a1b2) - (a2b1 - a2b2) against zero; the repeated-measures
% interaction F equals this t squared.
sc = S.rm_anova_interaction;
n_sub = max(sc.subject);
dd = zeros(n_sub, size(sc.data, 2), size(sc.data, 3));
sgn_map = containers.Map( ...
    {'a1_b1', 'a1_b2', 'a2_b1', 'a2_b2'}, {1, -1, -1, 1});
for r = 1:size(sc.data, 1)
    key = sprintf('%s_%s', sc.factor1{r}, sc.factor2{r});
    si = sc.subject(r);
    dd(si, :, :) = dd(si, :, :) + sgn_map(key) * sc.data(r, :, :);
end
scd = sc;
scd.data = dd;
scd.subject = 1:n_sub;
[tls, zero_tls] = deal(make_timelocks(scd), make_zero_timelocks(scd));
cfg = base_cfg(scd, 0);
cfg.statistic = 'depsamplesT';
cfg.design = [1:n_sub, 1:n_sub; ones(1, n_sub), 2 * ones(1, n_sub)];
cfg.uvar = 1;
cfg.ivar = 2;
cfg.numrandomization = 'all';
results.rm_anova_interaction = ...
    run_and_report('rm_anova_interaction (depsamplesT on double diff)', cfg, ...
                   [tls, zero_tls]);

% ------------------------------------------------------ spatiotemporal_paired
% paired t on 4 channels with chain neighbours A1-A2-A3-A4
sc = S.spatiotemporal_paired;
tls = make_timelocks(sc);
[~, ~, cond] = unique(sc.factor1);
cfg = base_cfg(sc, 0);
cfg.statistic = 'depsamplesT';
cfg.design = [sc.subject; cond'];
cfg.uvar = 1;
cfg.ivar = 2;
cfg.numrandomization = 'all';
results.spatiotemporal_paired = run_and_report('spatiotemporal_paired', cfg, tls);

% ----------------------------------------------------------------- rm_3level
% one within-subject factor, 3 levels: depsamplesFunivariate
sc = S.rm_3level;
tls = make_timelocks(sc);
[~, ~, cond] = unique(sc.factor1);
cfg = base_cfg(sc, 1);
cfg.statistic = 'depsamplesFunivariate';
cfg.design = [sc.subject; cond'];
cfg.uvar = 1;
cfg.ivar = 2;
cfg.numrandomization = 20000;
results.rm_3level = run_and_report('rm_3level (depsamplesFunivariate)', cfg, tls);

save(out_file, 'results', '-v7');
fprintf('saved %s\n', out_file);

% ---------------------------------------------------------- helper functions
function cfg = base_cfg(sc, tail)
    cfg = [];
    cfg.method = 'montecarlo';
    cfg.correctm = 'cluster';
    cfg.clusterstatistic = 'maxsum';
    cfg.clusteralpha = 0.05;
    cfg.tail = tail;
    cfg.clustertail = tail;
    cfg.alpha = 0.05;
    cfg.correcttail = 'no';
    cfg.randomseed = 42;
    cfg.feedback = 'no';
    cfg.channel = 'all';
    cfg.latency = 'all';
    cfg.neighbours = chain_neighbours(cellstr(sc.chan_labels));
end

function neigh = chain_neighbours(labels)
    % chain adjacency: channel i is neighbour of i-1 and i+1
    neigh = struct('label', {}, 'neighblabel', {});
    for i = 1:numel(labels)
        nb = {};
        if i > 1, nb{end + 1} = labels{i - 1}; end %#ok<*AGROW>
        if i < numel(labels), nb{end + 1} = labels{i + 1}; end
        neigh(i).label = labels{i};
        neigh(i).neighblabel = nb;
    end
end

function tls = make_timelocks(sc)
    n = size(sc.data, 1);
    tls = cell(1, n);
    for i = 1:n
        tl = [];
        tl.label = cellstr(sc.chan_labels);
        tl.time = sc.times;
        tl.avg = reshape(sc.data(i, :, :), numel(tl.label), numel(tl.time));
        tl.dimord = 'chan_time';
        tls{i} = tl;
    end
end

function tls = make_zero_timelocks(sc)
    scz = sc;
    scz.data = zeros(size(sc.data));
    tls = make_timelocks(scz);
end

function out = run_and_report(name, cfg, tls)
    stat = ft_timelockstatistics(cfg, tls{:});
    out = [];
    out.stat = stat.stat;
    if isfield(stat.cfg, 'clustercritval')
        out.clustercritval = stat.cfg.clustercritval;
    end
    for direction = {'pos', 'neg'}
        d = direction{1};
        cl_field = [d 'clusters'];
        lab_field = [d 'clusterslabelmat'];
        if isfield(stat, cl_field) && ~isempty(stat.(cl_field))
            out.([d '_prob']) = [stat.(cl_field).prob];
            out.([d '_clusterstat']) = [stat.(cl_field).clusterstat];
            out.([d '_labelmat']) = stat.(lab_field);
        else
            out.([d '_prob']) = [];
            out.([d '_clusterstat']) = [];
            out.([d '_labelmat']) = [];
        end
    end
    fprintf('\n=== %s ===\n', name);
    if isfield(out, 'clustercritval')
        fprintf('  clustercritval: %s\n', mat2str(out.clustercritval, 10));
    end
    fprintf('  stat: max=%.6f min=%.6f sum=%.6f\n', ...
            max(out.stat(:)), min(out.stat(:)), sum(out.stat(:)));
    for direction = {'pos', 'neg'}
        d = direction{1};
        probs = out.([d '_prob']);
        cstats = out.([d '_clusterstat']);
        lab = out.([d '_labelmat']);
        for k = 1:numel(probs)
            inds = find(lab == k) - 1;  % 0-based linear (chan fastest) indices
            fprintf('  %s cluster %d: clusterstat=%+.6f prob=%.6f n=%d inds0=%s\n', ...
                    d, k, cstats(k), probs(k), numel(inds), mat2str(inds'));
        end
    end
end

export_mat.py

"""Export cluster_test scenario data to .mat for the FieldTrip runs (out-of-PR)."""

from pathlib import Path
import sys

import numpy as np
from scipy.io import savemat

import mne

sys.path.insert(0, Path(mne.__file__).parent / "mne/stats/tests")
from test_cluster_equiv import SCENARIOS  # noqa: E402

import mne  # noqa: E402

out = {}
for name, builder in SCENARIOS.items():
    df, kwargs = builder()
    rows = []
    for _, row in df.iterrows():
        d = row["data"]
        if isinstance(d, mne.Evoked):
            d = d.get_data()  # (n_chan, n_times)
        d = np.atleast_2d(np.asarray(d, float))  # (n_chan, n_times)
        rows.append(d)
    data = np.stack(rows)  # (n_rows, n_chan, n_times)
    n_chan, n_times = data.shape[1:]
    sc = {
        "data": data,
        "n_chan": n_chan,
        "n_times": n_times,
        "chan_labels": np.array(
            [f"A{i + 1}" for i in range(n_chan)], dtype=object
        ),
        "times": np.arange(n_times) / 1000.0,
    }
    if "subject" in df.columns:
        sc["subject"] = df["subject"].to_numpy().astype(float) + 1  # 1-based
    else:
        sc["subject"] = np.zeros(0)
    factor_cols = [c for c in df.columns if c not in ("data", "subject")]
    sc["factor_names"] = np.array(factor_cols, dtype=object)
    for fi, col in enumerate(factor_cols):
        sc[f"factor{fi + 1}"] = np.array(df[col].astype(str).to_list(), dtype=object)
    out[name] = sc

savemat("cluster_equiv_data.mat", out, oned_as="row")
print("wrote cluster_equiv_data.mat")
for name, sc in out.items():
    print(f"  {name}: data{sc['data'].shape} factors={list(sc['factor_names'])}")