mne.stats.permutation_cluster_1samp_test

💬 Support & Discussions
1

Hello MNE community,

I am a bit confused about the input to mne.stats.permutation_cluster_1samp_test(X) (see related topic Data for permutation cluster 1samp test)

I have 40 participants that were tested in 2 conditions so a within subject design. I am interested in pre stimulus alpha power differences between my 2 conditions.
I calculated power using tfr_morlet(epochs_cond1) per participant and per condition and I stored it in a list. So I have two lists with length of number of participants.
Example output:

Condition 1 for 6 participants:

[<AverageTFR  |  time : [-1.000000, 1.499600], freq : [6.000000, 35.000000], nave : 215, channels : 61, ~7.9 MB>,
 <AverageTFR  |  time : [-1.000000, 1.499600], freq : [6.000000, 35.000000], nave : 154, channels : 61, ~7.9 MB>,
 <AverageTFR  |  time : [-1.000000, 1.499600], freq : [6.000000, 35.000000], nave : 185, channels : 61, ~7.9 MB>,
 <AverageTFR  |  time : [-1.000000, 1.499600], freq : [6.000000, 35.000000], nave : 195, channels : 61, ~7.9 MB>,
 <AverageTFR  |  time : [-1.000000, 1.499600], freq : [6.000000, 35.000000], nave : 228, channels : 61, ~7.9 MB>,
 <AverageTFR  |  time : [-1.000000, 1.499600], freq : [6.000000, 35.000000], nave : 102, channels : 61, ~7.9 MB>]

If I want to use the cluster permutation test, the documentation says that X should be a array with the first dimension reflecting the difference between conditions. I tried various ways but I am a bit lost.
I would be very grateful for an example.

Thank you for this awesome toolbox,

Carina

1 Like
2

Hi @CarinaFo,
if you want to compare two conditions, you want to use mne.stats.permutations_cluster_test, which as the first argument takes a list of conditions, where each condition is an array of (participants, ...) shape.

To compare your tfrs, I would:

  • first crop the tfrs to the desired time range (in your case - prestimulus) and frequency range:
     tfr_cond1_crop = [tfr.copy().crop(tmin=..., tmax=..., fmin=..., fmax=...) for tfr in tfr_cond1]
  • then stack all tfr arrays for each condition:
    tfr_array_cond1 = np.stack([tfr.data for tfr in tfr_cond1_crop], axis=0)
  • your X argument is now [tfr_array_cond1, tfr_array_cond2] (if you are going to use a t test the order of the conditions will affect the sign of the t values)
  • now you are almost ready to conduct the cluster permutation test, you just need to specify channels adjacency and then build the whole adjacency for the channels-time-frequency space using mne.stats.combine_adjacency. You also need to choose the correct stat_fun (that depends on the test you want to perform - I guess it is within subjects t test in your case?)
2 Likes
3

Hello @mmagnuski,

thank you for your help. Yes I have a within design, so I would specify mne.stats.ttest_1samp_no_p() in stats_fun, right? But I still have the same issue with the X array I provide as an input. It can’t be the same format as the X array I provide to mne.stats.permutation_cluster_test?

Update: [ mne.stats.combine_adjacency ] does not seem to be supported anymore.

Carina

4

If you have two with-subject conditions then you want to use ttest_rel as your stat_fun. Because mne does not have a specialized ttest_rel_nop you can create it yourself:

from scipy.stats import ttest_rel

def ttest_rel_nop(*args):
    tvals, _ = ttest_rel(*args)
    return tvals

Then, your X is just as in my examples: a list with two arrays (each array is one condition, all subjects stacked in the first dimension).

I’m not sure what you mean by combine_adjacency not being supported anymore. What problem do you experience? Maybe you have an outdated version of mne - combine_adjacency was introduced in mne 0.21.

3 Likes
5

Yes, you were right I used a old version, sorry for that.

Trying to run the permutation test using this command:

n_permutations = 1000
threshold = 2.015 #checked in table?

T_obs, clusters, cluster_p_values, H0 =
mne.stats.permutation_cluster_test([tfr_array_cond1, tfr_array_cond2], n_permutations=n_permutations,
threshold=threshold, tail=0,
out_type=‘mask’, verbose=True, stat_fun = ttest_rel_nop(tfr_array_cond1, tfr_array_cond2))

I encounter the following error:

847 # Step 1: Calculate t-stat for original data
848 # -------------------------------------------------------------
→ 849 t_obs = stat_fun(*X)
850 _validate_type(t_obs, np.ndarray, ‘return value of stat_fun’)
851 logger.info(‘stat_fun(H1): min=%f max=%f’ % (np.min(t_obs), np.max(t_obs)))

TypeError: ‘numpy.ndarray’ object is not callable

Also threshold = None doesnt’t work here, so I just took the respective t value from the table for n-1 degrees of freedom.

Thank you for your help.

Carina

6

The stat_fun argument should be the function itself, not the result of applying the function, so in your snippet it should be:

..., stat_fun=ttest_rel_nop)
1 Like