Cluster-based Permutation Analysis

💬 Support & Discussions
1

Hello MNE community,

I am in the process of doing some ERD/ERS analysis. After preprocessing steps, I compute time/frequency ERD/ERS maps using tfr_multitaper(). and can plot the maps with no issues. Following the instructions given here, I plot the corrected maps. However, for all events, the maps look exactly the same. I was wondering if you could help me with this issue. Below is a part of my code for ERDS computation and visualization. Thanks.

events, _ = mne.events_from_annotations(raw)
events = events[np.in1d(events[:, 2], (10002, 10004, 10006)), :]


tmin, tmax = -1, 3
epochs = mne.Epochs(raw, events, dict(down=10002, left=10004, right=10006), tmin, tmax,
                    picks=("C3", "Cz", "C4"), baseline=None, preload=True)

freqs = np.arange(1, 31)
baseline = -1.5, -0.5
cnorm = TwoSlopeNorm(vmin=-1, vcenter=0, vmax=1.5)  # min, center, and max ERDS
vmin, vmax = -0.5, 1.5
kwargs = dict(n_permutations=100, step_down_p=0.05, seed=1,
              buffer_size=None, out_type='mask')  # for cluster test
cmap = center_cmap(plt.get_cmap("RdBu"), vmin, vmax)

for event in epochs.event_id:  # for each condition
    tfr = tfr_multitaper(epochs, freqs=freqs, n_cycles=freqs, use_fft=True,
                     return_itc=False, average=False, decim=2)
    tfr.crop(-0.5, 2.5).apply_baseline(baseline, mode="percent")
    tfr_ev = tfr
    fig, axes = plt.subplots(1, 3, figsize=(15, 5))
    
    for ch, ax in enumerate(axes):
        # positive clusters
        _, c1, p1, _ = pcluster_test(tfr_ev.data[:, ch], tail=1, **kwargs)
        print(p1)
        # negative clusters
        _, c2, p2, _ = pcluster_test(tfr_ev.data[:, ch], tail=-1, **kwargs)
        print(p2)
        # note that we keep clusters with p <= 0.05 from the combined clusters
        # of two independent tests; in this example, we do not correct for
        # these two comparisons
        c = np.stack(c1 + c2, axis=2)  # combined clusters
        p = np.concatenate((p1, p2))  # combined p-values
        mask = c[..., p <= 0.05].any(axis=-1)

        # plot TFR (ERDS map with masking)
        tfr_ev.average().plot([ch], cmap="RdBu", axes=ax,
                              colorbar=False, show=False, mask=mask,
                              mask_style="mask")
        
        ax.set_title(epochs.info["ch_names"][ch])
plt.show()
2

Hi @yvaghei,

Questions like this, that ask to check a piece of code and see if there are any bugs are less likely to be answered. We usually expect the user to do a bit of work to narrow down the problem (if the problem is so broad as in your case). The short story is: I doubt many people on the forum will have more motivation than you yourself to carefully read that piece of code. :slight_smile: I would advice you to for example go through the code in debug mode to see if something unexpected is not happening at some point.

However, since I already looked into this code and found the potential problem, two more notes:

  1. Notice that your events loop goes like this:

    for event in epochs.event_id:
    tfr = tfr_multitaper(epochs, ...)

    you don’t select the epochs - so in each iteration you perform tfr on all the epochs.

  2. I’d advice to use more permutations than just 100 (I usually use at least 1000) - this affects the “resolution” of the cluster p values that you get.

1 Like