Compute and visualize ERDS maps

This example calculates and displays ERDS maps of event-related EEG data. ERDS (sometimes also written as ERD/ERS) is short for event-related desynchronization (ERD) and event-related synchronization (ERS) [1]. Conceptually, ERD corresponds to a decrease in power in a specific frequency band relative to a baseline. Similarly, ERS corresponds to an increase in power. An ERDS map is a time/frequency representation of ERD/ERS over a range of frequencies [2]. ERDS maps are also known as ERSP (event-related spectral perturbation) [3].

In this example, we use an EEG BCI data set containing two different motor imagery tasks (imagined hand and feet movement). Our goal is to generate ERDS maps for each of the two tasks.

First, we load the data and create epochs of 5s length. The data set contains multiple channels, but we will only consider C3, Cz, and C4. We compute maps containing frequencies ranging from 2 to 35Hz. We map ERD to red color and ERS to blue color, which is customary in many ERDS publications. Finally, we perform cluster-based permutation tests to estimate significant ERDS values (corrected for multiple comparisons within channels).

# Authors: Clemens Brunner <clemens.brunner@gmail.com>
#          Felix Klotzsche <klotzsche@cbs.mpg.de>
#
# License: BSD-3-Clause
# Copyright the MNE-Python contributors.

As usual, we import everything we need.

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
from matplotlib.colors import TwoSlopeNorm

import mne
from mne.datasets import eegbci
from mne.io import concatenate_raws, read_raw_edf
from mne.stats import permutation_cluster_1samp_test as pcluster_test

First, we load and preprocess the data. We use runs 6, 10, and 14 from subject 1 (these runs contains hand and feet motor imagery).

fnames = eegbci.load_data(subjects=1, runs=(6, 10, 14))
raw = concatenate_raws([read_raw_edf(f, preload=True) for f in fnames])

raw.rename_channels(lambda x: x.strip("."))  # remove dots from channel names
# rename descriptions to be more easily interpretable
raw.annotations.rename(dict(T1="hands", T2="feet"))

Extracting EDF parameters from /home/circleci/mne_data/MNE-eegbci-data/files/eegmmidb/1.0.0/S001/S001R06.edf...
EDF file detected
Setting channel info structure...
Creating raw.info structure...
Reading 0 ... 19999  =      0.000 ...   124.994 secs...
Extracting EDF parameters from /home/circleci/mne_data/MNE-eegbci-data/files/eegmmidb/1.0.0/S001/S001R10.edf...
EDF file detected
Setting channel info structure...
Creating raw.info structure...
Reading 0 ... 19999  =      0.000 ...   124.994 secs...
Extracting EDF parameters from /home/circleci/mne_data/MNE-eegbci-data/files/eegmmidb/1.0.0/S001/S001R14.edf...
EDF file detected
Setting channel info structure...
Creating raw.info structure...
Reading 0 ... 19999  =      0.000 ...   124.994 secs...

Now we can create 5-second epochs around events of interest.

tmin, tmax = -1, 4
event_ids = dict(hands=2, feet=3)  # map event IDs to tasks

epochs = mne.Epochs(
    raw,
    event_id=["hands", "feet"],
    tmin=tmin - 0.5,
    tmax=tmax + 0.5,
    picks=("C3", "Cz", "C4"),
    baseline=None,
    preload=True,
)

Used Annotations descriptions: [np.str_('T0'), np.str_('feet'), np.str_('hands')]
Ignoring annotation durations and creating fixed-duration epochs around annotation onsets.
Not setting metadata
45 matching events found
No baseline correction applied
0 projection items activated
Using data from preloaded Raw for 45 events and 961 original time points ...
0 bad epochs dropped

Here we set suitable values for computing ERDS maps. Note especially the cnorm variable, which sets up an asymmetric colormap where the middle color is mapped to zero, even though zero is not the middle value of the colormap range. This does two things: it ensures that zero values will be plotted in white (given that below we select the RdBu colormap), and it makes synchronization and desynchronization look equally prominent in the plots, even though their extreme values are of different magnitudes.

