ERDS Group analysis

💬 Support & Discussions
1

Hi all,

I am using the ERDS tutorial and trying to translate it to get a group statistic instead, though, I am not completely sure how to do it. I am guessing that just like other group ERP like analyses I want to create a group variable the stores each participant’s tfr that gets fed into the pcluster_test. Could anyone help me understand where to feed that into and if i need to do something to each participant first? starting from the fifth cell in the tutorial at Compute and visualize ERDS maps — MNE 1.4.2 documentation
I have added in a tfr_par that houses each participant’s tfr. I am just not sure how to plug that into the pcluster_test() and if I need to do something like first take the averages from each participant and then move forward with those as inputs to that test.

tfr_par=[  ]
for par in epochs:

    tfr = tfr_multitaper(
        par,
        freqs=freqs,
        n_cycles=freqs,
        use_fft=True,
        return_itc=False,
        average=False,
        decim=2,)
    tfr.crop(tmin, tmax).apply_baseline(baseline, mode="percent")

    event_ids=par.event_id

    tfr_par.append(tfr)

    for event in event_ids:
        # select desired epochs for visualization
        tfr_ev = tfr[event]
        fig, axes = plt.subplots(
            1, 5, figsize=(12, 5), gridspec_kw={"width_ratios": [10, 10, 10, 10, 1]}
        )
        for ch, ax in enumerate(axes[:-1]):  # for each channel
            # positive clusters
            _, c1, p1, _ = pcluster_test(tfr_ev.data[:, ch], tail=1, **kwargs)
            # negative clusters
            _, c2, p2, _ = pcluster_test(tfr_ev.data[:, ch], tail=-1, **kwargs)
            # 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)
2

In the end, it all comes down to building the data array to feed into pcluster_test. Let’s call that array X. How to build this array is very much dependent on exactly how you would like to form clusters.

What the tutorial is doing is making clusters inside a TFR image for a single channel, hence X has dimensions n_epochs x n_frequencies x n_times. In your case, I think you have instead of n_epochs participants who each contribute a single averaged TFR. Hence, the dimensions of your X are: n_participants x n_frequencies x n_times.

If you need other kinds of clusters, for example clusters also need to span across channels, then it becomes slightly more complicated.

I hope that helps a bit.

2 Likes
3

Ah, okay, that makes sense. Thanks for the clarification.
I have changed my code to

tfr_par={‘Cong’:,‘Incong’:,‘Neut’:}
for par in epochs:

epochs_power = list()
for condition in [par[k] for k in ("Cong", "Incong","Neut")]:
    tfr = tfr_multitaper(
        condition,
        freqs=freqs,
        n_cycles=freqs,
        use_fft=True,
        return_itc=False,
        average=False,
        decim=2,
    )
    tfr.crop(tmin, tmax).apply_baseline(baseline, mode="percent")
    tfr_par[condition._name].append(tfr.average().data)

for con in tfr_par:
tfr_par[con]=np.stack(tfr_par[con])

This gives me a dictionary of my conditions with each condition as a data array of Par_avg x ch x freqs x times
which allows me to run the tests:

for event in event_id:

# select desired epochs for visualization
tfr_ev = tfr_par[event]
fig, axes = plt.subplots(
    1, 5, figsize=(12, 5), gridspec_kw={"width_ratios": [10, 10, 10, 10, 1]}
)
for ch, ax in enumerate(axes[:-1]):  # for each channel
    # positive clusters
    _, c1, p1, _ = pcluster_test(tfr_ev[:, ch], tail=1, **kwargs)
    # negative clusters
    _, c2, p2, _ = pcluster_test(tfr_ev[:, ch], tail=-1, **kwargs)

    # 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)

    tfev=mne.time_frequency.EpochsTFR(info=tfr.info, data=tfr_ev, times=tfr.times, freqs=tfr.freqs)

    # plot TFR (ERDS map with masking)
    tfev.average().plot(
        [ch],
        cmap="RdBu",
        cnorm=cnorm,
        axes=ax,
        colorbar=False,
        show=False,
        mask=mask,
        mask_style="mask",
    )

    ax.set_title(par.ch_names[ch], fontsize=10)
    ax.axvline(0, linewidth=1, color="black", linestyle=":")  # event
    if ch != 0:
        ax.set_ylabel("")
        ax.set_yticklabels("")
fig.colorbar(axes[0].images[-1], cax=axes[-1]).ax.set_yscale("linear")
fig.suptitle(f"ERDS ({event})")
plt.show()

however the last part for plotting just gives me blank plots. it seems the test are ran correctly, though I will admit I am not 100% certain. Note I did have to remove the event in event_id loop out of the main for loop, as it is in the turorial, in order to create the conditions dictionary (tfr_par), so I am sure that is my problem. I am just not sure if there is a good workaround .

4

Update: The last code does work. I had originally tested with only 2 participants which made the mask cancel everything out. with all participants it works fine it seems. Thanks again @wmvanvliet for clarifying

5

Can you show what your clusters look like?

6

I unfortunately cannot share the actual output plots due to company standards until a publication comes out. But, this is how I clustered the data in the tfr_par variable with it being a dictionary of conditions and each condition being (par_avg x sensors x freqs x times)
image