TFR difference maps

[gmj3]

Hello,

I’m a beginner practicing analyzing data for a motor dataset in using MNE python. I tried analyzing data for all 109 participants in a dataset : EEG Motor Movement/Imagery Dataset v1.0.0 . For this project, I am working on run 3, which focuses on left and right hand movement, and the c3 and c4 mu band activity for each movement. According to the literature, there is supposed to be a mu band suppression in the contralateral hemisphere to the movement :

Introduction to Mu Waves for BCI | OpenBCI Community

So I first tried analyzing and plotting one participant, then grand average for left vs right hand across participants, and then did difference maps.

Please take a look at my difference maps, and let me know if it agrees with the expected results.

Thanks!

[ZyntZ (Anna Sok)]

Hi,

Your plots look qualitatively consistent with the expected lateralization, but there is an important distinction between a hemispheric difference and mu suppression itself

For the left-hand condition, the post-event C4-C3 values are mostly negative, meaning that power at C4 is lower than at C3. Since C4 is over the right sensorimotor cortex, this is compatible with contralateral mu suppression during left-hand movement. For the right hand condition, C4-C3 is mostly positive after the event, meaning that power at C3 is lower than at C4. That is compatible with contralateral suppression during right-hand movement

However a C4-C3 map alone cannot show that either channel was actually suppressed relative to baseline. For example, a positive C4-C3 value could result from a decrease at C3, an increase at C4, or both. I would therefore also plot the baseline-corrected TFRs for C3 and C4 separately. Apply baseline correction within each participant and condition before calculating the grand average, using the same baseline interval and preferably a relative or log-ratio measure. Then check whether C4 decreases from baseline for left-hand movement and C3 decreases for right-hand movement. A few other points may make the result easier to interpret:

1. The C3-C4 plot is exactly the negative of the C4-C3 plot, so showing both is redundant.
2. Use the same symmetric color limits for the left-and right-hand figures. Automatic limits can make the two conditions look more different or similar than they are.
3. Avoid interpreting the edges of the TFR too strongly because wavelet estimates there can be affected by edge effects.
4. Confirm that the Run 3 event labels were mapped correctly to left and right-hand movement.
5. To support a group-level conclusion, retain each participant’s contrast and test it across participants rather than interpreting only the grand-average image. A withinsubject permutation-cluster test across time and frequency would be suitable; report the significant clusters rather than judging significance from color alone. MNE provides examples for one-sample time-frequency cluster tests: Non-parametric 1 sample cluster statistic on single trial power — MNE 1.12.1 documentation

TL;DR: yes. the direction of the post-event C3/C4 asymmetry appears consistent with the expected contralateral pattern. But these difference maps by themselves do not yet demonstrate event-related mu suppression. Separate baseline corrected channel TFRs, consistent plotting limits, and a participant lvl statistical test are needed before making that conclusion

[gmj3]

Hello, thanks for your reply.

I already looped through all the subjects and applied baseline correction before averaging. Here is my code.

````power_list_right_mu = []

 power_list_left_mu = \[\] 

print(type(raw_resampled))
for subject, raw_resampled in all_data.items():
#raw_resampled = raws[0] # this calls the first element of raws, which is subject 1
print(subject, raw.info[“sfreq”])
print(type(raw_resampled))
events_resampled, event_id = mne.events_from_annotations(raw_resampled)raw resampled is a tuple. but events from annotations expects single object
epochs = mne.Epochs(raw_resampled, events_resampled, event_id, tmin=-1.5, tmax=2, baseline=None, preload=True) #we have to get a bigger range for alpha because it has longer wavelengths

left_epochs = epochs\["T1"\]
right_epochs = epochs\["T2"\]

power_mu_left = mne.time_frequency.tfr_morlet(
    left_epochs,
    freqs=np.arange(8, 13),
    n_cycles=7,
    return_itc=False
)
power_mu_right = mne.time_frequency.tfr_morlet(
    right_epochs,
    freqs=np.arange(8, 13),
    n_cycles=7,
    return_itc=False
)
power_mu_left.apply_baseline(baseline=(-1, 0), mode="percent") #they told me to do this before grand averaging
power_mu_right.apply_baseline(baseline=(-1, 0), mode="percent") #they told me to do this before grand averaging

power_list_left_mu.append(power_mu_left)
power_list_right_mu.append(power_mu_right)

