<jupyter_start><jupyter_text>MNE-Python: From raw data to epochs and evoked responses (ERF/ERP)`Authors:<NAME> <NAME> <NAME>`<jupyter_code>import os
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt<jupyter_output><empty_output><jupyter_text>First, load the mne package:<jupyter_code>import mne<jupyter_output><empty_output><jupyter_text>We set the logging level to 'warning' so the output will be less verbose and we can focus on the important bits 🧐<jupyter_code>mne.set_log_level('warning')<jupyter_output><empty_output><jupyter_text>Remember: If you need help, just ask... the machine!<jupyter_code>mne.pick_types?<jupyter_output><empty_output><jupyter_text>Access raw data Now we import the `sample` dataset.It will be downloaded automatically (approx. 2 GB)<jupyter_code>from mne.datasets import sample
data_path = sample.data_path()

raw_fname = os.path.join(data_path, 'MEG/sample/sample_audvis_filt-0-40_raw.fif')
print(raw_fname)<jupyter_output><empty_output><jupyter_text>Read data from file:<jupyter_code>mne.io.read_raw_fif?
raw = mne.io.read_raw_fif(raw_fname, preload=True)
print(raw)<jupyter_output><empty_output><jupyter_text>Note that by default, the data will actually not be loaded into memory automatically to preserve memory. To actually load the data, we have to pass `preload=True`.Now let's look at the measurement info. It will give details about:   - sampling rate   - filtering parameters   - available channel types   - bad channels   - etc.<jupyter_code>print(raw.info)<jupyter_output><empty_output><jupyter_text>raw.info is just a dictionary:<jupyter_code>isinstance(raw.info, dict)<jupyter_output><empty_output><jupyter_text>So we can access its elements this way:<jupyter_code>raw.info['sfreq']  # Sampling frequency
raw.info['bads']  # list of marked bad channels<jupyter_output><empty_output><jupyter_text>Next let's see what channels are present in the data. We simply have to take a loot at the `raw.ch_names` attribute.<jupyter_code>raw.ch_names<jupyter_output><empty_output><jupyter_text>You can index it as a list<jupyter_code>raw.ch_names[42]
raw.ch_names[:10]<jupyter_output><empty_output><jupyter_text>Channel type of a specific channel<jupyter_code>channel_type = mne.io.pick.channel_type(info=raw.info, idx=75)
print('Channel #75 is of type:', channel_type)

channel_type = mne.io.pick.channel_type(info=raw.info, idx=320)
print('Channel #320 is of type:', channel_type)<jupyter_output><empty_output><jupyter_text>`raw.info['chs']` contains all the details about the sensors (type, locations, coordinate frame etc.)<jupyter_code>len(raw.info['chs'])
type(raw.info['chs'])
raw.info['chs'][0]
raw.info['chs'][330]
raw.plot_sensors(kind='topomap', ch_type='grad');<jupyter_output><empty_output><jupyter_text>Accessing the dataTo access the data, just use the `[]` syntax as to access any element of a list, dict etc.<jupyter_code>start, stop = 0, 10
data, times = raw[:, start:stop]  # fetch all channels and the first 10 time points
print(data.shape)
print(times.shape)
times<jupyter_output><empty_output><jupyter_text>Note that it returns both the data and the times array.  Visualizing raw dataNote : we will use the QT backend from matplotlib that will open a separate window.<jupyter_code>%matplotlib qt
fig = raw.plot()
fig = raw.copy().pick_types(meg=False, eeg=True).plot()<jupyter_output><empty_output><jupyter_text>Filtering<jupyter_code>raw = mne.io.read_raw_fif(raw_fname, preload=True)
raw_beta = raw.copy().filter(l_freq=13, h_freq=30, verbose=True)
print(raw_beta.info)  # note the update of raw.info['lowpass'] and raw.info['highpass']
raw_beta.filter?<jupyter_output><empty_output><jupyter_text>ExercisePlot the 10 first seconds of the stimutation channel `STI 014` just using matplotlib.Tips:- Find the channel index using `raw.ch_names.index('STI 014')`- Get the data for this channel- Plot it using `plt.plot`  Define and read epochs  First, extract events.<jupyter_code>events = mne.find_events(raw, stim_channel='STI 014', verbose=True)
events.shape
type(events)
print(events[:5])  # events is a 2d array, (time, previous, trigger)
len(events[events[:, 2] == 2])
len(events)<jupyter_output><empty_output><jupyter_text>Let's visualize the paradigm:<jupyter_code>%matplotlib inline
fig = mne.viz.plot_events(events, sfreq=raw.info['sfreq'])<jupyter_output><empty_output><jupyter_text>For describing which event IDs code for which experimental events or conditions, we use a Python dictionary to store the mapping. The dictionary keys can contain `/` for grouping of sub-conditions.<jupyter_code>event_id = {"visual/left": 3, "visual/right": 4,
            "auditory/left": 1, "auditory/right": 2}

fig = mne.viz.plot_events(events, sfreq=raw.info['sfreq'], event_id=event_id)<jupyter_output><empty_output><jupyter_text>The events can be visualized together with the raw data:<jupyter_code>fig = raw.plot(event_id=event_id, events=events)<jupyter_output><empty_output><jupyter_text>Create epochsFirst, define parameters: start, stop, and baseline period of the epochs.<jupyter_code>tmin = -0.2  # start of each epoch (200ms before the trigger)
tmax = 0.5   # end of each epoch (500ms after the trigger)
baseline = (None, 0)  # from the first time instant to the trigger pulse<jupyter_output><empty_output><jupyter_text>Define peak-to-peak (amplitude range) rejection parameters for gradiometers, magnetometers, and EOG:<jupyter_code>reject = dict(grad=4000e-13, mag=4e-12, eog=150e-6)  # this can be highly data dependent
# we are picky again, this time we select MEG and EOG channels
picks_meg = mne.pick_types(raw.info, meg=True, eeg=False, eog=True,
                           stim=False, exclude='bads')<jupyter_output><empty_output><jupyter_text>Extract epochs:<jupyter_code>epochs = mne.Epochs(raw, events, event_id, tmin, tmax, proj=True,
                    picks=picks_meg, baseline=baseline,
                    reject=reject)
