diff --git a/.circleci/config.yml b/.circleci/config.yml index 46b55f2ea3c..ec5a3cca56f 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -193,7 +193,7 @@ jobs: python -c "import mne; print(mne.datasets.phantom_4dbti.data_path(update_path=True))"; fi; if [[ $(cat $FNAME | grep -x ".*datasets.*limo.*" | wc -l) -gt 0 ]]; then - python -c "import mne; print(mne.datasets.limo.data_path(subject=2, update_path=True))"; + python -c "import mne; print(mne.datasets.limo.data_path(subject=1, update_path=True))"; fi; fi; done; diff --git a/doc/changes/latest.inc b/doc/changes/latest.inc index a111bf2c0d7..2b10d716e71 100644 --- a/doc/changes/latest.inc +++ b/doc/changes/latest.inc @@ -42,6 +42,8 @@ Changelog - Add reader for ``*.dat`` electrode position files :func:`mne.channels.read_dig_dat` by `Christian Brodbeck`_ +- Improved :ref:`limo-dataset` usage and :ref:`example ` for usage of :func:`mne.stats.linear_regression` by `Jose Alanis`_ + - For KIT systems without built-in layout, :func:`mne.channels.find_layout` now falls back on an automatically generated layout, by `Christian Brodbeck`_ - :meth:`mne.Epochs.plot` now takes a ``epochs_colors`` parameter to color specific epoch segments by `Mainak Jas`_ diff --git a/doc/changes/names.inc b/doc/changes/names.inc index 9c61f099fa6..0a2258a45a0 100644 --- a/doc/changes/names.inc +++ b/doc/changes/names.inc @@ -270,6 +270,8 @@ .. _Fahimeh Mamashli: https://github.com/fmamashli +.. _Jose Alanis: https://github.com/JoseAlanis + .. _Adonay Nunes: https://github.com/AdoNunes .. _Victor Ferat: https://github.com/vferat diff --git a/doc/overview/datasets_index.rst b/doc/overview/datasets_index.rst index f04b5bfea6c..3f9bd56d2b0 100644 --- a/doc/overview/datasets_index.rst +++ b/doc/overview/datasets_index.rst @@ -345,6 +345,32 @@ functions in MNE and does not contain useful metadata for analysis. how to project a 3D electrode location onto a 2D image, a common procedure in electrocorticography. +.. _limo-dataset: + +LIMO Dataset +^^^^^^^^^^^^ +:func:`mne.datasets.limo.load_data`. + +In the original LIMO experiment (see [9]_), participants performed a +two-alternative forced choice task, discriminating between two face stimuli. +Subjects discriminated the same two faces during the whole experiment. +The critical manipulation consisted of the level of noise added to the +face-stimuli during the task, making the faces more or less discernible to the +observer. + +The presented faces varied across a noise-signal (or phase-coherence) continuum +spanning from 0 to 100% in increasing steps of 10%. In other words, faces with +high phase-coherence (e.g., 90%) were easy to identify, while faces with low +phase-coherence (e.g., 10%) were hard to identify and by extension hard to +discriminate. + +.. topic:: Examples + + * :ref:`Single trial linear regression analysis with the LIMO dataset + `: Explores data from a single subject of the LIMO dataset + and demonstrates how to fit a single trial linear regression using the + information contained in the metadata of the individual datasets. + References ========== @@ -362,6 +388,7 @@ References .. [8] Goldberger AL, Amaral LAN, Glass L, Hausdorff JM, Ivanov PCh, Mark RG, Mietus JE, Moody GB, Peng C-K, Stanley HE. PhysioBank, PhysioToolkit, and PhysioNet: Components of a New Research Resource for Complex Physiologic Signals. Circulation 101(23):e215-e220 [Circulation Electronic Pages; http://circ.ahajournals.org/cgi/content/full/101/23/e215]; 2000 (June 13). +.. [9] Rousselet, G. A., Gaspar, C. M., Pernet, C. R., Husk, J. S., Bennett, P. J., & Sekuler, A. B. (2010). Healthy aging delays scalp EEG sensitivity to noise in a face discrimination task. Frontiers in psychology, 1, 19. https://doi.org/10.3389/fpsyg.2010.00019 .. _auditory dataset tutorial: https://neuroimage.usc.edu/brainstorm/DatasetAuditory .. _resting state dataset tutorial: https://neuroimage.usc.edu/brainstorm/DatasetResting diff --git a/examples/datasets/plot_limo_data.py b/examples/datasets/plot_limo_data.py index d5c53f86521..038b140fab3 100644 --- a/examples/datasets/plot_limo_data.py +++ b/examples/datasets/plot_limo_data.py @@ -1,28 +1,30 @@ """ -======================================================== -Import single subject from LIMO data set into MNE-Python -======================================================== +.. _ex-limo-data: -Here we a define function to extract the eeg signal data -from the LIMO structures in the LIMO dataset, see [1]_, [2]_, and: +============================================================= +Single trial linear regression analysis with the LIMO dataset +============================================================= - https://datashare.is.ed.ac.uk/handle/10283/2189?show=full +Here we explore the structure of the data contained in the +`LIMO dataset`_. +This example replicates and extends some of the main analysis +and tools integrated in `LIMO MEEG`_, a MATLAB toolbox originally designed +to interface with EEGLAB_. - https://github.com/LIMO-EEG-Toolbox +In summary, the example: -The code allows to: +- Fetches epoched data files for a single subject of the LIMO dataset [1]_. + If the LIMO files are not found on disk, the + fetcher :func:`mne.datasets.limo.load_data()` will automatically download + the files from a remote repository. -Fetch single subjects epochs data for the LIMO data set. -Epochs information (i.e., sampling rate, number of epochs per condition, -number and name of EEG channels per subject, etc.) is extracted from -the LIMO .mat files stored on disk. -If files are not found, the function mne.datasets.limo.load_data() will -automatically download the data from a remote repository. - -:func:`mne.datasets.limo.load_data` creates a custom info and -epochs structure in MNE-Python. -Missing channels can be interpolated if desired. +- During import, information about the data (i.e., sampling rate, number of + epochs per condition, number and name of EEG channels per subject, etc.) is + extracted from the LIMO :file:`.mat` files stored on disk and added to the + epochs structure as metadata. +- Fits linear models on the single subject's data and visualizes inferential + measures to evaluate the significance of the estimated effects. References ---------- @@ -33,8 +35,16 @@ Bennett, P. J., & Sekuler, A. B. (2010). Healthy aging delays scalp EEG sensitivity to noise in a face discrimination task. Frontiers in psychology, 1, 19. https://doi.org/10.3389/fpsyg.2010.00019 -""" - +.. [3] Rousselet, G. A., Pernet, C. R., Bennett, P. J., & Sekuler, A. B. + (2008). Parametric study of EEG sensitivity to phase noise during face + processing. BMC neuroscience, 9(1), 98. + https://doi.org/10.1186/1471-2202-9-98 +.. _LIMO dataset: https://datashare.is.ed.ac.uk/handle/10283/2189?show=full +.. _LIMO MEEG: https://github.com/LIMO-EEG-Toolbox +.. _EEGLAB: https://sccn.ucsd.edu/eeglab/index.php +.. _Fig 1: https://bmcneurosci.biomedcentral.com/articles/10.1186/1471-2202-9-98/figures/1 +.. _least squares: https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.lstsq.html +""" # noqa: E501 # Authors: Jose C. Garcia Alanis # # License: BSD (3-clause) @@ -42,68 +52,105 @@ import numpy as np import matplotlib.pyplot as plt -import mne -from mne.datasets import limo +from mne.datasets.limo import load_data from mne.stats import linear_regression +from mne.viz import plot_events, plot_compare_evokeds +from mne import combine_evoked + print(__doc__) -# fetch data from subject 2 and interpolate missing channels -limo_epochs = limo.load_data(subject=2) +# subject to use +subj = 1 + +############################################################################### +# About the data +# -------------- +# +# In the original LIMO experiment (see [2]_), participants performed a +# two-alternative forced choice task, discriminating between two face stimuli. +# The same two faces were used during the whole experiment, +# with varying levels of noise added, making the faces more or less +# discernible to the observer (see `Fig 1`_ in [3]_ for a similar approach). +# +# The presented faces varied across a noise-signal (or phase-coherence) +# continuum spanning from 0 to 85% in increasing steps of 5%. +# In other words, faces with high phase-coherence (e.g., 85%) were easy to +# identify, while faces with low phase-coherence (e.g., 5%) were hard to +# identify and by extension very hard to discriminate. +# +# +# Load the data +# ------------- +# +# We'll begin by loading the data from subject 1 of the LIMO dataset. + +# This step can take a little while if you're loading the data for the +# first time. +limo_epochs = load_data(subject=subj) + +############################################################################### +# Note that the result of the loading process is an +# :class:`mne.EpochsArray` containing the data ready to interface +# with MNE-Python. + +print(limo_epochs) ############################################################################### -# In the original LIMO experiment, participants performed a two-alternative -# forced choice task discriminating between the same two faces. -# The critical manipulation in the experiment was that the phase-coherence of -# the presented face-stimuli was varied across a noise-signal continuum -# spanning from 0 to 100 %. In other words, faces with high phase coherence -# were easily discernible, while faces with low phase-coherence were hard to -# identify as such). -# The events coding the presentation of each of these two faces are stored in -# ``limo_epochs.events``. +# Visualize events +# ---------------- # # We can visualise the distribution of the face events contained in the -# epochs structure. Events should appear clearly grouped, as they are ordered -# during the import process. -fig, ax = plt.subplots(figsize=(7, 5)) -mne.viz.plot_events(limo_epochs.events, event_id=limo_epochs.event_id, axes=ax) -ax.set(title="Distribution of events") -plt.legend(loc='lower left', borderaxespad=1.) -plt.tight_layout() -plt.show() - -############################################################################### -# As it can be seen above, events are coded as ``Face/A`` and ``Face/B``. -# Information about the phase-coherence of the presented faces is stored in -# ``limo_epochs.metadata``, which also contains information about the presented -# faces for convenience. -# -# Check epochs metadata +# ``limo_epochs`` structure. Events should appear clearly grouped, as the +# epochs are ordered by condition. + +fig = plot_events(limo_epochs.events, event_id=limo_epochs.event_id) +fig.suptitle("Distribution of events in LIMO epochs") + +############################################################################### +# As it can be seen above, conditions are coded as ``Face/A`` and ``Face/B``. +# Information about the phase-coherence of the presented faces is stored in the +# epochs metadata. These information can be easily accessed by calling +# ``limo_epochs.metadata``. As shown below, the epochs metadata also contains +# information about the presented faces for convenience. + print(limo_epochs.metadata.head()) ############################################################################### -# Before going on, we'll drop the EOG channels present in the LIMO epochs -# (coded with EXG-prefix in ``limo_epochs.info['ch_names']``, as data has -# already been cleaned. -limo_epochs.drop_channels(['EXG1', 'EXG2', 'EXG3', 'EXG4']) +# Now let's take a closer look at the information in the epochs +# metadata. -# Furthermore, some datasets contain missing channels (stored in -# ``limo_epochs.info[‘bads’]``), which were dropped during preprocessing of the -# data. We’ll interpolate this channels for convenience. -limo_epochs.interpolate_bads(reset_bads=True) +# We want include all columns in the summary table +epochs_summary = limo_epochs.metadata.describe(include='all').round(3) +print(epochs_summary) ############################################################################### -# Now we can go ahead and plot the ERPs evoked by Face A and Face B +# The first column of the summary table above provides more or less the same +# information as the ``print(limo_epochs)`` command we ran before. There are +# 1055 faces (i.e., epochs), subdivided in 2 conditions (i.e., Face A and +# Face B) and, for this particular subject, there are more epochs for the +# condition Face B. +# +# In addition, we can see in the second column that the values for the +# phase-coherence variable range from -1.619 to 1.642. This is because the +# phase-coherence values are provided as a z-scored variable in the LIMO +# dataset. Note that they have a mean of zero and a standard deviation of 1. +# +# +# Visualize condition ERPs +# ------------------------ +# +# Let's plot the ERPs evoked by Face A and Face B, to see how similar they are. # only show -250 to 500 ms -ts_args = dict(xlim=(-.25, 0.5)) +ts_args = dict(xlim=(-0.25, 0.5)) -# plot evoked response for faces A & B -limo_epochs['Face/A'].average().plot_joint(times=[.15], +# plot evoked response for face A +limo_epochs['Face/A'].average().plot_joint(times=[0.15], title='Evoked response: Face A', ts_args=ts_args) - -limo_epochs['Face/B'].average().plot_joint(times=[.15], +# and face B +limo_epochs['Face/B'].average().plot_joint(times=[0.15], title='Evoked response: Face B', ts_args=ts_args) @@ -111,68 +158,162 @@ # We can also compute the difference wave contrasting Face A and Face B. # Although, looking at the evoked responses above, we shouldn't expect great # differences among these face-stimuli. -# -# Compute difference wave (Face A minus Face B) -difference_wave = mne.combine_evoked([limo_epochs['Face/A'].average(), - -limo_epochs['Face/B'].average()], - weights='equal') -# Plot difference between Face A and Face B -difference_wave.plot_joint(times=[.15], title='Difference Face A - Face B') +# Face A minus Face B +difference_wave = combine_evoked([limo_epochs['Face/A'].average(), + -limo_epochs['Face/B'].average()], + weights='equal') + +# plot difference wave +difference_wave.plot_joint(times=[0.15], title='Difference Face A - Face B') ############################################################################### -# As expected, no see clear differential patterns appears when contrasting -# Face A and Face B. However, we could narrow our search to -# since this is a "visual paradigm" it might be best to electrodes located over -# the occipital lobe. After all this is "visual paradigm". Thus, differences -# between stimuli (if any) might easier to spot over "more visual areas". -# +# As expected, no clear pattern appears when contrasting +# Face A and Face B. However, we could narrow our search a little bit more. +# Since this is a "visual paradigm" it might be best to look at electrodes +# located over the occipital lobe, as differences between stimuli (if any) +# might easier to spot over visual areas. + # Create a dictionary containing the evoked responses conditions = ["Face/A", "Face/B"] -evoked_dict = dict() -for condition in conditions: - evoked_dict[condition] = limo_epochs[condition].average() -print(evoked_dict) +evokeds = {condition: limo_epochs[condition].average() + for condition in conditions} -# concentrate on an occipital electrode -pick = evoked_dict["Face/A"].ch_names.index('B11') +# concentrate analysis an occipital electrodes (e.g. B11) +pick = evokeds["Face/A"].ch_names.index('B11') # compare evoked responses -mne.viz.plot_compare_evokeds(evoked_dict, picks=pick) +plot_compare_evokeds(evokeds, picks=pick, ylim=dict(eeg=(-15, 7.5))) ############################################################################### -# Next, we can inspect the effect of phase-coherence on the activation -# patterns evoked by the presented face-stimuli. -# Here, one would expect that faces with high phase-coherence evoke a stronger -# response, as participants should be better at identifying these faces. +# We do see a difference between Face A and B, but it is pretty small. # -# Create design matrix for linear regression. We'll use the information -# contained in the ``limo_epochs.metadata``. -design = limo_epochs.metadata.copy() -design = design.assign(intercept=1) # add intercept +# +# Visualize effect of stimulus phase-coherence +# -------------------------------------------- +# +# Since phase-coherence +# determined whether a face stimulus could be easily identified, +# one could expect that faces with high phase-coherence should evoke stronger +# activation patterns along occipital electrodes. + +phase_coh = limo_epochs.metadata['phase-coherence'] +# get levels of phase coherence +levels = sorted(phase_coh.unique()) +# create labels for levels of phase coherence (i.e., 0 - 85%) +labels = ["{0:.2f}".format(i) for i in np.arange(0., 0.90, 0.05)] + +# create dict of evokeds for each level of phase-coherence +evokeds = {label: limo_epochs[phase_coh == level].average() + for level, label in zip(levels, labels)} + +# pick channel to plot +electrodes = ['C22', 'B11'] +# create figures +for electrode in electrodes: + fig, ax = plt.subplots(figsize=(8, 4)) + plot_compare_evokeds(evokeds, + axes=ax, + ylim=dict(eeg=(-20, 15)), + picks=electrode, + cmap=("Phase coherence", "magma")) + +############################################################################### +# As shown above, there are some considerable differences between the +# activation patterns evoked by stimuli with low vs. high phase-coherence at +# the chosen electrodes. +# +# +# Prepare data for linear regression analysis +# -------------------------------------------- +# +# Before we test the significance of these differences using linear +# regression, we'll interpolate missing channels that were +# dropped during preprocessing of the data. +# Furthermore, we'll drop the EOG channels (marked by the "EXG" prefix) +# present in the data: + +limo_epochs.interpolate_bads(reset_bads=True) +limo_epochs.drop_channels(['EXG1', 'EXG2', 'EXG3', 'EXG4']) + +############################################################################### +# Define predictor variables and design matrix +# -------------------------------------------- +# +# To run the regression analysis, +# we need to create a design matrix containing information about the +# variables (i.e., predictors) we want to use for prediction of brain +# activity patterns. For this purpose, we'll use the information we have in +# ``limo_epochs.metadata``: phase-coherence and Face A vs. Face B. + +# name of predictors + intercept +predictor_vars = ['face a - face b', 'phase-coherence', 'intercept'] + +# create design matrix +design = limo_epochs.metadata[['phase-coherence', 'face']].copy() design['face a - face b'] = np.where(design['face'] == 'A', 1, -1) -names = ['intercept', 'face a - face b', 'phase-coherence'] +design['intercept'] = 1 +design = design[predictor_vars] + +############################################################################### +# Now we can set up the linear model to be used in the analysis using +# MNE-Python's func:`~mne.stats.linear_regression` function. -# fit linear model -reg = linear_regression(limo_epochs, design[names], names=names) +reg = linear_regression(limo_epochs, + design_matrix=design, + names=predictor_vars) ############################################################################### -# Visualise effect of phase-coherence. +# Extract regression coefficients +# ------------------------------- +# +# The results are stored within the object ``reg``, +# which is a dictionary of evoked objects containing +# multiple inferential measures for each predictor in the design matrix. + +print('predictors are:', list(reg)) +print('fields are:', [field for field in getattr(reg['intercept'], '_fields')]) + +############################################################################### +# Plot model results +# ------------------ +# +# Now we can access and plot the results of the linear regression analysis by +# calling :samp:`reg['{}'].{}` and +# using the +# :meth:`~mne.Evoked.plot_joint` method just as we would do with any other +# evoked object. +# Below we can see a clear effect of phase-coherence, with higher +# phase-coherence (i.e., better "face visibility") having a negative effect on +# the activity measured at occipital electrodes around 200 to 250 ms following +# stimulus onset. + reg['phase-coherence'].beta.plot_joint(ts_args=ts_args, - title='Effect of phase-coherence', - times=[.23]) + title='Effect of Phase-coherence', + times=[0.23]) ############################################################################### -# Here we can see a clear effect of phase-coherence, with higher -# phase-coherence (i.e., better "face visibility") being associated with -# stronger activity patterns. +# We can also plot the corresponding T values. + +# use unit=False and scale=1 to keep values at their original +# scale (i.e., avoid conversion to micro-volt). +ts_args = dict(xlim=(-0.25, 0.5), + unit=False) +topomap_args = dict(scalings=dict(eeg=1), + average=0.05) + +# sphinx_gallery_thumbnail_number = 9 +fig = reg['phase-coherence'].t_val.plot_joint(ts_args=ts_args, + topomap_args=topomap_args, + times=[0.23]) +fig.axes[0].set_ylabel('T-value') ############################################################################### # Conversely, there appears to be no (or very small) systematic effects when -# constraining Face A and Face B. This is largely consistent with the +# comparing Face A and Face B stimuli. This is largely consistent with the # difference wave approach presented above. -# -# Visualise effect of face condition (Face A vs. Face B). -reg['face a - face b'].beta.plot_joint(title='Contrast: Face A - Face B', - ts_args=ts_args, - times=[.15]) +ts_args = dict(xlim=(-0.25, 0.5)) + +reg['face a - face b'].beta.plot_joint(ts_args=ts_args, + title='Effect of Face A vs. Face B', + times=[0.23]) diff --git a/mne/datasets/__init__.py b/mne/datasets/__init__.py index 80a43996660..dacdfb8e1ff 100644 --- a/mne/datasets/__init__.py +++ b/mne/datasets/__init__.py @@ -31,4 +31,5 @@ 'fetch_hcp_mmp_parcellation', 'fieldtrip_cmc', 'hf_sef', 'kiloword', 'misc', 'mtrf', 'multimodal', 'opm', 'phantom_4dbti', 'sample', 'sleep_physionet', 'somato', 'spm_face', 'testing', 'visual_92_categories', + 'limo', ] diff --git a/mne/datasets/utils.py b/mne/datasets/utils.py index 01fa31eede8..de7947a8d5a 100644 --- a/mne/datasets/utils.py +++ b/mne/datasets/utils.py @@ -231,7 +231,8 @@ def _data_path(path=None, force_update=False, update_path=True, download=True, 'kiloword': 'MNE_DATASETS_KILOWORD_PATH', 'mtrf': 'MNE_DATASETS_MTRF_PATH', 'fieldtrip_cmc': 'MNE_DATASETS_FIELDTRIP_CMC_PATH', - 'phantom_4dbti': 'MNE_DATASETS_PHANTOM_4DBTI_PATH' + 'phantom_4dbti': 'MNE_DATASETS_PHANTOM_4DBTI_PATH', + 'limo': 'MNE_DATASETS_LIMO_PATH', }[name] path = _get_path(path, key, name) @@ -598,7 +599,7 @@ def _download_all_example_data(verbose=True): fetch_hcp_mmp_parcellation() finally: sys.argv.pop(-1) - limo.load_data(subject=2, update_path=True) + limo.load_data(subject=1, update_path=True) @verbose diff --git a/mne/utils/config.py b/mne/utils/config.py index 8d039314fbc..d1fe9c8c617 100644 --- a/mne/utils/config.py +++ b/mne/utils/config.py @@ -103,6 +103,7 @@ def set_memmap_min_size(memmap_min_size): 'MNE_DATASETS_KILOWORD_PATH', 'MNE_DATASETS_FIELDTRIP_CMC_PATH', 'MNE_DATASETS_PHANTOM_4DBTI_PATH', + 'MNE_DATASETS_LIMO_PATH', 'MNE_FORCE_SERIAL', 'MNE_KIT2FIFF_STIM_CHANNELS', 'MNE_KIT2FIFF_STIM_CHANNEL_CODING',