grand_power_left_mu = mne.grand_average(power_list_left_mu)
grand_power_right_mu = mne.grand_average(power_list_right_mu)
print(‘done’)```

I already did the grand average left and right hand T3 and T4 before, and plotted them separately. I see a slightly more negative value for the contralateral hemisphere after time 0, for both hands. Here is the code I wrote for one of them (I did the same for both hands), and the two plots:

```motor_data2 = grand_power_right_mu.copy().pick([“C3”, “C4”]).data
limit = np.nanpercentile(np.abs(motor_data2), 99)
for channel in (“C3”, “C4”):
grand_power_right_mu.plot(
picks=channel,
vlim=(-limit, limit),
title=f"Grand average T2 {channel}",
)```

Is this what you mean by identical colour limits? I’m not sure how to do that for the different maps too.

Here is proof that T1 and T2 in run 3 corresponds to left and right fist:

EEG Motor Movement/Imagery Dataset

For the statistical test, I have to loop it over all participants, right? Before the grand average?

Thank you.

[ZyntZ (Anna Sok)]

Thanks, that clears it up. Your baseline correction is already in the right place: you apply it to each participant before computing the grand average. And yes, for Run 3, T1 is left-fist movement and T2 is right-fist movement

The new plots look qualitatively consistent with contralateral mu suppression: C4 is slightly more negative for left-hand movement, while C3 is slightly more negative for righthand movement. I would still avoid claiming a reliable group effect until you run the statistics.

Your colour limits are identical for C3 and C4 within each condition!, which is good. To make T1 and T2 directly comparable as well, calculate one limit from both datasets:

data = np.concatenate([
    grand_power_left_mu.copy().pick(["C3", "C4"]).data.ravel(),
    grand_power_right_mu.copy().pick(["C3", "C4"]).data.ravel(),
])
limit = np.nanpercentile(np.abs(data), 99)

Then use vlim=(-limit, limit) for all four plots. The same idea applies to the difference maps: calculate one limit_diff from both difference arrays and use vmin=-limit_diff, vmax=limit_diff for both figures

For the statistical test, you do need the dta from every participant, but you do not run a separate test for each person. Create one contrast map per participant, for example C4-C3 for left-hand movement and C3-C4 for right-hand movement, then stack those maps and run one group-level one-sample cluster permutation test. Do this using the subject-level TFRs, not the grand average

One thing I would check first is the strong narrow increase around time zero. With seven-cycle wavelets, activity at the cue can spread into nearby time points and even into the end of your baseline. Longer input epochs and a baseline that ends before zero, such as -1.5 to-0.5 s, would be a useful sensitivity check

Also print(subject, raw.info["sfreq"]) probably should be raw_resampled.info["sfreq"] and I assume “T3 and T4” was just a typo for C3 and C4

[gmj3]

Here is the extracting epochs step, adapted to my project. I have marked what I think I need to add, and what I don’t need (marked with ***). I have done most of the steps when calculating individual TFRs before doing the grand average. Do I have to run this separate from the process I used to calculate grand average? I noticed some differences from the grand average pathway: for example, there is tfr_epochs = epochs.compute_tfr( instead of mne.time_frequency.tfr_morlet( and there is logratio baseline instead of percentage. Please take a look and let me know anything I need to change:

#already have with loading EEG datasets, no need to add

data_path = sample.data_path()

meg_path = data_path / “MEG” / “sample”

raw_fname = meg_path / “sample_audvis_raw.fif” # did this with EDF

tmin, tmax, event_id = -0.3, 0.6, 1

raw = mne.io.read_raw_fif(raw_fname)

events = mne.find_events(raw, stim_channel=“STI 014”)

#don’t need

***include = []

raw.info[“bads”] += [“MEG 2443”, “EEG 053”] # bads + 2 more***

#don’t need

*** for speed, we’ll only look at right-temporal gradiometers (and EOG)

picks_eog = mne.pick_types(raw.info, eog=True)

picks_grad = mne.pick_types(raw.info, meg=“grad”, exclude=“bads”)

picks_rtemp = mne.pick_channels(

raw.info[“ch_names”], mne.read_vectorview_selection(“Right-temporal”), ordered=True

)

picks = list((set(picks_rtemp) & set(picks_grad)) | set(picks_eog))***

# Load condition 1

event_id = 1

epochs = mne.Epochs(

raw,

events,

event_id,

tmin,

tmax,

picks=picks,

baseline=(None, 0),

preload=True,

***reject=dict(grad=4000e-13, eog=150e-6),*** #don’t need

)

#do same with condition 2

evoked = epochs.average()

#do the same with condition 2

freqs = np.arange(8, 40, 2)

#do the same with condition 2

tfr_epochs = epochs.compute_tfr(

“morlet”,

freqs,

n_cycles=4.0,

decim=decim,

average=False,

return_itc=False,

n_jobs=None,

)

#do the same with condition 2```

