TFR Analysis and Cluster-Based Permutation with 1 sample and two conditions

💬 Support & Discussions
,
1
  • MNE version: 1.8.0
  • operating system: macOS Sonoma 14.3.1

Hi everyone,
I’m a beginner with MNE and I’ve recently started working on time-frequency analysis with my EEG data. I’m analyzing data from a within-subjects design with 36 participants and two conditions (metaphor and literal) that I want to compare. The event of interest spans from 0 ms to 2000 ms, so I set up my epochs as follows::

epochs_images = mne.Epochs(raw, events, event_id=event_id, tmin=-1.000, tmax=2.500, baseline=(-0.750, -0.250), preload=True)

After that, I ran a time-frequency analysis using Morlet wavelets. Starting from prior literature, I tried two approaches for the n_cycles parameter: firstly using a fixed value of 3 cycles, secondly half the frequency:

freqs_low = np.arange(5, 41, 2)
n_cycles = 3

freqs_low = np.arange(5, 41, 2)
n_cycles = freqs_low / 2

I then ran the TFR analysis like this:
power_morlet = mne.time_frequency.tfr_morlet(epochs_condition, freqs=freqs_low, n_cycles=n_cycles, use_fft=True, return_itc=False, decim=1, n_jobs=-1)

Next, I organized my data for the two conditions into separate lists and cropped the time range of interest tmin, tmax = -0.7, 2.0. Data were then transposed as (0, 2, 3, 1) and an adjacency matrix was created combined_adjacency = combine_adjacency(len(metaphor_info.freqs), len(metaphor_info.times), filtered_sensor_adjacency)

I tried a cluster permutation test with the following code:

from mne.stats import spatio_temporal_cluster_test
import scipy.stats

n_permutations = 1040
n_observations = len(metaphor_data_transposed)

pval = 0.05
df = n_observations - 1
threshold = scipy.stats.t.ppf(1 - pval / 2, df)

T_obs, clusters, cluster_p_values, H0 = mne.stats.permutation_cluster_test(
    [metaphor_data_transposed, literal_data_transposed],
    n_permutations=n_permutations,
    threshold=threshold, 
    tail=0,
    out_type='mask', 
    verbose=True,
    adjacency=combined_adjacency, 
    seed=42,
    stat_fun=ttest_rel_nop
)

I also tried an alternative approach using permutation_cluster_1samp_test


X = metaphor_data_transposed - literal_data_transposed
n_observations = len(metaphor_data_transposed)
pval = 0.05
df = n_observations - 1
thresh = scipy.stats.t.ppf(1 - pval / 2, df)

n_permutations = 1024

T_obs, clusters, cluster_p_values, H0 = permutation_cluster_1samp_test(
    X,
    n_permutations=n_permutations,
    threshold=thresh,
    tail=0,
    adjacency=combined_adjacency,
    out_type="mask",
    verbose=True,
    n_jobs=-1,
    seed=42
)

The tests returned 10 clusters, but none of them were significant. I think it is quite strange that there are no significant differences between the conditions, so I’m wondering if I might be doing something incorrectly.
Any guidance on my approach or suggestions for troubleshooting would be greatly appreciated. Thanks in advance for your help!

2

Hello,

that is a tricky one. You run a cluster test over 3D data (time, frequencies and channels). Maybe the clusters found are rather small or very large? You could change the cluster-forming threshold. If you plot the t-values do you see the strongest differences between conditions in channels that you would expect (sanity check?). The frequency resolution might also make a difference (but there is no right or wrong). Regarding the test itself, you want to use the 1-samp test, as it runs a paired t-test which is the one you want for your design. Also make sure that the dimensions of X (the contrast between conditions) matches the adjacency matrix (which I think it does).

Hope this is helpful.

Best,

Carina

3

Dear Carina, if I may,
thanks for the reply! Actually, if I understand correctly, the data I am using for the cluster permutation are 4-dimensional: X.shape (36, 18, 1551, 58) . When you talk about the cluster forming threshold you refer to the parameter threshold and not pval, right? Anyway, since I am running exploratory analyses, I have no strong hypotheses about the channles in which I expect to observe the strongest differences between conditions.

Here is the shape of adjacency matrix: (1671444, 1671444) combined_adjacency

COOrdinate sparse array of dtype 'float64'
	with 19062620 stored elements and shape (1671444, 1671444)

adjacency_matrix
adjacency_matrix651×651 44.8 KB

To create the adjacency matrix, I uploaded a standard montage, i.e., easycapM1 and then removed some electrodes that I don’t have, since my data refer to 58 electrodes:

sensor_adjacency, ch_names = mne.channels.read_ch_adjacency("easycapM1")
common_channels = [ch for ch in metaphor_info.info['ch_names'] if ch in ch_names]
common_indices = [ch_names.index(ch) for ch in common_channels]

# Aggiorna sensor_adjacency per i canali comuni
filtered_sensor_adjacency = sensor_adjacency[common_indices][:, common_indices]

# Combina la matrice di adiacenza con quella temporale e frequenziale
combined_adjacency = combine_adjacency(len(metaphor_info.freqs), len(metaphor_info.times),filtered_sensor_adjacency)