freqs = np.arange(2, 36)  # frequencies from 2-35Hz
vmin, vmax = -1, 1.5  # set min and max ERDS values in plot
baseline = (-1, 0)  # baseline interval (in s)
cnorm = TwoSlopeNorm(vmin=vmin, vcenter=0, vmax=vmax)  # min, center & max ERDS

kwargs = dict(
    n_permutations=100, step_down_p=0.05, seed=1, buffer_size=None, out_type="mask"
)  # for cluster test

Finally, we perform time/frequency decomposition over all epochs.

tfr = epochs.compute_tfr(
    method="multitaper",
    freqs=freqs,
    n_cycles=freqs,
    use_fft=True,
    return_itc=False,
    average=False,
    decim=2,
)
tfr.crop(tmin, tmax).apply_baseline(baseline, mode="percent")

for event in event_ids:
    # select desired epochs for visualization
    tfr_ev = tfr[event]
    fig, axes = plt.subplots(
        1, 4, figsize=(12, 4), gridspec_kw={"width_ratios": [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)

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

        ax.set_title(epochs.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()

• 
• 

Applying baseline correction (mode: percent)
Using a threshold of 1.724718
stat_fun(H1): min=-8.55207615547708 max=3.1832308002458687
Running initial clustering …
Found 80 clusters

  0%|          | Permuting : 0/99 [00:00<?,       ?it/s]
  9%|▉         | Permuting : 9/99 [00:00<00:00,  255.21it/s]
 18%|█▊        | Permuting : 18/99 [00:00<00:00,  261.18it/s]
 27%|██▋       | Permuting : 27/99 [00:00<00:00,  263.30it/s]
 35%|███▌      | Permuting : 35/99 [00:00<00:00,  256.51it/s]
 44%|████▍     | Permuting : 44/99 [00:00<00:00,  258.84it/s]
 55%|█████▍    | Permuting : 54/99 [00:00<00:00,  263.55it/s]
 64%|██████▎   | Permuting : 63/99 [00:00<00:00,  263.74it/s]
 73%|███████▎  | Permuting : 72/99 [00:00<00:00,  262.77it/s]
 83%|████████▎ | Permuting : 82/99 [00:00<00:00,  267.31it/s]
 93%|█████████▎| Permuting : 92/99 [00:00<00:00,  270.96it/s]
100%|██████████| Permuting : 99/99 [00:00<00:00,  272.90it/s]
100%|██████████| Permuting : 99/99 [00:00<00:00,  270.71it/s]
Step-down-in-jumps iteration #1 found 0 clusters to exclude from subsequent iterations
Using a threshold of -1.724718
stat_fun(H1): min=-8.55207615547708 max=3.1832308002458687
Running initial clustering …
Found 67 clusters

  0%|          | Permuting : 0/99 [00:00<?,       ?it/s]
  9%|▉         | Permuting : 9/99 [00:00<00:00,  249.95it/s]
 18%|█▊        | Permuting : 18/99 [00:00<00:00,  258.27it/s]
 27%|██▋       | Permuting : 27/99 [00:00<00:00,  261.31it/s]
 36%|███▋      | Permuting : 36/99 [00:00<00:00,  262.81it/s]
 45%|████▌     | Permuting : 45/99 [00:00<00:00,  263.75it/s]
 55%|█████▍    | Permuting : 54/99 [00:00<00:00,  264.37it/s]
 64%|██████▎   | Permuting : 63/99 [00:00<00:00,  263.09it/s]
 72%|███████▏  | Permuting : 71/99 [00:00<00:00,  259.32it/s]
 81%|████████  | Permuting : 80/99 [00:00<00:00,  260.37it/s]
 90%|████████▉ | Permuting : 89/99 [00:00<00:00,  261.22it/s]
 99%|█████████▉| Permuting : 98/99 [00:00<00:00,  259.92it/s]
100%|██████████| Permuting : 99/99 [00:00<00:00,  261.92it/s]
Step-down-in-jumps iteration #1 found 1 cluster to exclude from subsequent iterations

  0%|          | Permuting : 0/99 [00:00<?,       ?it/s]
  9%|▉         | Permuting : 9/99 [00:00<00:00,  248.07it/s]
 18%|█▊        | Permuting : 18/99 [00:00<00:00,  257.10it/s]
 27%|██▋       | Permuting : 27/99 [00:00<00:00,  260.53it/s]
 37%|███▋      | Permuting : 37/99 [00:00<00:00,  270.13it/s]
 46%|████▋     | Permuting : 46/99 [00:00<00:00,  269.47it/s]
 55%|█████▍    | Permuting : 54/99 [00:00<00:00,  263.48it/s]
 64%|██████▎   | Permuting : 63/99 [00:00<00:00,  264.10it/s]
 73%|███████▎  | Permuting : 72/99 [00:00<00:00,  264.55it/s]
 82%|████████▏ | Permuting : 81/99 [00:00<00:00,  263.47it/s]
 91%|█████████ | Permuting : 90/99 [00:00<00:00,  263.93it/s]
 99%|█████████▉| Permuting : 98/99 [00:00<00:00,  260.87it/s]
100%|██████████| Permuting : 99/99 [00:00<00:00,  263.31it/s]
Step-down-in-jumps iteration #2 found 0 additional clusters to exclude from subsequent iterations
No baseline correction applied
Using a threshold of 1.724718
stat_fun(H1): min=-4.528367198129923 max=3.7064219374431917
Running initial clustering …
Found 88 clusters

  0%|          | Permuting : 0/99 [00:00<?,       ?it/s]
  9%|▉         | Permuting : 9/99 [00:00<00:00,  248.82it/s]
 18%|█▊        | Permuting : 18/99 [00:00<00:00,  257.94it/s]
 28%|██▊       | Permuting : 28/99 [00:00<00:00,  271.35it/s]
 36%|███▋      | Permuting : 36/99 [00:00<00:00,  262.39it/s]
 46%|████▋     | Permuting : 46/99 [00:00<00:00,  269.85it/s]
 56%|█████▌    | Permuting : 55/99 [00:00<00:00,  269.32it/s]
 66%|██████▌   | Permuting : 65/99 [00:00<00:00,  273.88it/s]
 75%|███████▍  | Permuting : 74/99 [00:00<00:00,  272.90it/s]
 85%|████████▍ | Permuting : 84/99 [00:00<00:00,  274.44it/s]
 95%|█████████▍| Permuting : 94/99 [00:00<00:00,  277.23it/s]
100%|██████████| Permuting : 99/99 [00:00<00:00,  277.55it/s]
100%|██████████| Permuting : 99/99 [00:00<00:00,  275.91it/s]
Step-down-in-jumps iteration #1 found 1 cluster to exclude from subsequent iterations

  0%|          | Permuting : 0/99 [00:00<?,       ?it/s]
  9%|▉         | Permuting : 9/99 [00:00<00:00,  247.63it/s]
 18%|█▊        | Permuting : 18/99 [00:00<00:00,  257.19it/s]
 28%|██▊       | Permuting : 28/99 [00:00<00:00,  270.84it/s]
 36%|███▋      | Permuting : 36/99 [00:00<00:00,  262.08it/s]
 46%|████▋     | Permuting : 46/99 [00:00<00:00,  269.69it/s]
 56%|█████▌    | Permuting : 55/99 [00:00<00:00,  269.20it/s]
 65%|██████▍   | Permuting : 64/99 [00:00<00:00,  268.80it/s]
 74%|███████▎  | Permuting : 73/99 [00:00<00:00,  268.58it/s]
 84%|████████▍ | Permuting : 83/99 [00:00<00:00,  272.38it/s]
 94%|█████████▍| Permuting : 93/99 [00:00<00:00,  274.58it/s]
100%|██████████| Permuting : 99/99 [00:00<00:00,  277.66it/s]
100%|██████████| Permuting : 99/99 [00:00<00:00,  275.48it/s]
Step-down-in-jumps iteration #2 found 0 additional clusters to exclude from subsequent iterations
Using a threshold of -1.724718
stat_fun(H1): min=-4.528367198129923 max=3.7064219374431917
Running initial clustering …
Found 58 clusters

  0%|          | Permuting : 0/99 [00:00<?,       ?it/s]
  9%|▉         | Permuting : 9/99 [00:00<00:00,  252.84it/s]
 19%|█▉        | Permuting : 19/99 [00:00<00:00,  274.86it/s]
 28%|██▊       | Permuting : 28/99 [00:00<00:00,  272.28it/s]
 38%|███▊      | Permuting : 38/99 [00:00<00:00,  278.85it/s]
 47%|████▋     | Permuting : 47/99 [00:00<00:00,  276.34it/s]
 57%|█████▋    | Permuting : 56/99 [00:00<00:00,  274.46it/s]
 66%|██████▌   | Permuting : 65/99 [00:00<00:00,  273.29it/s]
 75%|███████▍  | Permuting : 74/99 [00:00<00:00,  269.79it/s]
 84%|████████▍ | Permuting : 83/99 [00:00<00:00,  269.45it/s]
 93%|█████████▎| Permuting : 92/99 [00:00<00:00,  269.18it/s]
100%|██████████| Permuting : 99/99 [00:00<00:00,  272.10it/s]
100%|██████████| Permuting : 99/99 [00:00<00:00,  271.85it/s]
Step-down-in-jumps iteration #1 found 0 clusters to exclude from subsequent iterations
No baseline correction applied
Using a threshold of 1.724718
stat_fun(H1): min=-6.581589289696407 max=3.346448151996899
Running initial clustering …
Found 67 clusters

  0%|          | Permuting : 0/99 [00:00<?,       ?it/s]
  9%|▉         | Permuting : 9/99 [00:00<00:00,  255.88it/s]
 18%|█▊        | Permuting : 18/99 [00:00<00:00,  261.47it/s]
 27%|██▋       | Permuting : 27/99 [00:00<00:00,  263.38it/s]
 36%|███▋      | Permuting : 36/99 [00:00<00:00,  264.40it/s]
 45%|████▌     | Permuting : 45/99 [00:00<00:00,  260.85it/s]
 56%|█████▌    | Permuting : 55/99 [00:00<00:00,  267.56it/s]
 66%|██████▌   | Permuting : 65/99 [00:00<00:00,  272.36it/s]
 75%|███████▍  | Permuting : 74/99 [00:00<00:00,  271.63it/s]
 84%|████████▍ | Permuting : 83/99 [00:00<00:00,  271.05it/s]
 94%|█████████▍| Permuting : 93/99 [00:00<00:00,  274.25it/s]
100%|██████████| Permuting : 99/99 [00:00<00:00,  274.96it/s]
100%|██████████| Permuting : 99/99 [00:00<00:00,  273.35it/s]
Step-down-in-jumps iteration #1 found 0 clusters to exclude from subsequent iterations
Using a threshold of -1.724718
stat_fun(H1): min=-6.581589289696407 max=3.346448151996899
Running initial clustering …
Found 69 clusters

  0%|          | Permuting : 0/99 [00:00<?,       ?it/s]
  9%|▉         | Permuting : 9/99 [00:00<00:00,  259.07it/s]
 19%|█▉        | Permuting : 19/99 [00:00<00:00,  278.13it/s]
 28%|██▊       | Permuting : 28/99 [00:00<00:00,  274.38it/s]
 37%|███▋      | Permuting : 37/99 [00:00<00:00,  272.53it/s]
 47%|████▋     | Permuting : 47/99 [00:00<00:00,  273.71it/s]
 57%|█████▋    | Permuting : 56/99 [00:00<00:00,  272.53it/s]
 66%|██████▌   | Permuting : 65/99 [00:00<00:00,  271.66it/s]
 75%|███████▍  | Permuting : 74/99 [00:00<00:00,  271.08it/s]
 85%|████████▍ | Permuting : 84/99 [00:00<00:00,  274.56it/s]
 94%|█████████▍| Permuting : 93/99 [00:00<00:00,  273.67it/s]
100%|██████████| Permuting : 99/99 [00:00<00:00,  275.52it/s]
100%|██████████| Permuting : 99/99 [00:00<00:00,  274.71it/s]
Step-down-in-jumps iteration #1 found 0 clusters to exclude from subsequent iterations
No baseline correction applied
Using a threshold of 1.713872
stat_fun(H1): min=-3.754759498497248 max=3.3607039428313787
Running initial clustering …
Found 71 clusters

  0%|          | Permuting : 0/99 [00:00<?,       ?it/s]
  9%|▉         | Permuting : 9/99 [00:00<00:00,  252.17it/s]
 18%|█▊        | Permuting : 18/99 [00:00<00:00,  259.53it/s]
 26%|██▋       | Permuting : 26/99 [00:00<00:00,  251.85it/s]
 34%|███▍      | Permuting : 34/99 [00:00<00:00,  248.08it/s]
 43%|████▎     | Permuting : 43/99 [00:00<00:00,  252.30it/s]
 53%|█████▎    | Permuting : 52/99 [00:00<00:00,  255.07it/s]
 62%|██████▏   | Permuting : 61/99 [00:00<00:00,  257.09it/s]
 71%|███████   | Permuting : 70/99 [00:00<00:00,  257.25it/s]
 79%|███████▉  | Permuting : 78/99 [00:00<00:00,  254.61it/s]
 88%|████████▊ | Permuting : 87/99 [00:00<00:00,  256.18it/s]
 97%|█████████▋| Permuting : 96/99 [00:00<00:00,  257.47it/s]
100%|██████████| Permuting : 99/99 [00:00<00:00,  257.29it/s]
Step-down-in-jumps iteration #1 found 0 clusters to exclude from subsequent iterations
Using a threshold of -1.713872
stat_fun(H1): min=-3.754759498497248 max=3.3607039428313787
Running initial clustering …
Found 80 clusters

  0%|          | Permuting : 0/99 [00:00<?,       ?it/s]
  8%|▊         | Permuting : 8/99 [00:00<00:00,  224.84it/s]
 17%|█▋        | Permuting : 17/99 [00:00<00:00,  246.03it/s]
 25%|██▌       | Permuting : 25/99 [00:00<00:00,  243.05it/s]
 34%|███▍      | Permuting : 34/99 [00:00<00:00,  249.16it/s]
 44%|████▍     | Permuting : 44/99 [00:00<00:00,  255.11it/s]
 54%|█████▎    | Permuting : 53/99 [00:00<00:00,  257.35it/s]
 63%|██████▎   | Permuting : 62/99 [00:00<00:00,  258.96it/s]
 73%|███████▎  | Permuting : 72/99 [00:00<00:00,  264.28it/s]
 82%|████████▏ | Permuting : 81/99 [00:00<00:00,  264.68it/s]
 92%|█████████▏| Permuting : 91/99 [00:00<00:00,  266.29it/s]
100%|██████████| Permuting : 99/99 [00:00<00:00,  267.51it/s]
100%|██████████| Permuting : 99/99 [00:00<00:00,  265.27it/s]
Step-down-in-jumps iteration #1 found 0 clusters to exclude from subsequent iterations
No baseline correction applied
Using a threshold of 1.713872
stat_fun(H1): min=-4.992503139561789 max=5.416450269990024
Running initial clustering …
Found 103 clusters

  0%|          | Permuting : 0/99 [00:00<?,       ?it/s]
  9%|▉         | Permuting : 9/99 [00:00<00:00,  248.83it/s]
 19%|█▉        | Permuting : 19/99 [00:00<00:00,  272.59it/s]
 28%|██▊       | Permuting : 28/99 [00:00<00:00,  270.80it/s]
 37%|███▋      | Permuting : 37/99 [00:00<00:00,  269.91it/s]
 45%|████▌     | Permuting : 45/99 [00:00<00:00,  262.82it/s]
 56%|█████▌    | Permuting : 55/99 [00:00<00:00,  269.16it/s]
 65%|██████▍   | Permuting : 64/99 [00:00<00:00,  268.86it/s]
 74%|███████▎  | Permuting : 73/99 [00:00<00:00,  268.65it/s]
 83%|████████▎ | Permuting : 82/99 [00:00<00:00,  268.45it/s]
 92%|█████████▏| Permuting : 91/99 [00:00<00:00,  266.44it/s]
100%|██████████| Permuting : 99/99 [00:00<00:00,  269.09it/s]
100%|██████████| Permuting : 99/99 [00:00<00:00,  268.59it/s]
Step-down-in-jumps iteration #1 found 1 cluster to exclude from subsequent iterations

  0%|          | Permuting : 0/99 [00:00<?,       ?it/s]
  9%|▉         | Permuting : 9/99 [00:00<00:00,  249.83it/s]
 19%|█▉        | Permuting : 19/99 [00:00<00:00,  273.00it/s]
 27%|██▋       | Permuting : 27/99 [00:00<00:00,  260.86it/s]
 37%|███▋      | Permuting : 37/99 [00:00<00:00,  270.50it/s]
 46%|████▋     | Permuting : 46/99 [00:00<00:00,  265.43it/s]
 56%|█████▌    | Permuting : 55/99 [00:00<00:00,  265.78it/s]
 66%|██████▌   | Permuting : 65/99 [00:00<00:00,  270.85it/s]
 75%|███████▍  | Permuting : 74/99 [00:00<00:00,  270.34it/s]
 83%|████████▎ | Permuting : 82/99 [00:00<00:00,  265.94it/s]
 92%|█████████▏| Permuting : 91/99 [00:00<00:00,  266.09it/s]
100%|██████████| Permuting : 99/99 [00:00<00:00,  270.86it/s]
100%|██████████| Permuting : 99/99 [00:00<00:00,  269.83it/s]
Step-down-in-jumps iteration #2 found 0 additional clusters to exclude from subsequent iterations
Using a threshold of -1.713872
stat_fun(H1): min=-4.992503139561789 max=5.416450269990024
Running initial clustering …
Found 67 clusters

  0%|          | Permuting : 0/99 [00:00<?,       ?it/s]
  9%|▉         | Permuting : 9/99 [00:00<00:00,  251.01it/s]
 18%|█▊        | Permuting : 18/99 [00:00<00:00,  259.03it/s]
 27%|██▋       | Permuting : 27/99 [00:00<00:00,  261.88it/s]
 36%|███▋      | Permuting : 36/99 [00:00<00:00,  263.32it/s]
 45%|████▌     | Permuting : 45/99 [00:00<00:00,  264.20it/s]
 56%|█████▌    | Permuting : 55/99 [00:00<00:00,  270.34it/s]
 64%|██████▎   | Permuting : 63/99 [00:00<00:00,  264.94it/s]
 73%|███████▎  | Permuting : 72/99 [00:00<00:00,  265.29it/s]
 83%|████████▎ | Permuting : 82/99 [00:00<00:00,  268.24it/s]
 92%|█████████▏| Permuting : 91/99 [00:00<00:00,  268.14it/s]
100%|██████████| Permuting : 99/99 [00:00<00:00,  267.49it/s]
100%|██████████| Permuting : 99/99 [00:00<00:00,  266.91it/s]
Step-down-in-jumps iteration #1 found 0 clusters to exclude from subsequent iterations
No baseline correction applied
Using a threshold of 1.713872
stat_fun(H1): min=-6.044339991978302 max=4.070443782271568
Running initial clustering …
Found 92 clusters

  0%|          | Permuting : 0/99 [00:00<?,       ?it/s]
  9%|▉         | Permuting : 9/99 [00:00<00:00,  248.84it/s]
 18%|█▊        | Permuting : 18/99 [00:00<00:00,  257.85it/s]
 27%|██▋       | Permuting : 27/99 [00:00<00:00,  261.11it/s]
 35%|███▌      | Permuting : 35/99 [00:00<00:00,  254.90it/s]
 45%|████▌     | Permuting : 45/99 [00:00<00:00,  264.09it/s]
 55%|█████▍    | Permuting : 54/99 [00:00<00:00,  264.70it/s]
 64%|██████▎   | Permuting : 63/99 [00:00<00:00,  265.14it/s]
 73%|███████▎  | Permuting : 72/99 [00:00<00:00,  265.48it/s]
 82%|████████▏ | Permuting : 81/99 [00:00<00:00,  265.71it/s]
 92%|█████████▏| Permuting : 91/99 [00:00<00:00,  268.63it/s]
100%|██████████| Permuting : 99/99 [00:00<00:00,  267.31it/s]
100%|██████████| Permuting : 99/99 [00:00<00:00,  266.39it/s]
Step-down-in-jumps iteration #1 found 1 cluster to exclude from subsequent iterations

  0%|          | Permuting : 0/99 [00:00<?,       ?it/s]
  9%|▉         | Permuting : 9/99 [00:00<00:00,  248.00it/s]
 18%|█▊        | Permuting : 18/99 [00:00<00:00,  257.56it/s]
 27%|██▋       | Permuting : 27/99 [00:00<00:00,  260.93it/s]
 35%|███▌      | Permuting : 35/99 [00:00<00:00,  254.64it/s]
 45%|████▌     | Permuting : 45/99 [00:00<00:00,  263.90it/s]
 55%|█████▍    | Permuting : 54/99 [00:00<00:00,  264.54it/s]
 64%|██████▎   | Permuting : 63/99 [00:00<00:00,  265.05it/s]
 73%|███████▎  | Permuting : 72/99 [00:00<00:00,  265.38it/s]
 82%|████████▏ | Permuting : 81/99 [00:00<00:00,  265.64it/s]
 91%|█████████ | Permuting : 90/99 [00:00<00:00,  265.86it/s]
100%|██████████| Permuting : 99/99 [00:00<00:00,  268.60it/s]
100%|██████████| Permuting : 99/99 [00:00<00:00,  267.25it/s]
Step-down-in-jumps iteration #2 found 0 additional clusters to exclude from subsequent iterations
Using a threshold of -1.713872
stat_fun(H1): min=-6.044339991978302 max=4.070443782271568
Running initial clustering …
Found 51 clusters

  0%|          | Permuting : 0/99 [00:00<?,       ?it/s]
  9%|▉         | Permuting : 9/99 [00:00<00:00,  250.13it/s]
 18%|█▊        | Permuting : 18/99 [00:00<00:00,  258.58it/s]
 28%|██▊       | Permuting : 28/99 [00:00<00:00,  271.77it/s]
 37%|███▋      | Permuting : 37/99 [00:00<00:00,  270.63it/s]
 46%|████▋     | Permuting : 46/99 [00:00<00:00,  269.91it/s]
 56%|█████▌    | Permuting : 55/99 [00:00<00:00,  269.46it/s]
 64%|██████▎   | Permuting : 63/99 [00:00<00:00,  264.21it/s]
 73%|███████▎  | Permuting : 72/99 [00:00<00:00,  264.68it/s]
 82%|████████▏ | Permuting : 81/99 [00:00<00:00,  264.31it/s]
 91%|█████████ | Permuting : 90/99 [00:00<00:00,  264.57it/s]
100%|██████████| Permuting : 99/99 [00:00<00:00,  267.03it/s]
100%|██████████| Permuting : 99/99 [00:00<00:00,  266.59it/s]
Step-down-in-jumps iteration #1 found 0 clusters to exclude from subsequent iterations
No baseline correction applied

Similar to Epochs objects, we can also export data from EpochsTFR and AverageTFR objects to a Pandas DataFrame. By default, the time column of the exported data frame is in milliseconds. Here, to be consistent with the time-frequency plots, we want to keep it in seconds, which we can achieve by setting time_format=None:

df = tfr.to_data_frame(time_format=None)
df.head()

	time	freq	epoch	condition	C3	Cz	C4
0	-1.0000	2.0	1	feet	-0.473093	-0.474571	-0.764021
1	-0.9875	2.0	1	feet	-0.491659	-0.480088	-0.762445
2	-0.9750	2.0	1	feet	-0.500872	-0.469376	-0.749421
3	-0.9625	2.0	1	feet	-0.500605	-0.449373	-0.728686
4	-0.9500	2.0	1	feet	-0.515911	-0.443439	-0.719536

This allows us to use additional plotting functions like seaborn.lineplot() to plot confidence bands:

df = tfr.to_data_frame(time_format=None, long_format=True)

# Map to frequency bands:
freq_bounds = {"_": 0, "delta": 3, "theta": 7, "alpha": 13, "beta": 35, "gamma": 140}
df["band"] = pd.cut(
    df["freq"], list(freq_bounds.values()), labels=list(freq_bounds)[1:]
)

# Filter to retain only relevant frequency bands:
freq_bands_of_interest = ["delta", "theta", "alpha", "beta"]
df = df[df.band.isin(freq_bands_of_interest)]
df["band"] = df["band"].cat.remove_unused_categories()

# Order channels for plotting:
df["channel"] = df["channel"].cat.reorder_categories(("C3", "Cz", "C4"), ordered=True)

g = sns.FacetGrid(df, row="band", col="channel", margin_titles=True)
g.map(sns.lineplot, "time", "value", "condition", n_boot=10)
axline_kw = dict(color="black", linestyle="dashed", linewidth=0.5, alpha=0.5)
g.map(plt.axhline, y=0, **axline_kw)
g.map(plt.axvline, x=0, **axline_kw)
g.set(ylim=(None, 1.5))
g.set_axis_labels("Time (s)", "ERDS")
g.set_titles(col_template="{col_name}", row_template="{row_name}")
g.add_legend(ncol=2, loc="lower center")
g.fig.subplots_adjust(left=0.1, right=0.9, top=0.9, bottom=0.08)

Converting "condition" to "category"...
Converting "epoch" to "category"...
Converting "channel" to "category"...
Converting "ch_type" to "category"...

Having the data as a DataFrame also facilitates subsetting, grouping, and other transforms. Here, we use seaborn to plot the average ERDS in the motor imagery interval as a function of frequency band and imagery condition:

df_mean = (
    df.query("time > 1")
    .groupby(["condition", "epoch", "band", "channel"], observed=False)[["value"]]
    .mean()
    .reset_index()
)

g = sns.FacetGrid(
    df_mean, col="condition", col_order=["hands", "feet"], margin_titles=True
)
g = g.map(
    sns.violinplot,
    "channel",
    "value",
    "band",
    cut=0,
    palette="deep",
    order=["C3", "Cz", "C4"],
    hue_order=freq_bands_of_interest,
    linewidth=0.5,
).add_legend(ncol=4, loc="lower center")

g.map(plt.axhline, **axline_kw)
g.set_axis_labels("", "ERDS")
g.set_titles(col_template="{col_name}", row_template="{row_name}")
g.fig.subplots_adjust(left=0.1, right=0.9, top=0.9, bottom=0.3)

References

[1]

Gert Pfurtscheller and Fernando H. Lopes da Silva. Event-related EEG/MEG synchronization and desynchronization: basic principles. Clinical Neurophysiology, 110(11):1842–1857, 1999. doi:10.1016/S1388-2457(99)00141-8.

[2]

Bernhard Graimann, Jane E. Huggins, Simon P. Levine, and Gert Pfurtscheller. Visualization of significant ERD/ERS patterns in multichannel EEG and ECoG data. Clinical Neurophysiology, 113(1):43–47, 2002. doi:10.1016/S1388-2457(01)00697-6.

[3]

Scott Makeig. Auditory event-related dynamics of the EEG spectrum and effects of exposure to tones. Electroencephalography and Clinical Neurophysiology, 86(4):283–293, 1993. doi:10.1016/0013-4694(93)90110-H.