I already did similar steps as the above code when looping through participants before doing the grand average. I just want to know if I should do it new for statistical test.

tfr_epochs.apply_baseline(mode="logratio", baseline=(-0.100, 0)) - #why are we doing logratio? Do I need to add this?

evoked.crop(-0.1, 0.4)

tfr_epochs.crop(-0.1, 0.4)

#do the same with condition 2

epochs_power = tfr_epochs.data

#do the same with condition 2

power_list_left_mu.append(power_mu_left)

power_list_right_mu.append(power_mu_right)```

[ZyntZ (Anna Sok)]

No, you do not need to repeat the epoching and TFR calculation. The tutorial uses the sample MEG dataset, so most of that code is specific to that example

Your existing lists already contain one baseline-corrected TFR per participant:

power_list_left_mu
power_list_right_mu

Use those lists for the group test, before reducing them to a grand average. epochs.compute_tfr(method="morlet", ...) is the newer form of mne.time_frequency.tfr_morlet(...); it is not a different analysis. You can also keep mode="percent". logratio is anther valid baseline method, but you should not mix the two within one analysis

For example, you can create one lateralization map per participant like this:

contrasts = []

for power_left, power_right in zip(
    power_list_left_mu, power_list_right_mu
):
    left = power_left.copy().pick(["C3", "C4"]).data
    right = power_right.copy().pick(["C3", "C4"]).data

    # Positive means lower power over the contralateral hemisphere
    left_contrast = left[0] - left[1]     # C3 - C4
    right_contrast = right[1] - right[0]  # C4 - C3

    contrasts.append((left_contrast + right_contrast) / 2)

X = np.stack(contrasts)  # participants × frequencies × times

Then run one group-level test on X, not one test for each participant and not a test on the grand average:

from scipy import stats

threshold = stats.t.ppf(1 - 0.05 / 2, df=X.shape[0] - 1)

T_obs, clusters, cluster_p, H0 = \
    mne.stats.permutation_cluster_1samp_test(
        X,
        threshold=threshold,
        tail=0,
        n_permutations=5000,
        seed=42,
        out_type="mask",
    )

significant = [
    (cluster, p)
    for cluster, p in zip(clusters, cluster_p)
    if p < 0.05
]

You do not need evoked = epochs.average(), the MEG channel selections or the MEG rejection values from the tutorial. Also, average=False keeps every trial, whereas your current tfr_morlet call returns a trialaveraged TFR for each participant. That is suitable for this participant-level group test

One separate issue is artifact handling.You should not copy the tutorial’s MEG rejection threshold, but the EEG data should still be checked for bad channels and artifacts before calculating the TFR

If you want to make separate claims about left- and right-hand movement, test left_contrast and right_contrast separately instead of averaging them, and correct for running two tests

[gmj3]

I did the above two steps and got these values:

Number of clusters: 5
Cluster p-values: [0.7358 0.2864 0.7646 0.5208 0.9394]
Number of significant clusters: 0
Significant clusters: []

This means that the effect is not enough to be considered significant, right? If so, is there a simple, beginner friendly way to make the results a little bit more robust? For example you mentioned excluding bad channels. So far, the only thing I learned in my EEG online course was artifact removal with ICA.. I know there are many ways to make the results more robust, but I don’t want to expand into may preprocessing methods indefinitely. I would only like to know what is a reasonable, sufficient QC for this analysis. Please let me know a reasonable, simple next step.

[viranovskaya (Daria Agafonova)]

Yes. With alpha = 0.05, none of the five clusters reached cluster-level significance; the smallest p-value is 0.2864. I would phrase this as “this analysis did not detect a significant cluster,” rather than “there is no effect.”

I would not add or tune preprocessing steps after seeing these p-values in order to obtain significance. A simple, defensible QC would be:

1. Confirm the Run 3 event mapping and use the same epoch and baseline settings for every participant.
2. Inspect the raw traces and PSD for persistently flat or clearly noisy channels. Mark those channels as bad before average referencing, ICA, and TFR calculation.
3. Use ICA only for components that are clearly supported as artifacts by their topography and time course, plus EOG/ECG signals if available. Do not remove a component because doing so improves the final statistical result.
4. Reject only epochs with clear gross artifacts, using the same predefined rule for all participants. Do not copy the MEG rejection threshold from the tutorial; there is no universal EEG threshold that is appropriate for every dataset.
5. Record the retained T1/T2 epoch counts, bad channels, excluded ICA components, and drop reasons for each participant.

Then freeze those QC decisions, recompute the participant-level TFR contrasts once, and rerun the same group test. If it remains non-significant, that is a valid result to report; it does not mean the code is wrong.

Because you also saw a strong feature near time zero, one predefined sensitivity check with a baseline ending before the cue (for example, -1.5 to -0.5 s) would be reasonable. I would not try several baselines and keep the one with the smallest p-value.

[gmj3]

Before doing the QC, someone else said I was missing one more step which is re-referencing. I got the information from here on EEG lab about rereferencing Re-referencing - EEGLAB Wiki , and then I looked at the MNE specific documentation: Setting the EEG reference — MNE 1.12.1 documentation . However, in MNE there are many ways to set the reference beyond picking non-scalp regions, such as using the average as a reference. Which approach is the best and why?

For the second step, I realized I should make the raw plot more interactive, so I used an interactive backend to do so. When looking at each channel, I can’t seem to find one that stands out; which is whether one is more noisy or flat than the others.

This pattern remains as I scroll vertically and horizontally. How can I inspect bad channels from here? If the data is generally noisy, do I have to change my filter? And does rereferencing have an effect?

[gmj3]

Hello, did you see my message? Please let me know.

[wmvanvliet (Marijn van Vliet)]

Could you post/link your analysis code so we could see what analysis steps you are and are not performing?

[wmvanvliet (Marijn van Vliet)]

By the way, have you seen: Motor imagery decoding from EEG data using the Common Spatial Pattern (CSP) — MNE 1.12.1 documentation?

[wmvanvliet (Marijn van Vliet)]

I tried looking at the data myself. It looks like there is already an average mastoid reference applied to it. For some subjects, the data has a different sample rate or something, so I just skipped over them. There are notable eye-blink artifacts, but the C3 and C4 channels are not that much affected by them. Here is an example script to show the ERD (a drop in signal power between 8-12 Hz after imagined movement onset at t=0):

import matplotlib.pyplot as plt
import mne
import numpy as np
from tqdm import tqdm

subjects = list(range(1, 81))
fnames = mne.datasets.eegbci.load_data(subjects=subjects, runs=[3])
tfr_left = list()
tfr_right = list()
for fname in tqdm(fnames, unit="subjects"):
    raw = mne.io.read_raw(fname, preload=True)
    mne.datasets.eegbci.standardize(raw)  # standardize channels names: "Cz.." -> "Cz"
    raw.set_montage("spherical_1005")
    raw.annotations.rename(dict(T1="left", T2="right"))  # as documented on PhysioNet

    epochs = mne.Epochs(raw, event_id=["left", "right"], tmin=-2.0, tmax=3.0)
    epochs_tfr = epochs.compute_tfr(method="multitaper", freqs=np.linspace(5, 20, 20))
    tfr_left.append(epochs_tfr["left"].average())
    tfr_right.append(epochs_tfr["right"].average())

mean_tfr_left = np.mean(tfr_left)
mean_tfr_right = np.mean(tfr_right)

fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(15, 10))
mean_tfr_left.plot("C3", baseline=(-1.0, 0), axes=axes[0, 0])
axes[0, 0].set_title("Left hand movement, channel C3")
mean_tfr_left.plot("C4", baseline=(-1.0, 0), axes=axes[0, 1])
axes[0, 1].set_title("Left hand movement, channel C4")
mean_tfr_right.plot("C3", baseline=(-1.0, 0), axes=axes[1, 0])
axes[1, 0].set_title("Right hand movement, channel C3")
mean_tfr_right.plot("C4", baseline=(-1.0, 0), axes=axes[1, 1])
axes[1, 1].set_title("Right hand movement, channel C4")
plt.tight_layout()

[gmj3]

Here is my code, from loading all the participants, checking a single subject data, and then plotting grand averages for left hand both electrodes and creating a difference map. I just need to add more preprocessing steps to clean the data, so any guidance would be appreciated.

[gmj3]

I sent you my steps. Please let me know if there is anything to add, expecially regarding preprocessing

[gmj3]

Not yet, but I will look at it.