print(epochs)<jupyter_output><empty_output><jupyter_text>Remove bad epochs based on the `reject` parameter we passed to `Epochs`.<jupyter_code>epochs.drop_bad()<jupyter_output><empty_output><jupyter_text>Explore the epochs namespaceHit ``epochs.`` See how epochs were dropped<jupyter_code>fig = epochs.plot_drop_log()<jupyter_output><empty_output><jupyter_text>Visualization EpochsSee [this page](https://mne.tools/stable/auto_tutorials/epochs/plot_visualize_epochs.html) for options on how to visualize epochs.Here is just an illustration to make a so-called ERP/ERF image:<jupyter_code>fig = raw.plot_psd(fmax=40)
figs = epochs.plot_image()
%matplotlib qt
fig = epochs.plot()<jupyter_output><empty_output><jupyter_text>Average the epochs to get the evoked response (ERF/ERP), and plot it!<jupyter_code>evoked = epochs.average()
print(evoked)
%matplotlib inline
fig = evoked.plot(spatial_colors=True)<jupyter_output><empty_output><jupyter_text>This created an average across **all** conditions. Let's now estimate evoked responses for **individual** conditions.<jupyter_code>print(event_id)
fig = epochs['auditory/left'].average().plot(spatial_colors=True)<jupyter_output><empty_output><jupyter_text>Accessing and indexing epochs by conditionEpochs can be indexed by integers or slices to select a subset of epochs but also with strings to select by conditions `epochs[condition]` Remember `/` serves as a grouping operator. To calculate the evoked response across **all** "left" stimulations, do the following:<jupyter_code>fig = epochs['left'].average().plot(spatial_colors=True);  # note the legend
# remember ...
event_id
epochs[0]  # first epoch
epochs[:10]  # first 10 epochs
epochs['visual/left']  # epochs for the left visual condition<jupyter_output><empty_output><jupyter_text>In event_id, `/` selects conditions in a hierarchical way, e.g. here, "auditory" vs. "visual", "left" vs. "right", and MNE can select them individually.<jupyter_code>evoked_auditory_left = epochs['auditory/left'].average().pick_types(meg='grad')
evoked_auditory_left.crop(None, 0.2) # Beginning of evoked until 0.2s after stimulus onset.
fig = evoked_auditory_left.plot(spatial_colors=True)
epochs['visual']  # epochs for the visual condition (either left or right)
epochs['left']<jupyter_output><empty_output><jupyter_text>To access the data of some epochs use the `get_data` method.<jupyter_code>epochs_data = epochs.get_data()
type(epochs_data), epochs_data.shape<jupyter_output><empty_output><jupyter_text>`epochs_data` is a 3D array of dimension (239 epochs, 306 channels, 106 time instants).  Visualize Topographies First, let's plot an illustration of the sensor locations on the scalp<jupyter_code>fig = raw.plot_sensors(ch_type='mag')
fig = raw.plot_sensors(ch_type='grad')<jupyter_output><empty_output><jupyter_text>Now, let's plot the scalp topography at different time points. These topography plots called `topomap` in the MNE jargon.<jupyter_code>fig = evoked.plot_topomap(ch_type='mag', times=[0.05, 0.1, 0.15])
fig = evoked.plot_topomap(ch_type='grad', times=[0.05, 0.1, 0.15])
import numpy as np


times = np.linspace(0.05, 0.15, 10)

for ch_type in ('mag', 'grad'):
    fig = evoked.plot_topomap(times=times, ch_type=ch_type)<jupyter_output><empty_output><jupyter_text>Topoplot and time series can also be shown in one single "joint" plot. By default, MNE will plot the topographies at the peaks of the time series signal.<jupyter_code>figs = evoked.plot_joint()<jupyter_output><empty_output><jupyter_text>But of course, you can also specify custom time points for the topomaps.<jupyter_code>figs = evoked.plot_joint(times=[0.1, 0.3])<jupyter_output><empty_output><jupyter_text>Let's visualize topomaps for all experimental conditions.<jupyter_code>for condition in event_id:
    fig = epochs[condition].average().plot_topomap(times=[0.1, 0.15], title=condition)<jupyter_output><empty_output><jupyter_text>Compute a contrast:<jupyter_code>evoked1 = epochs['left'].average()
evoked2 = epochs['right'].average()

contrast = mne.combine_evoked([evoked1, evoked2], weights=[1, -1])<jupyter_output><empty_output><jupyter_text>Note that this combines evokeds taking into account the number of averaged epochs (to scale the noise variance)<jupyter_code>print(evoked1.nave)  # average of 55 epochs
print(contrast.nave)  # average of 116 epochs
print(contrast)
fig = contrast.plot_joint()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11364735467646775, 'fraction_numerical': 0.013540961408259987, 'mean_word_length': 7.027950310559007, 'pattern_counts': {'":': 4, '<': 156, '<?xml version=': 0, '>': 156, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '38', 'dir': 'jupyter-structured-clean-dedup', 'id': '106648', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 3710, 'n_tokens_neox': 3494, 'n_words': 887, 'path': '1-From_raw_to_epochs_evoked_ERF_ERP.ipynb', 'repo_name': 'agramfort/AI4Health_ml_eeg', 'size': '10390'}