Observations dimension for permutation cluster test with time-frequency data

💬 Support & Discussions
,
1

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

2

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

2 Likes
3

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?

1 Like
4

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

1 Like
5

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.

Figure_TFR_none
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.
Figure_TFR_2.3

Thanks again!

6

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

2 Likes
7

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!

1 Like
8

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

9

Hi guys, currently I am trying to analyze Time frequency analysis. sorry for joining with a new problem. but seems the problem is quite straightforward from here. From this discussion I got that for paired t test, one sample permutation is the correct one which is I agree with both of you. The problem I have is, seems there is a miss match with the visualization before permutation and the result of permutation. In the contrast image before applying permutation I saw the contrast is quite obvious when I have subtracted two power conditions. However I put it in permutation, there is nothing significant. I mean this is unlikely because we saw it and PI said the contrast value amplitude is still quite big for contrast value it’s (score 1e*9).

I am figuring out, does the threshold I use, is too high because the p value cluster show everything almost 1. I feel like really stuck, where to debug the problem. I have tried to modify the data smaller frequency, modify time duration, pick signficant channels from ERP. but nothing work to me and currently I don’t have someone to discuss in python or mne point of view what to change.

I am sorry for disturbing your notification but I am somehow stuck.

20240721_142213
20240721_1422131089×465 88.2 KB

10

Hello,

Could you provide the code you use for calculating the TFR you feed to the cluster permutation test as well as the code for the permutation test?

Best,

Carina

11

Dear Carina,

Thank you for replying me.

Yes, from the experiment we expect to see the low frequency around theta and alpha. According to the guide FLUX, I can successfully extract the power with multitapper.

def power_function(data):
freqs = np.arange(2, 31, 1)
n_cycles = freqs / 2
time_bandwidth = 2.0

power_data = mne.time_frequency.tfr_multitaper(
data, 
freqs=freqs, 
n_cycles=n_cycles,
time_bandwidth=time_bandwidth, 
picks = 'eeg',
use_fft=True, 
return_itc=False,
average=True, 
decim=2,
n_jobs= -1, 
verbose=True)

return power_data

then I can extract the power from each condition by this code

power_stay_similarity= power_function(epochs_st_sim)
power_switch_similarity= power_function(epochs_sw_sim)

then I subtracted to get the contrast of the condition

#contrast
contrast_power_simstsw= power_stay_similarity.copy()
contrast_power_simstsw_data = power_switch_similarity.get_data() - power_stay_similarity.get_data()
contrast_power_simstsw._data= contrast_power_simstsw_data

Then after I get the contrast power, I save the contrast power file with this format

contrast_power_simstsw.save(power_stsw_sim + sub + ‘_contrast_power_stsw_sim_tapper12_b-ave-tfr.h5’,overwrite=True)

Then I apply permutation to the evoked contrast of all participants total: 35 participants by this code,

from scipy import stats

n_subjects=len(subjects)

p_threshold = 0.05
df = n_subjects - 1 # degrees of freedom for the test
t_threshold = stats.distributions.t.ppf(1 - p_threshold / 2, df=df)

tail = 0 # for two sided test
n_permutations = 50

Now use the combined adjacency in your cluster analysis

cluster_stats = mne.stats.spatio_temporal_cluster_1samp_test( sub_data, threshold=t_threshold,n_permutations=n_permutations, n_jobs=2, verbose=True, tail=tail,adjacency=tfr_adjacency,step_down_p=0.05, out_type=‘indices’,check_disjoint=True, seed=None)

T_obs, clusters, p_values, _ = cluster_stats

The original epochs was without baseline then I realized all TFR tutorial apply baseline (None,0) with mode mean by default.

rightnow, by applying the baseline I can get the significant, but it has so many cluster I get. because I expect to sync the result my ERP result but now I got so many cluster, does’t know which one reliable. I will play with the threshold.

my question now, is there any guideline to apply the correct baseline, I am worry that this is just I got by chance cause I am trying to get solution.

I am appreciate your reply. Thank you very much. It’s really kind of you.

here I attached the result of my cluster permutation

Best,
Risa

12

Screenshot from 2024-07-23 22-22-05
Screenshot from 2024-07-23 22-22-05826×816 301 KB

Screenshot from 2024-07-23 22-22-19
Screenshot from 2024-07-23 22-22-19831×830 311 KB

13

Hello Risa,

baseline correction depends on your research question and the trial design, it should definitely not be (ab)used to get significant results. Do you expect differences in power between conditions before stimulation onset? First thing I noticed: you should never run a permutation test on 50 permutations only, 1000 permutations is minimum to get a stable result. I would try that first. If you still get so many small clusters, you can adjust the cluster threshold t which is not p-hacking as it only controls the cluster extent, not significance.

Hope this helps,

Carina

16

Hi Carina, Yes thank you, I have applied with 1024 n permutation. I think, I found the problem. I think I didn’t set the adjacency correctly. I need to create a connectivity via adjacency matrix, but the problem is the adjacency now wasn’t matched with the shape of the data.

I need to make it match.

1 Like
17

Hello! Did you ever manage to solve this by fixing the adjacency code? I am in a similar situation, running cluster based permutation tests (1-sample) and finding no significant clusters (which makes no sense, since these are event-locked epochs with ERPs, significant clusters should be a given). I assume it has something to do with the adjacency matrix not matching the data, but I am not well versed enough in the underlying maths to figure out what could be going wrong.