Observations dimension for permutation cluster test with time-frequency data

[Sebasdevrs]

Hi MNE team!

I’m doing a time-frequency analysis on the difference between two conditions with regards to activity in one EEG channel. Everything has been going well so far, but the finding (and visualizing) of significant clusters has been a bit confusing to me as I’m unsure in what format the permutation_cluster_test wants me to present the data.

Here is a quick rundown of my data:
Basically, every participant has a number of recordings in each of the two conditions. I’ve made a list of raws for each of the conditions, containing every recording that was made within a condition. I’ve then extracted the epochs for each of the recordings in each of the conditions, and averaged these out to end up with one evoked per recording per condition. I’ve then combined these evoked objects into one evoked object per condition, after which I combined these combined evokeds into a evoked for the contrast between the two. Finally, I used tfr_morlet to compute the time-frequency representation of this contrast.

This all works fine, and I can plot the contrast TFR on which clusters appear to be clearly visible. However, I am having trouble determining and visualizing the signifcant clusters. The TFR.data for each condition that I use as input for the permutation cluster test is of the shape [n_channels(1), n_frequencies(35), n_times(449)], but the documentation for the mne.stats.permutation_cluster_test function states that it wants the data for each condition to be of shape [n_observations, n_frequencies, n_times] as I understand it.

What does the observations dimension here refer to? Is it the different epochs within the condition? the evokeds of single recordings in a condition? Any help would be much appreciated!

Thanks in advance. I feel like I’m probably missing something really basic here, so apologies if the answer is very obvious!

Kind regards,

Sebastiaan

[CarinaFo (Carina)]

Hello Sebastian,

The observations dimension refers to the number of participants you want to run your cluster based permutation test, e.g. compare whether there is a statisticially significant difference between two conditions. So you would have to create a list with the average TFR per condition and per participant (length of list should be n of participants). If you have two conditions and it is a within subject design you can subtract the conditions and you end up with one list that contains the difference (this would be a paired t-test).

Hope this helps?

Best,

Carina

[Sebasdevrs]

That makes perfect sense, Carina. Thank you so much! For the within-subject design, subtracting the conditions, I would use the permutation_cluster_1samp_test instead of permutation_cluster_test then, right?

[CarinaFo (Carina)]

Yes this is exactly what I did for my within-participants design.

[Sebasdevrs]

Thanks again for your reply! I’ve tried implementing this method and it results in a mask with the right size, but the mask does not really seem to fit the data yet (I added a figure of the TFA with corresponding mask below, so you can see what I mean), and am wondering whether there is something I did wrong.

Long read, I know. I understand if you do not have the time to go through my entire design, but if you could help me out that would be a lifesaver!

For the array to be used for the cluster test, I first concatenated all the epochs from the different raw files with their corresponding condition and subject. This resulted in 2 sets with 7 epoch objects in each, corresponding to the 7 subjects that went through both conditions. I then made an evoked object from each of these 2x7 epoch objects.

I then used the mne.combine_evoked function with weights = [1,-1] to subtract one condition’s evoked from the other’s, for each participant separately. I created a TFR of each participant’s contrast evoked separately (resulting in seven [n_frequencies, n_timepoints] arrays. I then put these together in one [7, n_frequencies, n_timepoints] array, and used this as input for the mask.

For the actual TFR plot, I concatenated all epochs corresponding to a condition into one epochs object, created an evoked object out of both, and combined these into a contrast (again with mne.combine_evoked with weights = [1,-1]. I also created the TFR of this contrast, which I then plotted along with the mask that I had created.

This results in the figure, which doesn’t seem to fully make sense. Furthermore, when I increase the threshold (which I tried to cut off some of the area’s with lower power), the mask doesn’t become tighter around the cluster, but instead cuts into the cluster, while leaving some of the less extreme clusters in there.

Thanks again!

[CarinaFo (Carina)]

Could you provide us with a short example code that implemented the cluster test as well as the plot?

[Sebasdevrs]

For sure. Here it is:
Using python 3.6.6 on Windows 10

# Parameters:
freqs = np.linspace(0.6, 2, 32)
n_cycles = 3
tmin = -1.25 #seconds
tmax = 8 #seconds
baseline = (None,0)
channel='EEG right (uV)'
n_permutations = 'all'
threshold = None
n_subj = 7
n_freq = 32
n_timepoints = 790
# TFR for the contrast per subject
tfr_contrast_mask = []
for contrast in evoked_combined_sub:
    tfr_contrast_mask.append(tfr_morlet(contrast, freqs=freqs,
                                        n_cycles=n_cycles, decim=3, n_jobs=4,
                                        return_itc=False, zero_mean=True))
tfr_contrast_array = np.zeros([n_subj,n_freq ,n_timepoints])

for i in range(len(tfr_contrast_mask)):
    tfr_contrast_array[i,:,:] = np.squeeze(tfr_contrast_mask[i].data)
# Cluster permutation test:
t_obs, clusters, cluster_pv, H0 = \
    permutation_cluster_1samp_test(tfr_contrast_array, out_type='mask',
                            n_permutations=n_permutations, threshold=threshold, tail=1)
# Creating a mask of all the clusters:

mask = np.zeros_like(t_obs,dtype=bool)

for cluster in clusters:
    mask[cluster] = 1

# TFR for the contrast over all subjects combined:
tfr_combined = tfr_morlet(evoked_combined, freqs=freqs,
                          n_cycles=n_cycles, decim=3, n_jobs=4,
                          return_itc=False, zero_mean=True)
# Plotting the TFR with the mask:
fig_tfr_combined = tfr_combined.plot(baseline = baseline,
    tmin=tmin, 
    tmax=tmax, 
    cmap=cmap, # 'jet'
    title='TFR: Contrast, Threshold = None',
    mask = mask,
    mask_style = mask_style, # 'both'
    mask_cmap = cmap, # 'jet'
    mask_alpha = 0.5)

Thank you so much!

[CarinaFo (Carina)]

Hi Sebas,

Can you calculate separate TF estimates similar to this example, I am not sure if the combine_evokeds messes things up here, maybe someone else knows?

I suggest to creat the contrast like shown in this example Non-parametric between conditions cluster statistic on single trial power — MNE 1.6.0 documentation

decim = 2
freqs = np.arange(7, 30, 3) # define frequencies of interest
n_cycles = 1.5
tfr_epochs_1 = tfr_morlet(epochs_condition_1, freqs,
n_cycles=n_cycles, decim=decim,
return_itc=False, average=False)
tfr_epochs_2 = tfr_morlet(epochs_condition_2, freqs,
n_cycles=n_cycles, decim=decim,
return_itc=False, average=False)
tfr_epochs_1.apply_baseline(mode=‘ratio’, baseline=(None, 0))
tfr_epochs_2.apply_baseline(mode=‘ratio’, baseline=(None, 0))
epochs_power_1 = tfr_epochs_1.data[:, 0, :, :] # only 1 channel as 3D matrix
epochs_power_2 = tfr_epochs_2.data[:, 0, :, :] # only 1 channel as 3D matrix
contrast = epochs_power_1 - epochs_power_2

and then feed the conrast to your permutation test.

Then try to plot the output using the example in this tutorial Non-parametric 1 sample cluster statistic on single trial power — MNE 1.6.0 documentation

Let me know if you get a meaningful cluster this way?

Best,

Carina