Next, I created X: X = metaphor_data_transposed - literal_data_transposed

Since I don’t know which threshold to use for the permutation test I left None:

n_observations= len(X)
tail= 0

pval = 0.05  # arbitrary
df = n_observations - 1  # degrees of freedom for the test
thresh = None  # two-tailed, t distribution

# Set the number of permutations to run
n_permutations = 1024  # Increase this number for real analysis

# Run the cluster-based permutation test
T_obs, clusters, cluster_p_values, H0 = spatio_temporal_cluster_1samp_test(
    X,
    n_permutations=n_permutations,
    threshold=thresh,
    tail=tail,
    adjacency=combined_adjacency,
    out_type="mask",
    verbose=True,
    n_jobs= -1,
    seed=42
)

This is the output of the cluster permutation:

Using a threshold of 2.030108
stat_fun(H1): min=-3.882895189270287 max=2.4152846969668698
Running initial clustering …
Found 11 clusters

However, when I inspect cluster_p_values I get:

array([1.        , 1.        , 1.        , 0.94628906, 1.        ,
       1.        , 1.        , 1.        , 1.        , 1.        ,
       1.        ])

I have the feeling there is a discrepancy between the adjacency matrix and the X array I use for the permutation, but I cannot understand how to solve it.

I hope I provided all useful information. Thank you in advance!

5

Hello MNE community!
This morning, I am trying to perform cluster permutation using TFCE method described in the documentation: Statistical inference — MNE 1.8.0 documentation

Here is what the operation is printing:

from mne.stats import spatio_temporal_cluster_1samp_test
titles = ["t"]
n_permutations = 1000  # Increase this number for real analysis
titles.append(r"$\mathbf{C_{TFCE}}$")
threshold_tfce = dict(start=3, step=1)#consigliato qui: https://mail.nmr.mgh.harvard.edu/pipermail/mne_analysis/2014-August/002306.html, faccio 3 visto che nei miei dati il valore statistico più estremo è 3.88
t_tfce, _, p_tfce, H0 = spatio_temporal_cluster_1samp_test(
    X,
    n_jobs=4,
    threshold=threshold_tfce,
    adjacency=combined_adjacency,
    n_permutations=n_permutations,
    verbose=True,
    seed= 42,
    check_disjoint=True,
    out_type="mask",
)
ts.append(t_tfce)
ps.append(p_tfce)
mccs.append(True)
plot_t_p(ts[-1], ps[-1], titles[-1], mccs[-1])

However, time processing is too long. I left my computer working yesterday night and this morning it still was! Also, not having progression bar is not helpful for me.

Here is the output of the code:

stat_fun(H1): min=-3.882895189270287 max=2.4152846969668698
No disjoint adjacency sets found
Running initial clustering …
Using 1 thresholds from 3.00 to 3.00 for TFCE computation (h_power=2.00, e_power=0.50)
Found 1514844 clusters

I understand that the main reason for such long prossing time might be the great number of clusters found. Anyway, even if I tried to play around with the parameters, I still get 1514844 clusters. Since these analyses are exploratory I don’t want to focus on a specific frequency range and crop the signal according to.
Do you have any recommendations on additional parameter adjustments or approaches to reduce processing time without focusing on a specific frequency range?

I think I am stuck :confused:

6

This looks fine. Can you check the shape of T_obs from the cluster permutation test? Also if the adjacency matrix would not match the X shape you would get an error, which you apparently don’t get, so this should be fine.

1 Like
7

TFCE does take a lot of time, especially with your 3D input. I noticed that the time dimension has a lot of samples, maybe you could downsample the data for TFCE?

1 Like
8

Hi Carina, thank you for the reply!

In these days I made several attempts but could not find any difference. Last code I used for the permutation is the following:

from mne.stats import spatio_temporal_cluster_1samp_test
import scipy.stats
# Compute the difference data
X = metaphor_data_transposed-literal_data_transposed
n_observations= len(metaphor_data_transposed)
tail= 0

pval = 0.05  # arbitrary
df = n_observations - 1  # degrees of freedom for the test
thresh = scipy.stats.t.ppf(1 - pval / 2, df)  # two-tailed, t distribution

# Set the number of permutations to run
n_permutations = 1024  # Increase this number for real analysis

# Run the cluster-based permutation test
T_obs, clusters, cluster_p_values, H0 = clu = spatio_temporal_cluster_1samp_test(
    X,
    n_permutations=n_permutations,
    threshold=thresh,
    tail=tail,
    adjacency=combined_adjacency,
    out_type="indices",
    verbose=True,
    check_disjoint=True,#PROVO AD INSERIRE QUESTO PARAMETRO VISTA LA COMPLESSITà DELLA MATRICE
    n_jobs= -1,
    max_step=0,
    seed=42
)
print("Eseguito il test di permutazione cluster-based.")

I get many clusters but only one is marginally significant (p = 0.57)

This is the shape of T I get from cluster permutation: (46, 1501, 58)
Also, I attach a plot of the significant cluster I get:

Cluster_Permutation_18_November
Cluster_Permutation_18_November1046×276 140 KB

I don’t know if there are any error but I am thinking that there are no differences between the two conditions.