From a348c7958b0b95dbf9f73d447f034d5c79cc8204 Mon Sep 17 00:00:00 2001 From: "Jose C. G. Alanis" <12409129+JoseAlanis@users.noreply.github.com> Date: Sat, 12 Oct 2019 22:53:24 +0200 Subject: [PATCH 01/10] add info / refs on limo dataset add sklearn regression add class linear regression update limo doc fix limo fetcher docstring fix merging issues --- doc/overview/datasets_index.rst | 25 ++ doc/python_reference.rst | 1 + examples/datasets/plot_limo_data.py | 428 ++++++++++++++++++++++------ mne/datasets/__init__.py | 1 + mne/datasets/limo/limo.py | 3 +- mne/datasets/utils.py | 5 +- 6 files changed, 365 insertions(+), 98 deletions(-) diff --git a/doc/overview/datasets_index.rst b/doc/overview/datasets_index.rst index f04b5bfea6c..b9e544b7c4f 100644 --- a/doc/overview/datasets_index.rst +++ b/doc/overview/datasets_index.rst @@ -345,6 +345,30 @@ 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 +^^^^^^^^^^^^ +: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 in 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 +386,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/doc/python_reference.rst b/doc/python_reference.rst index dcf7db56638..c0764fd7494 100644 --- a/doc/python_reference.rst +++ b/doc/python_reference.rst @@ -202,6 +202,7 @@ Datasets spm_face.data_path visual_92_categories.data_path phantom_4dbti.data_path + limo.load_data Visualization diff --git a/examples/datasets/plot_limo_data.py b/examples/datasets/plot_limo_data.py index d5c53f86521..a3efc5a0454 100644 --- a/examples/datasets/plot_limo_data.py +++ b/examples/datasets/plot_limo_data.py @@ -1,28 +1,31 @@ """ -======================================================== -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`_. +Furthermore, the code 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 code allows to: -The code allows to: +Fetch epoched data from single subject files 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, the 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 and added to the epochs structure as +metadata. +In addition, the code shows how to to fit linear models on single subject data +and derive inferential measures to evaluate the significance of the estimated +effects using bootstrap and spatio-temporal clustering techniques. References ---------- @@ -33,76 +36,131 @@ 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) import numpy as np +import pandas as pd import matplotlib.pyplot as plt -import mne -from mne.datasets import limo +from sklearn.linear_model import LinearRegression + +from mne.datasets.limo import load_data +from mne.decoding import Vectorizer, get_coef +from mne.evoked import EvokedArray 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 ############################################################################### -# 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``. +# About the data +# -------------- +# +# In the original LIMO experiment (see [2]_), participants performed a +# two-alternative forced choice task, discriminating between two face stimuli. +# Subjects discriminated the same two faces during the whole experiment. +# However, the critical manipulation consisted in the level of noise +# added to the face-stimuli during the task, making the faces more or less +# discernible to the observer (see `Fig 1`_ in [3]_ for instance). +# +# 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 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) + +############################################################################### +# 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") +# ``limo_epochs`` structure. Events should appear clearly grouped, as they +# were ordered during the import process. + +fig, ax = plt.subplots(figsize=(7, 4)) +plot_events(limo_epochs.events, event_id=limo_epochs.event_id, axes=ax) +ax.set(title="Distribution of events in LIMO epochs") 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 +# 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 we can take a closer look at the information entailed 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 +# ------------------------ +# +# We can go ahead and plot the ERPs evoked by Face A and Face B # only show -250 to 500 ms ts_args = dict(xlim=(-.25, 0.5)) -# plot evoked response for faces A & B +# plot evoked response for face A limo_epochs['Face/A'].average().plot_joint(times=[.15], title='Evoked response: Face A', ts_args=ts_args) - +# and face B limo_epochs['Face/B'].average().plot_joint(times=[.15], title='Evoked response: Face B', ts_args=ts_args) @@ -111,68 +169,250 @@ # 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 +# 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=[.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 see clear differential 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) -# concentrate on an occipital electrode +# concentrate analysis an occipital electrodes (e.g. B11) pick = evoked_dict["Face/A"].ch_names.index('B11') # compare evoked responses -mne.viz.plot_compare_evokeds(evoked_dict, picks=pick) +plot_compare_evokeds(evoked_dict, picks=pick, ylim=dict(eeg=[-15, 5])) + +############################################################################### +# As expected, the difference between Face A and Face B are very small. ############################################################################### -# 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. +# Visualize effect of stimulus phase-coherence +# -------------------------------------------- # -# 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 +# For visualization purposes, we can transform the phase-coherence +# variable back to fit roughly to it's original scale. Phase-coherence +# determined whether a face stimulus could be identified as such. Thus, +# one could expect that faces with high phase-coherence should evoke stronger +# activation patterns along occipital electrodes. +# +# As phase-coherence variable is a continuous variable, we'll need to +# split the values of face coherence to percentiles to match roughly their +# original scale (i.e, 0 to 100 % signal-to-noise ratio in 10% steps). +# **Note:** We'll only do this for visualization purposes, + +name = "phase-coherence" +factor = 'factor-' + name + +# create phase-coherence percentiles +df = limo_epochs.metadata +df[factor] = pd.cut(df[name], 11, labels=False) / 10 +# overwrite metadata +limo_epochs.metadata = df + +# color scheme for percentile plot +colors = {str(val): val for val in np.sort(df[factor].unique())} +# compute evoked for each phase-coherence percentile +evokeds = {str(val): limo_epochs[limo_epochs.metadata[factor] == val].average() + for val in colors.values()} + +# 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, 20]), + colors=colors, + split_legend=True, + picks=electrode, + cmap=(name + " Percentile", "magma")) + fig.show() + +############################################################################### +# As shown above, there are some considerable differences between the +# activation patterns evoked by stimuli with low vs. high phase-coherence. +# +# In particular, differences appear to be most pronounced along fronto-central +# and occipital electrodes. + +############################################################################### +# Prepare data for linear regression analysis +# -------------------------------------------- +# +# Next, we can test the significance of these differences using linear +# regression.But before we go on, we'll start by interpolating any missing +# channels in the LIMO epochs structure. Some subjects in the datasets contain +# missing channels (stored in ``limo_epochs.info[‘bads’]``) as these were +# dropped during preprocessing of the data. + +limo_epochs.interpolate_bads(reset_bads=True) + +# Furthermore, we'll drop the EOG channels (marked by the "EXG" prefix) +# present in the LIMO epochs structure. +limo_epochs.drop_channels(['EXG1', 'EXG2', 'EXG3', 'EXG4']) + +############################################################################### +# Define predictor variables and design matrix +# -------------------------------------------- +# +# First, 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 contained in +# the ``limo_epochs.metadata``. Here, we'll explore the effect of +# phase-coherence as well as the effect of 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']] design['face a - face b'] = np.where(design['face'] == 'A', 1, -1) -names = ['intercept', 'face a - face b', 'phase-coherence'] +design = design.assign(intercept=1) # add intercept +design = design[predictor_vars] + +############################################################################### +# Now we can set up the linear model to be used in the analysis. +# For the purpose of demonstration, we'll use scikit-learn's LinearRegression +# (see :class:`sklearn.linear_model.LinearRegression`) in this example. +# +# Basically, we're creating a wrapper, ``STLinearRegression``, for +# sklearn's ``LinearRegression``, which computes a +# `least squares`_ solution for our data given the provided design matrix. The +# and results are stored within ``LinearRegression`` object for convenience. + + +class STLinearRegression(LinearRegression): + """ + Create linear model object. + + Notes + ----- + Currently, the input data has to have a shape of samples by channels by + time points. The data will be automatically vectorized for easier / faster + handling. Thus the vectorized data (Y) within STLinearRegression + has shape of samples by channels * time points. + """ + def __init__(self, predictors, design_matrix, data, weights=None, + fit_intercept=True, normalize=False, + n_jobs=None): + + # store model parameters + super().__init__(fit_intercept=fit_intercept, normalize=normalize, + n_jobs=n_jobs) + self.predictors = predictors + self.design = design_matrix + self.orig_shape = data.shape[1:] + self.Y = Vectorizer().fit_transform(data) + self.weights = weights + + # automatically fit the linear model + self.fit(X=self.design, y=self.Y, sample_weight=self.weights) + + # compute beta coefficients + def compute_beta_coefs(self, predictor, output='evoked'): + # extract coefficients from linear model estimator + beta_coefs = get_coef(self, 'coef_') + # select predictor in question + pred_col = predictor_vars.index(predictor) -# fit linear model -reg = linear_regression(limo_epochs, design[names], names=names) + if output == 'evoked': + # coefficients are projected back to a channels x time points + betas = beta_coefs[:, pred_col] + betas = betas.reshape(self.orig_shape) + # create evoked object containing the back projected coefficients + betas = EvokedArray(betas, + comment=predictor, + info=limo_epochs.info, + tmin=limo_epochs.tmin) + else: + # return raw values + betas = beta_coefs[:, pred_col] + + # return beta coefficients + return betas + + # compute model predictions + def compute_predictions(self): + # compute predicted values + predictions = self.predict(X=self.design) + # return beta predictions + return predictions + + +############################################################################### +# Set up the model +# ---------------- +# +# We already have an intercept column in the design matrix, +# thus we'll call STLinearRegression with fit_intercept=False + +linear_model = STLinearRegression(fit_intercept=False, + predictors=predictor_vars, + design_matrix=design, + data=limo_epochs.get_data()) ############################################################################### -# Visualise effect of phase-coherence. -reg['phase-coherence'].beta.plot_joint(ts_args=ts_args, - title='Effect of phase-coherence', - times=[.23]) +# Extract regression coefficients +# ------------------------------- +# +# As described above, the results stored within the object ``linear_model``. +# We can extract the coefficients (i.e., the betas) from the +# linear model estimator by calling ``linear_model.compute_beta_coefs()``. +# This will automatically create an evoked object that can used later + +pc_betas = linear_model.compute_beta_coefs(predictor='phase-coherence') +face_betas = linear_model.compute_beta_coefs(predictor='face a - face b') ############################################################################### -# 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. +# Plot model results +# ------------------ +# +# Now we can plot results of the linear regression analysis. +# 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. + +# only show -250 to 500 ms +ts_args = dict(xlim=(-.25, 0.5)) + +pc_betas.plot_joint(ts_args=ts_args, + title='Phase-coherence - sklearn betas', + times=[.23]) ############################################################################### # Conversely, there appears to be no (or very small) systematic effects when -# constraining Face A and Face B. This is largely consistent with the +# constraining 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]) + +face_betas.plot_joint(title='Face A - Face B - sklearn betas', + ts_args=ts_args, + times=[.15]) + +############################################################################### +# Finally we can compare the output to MNE's linear_regression function +# see func:`mne.stats.linear_regression`. + +mne_reg = linear_regression(limo_epochs, + design_matrix=design, + names=predictor_vars) + +mne_reg['phase-coherence'].beta.plot_joint(ts_args=ts_args, + title='Phase-coherence - MNE', + times=[.23]) diff --git a/mne/datasets/__init__.py b/mne/datasets/__init__.py index 80a43996660..86bc5796f5c 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/limo/limo.py b/mne/datasets/limo/limo.py index 8450c1d083b..eb01053be09 100644 --- a/mne/datasets/limo/limo.py +++ b/mne/datasets/limo/limo.py @@ -158,8 +158,7 @@ def load_data(subject, path=None, force_update=False, update_path=None, Returns ------- - epochs : instance of Epochs - The epochs. + epochs : Epochs """ # noqa: E501 pd = _check_pandas_installed() from scipy.io import loadmat diff --git a/mne/datasets/utils.py b/mne/datasets/utils.py index 01fa31eede8..d3ade7ef09c 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 From 387cd30e85e88152503d495f8ceabde96d2c8cb3 Mon Sep 17 00:00:00 2001 From: Eric Larson Date: Mon, 14 Oct 2019 09:29:25 -0400 Subject: [PATCH 02/10] FIX: Indentation --- examples/datasets/plot_limo_data.py | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/examples/datasets/plot_limo_data.py b/examples/datasets/plot_limo_data.py index a3efc5a0454..237f8533860 100644 --- a/examples/datasets/plot_limo_data.py +++ b/examples/datasets/plot_limo_data.py @@ -13,19 +13,19 @@ In summary, the code allows to: -Fetch epoched data from single subject files 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. - -During import, the 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 and added to the epochs structure as -metadata. - -In addition, the code shows how to to fit linear models on single subject data -and derive inferential measures to evaluate the significance of the estimated -effects using bootstrap and spatio-temporal clustering techniques. +- Fetch epoched data from single subject files 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. + +- During import, the 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 and added to the epochs + structure as metadata. + +- In addition, the code shows how to to fit linear models on single subject + data and derive inferential measures to evaluate the significance of the + estimated effects using bootstrap and spatio-temporal clustering techniques. References ---------- From 9ea120bf958e90ceabe057c617d07b4ad287ed93 Mon Sep 17 00:00:00 2001 From: Eric Larson Date: Mon, 14 Oct 2019 09:31:29 -0400 Subject: [PATCH 03/10] FIX: A couple more --- .circleci/config.yml | 2 +- mne/datasets/utils.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) 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/mne/datasets/utils.py b/mne/datasets/utils.py index d3ade7ef09c..de7947a8d5a 100644 --- a/mne/datasets/utils.py +++ b/mne/datasets/utils.py @@ -232,7 +232,7 @@ def _data_path(path=None, force_update=False, update_path=True, download=True, 'mtrf': 'MNE_DATASETS_MTRF_PATH', 'fieldtrip_cmc': 'MNE_DATASETS_FIELDTRIP_CMC_PATH', 'phantom_4dbti': 'MNE_DATASETS_PHANTOM_4DBTI_PATH', - 'limo': 'MNE_DATASETS_LIMO_PATH' + 'limo': 'MNE_DATASETS_LIMO_PATH', }[name] path = _get_path(path, key, name) From 6cb374eb50bbfd16c225852ffcd6cdc5ad979b02 Mon Sep 17 00:00:00 2001 From: Eric Larson Date: Mon, 14 Oct 2019 09:36:15 -0400 Subject: [PATCH 04/10] STY: Nitpicks [ci skip] --- doc/python_reference.rst | 1 - mne/datasets/__init__.py | 2 +- mne/datasets/limo/limo.py | 3 ++- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/doc/python_reference.rst b/doc/python_reference.rst index c0764fd7494..dcf7db56638 100644 --- a/doc/python_reference.rst +++ b/doc/python_reference.rst @@ -202,7 +202,6 @@ Datasets spm_face.data_path visual_92_categories.data_path phantom_4dbti.data_path - limo.load_data Visualization diff --git a/mne/datasets/__init__.py b/mne/datasets/__init__.py index 86bc5796f5c..dacdfb8e1ff 100644 --- a/mne/datasets/__init__.py +++ b/mne/datasets/__init__.py @@ -31,5 +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' + 'limo', ] diff --git a/mne/datasets/limo/limo.py b/mne/datasets/limo/limo.py index eb01053be09..8450c1d083b 100644 --- a/mne/datasets/limo/limo.py +++ b/mne/datasets/limo/limo.py @@ -158,7 +158,8 @@ def load_data(subject, path=None, force_update=False, update_path=None, Returns ------- - epochs : Epochs + epochs : instance of Epochs + The epochs. """ # noqa: E501 pd = _check_pandas_installed() from scipy.io import loadmat From 6d6d86d7de47d018e088bc0c3a5b7c1da3de82be Mon Sep 17 00:00:00 2001 From: Eric Larson Date: Mon, 14 Oct 2019 09:49:38 -0400 Subject: [PATCH 05/10] ENH: Links --- doc/changes/latest.inc | 2 ++ doc/changes/names.inc | 2 ++ doc/overview/datasets_index.rst | 2 ++ 3 files changed, 6 insertions(+) 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 b9e544b7c4f..966ca8fe58a 100644 --- a/doc/overview/datasets_index.rst +++ b/doc/overview/datasets_index.rst @@ -345,6 +345,8 @@ 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`. From 80de925a4f80572d79e63f825eb8867dc1aa5e00 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jos=C3=A9=20C=2E=20Garc=C3=ADa=20Alanis?= <12409129+JoseAlanis@users.noreply.github.com> Date: Tue, 22 Oct 2019 11:25:43 +0200 Subject: [PATCH 06/10] Apply suggestions from code review add jonas suggestions Co-Authored-By: jona-sassenhagen --- doc/overview/datasets_index.rst | 2 +- examples/datasets/plot_limo_data.py | 10 +++++----- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/doc/overview/datasets_index.rst b/doc/overview/datasets_index.rst index 966ca8fe58a..3f9bd56d2b0 100644 --- a/doc/overview/datasets_index.rst +++ b/doc/overview/datasets_index.rst @@ -354,7 +354,7 @@ LIMO Dataset 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 in the level of noise added to the +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. diff --git a/examples/datasets/plot_limo_data.py b/examples/datasets/plot_limo_data.py index 237f8533860..774f66f694c 100644 --- a/examples/datasets/plot_limo_data.py +++ b/examples/datasets/plot_limo_data.py @@ -195,7 +195,7 @@ pick = evoked_dict["Face/A"].ch_names.index('B11') # compare evoked responses -plot_compare_evokeds(evoked_dict, picks=pick, ylim=dict(eeg=[-15, 5])) +plot_compare_evokeds(evoked_dict, picks=pick, ylim=dict(eeg=(-15, 5))) ############################################################################### # As expected, the difference between Face A and Face B are very small. @@ -216,7 +216,7 @@ # **Note:** We'll only do this for visualization purposes, name = "phase-coherence" -factor = 'factor-' + name +factor = 'factor_' + name # create phase-coherence percentiles df = limo_epochs.metadata @@ -227,7 +227,7 @@ # color scheme for percentile plot colors = {str(val): val for val in np.sort(df[factor].unique())} # compute evoked for each phase-coherence percentile -evokeds = {str(val): limo_epochs[limo_epochs.metadata[factor] == val].average() +evokeds = {str(val): limo_epochs[factor + ' == ' str(val)].average() for val in colors.values()} # pick channel to plot @@ -237,7 +237,7 @@ fig, ax = plt.subplots(figsize=(8, 4)) plot_compare_evokeds(evokeds, axes=ax, - ylim=dict(eeg=[-20, 20]), + ylim=dict(eeg=(-20, 20)), colors=colors, split_legend=True, picks=electrode, @@ -283,7 +283,7 @@ # create design matrix design = limo_epochs.metadata[['phase-coherence', 'face']] design['face a - face b'] = np.where(design['face'] == 'A', 1, -1) -design = design.assign(intercept=1) # add intercept +design['intercept'] = 1 design = design[predictor_vars] ############################################################################### From ec85b913391fc51874fc1498a47dc3d9e30842d6 Mon Sep 17 00:00:00 2001 From: "Jose C. G. Alanis" <12409129+JoseAlanis@users.noreply.github.com> Date: Tue, 22 Oct 2019 23:38:31 +0200 Subject: [PATCH 07/10] fix syntax + remaining issues --- examples/datasets/plot_limo_data.py | 24 +++++++++++------------- 1 file changed, 11 insertions(+), 13 deletions(-) diff --git a/examples/datasets/plot_limo_data.py b/examples/datasets/plot_limo_data.py index 774f66f694c..992c04c0194 100644 --- a/examples/datasets/plot_limo_data.py +++ b/examples/datasets/plot_limo_data.py @@ -187,15 +187,14 @@ # 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() +evokeds = {condition: limo_epochs[condition].average() + for condition in conditions} # concentrate analysis an occipital electrodes (e.g. B11) -pick = evoked_dict["Face/A"].ch_names.index('B11') +pick = evokeds["Face/A"].ch_names.index('B11') # compare evoked responses -plot_compare_evokeds(evoked_dict, picks=pick, ylim=dict(eeg=(-15, 5))) +plot_compare_evokeds(evokeds, picks=pick, ylim=dict(eeg=(-15, 5))) ############################################################################### # As expected, the difference between Face A and Face B are very small. @@ -218,16 +217,15 @@ name = "phase-coherence" factor = 'factor_' + name -# create phase-coherence percentiles -df = limo_epochs.metadata -df[factor] = pd.cut(df[name], 11, labels=False) / 10 -# overwrite metadata -limo_epochs.metadata = df +# color scheme for percentile plot +limo_epochs.metadata[factor] = pd.cut(limo_epochs.metadata[name], 11, + labels=False) / 10 # color scheme for percentile plot -colors = {str(val): val for val in np.sort(df[factor].unique())} +colors = {str(val): val + for val in np.sort(limo_epochs.metadata[factor].unique())} # compute evoked for each phase-coherence percentile -evokeds = {str(val): limo_epochs[factor + ' == ' str(val)].average() +evokeds = {str(val): limo_epochs[limo_epochs.metadata[factor] == val].average() for val in colors.values()} # pick channel to plot @@ -351,7 +349,7 @@ def compute_beta_coefs(self, predictor, output='evoked'): def compute_predictions(self): # compute predicted values predictions = self.predict(X=self.design) - # return beta predictions + # return beta coefficients return predictions From ba696f7dc550efd7e3d98155c5445f43c701c61c Mon Sep 17 00:00:00 2001 From: Daniel McCloy Date: Fri, 25 Oct 2019 03:45:06 -0700 Subject: [PATCH 08/10] tighten up prose and code (#1) --- examples/datasets/plot_limo_data.py | 141 ++++++++++------------------ 1 file changed, 52 insertions(+), 89 deletions(-) diff --git a/examples/datasets/plot_limo_data.py b/examples/datasets/plot_limo_data.py index 992c04c0194..3b4877b4548 100644 --- a/examples/datasets/plot_limo_data.py +++ b/examples/datasets/plot_limo_data.py @@ -7,25 +7,25 @@ Here we explore the structure of the data contained in the `LIMO dataset`_. -Furthermore, the code replicates and extends some of the main analysis +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_. -In summary, the code allows to: +In summary, the example: -- Fetch epoched data from single subject files of the LIMO dataset [1]_. +- 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. -- During import, the 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 and added to the epochs - structure as metadata. +- 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. -- In addition, the code shows how to to fit linear models on single subject - data and derive inferential measures to evaluate the significance of the - estimated effects using bootstrap and spatio-temporal clustering techniques. +- Fits linear models on the single subject's data and derives inferential + measures to evaluate the significance of the estimated effects using + bootstrapping and spatio-temporal clustering techniques. References ---------- @@ -51,7 +51,6 @@ # License: BSD (3-clause) import numpy as np -import pandas as pd import matplotlib.pyplot as plt from sklearn.linear_model import LinearRegression @@ -75,9 +74,8 @@ # # In the original LIMO experiment (see [2]_), participants performed a # two-alternative forced choice task, discriminating between two face stimuli. -# Subjects discriminated the same two faces during the whole experiment. -# However, the critical manipulation consisted in the level of noise -# added to the face-stimuli during the task, making the faces more or less +# 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 instance). # # The presented faces varied across a noise-signal (or phase-coherence) @@ -85,8 +83,8 @@ # 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 very hard to discriminate. - -############################################################################### +# +# # Load the data # ------------- # @@ -108,15 +106,11 @@ # ---------------- # # We can visualise the distribution of the face events contained in the -# ``limo_epochs`` structure. Events should appear clearly grouped, as they -# were ordered during the import process. +# ``limo_epochs`` structure. Events should appear clearly grouped, as the +# epochs are ordered by condition. -fig, ax = plt.subplots(figsize=(7, 4)) -plot_events(limo_epochs.events, event_id=limo_epochs.event_id, axes=ax) -ax.set(title="Distribution of events in LIMO epochs") -plt.legend(loc='lower left', borderaxespad=1.) -plt.tight_layout() -plt.show() +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``. @@ -146,22 +140,22 @@ # 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 # ------------------------ # # We can go ahead and plot the ERPs evoked by Face A and Face B # 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 face A -limo_epochs['Face/A'].average().plot_joint(times=[.15], +limo_epochs['Face/A'].average().plot_joint(times=[0.15], title='Evoked response: Face A', ts_args=ts_args) # and face B -limo_epochs['Face/B'].average().plot_joint(times=[.15], +limo_epochs['Face/B'].average().plot_joint(times=[0.15], title='Evoked response: Face B', ts_args=ts_args) @@ -176,10 +170,10 @@ weights='equal') # plot difference wave -difference_wave.plot_joint(times=[.15], title='Difference Face A - Face B') +difference_wave.plot_joint(times=[0.15], title='Difference Face A - Face B') ############################################################################### -# As expected, no see clear differential pattern appears when contrasting +# As expected, no see 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) @@ -198,35 +192,20 @@ ############################################################################### # As expected, the difference between Face A and Face B are very small. - -############################################################################### +# +# # Visualize effect of stimulus phase-coherence # -------------------------------------------- # -# For visualization purposes, we can transform the phase-coherence -# variable back to fit roughly to it's original scale. Phase-coherence -# determined whether a face stimulus could be identified as such. Thus, +# 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. -# -# As phase-coherence variable is a continuous variable, we'll need to -# split the values of face coherence to percentiles to match roughly their -# original scale (i.e, 0 to 100 % signal-to-noise ratio in 10% steps). -# **Note:** We'll only do this for visualization purposes, -name = "phase-coherence" -factor = 'factor_' + name - -# color scheme for percentile plot -limo_epochs.metadata[factor] = pd.cut(limo_epochs.metadata[name], 11, - labels=False) / 10 - -# color scheme for percentile plot -colors = {str(val): val - for val in np.sort(limo_epochs.metadata[factor].unique())} -# compute evoked for each phase-coherence percentile -evokeds = {str(val): limo_epochs[limo_epochs.metadata[factor] == val].average() - for val in colors.values()} +phase_coh = limo_epochs.metadata['phase-coherence'] +levels = sorted(phase_coh.unique()) +evokeds = {str(round(level, 2)): limo_epochs[phase_coh == level].average() + for level in levels} # pick channel to plot electrodes = ['C22', 'B11'] @@ -236,63 +215,52 @@ plot_compare_evokeds(evokeds, axes=ax, ylim=dict(eeg=(-20, 20)), - colors=colors, - split_legend=True, picks=electrode, - cmap=(name + " Percentile", "magma")) - fig.show() + cmap=("Standardized phase coherence", "magma")) ############################################################################### # As shown above, there are some considerable differences between the -# activation patterns evoked by stimuli with low vs. high phase-coherence. +# activation patterns evoked by stimuli with low vs. high phase-coherence at +# the chosen electrodes. +# # -# In particular, differences appear to be most pronounced along fronto-central -# and occipital electrodes. - -############################################################################### # Prepare data for linear regression analysis # -------------------------------------------- # -# Next, we can test the significance of these differences using linear -# regression.But before we go on, we'll start by interpolating any missing -# channels in the LIMO epochs structure. Some subjects in the datasets contain -# missing channels (stored in ``limo_epochs.info[‘bads’]``) as these were +# 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) -# Furthermore, we'll drop the EOG channels (marked by the "EXG" prefix) -# present in the LIMO epochs structure. limo_epochs.drop_channels(['EXG1', 'EXG2', 'EXG3', 'EXG4']) ############################################################################### # Define predictor variables and design matrix # -------------------------------------------- # -# First, we need to create a design matrix containing information about the +# 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 contained in -# the ``limo_epochs.metadata``. Here, we'll explore the effect of -# phase-coherence as well as the effect of Face A vs. Face B. +# 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']] +design = limo_epochs.metadata[['phase-coherence', 'face']].copy() design['face a - face b'] = np.where(design['face'] == 'A', 1, -1) design['intercept'] = 1 design = design[predictor_vars] ############################################################################### # Now we can set up the linear model to be used in the analysis. -# For the purpose of demonstration, we'll use scikit-learn's LinearRegression -# (see :class:`sklearn.linear_model.LinearRegression`) in this example. -# -# Basically, we're creating a wrapper, ``STLinearRegression``, for -# sklearn's ``LinearRegression``, which computes a -# `least squares`_ solution for our data given the provided design matrix. The -# and results are stored within ``LinearRegression`` object for convenience. +# We'll use :class:`sklearn.linear_model.LinearRegression`, and create a +# wrapper ``STLinearRegression`` that computes a +# `least squares`_ solution for our data given the provided design matrix. class STLinearRegression(LinearRegression): @@ -341,8 +309,6 @@ def compute_beta_coefs(self, predictor, output='evoked'): else: # return raw values betas = beta_coefs[:, pred_col] - - # return beta coefficients return betas # compute model predictions @@ -387,12 +353,9 @@ def compute_predictions(self): # the activity measured at occipital electrodes around 200 to 250 ms following # stimulus onset. -# only show -250 to 500 ms -ts_args = dict(xlim=(-.25, 0.5)) - pc_betas.plot_joint(ts_args=ts_args, title='Phase-coherence - sklearn betas', - times=[.23]) + times=[0.23]) ############################################################################### # Conversely, there appears to be no (or very small) systematic effects when @@ -401,7 +364,7 @@ def compute_predictions(self): face_betas.plot_joint(title='Face A - Face B - sklearn betas', ts_args=ts_args, - times=[.15]) + times=[0.15]) ############################################################################### # Finally we can compare the output to MNE's linear_regression function @@ -413,4 +376,4 @@ def compute_predictions(self): mne_reg['phase-coherence'].beta.plot_joint(ts_args=ts_args, title='Phase-coherence - MNE', - times=[.23]) + times=[0.23]) From f1553fa5192dd5ea81fdcaf573897bfd000de685 Mon Sep 17 00:00:00 2001 From: "Jose C. Garcia Alanis" <12409129+JoseAlanis@users.noreply.github.com> Date: Fri, 25 Oct 2019 14:34:40 +0200 Subject: [PATCH 09/10] drop class approach / add LIMO_PATH to config --- examples/datasets/plot_limo_data.py | 163 +++++++++------------------- mne/utils/config.py | 1 + 2 files changed, 51 insertions(+), 113 deletions(-) diff --git a/examples/datasets/plot_limo_data.py b/examples/datasets/plot_limo_data.py index 3b4877b4548..d0a2f7d8768 100644 --- a/examples/datasets/plot_limo_data.py +++ b/examples/datasets/plot_limo_data.py @@ -23,9 +23,8 @@ 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 derives inferential - measures to evaluate the significance of the estimated effects using - bootstrapping and spatio-temporal clustering techniques. +- Fits linear models on the single subject's data and visualizes inferential + measures to evaluate the significance of the estimated effects. References ---------- @@ -53,11 +52,7 @@ import numpy as np import matplotlib.pyplot as plt -from sklearn.linear_model import LinearRegression - from mne.datasets.limo import load_data -from mne.decoding import Vectorizer, get_coef -from mne.evoked import EvokedArray from mne.stats import linear_regression from mne.viz import plot_events, plot_compare_evokeds from mne import combine_evoked @@ -76,12 +71,12 @@ # 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 instance). +# 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 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 +# 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. # # @@ -188,7 +183,7 @@ pick = evokeds["Face/A"].ch_names.index('B11') # compare evoked responses -plot_compare_evokeds(evokeds, picks=pick, ylim=dict(eeg=(-15, 5))) +plot_compare_evokeds(evokeds, picks=pick, ylim=dict(eeg=(-15, 7.5))) ############################################################################### # As expected, the difference between Face A and Face B are very small. @@ -203,9 +198,14 @@ # activation patterns along occipital electrodes. phase_coh = limo_epochs.metadata['phase-coherence'] +# get levels of phase coherence levels = sorted(phase_coh.unique()) -evokeds = {str(round(level, 2)): limo_epochs[phase_coh == level].average() - for level in levels} +# create labels for levels of phase coherence (i.e., 0 - 85%) +labels = ["{0:.2f}".format(i) for i in np.arange(0., .90, .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'] @@ -214,9 +214,9 @@ fig, ax = plt.subplots(figsize=(8, 4)) plot_compare_evokeds(evokeds, axes=ax, - ylim=dict(eeg=(-20, 20)), + ylim=dict(eeg=(-15, 15)), picks=electrode, - cmap=("Standardized phase coherence", "magma")) + cmap=("Phase coherence", "magma")) ############################################################################### # As shown above, there are some considerable differences between the @@ -257,123 +257,60 @@ design = design[predictor_vars] ############################################################################### -# Now we can set up the linear model to be used in the analysis. -# We'll use :class:`sklearn.linear_model.LinearRegression`, and create a -# wrapper ``STLinearRegression`` that computes a -# `least squares`_ solution for our data given the provided design matrix. - - -class STLinearRegression(LinearRegression): - """ - Create linear model object. - - Notes - ----- - Currently, the input data has to have a shape of samples by channels by - time points. The data will be automatically vectorized for easier / faster - handling. Thus the vectorized data (Y) within STLinearRegression - has shape of samples by channels * time points. - """ - def __init__(self, predictors, design_matrix, data, weights=None, - fit_intercept=True, normalize=False, - n_jobs=None): - - # store model parameters - super().__init__(fit_intercept=fit_intercept, normalize=normalize, - n_jobs=n_jobs) - self.predictors = predictors - self.design = design_matrix - self.orig_shape = data.shape[1:] - self.Y = Vectorizer().fit_transform(data) - self.weights = weights - - # automatically fit the linear model - self.fit(X=self.design, y=self.Y, sample_weight=self.weights) - - # compute beta coefficients - def compute_beta_coefs(self, predictor, output='evoked'): - # extract coefficients from linear model estimator - beta_coefs = get_coef(self, 'coef_') - # select predictor in question - pred_col = predictor_vars.index(predictor) - - if output == 'evoked': - # coefficients are projected back to a channels x time points - betas = beta_coefs[:, pred_col] - betas = betas.reshape(self.orig_shape) - # create evoked object containing the back projected coefficients - betas = EvokedArray(betas, - comment=predictor, - info=limo_epochs.info, - tmin=limo_epochs.tmin) - else: - # return raw values - betas = beta_coefs[:, pred_col] - return betas - - # compute model predictions - def compute_predictions(self): - # compute predicted values - predictions = self.predict(X=self.design) - # return beta coefficients - return predictions - - -############################################################################### -# Set up the model -# ---------------- -# -# We already have an intercept column in the design matrix, -# thus we'll call STLinearRegression with fit_intercept=False +# Now we can set up the linear model to be used in the analysis using +# MNE's linear_regression function (see func:`mne.stats.linear_regression`). -linear_model = STLinearRegression(fit_intercept=False, - predictors=predictor_vars, - design_matrix=design, - data=limo_epochs.get_data()) +reg = linear_regression(limo_epochs, + design_matrix=design, + names=predictor_vars) ############################################################################### # Extract regression coefficients # ------------------------------- # -# As described above, the results stored within the object ``linear_model``. -# We can extract the coefficients (i.e., the betas) from the -# linear model estimator by calling ``linear_model.compute_beta_coefs()``. -# This will automatically create an evoked object that can used later +# The results stored within the object ``reg``. +# It basically consists of a dictionary of evoked objects containing +# multiple inferential measures for each predictor in the design matrix. -pc_betas = linear_model.compute_beta_coefs(predictor='phase-coherence') -face_betas = linear_model.compute_beta_coefs(predictor='face a - face b') +print('predictors are:', [key for key in reg.keys()]) +print('fields are:', [field for field in getattr(reg['intercept'], '_fields')]) ############################################################################### # Plot model results # ------------------ # -# Now we can plot results of the linear regression analysis. +# Now we can access and plot the results of the linear regression analysis by +# calling `reg[''].` and using the +# `.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. -pc_betas.plot_joint(ts_args=ts_args, - title='Phase-coherence - sklearn betas', - times=[0.23]) +reg['phase-coherence'].beta.plot_joint(ts_args=ts_args, + title='Effect of Phase-coherence', + times=[0.23]) ############################################################################### -# Conversely, there appears to be no (or very small) systematic effects when -# constraining Face A and Face B stimuli. This is largely consistent with the -# difference wave approach presented above. +# We can also plot the corresponding T values. -face_betas.plot_joint(title='Face A - Face B - sklearn betas', - ts_args=ts_args, - times=[0.15]) +# use unit=False and scale=1 to avoid conversion 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) -############################################################################### -# Finally we can compare the output to MNE's linear_regression function -# see func:`mne.stats.linear_regression`. +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') -mne_reg = linear_regression(limo_epochs, - design_matrix=design, - names=predictor_vars) +############################################################################### +# Conversely, there appears to be no (or very small) systematic effects when +# constraining Face A and Face B stimuli. This is largely consistent with the +# difference wave approach presented above. -mne_reg['phase-coherence'].beta.plot_joint(ts_args=ts_args, - title='Phase-coherence - MNE', - times=[0.23]) +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/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', From 6ee563083ed9efa5da6f15baadff08c21bef0873 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jos=C3=A9=20C=2E=20Garc=C3=ADa=20Alanis?= <12409129+JoseAlanis@users.noreply.github.com> Date: Fri, 25 Oct 2019 20:36:01 +0200 Subject: [PATCH 10/10] Apply suggestions from code review Co-Authored-By: Daniel McCloy --- examples/datasets/plot_limo_data.py | 35 ++++++++++++++++------------- 1 file changed, 19 insertions(+), 16 deletions(-) diff --git a/examples/datasets/plot_limo_data.py b/examples/datasets/plot_limo_data.py index d0a2f7d8768..038b140fab3 100644 --- a/examples/datasets/plot_limo_data.py +++ b/examples/datasets/plot_limo_data.py @@ -117,7 +117,7 @@ print(limo_epochs.metadata.head()) ############################################################################### -# Now we can take a closer look at the information entailed in the epochs +# Now let's take a closer look at the information in the epochs # metadata. # We want include all columns in the summary table @@ -126,7 +126,7 @@ ############################################################################### # 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 +# 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. @@ -140,7 +140,7 @@ # Visualize condition ERPs # ------------------------ # -# We can go ahead and plot the ERPs evoked by Face A and Face B +# 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=(-0.25, 0.5)) @@ -168,7 +168,7 @@ difference_wave.plot_joint(times=[0.15], title='Difference Face A - Face B') ############################################################################### -# As expected, no see clear pattern appears when contrasting +# 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) @@ -186,7 +186,7 @@ plot_compare_evokeds(evokeds, picks=pick, ylim=dict(eeg=(-15, 7.5))) ############################################################################### -# As expected, the difference between Face A and Face B are very small. +# We do see a difference between Face A and B, but it is pretty small. # # # Visualize effect of stimulus phase-coherence @@ -201,7 +201,7 @@ # 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., .90, .05)] +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() @@ -214,7 +214,7 @@ fig, ax = plt.subplots(figsize=(8, 4)) plot_compare_evokeds(evokeds, axes=ax, - ylim=dict(eeg=(-15, 15)), + ylim=dict(eeg=(-20, 15)), picks=electrode, cmap=("Phase coherence", "magma")) @@ -234,7 +234,6 @@ # present in the data: limo_epochs.interpolate_bads(reset_bads=True) - limo_epochs.drop_channels(['EXG1', 'EXG2', 'EXG3', 'EXG4']) ############################################################################### @@ -258,7 +257,7 @@ ############################################################################### # Now we can set up the linear model to be used in the analysis using -# MNE's linear_regression function (see func:`mne.stats.linear_regression`). +# MNE-Python's func:`~mne.stats.linear_regression` function. reg = linear_regression(limo_epochs, design_matrix=design, @@ -268,11 +267,11 @@ # Extract regression coefficients # ------------------------------- # -# The results stored within the object ``reg``. -# It basically consists of a dictionary of evoked objects containing +# 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:', [key for key in reg.keys()]) +print('predictors are:', list(reg)) print('fields are:', [field for field in getattr(reg['intercept'], '_fields')]) ############################################################################### @@ -280,8 +279,10 @@ # ------------------ # # Now we can access and plot the results of the linear regression analysis by -# calling `reg[''].` and using the -# `.plot_joint()` method just as we would do with any other evoked object. +# 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 @@ -294,13 +295,14 @@ ############################################################################### # We can also plot the corresponding T values. -# use unit=False and scale=1 to avoid conversion keep values at their original +# 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]) @@ -308,8 +310,9 @@ ############################################################################### # Conversely, there appears to be no (or very small) systematic effects when -# constraining Face A and Face B stimuli. This is largely consistent with the +# comparing Face A and Face B stimuli. This is largely consistent with the # difference wave approach presented above. +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',