Cluster-based permutation test for a 2x3 Design

💬 Support & Discussions
, ,
1

Hi everyone,

I am trying to compute statistics on ERP-Data (32 Channels x 720 time-points) for a 2-by-3 design. Following this tutorial (Mass-univariate twoway repeated measures ANOVA on single trial power — MNE 1.6.0 documentation) I was able to calculate clusters for the effects separately (2 main effects plus 1 interaction: “A”, “B”, then “A:B”).
My question: Is this valid (statistically speaking) or is there a way to somehow include all effects at once in a cluster-based permutation test?

Background:
If I am not mistaken, it is not the same to calculate an ANOVA with all effects as having to calculate it for every effect separately - specially when it comes to degrees of freedom.

Computing all effects together (effects = “A*B”), does throw an error stemming from line 982/983 in cluster_level.py when calling _permutation_cluster_test. This probably occurs because stat_fun ( f_mway_rm modified in the tutorial) outputs F-values for all data-points (time-points x channels) and for every effect (3 in my case). Line 982 compares the shape of the stats (t.obs_shape = 2, 23040 (32x720)) with the product of the sample_shape (32x720). It probably expects only one statistical value per data-point instead of several.

This is my error message: ValueError: t_obs.shape (2, 23040) provided by stat_fun <function stat_fun at 0x0000018392C144A0> is not compatible with the sample shape (32, 720)

FYI: my X for the permutation_cluster_test has the shape (30, 6, 32, 720)

#define function: swapaxes not necessary as performed when concatenating arrays
def stat_fun(*args):
    return f_mway_rm(
        #np.swapaxes(args, 1, 0),
        np.asarray(args),
        factor_levels=factor_levels,
        effects=effects,
        return_pvals=False,
    )[0]

#mock-data
X = np.random.rand(30, 6, 32, 720)

#define effects
factor_levels = [3,2]  # number of levels in each factor
effects = "A:B"  # R-Language for effects (or "A" or "B")

pthresh = 0.05  # set threshold rather high to save some time
f_thresh = f_threshold_mway_rm(n_subjects, factor_levels, effects, pthresh)
tail = 1  # f-test, so tail > 0

F_obs_int, clusters_int, cluster_p_values_int, _ = mne.stats.permutation_cluster_test(
    X,
    stat_fun=stat_fun,
    threshold=f_thresh,
    tail=tail,
    n_jobs=None,
    n_permutations=1000,
    buffer_size=None
)

Best regards,
Brais

  • MNE version: 1.4.2
  • operating system: Windows 11