Eelbrain

https://zenodo.org/badge/3651023.svg https://img.shields.io/conda/vn/conda-forge/eelbrain.svg https://img.shields.io/conda/pn/conda-forge/eelbrain.svg

Eelbrain is an open-source Python package for accessible statistical analysis of MEG and EEG data. It is maintained by Christian Brodbeck at the Computational sensorimotor systems lab at University of Maryland, College Park.

If you use Eelbrain in work that is published, please acknowledge it by citing it with the appropriate version and DOI.

Manual

Installing

Note

Because of the fluid nature of Python development, the recommended way of installing Eelbrain changes occasionally. For up-to-date information, see the corresponding Eelbrain wiki page.

Getting Started

Documentation

For an introduction to Eelbrain, see Introduction and the other Examples. For details on each functionality see the API Reference.

MacOS: Framework Build

On macOS, the GUI tool Eelbrain uses requires a special build of Python called a “Framework build”. You might see this error when trying to create a plot:

SystemExit: This program needs access to the screen.
Please run with a Framework build of python, and only when you are
logged in on the main display of your Mac.

In order to avoid this, Eelbrain installs a shortcut to start IPython with a Framework build:

$ eelbrain

This automatically launches IPython with the “eelbrain” profile. A default startup script that executes from eelbrain import * is created, and can be changed in the corresponding IPython profile.

Quitting iPython

Sometimes iPython seems to get stuck after this line:

Do you really want to exit ([y]/n)? y

In those instances, pressing ctrl-c usually terminates iPython immediately.

Windows: Scrolling

Scrolling inside a plot axes normally uses arrow keys, but this is currently not possible on Windows (due to an issue in Matplotlib). Instead, the following keys can be used:

↑ i

← j

→ l

↓ k

Changes

New in 0.34

  • API:

  • plot.Correlation renamed to plot.Scatter with some parameter changes for improved functionality.

  • New:

New in 0.33

  • API load.fiff.events(): The merge parameter is now by default inferred based on the raw data.

  • Boosting: plot data partitioning scheme (BoostingResult.splits.plot()).

  • Experiment pipeline:

New in 0.32

  • Requires at least Python 3.7

  • API changes:

    • Consistent class names for tests in test, testnd and pipeline.

    • plot.Timeplot argument order: second and third argument switched to facilitate plotting single category.

    • Significance markers for trends (.05 < p ≤ .1) are disabled by default.

  • boosting():

    • When using a basis, the function now considers the effect of the basis when normalizing predictors. This change leads to slightly different results, so TRFs should not be compared between this and previous versions. To facilitate keeping track of this change, it corresponds to a change in the BoostingResult.algorithm_version attribute from -1 to 0.

    • Different tstart/tstop for different predictors (contributed by Joshua Kulasingham)

    • Cross-validation of model fit (test parameter)

  • plot.Style to control advanced plotting options by category (see Boxplot example).

  • New functions/methods:

  • Experiment pipeline:

    • Methods with decim parameter now also have samplingrate parameter

    • More control over Events

New in 0.31

New in 0.30

New in 0.29

New in 0.28

  • Transition to Python 3.6

  • API changes:

    • testnd.ANOVA: The match parameter is now determined automatically and does not need to be specified anymore in most cases.

    • testnd.TTestOneSample.diff renamed to testnd.TTestOneSample.difference.

    • plot.Histogram: following matplotlib, the normed parameter was renamed to density.

    • Previously capitalized argument and attribute names Y, X and Xax are now lowercase.

    • Topomap-plot argument order changed to provide consistency between different plots.

  • NDVar and Var support for round(x)

  • MneExperiment pipeline:

    • Independent measures t-test

New in 0.27

New in 0.26

New in 0.25

  • Installation with conda (see Installing) and $ eelbrain launch script (see Getting Started).

  • API:

    • NDVar objects now inherit names through operations.

    • Assignment to a Dataset overwrites variable .name attributes, unless the Dataset key is a pythonified version of .name.

  • GUI/plotting:

    • When using iPython 5 or later, GUI start and stop is now automatic. It is possible to revert to the old behavior with plot.configure().

    • There are new hotkeys for most plots (see the individual plots’ help for details).

    • Plots automatically rescale when the window is resized.

  • MneExperiment:

    • A new MneExperiment.sessions attribute replaces defaults['experiment'], with support for multiple sessions in one experiment (see Setting up the file structure).

    • The MneExperiment.epochs parameter sel_epoch has been removed, use base instead.

    • The setting raw='clm' has been renamed to raw='raw'.

    • Custom preprocessing pipelines (see MneExperiment.raw).

    • The model parameter for ANOVA tests is now optional (see MneExperiment.tests).

  • Reverse correlation using boosting().

  • Loading and saving *.wav files (load.wav() and save.wav()).

New in 0.24

  • API:

    • MneExperiment: For data from the NYU New York system converted with mne < 0.13, the MneExperiment.meg_system attribute needs to be set to "KIT-157" to distinguish it from data collected with the KIT UMD system.

    • masked_parameter_map() method of cluster-based test results: use of pmin=None is deprecated, use pmin=1 instead.

  • New test: test.TTestRelated.

  • MneExperiment.make_report_rois() includes corrected p-values in reports for tests in more than one ROI

  • MneExperiment.make_rej() now has a decim parameter to improve display performance.

  • MneExperiment: BEM-solution files are now created dynamically with mne and are not cached any more. This can lead to small changes in results due to improved numerical precision. Delete old files to free up space with mne_experiment.rm('bem-sol-file', subject='*').

  • New MneExperiment.make_report_coreg() method.

  • New MneExperiment: analysis parameter connectivity

  • plot.TopoButterfly: press Shift-T for a large topo-map with sensor names.

New in 0.23

New in 0.22

  • Epoch Rejection GUI:

    • New “Tools” menu.

    • New “Info” tool to show summary info on the rejection.

    • New “Find Bad Channels” tool to automatically find bad channels.

    • Set marked channels by clicking on topo-map.

    • Faster page redraw.

  • plot.Barplot and plot.Boxplot: new cells argument to customize the order of bars/boxes.

  • MneExperiment: new method MneExperiment.show_rej_info().

  • NDVar: new method NDVar.label_clusters().

  • plot.configure(): option to revert to wxPython backend for plot.brain.

New in 0.21

  • MneExperiment:

    • New epoch parameters: trigger_shift and vars (see MneExperiment.epochs).

    • load_selected_events(): new vardef parameter to load variables from a test definition.

    • Log files stored in the root directory.

    • Parcellations (MneExperiment.parcs) based on combinations can also include split commands.

  • New Factor method: Factor.floodfill().

  • Model methods: get_table() replaced with as_table(), new head() and tail().

  • API: .sort_idx methods renamed to .sort_index.

New in 0.20

  • MneExperiment: new analysis parameter select_clusters='all' to keep all clusters in cluster tests (see select_clusters (cluster selection criteria)).

  • Use testnd.configure() to limit the number of CPUs that are used in permutation cluster tests.

New in 0.19

  • Two-stage tests (see MneExperiment.tests).

  • Safer cache-handling. See note at Analysis.

  • Dataset.head() and Dataset.tail() methods for more efficiently inspecting partial Datasets.

  • The default format for plots in reports is now SVG since they are displayed correctly in Safari 9.0. Use plot.configure() to change the default format.

  • API: Improvements in plot.Topomap with concomitant changes in the constructor signature. For examples see the meg/topographic plotting example.

  • API: plot.ColorList has a new argument called labels.

  • API: testnd.ANOVA attribute probability_maps renamed to p analogous to other testnd results.

  • Rejection-GUI: The option to plot the data range only has been removed.

New in 0.18

  • API: The first argument for MneExperiment.plot_annot() is now parc.

  • API: the fill_in_missing parameter to combine() has been deprecated and replaced with a new parameter called incomplete.

  • API: Several plotting functions have a new xticklabels parameter to suppress x-axis tick labels (e.g. plot.UTSStat).

  • The objects returned by plot.brain plotting functions now contain a plot_colorbar() method to create a corresponding plot.ColorBar plot.

  • New function choose() to combine data in different NDVars on a case by case basis.

  • Rejection-GUI (gui.select_epochs()): Press Shift-i when hovering over an epoch to enter channels for interpolation manually.

  • MneExperiment.show_file_status() now shows the last modification date of each file.

  • Under OS X 10.8 and newer running code under a notifier statement now automatically prevents the computer from going to sleep.

New in 0.17

  • MneExperiment.brain_plot_defaults can be used to customize PySurfer plots in movies and reports.

  • MneExperiment.trigger_shift can now also be a dictionary mapping subject name to shift value.

  • The rejection GUI now allows selecting individual channels for interpolation using the ‘i’ key.

  • Parcellations based on combinations of existing labels, as well as parcellations based on regions around points specified in MNI coordinates can now be defined in MneExperiment.parcs.

  • Source space NDVar can be indexed with lists of region names, e.g., ndvar.sub(source=['cuneus-lh', 'lingual-lh']).

  • API: plot.brain.bin_table() function signature changed slightly (more parameters, new hemi parameter inserted to match other functions’ argument order).

  • API: combine() now raises KeyError when trying to combine Dataset objects with unequal keys; set fill_in_missing=True to reproduce previous behavior.

  • API: Previously, Var.as_factor() mapped unspecified values to str(value). Now they are mapped to ''. This also applies to MneExperiment.variables entries with unspecified values.

New in 0.16

  • New function for plotting a legend for annot-files: plot.brain.annot_legend() (automatically used in reports).

  • Epoch definitions in MneExperiment.epochs can now include a 'base' parameter, which will copy the given “base” epoch and modify it with the current definition.

  • MneExperiment.make_mov_ttest() and MneExperiment.make_mov_ga_dspm() are fixed but require PySurfer 0.6.

  • New function: table.melt_ndvar().

  • API: plot.brain function signatures changed slightly to accommodate more layout-related arguments.

  • API: use Brain.image() instead of plot.brain.image().

New in 0.15

  • The Eelbrain package on the PYPI is now compiled with Anaconda. This means that the package can now be installed into an Anaconda distribution with pip, whereas easy_install has to be used for the Canopy distribution.

  • GUI gui.select_epochs(): Set marked channels through menu (View > Mark Channels)

  • Datasets can be saved as tables in RTF format (Dataset.save_rtf()).

  • API plot.Timeplot: the default spread indicator changed to SEM, and there is a new argument for timelabels.

  • API: test.anova is now a function with a slightly changed signature. The old class has been renamed to test.ANOVA.

  • API: test.oneway was removed. Use test.ANOVA.

  • API: the default value of the plot.Timeplot parameter bottom changed from 0 to None (determined by the data).

  • API: Factor.relabel() renamed to Factor.update_labels().

  • Plotting: New option for the figure legend 'draggable' (drag the legend with the mouse pointer).

New in 0.14

  • API: the plot.Topomap argument sensors changed to sensorlabels.

  • GUI: The python/Quit Eelbrain menu command now closes all windows to ensure that unsaved documents are handled properly. In order to yield to the terminal without closing windows, use the Go/Yield to Terminal command (Command-Alt-Q).

  • testnd.TContrastRelated: support for unary operation abs.

New in 0.13

  • The gui.select_epochs() GUI can now also be used to set bad channels. MneExperiment subclasses will combine bad channel information from rejection files with bad channel information from bad channel files. Note that while bad channel files set bad channels for a given raw file globally, rejection files set bad channels only for the given epoch.

  • Factor objects can now remember a custom cell order which determines the order in tables and plots.

  • The Var.as_factor() method can now assign all unmentioned codes to a default value.

  • MneExperiment:

    • API: Subclasses should remove the subject and experiment parameters from MneExperiment.label_events().

    • API: MneExperiment can now be imported directly from eelbrain.

    • API: The MneExperiment._defaults attribute should be renamed to MneExperiment.defaults.

    • A draft for a guide at The MneExperiment Pipeline.

    • Cached files are now saved in a separate folder at root/eelbrain-cache. The cache can be cleared using MneExperiment.clear_cache(). To preserve cached test results, move the root/test folder into the root/eelbrain-cache folder.

New in 0.12

  • API: Dataset construction changed, allows setting the number of cases in the Dataset.

  • API: plot.SensorMap2d was renamed to plot.SensorMap.

  • MneExperiment:

    • API: The default number of samples for reports is now 10’000.

    • New epoch parameter 'n_cases': raise an error if an epoch definition does not yield expected number of trials.

    • A custom baseline period for epochs can now be specified as a parameter in the epoch definition (e.g., 'baseline': (-0.2, -0.1)). When loading data, specifying baseline=True uses the epoch’s custom baseline.

New in 0.11

  • MneExperiment:

    • Change in the way the covariance matrix is defined: The epoch for the covariance matrix should be specified in MneExperiment.epochs['cov']. The regularization is no longer part of set_inv(), but is instead set with MneExperiment.set(cov='reg') or MneExperiment.set(cov='noreg').

    • New option cov='bestreg' automatically selects the regularization parameter for each subejct.

  • Var.as_factor() allows more efficient labeling when multiple values share the same label.

  • API: Previously plot.configure_backend() is now plot.configure()

New in 0.10

  • Tools for generating colors for categories (see Plotting).

  • Plots now all largely respect matplotlib rc-parameters (see Customizing Matplotlib).

  • Fixed an issue in the testnd module that could affect permutation based p-values when multiprocessing was used.

New in 0.9

  • Factor API change: The rep argument was renamed to repeat.

  • T-values for regression coefficients through NDVar.ols_t().

  • MneExperiment: subject name patterns and eog_sns are now handled automatically.

  • UTSStat and Barplot plots can use pooled error for variability estimates (on by default for related measures designs, can be turned off using the pool_error argument).

    • API: for consistency, the argument to specify the kind of error to plot changed to error in both plots.

New in 0.8

  • A new GUI application controls plots as well as the epoch selection GUI (see notes in the reference sections on Plotting and GUIs).

  • Randomization/Monte Carlo tests now seed the random state to make results replicable.

New in 0.6

  • New recipes (outdated).

New in 0.5

  • The eelbrain.lab and eelbrain.eellab modules are deprecated. Everything can now me imported from eelbrain directly.

New in 0.4

New in 0.3

  • Optimized clustering for cluster permutation tests.

New in 0.2

  • gui.SelectEpochs Epoch rejection GIU has a new “GA” button to plot the grand average of all accepted trials

  • Cluster permutation tests in testnd use multiple cores; To disable multiprocessing set eelbrain._stats.testnd.multiprocessing = False.

New in 0.1.7

  • gui.SelectEpochs can now be initialized with a single mne.Epochs instance (data needs to be preloaded).

  • Parameters that take NDVar objects now also accept mne.Epochs and mne.fiff.Evoked objects.

New in 0.1.5

  • plot.topo.TopoButterfly plot: new keyboard commands (t, left arrow, right arrow).

Publications using Eelbrain

Ordered alphabetically according to authors’ last name:

1

Esti Blanco-Elorrieta and Liina Pylkkänen. Bilingual language switching in the lab vs. in the wild: the spatio-temporal dynamics of adaptive language control. Journal of Neuroscience, pages 0553–17, 2017. URL: http://www.jneurosci.org/content/early/2017/08/16/JNEUROSCI.0553-17.2017.abstract.

2

Christian Brodbeck, Laura Gwilliams, and Liina Pylkkänen. EEG can track the time course of successful reference resolution in small visual worlds. Frontiers in psychology, 6:1787, 2015. URL: https://www.frontiersin.org/articles/10.3389/fpsyg.2015.01787.

3

Christian Brodbeck, Laura Gwilliams, and Liina Pylkkänen. Language in context: MEG evidence for modality general and specific responses to reference resolution. eNeuro, pages ENEURO–0145, 2016. URL: http://www.eneuro.org/content/early/2016/12/15/ENEURO.0145-16.2016.abstract.

4

Christian Brodbeck, L Elliot Hong, and Jonathan Z Simon. Rapid transformation from auditory to linguistic representations of continuous speech. Current Biology, 28(24):3976–3983, 2018. URL: https://www.sciencedirect.com/science/article/pii/S096098221831409X.

5

Christian Brodbeck, Alessandro Presacco, Samira Anderson, and Jonathan Z Simon. Over-representation of speech in older adults originates from early response in higher order auditory cortex. Acta Acustica united with Acustica, 104(5):774–777, 2018. URL: https://www.ingentaconnect.com/content/dav/aaua/2018/00000104/00000005/art00013.

6

Christian Brodbeck, Alessandro Presacco, and Jonathan Z Simon. Neural source dynamics of brain responses to continuous stimuli: speech processing from acoustics to comprehension. NeuroImage, 172:162–174, 2018. URL: https://www.sciencedirect.com/science/article/pii/S1053811918300429.

7

Christian Brodbeck and Liina Pylkkänen. Language in context: characterizing the comprehension of referential expressions with MEG. NeuroImage, 147:447–460, 2017. URL: https://www.sciencedirect.com/science/article/pii/S1053811916307169.

8

Teon L Brooks and Daniela Cid de Garcia. Evidence for morphological composition in compound words using MEG. Frontiers in human neuroscience, 9:215, 2015. URL: https://www.frontiersin.org/articles/10.3389/fnhum.2015.00215.

9

Julien Dirani and Liina Pylkkanen. Lexical access in comprehension vs. production: spatiotemporal localization of semantic facilitation and interference. bioRxiv, pages 449157, 2018. URL: https://www.biorxiv.org/content/early/2018/10/23/449157.1.abstract.

10

Graham Flick, Yohei Oseki, Amanda R Kaczmarek, Meera Al Kaabi, Alec Marantz, and Liina Pylkkänen. Building words and phrases in the left temporal lobe. Cortex, 2018. URL: https://www.sciencedirect.com/science/article/pii/S0010945218301904.

11

Graham Flick and Liina Pylkkänen. Isolating syntax in natural language: MEG evidence for an early contribution of left posterior temporal cortex. BioRxiv, pages 439158, 2018. URL: https://www.biorxiv.org/content/early/2018/10/09/439158.abstract.

12

Phoebe Gaston and Alec Marantz. The time course of contextual cohort effects in auditory processing of category-ambiguous words: MEG evidence for a single “clash” as noun or verb. Language, Cognition and Neuroscience, 33(4):402–423, 2018. URL: https://www.tandfonline.com/doi/abs/10.1080/23273798.2017.1395466.

13

Laura Gwilliams, GA Lewis, and Alec Marantz. Functional characterisation of letter-specific responses in time, space and current polarity using magnetoencephalography. NeuroImage, 132:320–333, 2016. URL: https://www.sciencedirect.com/science/article/pii/S105381191600166X.

14

Laura Gwilliams and Alec Marantz. Morphological representations are extrapolated from morpho-syntactic rules. Neuropsychologia, 114:77–87, 2018. URL: https://www.sciencedirect.com/science/article/pii/S0028393218301568.

15

Thomas Hartmann and Nathan Weisz. Auditory cortical generators of the frequency following response are modulated by intermodal attention. NeuroImage, 203:116185, 2019. doi:10.1016/j.neuroimage.2019.116185.

16

William Matchin, Christian Brodbeck, Christopher Hammerly, and Ellen Lau. The temporal dynamics of structure and content in sentence comprehension: evidence from fMRI-constrained MEG. Human brain mapping, 40(2):663–678, 2019. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/hbm.24403.

17

Kyriaki Neophytou, Christina Manouilidou, Linnaea Stockall, and Alec Marantz. Syntactic and semantic restrictions on morphological recomposition: MEG evidence from greek. Brain and language, 183:11–20, 2018. URL: https://www.sciencedirect.com/science/article/pii/S0093934X1730130X.

18

Krishna C Puvvada, Marisel Villafane-Delgado, Christian Brodbeck, and Jonathan Z Simon. Neural coding of noisy and reverberant speech in human auditory cortex. bioRxiv, pages 229153, 2017. URL: https://www.biorxiv.org/content/early/2017/12/04/229153.abstract.

19

Victoria Sharpe, Samir Reddigari, Liina Pylkkänen, and Alec Marantz. Automatic access to verb continuations on the lexical and categorical levels: evidence from MEG. Language, Cognition and Neuroscience, 34(2):137–150, 2019. URL: https://www.tandfonline.com/doi/abs/10.1080/23273798.2018.1531139.

20

Eline Verschueren, Jonas Vanthornhout, and Tom Francart. Semantic context enhances neural envelope tracking. bioRxiv, pages 421727, 2018. URL: https://www.biorxiv.org/content/early/2018/09/19/421727.abstract.

21

Eline Verschueren, Jonas Vanthornhout, and Tom Francart. The effect of stimulus choice on an EEG-based objective measure of speech intelligibility. bioRxiv, pages 421727, 2019. URL: https://www.biorxiv.org/content/10.1101/421727v2.full-text.

22

Adina Williams, Samir Reddigari, and Liina Pylkkänen. Early sensitivity of left perisylvian cortex to relationality in nouns and verbs. Neuropsychologia, 100:131–143, 2017. URL: https://www.sciencedirect.com/science/article/pii/S0028393217301586.

23

Linmin Zhang and Liina Pylkkänen. Composing lexical versus functional adjectives: evidence for uniformity in the left temporal lobe. Psychonomic bulletin & review, 25(6):2309–2322, 2018. URL: https://link.springer.com/article/10.3758/s13423-018-1469-y.

24

Linmin Zhang and Liina Pylkkänen. Semantic composition of sentences word by word: MEG evidence for shared processing of conceptual and logical elements. Neuropsychologia, 119:392–404, 2018. URL: https://www.sciencedirect.com/science/article/pii/S0028393218305037.

Development

Eelbrain is actively developed and maintained by Christian Brodbeck at the Computational sensorimotor systems lab at University of Maryland, College Park.

Eelbrain is fully open-source and new contributions are welcome on GitHub. Suggestions can be raised as issues, and modifications can be made as pull requests into the master branch.

The Development Version

The Eelbrain source code is hosted on GitHub. Development takes place on the master branch, while release versions are maintained on branches called r/0.26 etc. For further information on working with GitHub see GitHub’s instructions.

Installing the development version requires the presence of a compiler. On macOS, make sure Xcode is installed (open it once to accept the license agreement). Windows will indicate any needed files when the install command is run.

After cloning the repository, the development version can be installed by running, from the Eelbrain repository’s root directory:

$ python setup.py develop

On macOS, the $ eelbrain shell script to run iPython with the framework build is not installed properly by setup.py; in order to fix this, run:

$ ./fix-bin

In Python, you can make sure that you are working with the development version:

>>> import eelbrain
>>> eelbrain.__version__
'dev'

To switch back to the release version use $ pip uninstall eelbrain.

Building with Conda

To build Eelbrain with conda, make sure that conda-build is installed. Then, from Eelbrain/conda run:

$ conda build eelbrain

After building successfully, the build can be installed with:

$ conda install --use-local eelbrain

Contributing

Style guides:

Useful tools:

Testing

Tests for individual modules are included in folders called tests, usually on the same level as the module. To run all tests, run $ make test from the Eelbrain project directory. On macOS, tests needs to run with the framework build of Python; if you get a corresponding error, run $ ./fix-bin pytest from the Eelbrain repository root.

Reference

Data Classes

Primary data classes:

Dataset([items, name, caption, info, n_cases])

Store multiple variables pertaining to a common set of measurement cases

Factor(x[, name, random, repeat, tile, …])

Container for categorial data.

Var(x[, name, info, repeat, tile])

Container for scalar data.

NDVar(x, dims[, name, info])

Container for n-dimensional data.

Datalist([items, name, fmt])

list subclass for including lists in in a Dataset.

Model classes (not usually initialized by themselves but through operations on primary data-objects):

Interaction(base)

Represents an Interaction effect.

Model(x)

A list of effects.

NDVar dimensions (not usually initialized by themselves but through load functions):

Case(n[, connectivity])

Case dimension

Categorial(name, values[, connectivity])

Simple categorial dimension

Scalar(name, values[, unit, tick_format, …])

Scalar dimension

Sensor(locs[, names, sysname, proj2d, …])

Dimension class for representing sensor information

SourceSpace(vertices[, subject, src, …])

MNE surface-based source space

VolumeSourceSpace(vertices[, subject, src, …])

MNE volume source space

Space(directions[, name])

Represent multiple directions in space

UTS(tmin, tstep, nsamples[, unit])

Dimension object for representing uniform time series

File I/O

Eelbrain objects can be pickled. Eelbrain’s own pickle I/O functions provide backwards compatibility for pickled Eelbrain objects:

save.pickle(obj[, dest, protocol])

Pickle a Python object.

load.unpickle([path])

Load pickled Python objects from a file.

load.arrow([file_path])

Load object serialized with pyarrow.

save.arrow(obj[, dest])

Save a Python object with pyarrow.

load.update_subjects_dir(obj, subjects_dir)

Update NDVar SourceSpace.subjects_dir attributes

load.convert_pickle_protocol([root, …])

Re-save all pickle files with a specific protocol
Import

Functions and modules for loading specific file formats as Eelbrain object:

load.wav([filename, name, backend])

Load a wav file as NDVar

load.tsv([path, names, types, delimiter, …])

Load a Dataset from a text file.

load.eyelink

Tools for loading data form eyelink edf files.

load.fiff

Tools for importing data through mne.

load.txt

Tools for loading data from text files.

load.besa

Tools for loading data from the BESA-MN pipeline.
Export

Dataset with only univariate data can be saved as text using the save_txt() method. Additional export functions:

save.txt(iterator[, fmt, delim, dest])

Write any object that supports iteration to a text file

save.wav(ndvar[, filename, toint])

Save an NDVar as wav file

Sorting and Reordering

align(d1, d2[, i1, i2, out])

Align two data-objects based on index variables

align1(d, to[, by, out])

Align a data object to an index variable

Celltable(y[, x, match, sub, cat, ds, …])

Divide y into cells defined by x.

choose(choice, sources[, name])

Combine data-objects picking from a different object for each case

combine(items[, name, check_dims, …])

Combine a list of items of the same type into one item.

shuffled_index(n[, cells])

Return an index to shuffle a data-object

NDVar Initializers

gaussian(center, width, time)

Gaussian window NDVar

powerlaw_noise(dims, exponent[, seed])

Gaussian (1/f)^{exponent} noise.

NDVar Operations

NDVar objects support the native abs() and round() functions. See also NDVar methods.

Butterworth(low, high, order[, sfreq])

Butterworth filter

complete_source_space(ndvar[, fill, mask, to])

Fill in missing vertices on an NDVar with a partial source space

concatenate(ndvars[, dim, name, tmin, info, …])

Concatenate multiple NDVars

convolve(h, x[, ds, name])

Convolve h and x along the time dimension

correlation_coefficient(x, y[, dim, name])

Correlation between two NDVars along a specific dimension

cross_correlation(in1, in2[, name])

Cross-correlation between two NDVars along the time axis

cwt_morlet(y, freqs[, use_fft, n_cycles, …])

Time frequency decomposition with Morlet wavelets (mne-python)

dss(ndvar)

Denoising source separation (DSS)

filter_data(ndvar, l_freq, h_freq[, …])

Apply mne.filter.filter_data() to an NDVar

find_intervals(ndvar[, interpolate])

Find intervals from a boolean NDVar

find_peaks(ndvar)

Find local maxima in an NDVar

frequency_response(b[, frequencies])

Frequency response for a FIR filter

label_operator(labels[, operation, exclude, …])

Convert labeled NDVar into a matrix operation to extract label values

labels_from_clusters(clusters[, names])

Create Labels from source space clusters

maximum(ndvars[, name])

Element-wise maximum of multiple array-like objects

minimum(ndvars[, name])

Element-wise minimum of multiple array-like objects

morph_source_space(ndvar[, subject_to, …])

Morph source estimate to a different MRI subject

normalize_in_cells(y, for_dim[, in_cells, …])

Normalize data in cells to make it appropriate for ANOVA 1

neighbor_correlation(x[, dim, obs, name])

Calculate Neighbor correlation

psd_welch(ndvar[, fmin, fmax, n_fft, …])

Power spectral density with Welch’s method

resample(ndvar, sfreq[, npad, window, pad, name])

Resample an NDVar along the time dimension

segment(continuous, times, tstart, tstop[, …])

Segment a continuous NDVar

set_parc(data, parc[, dim])

Change the parcellation of an NDVar or SourceSpace dimension

set_time(ndvar, time[, mode])

Crop and/or pad an NDVar to match the time axis time

set_tmin(ndvar[, tmin])

Shift the time axis of an NDVar relative to its data

xhemi(ndvar[, mask, hemi, parc])

Project data from both hemispheres to hemi of fsaverage_sym

Reverse Correlation

boosting(y, x, tstart, tstop[, scale_data, …])

Estimate a linear filter with coordinate descent

BoostingResult(y, x, tstart, tstop, …[, …])

Result from boosting a temporal response function

epoch_impulse_predictor(shape[, value, …])

Time series with one impulse for each of n epochs

event_impulse_predictor(shape[, time, …])

Time series with multiple impulses

Tables

Manipulate data tables and compile information about data objects such as cell frequencies:

table.cast_to_ndvar(data, dim_values, match)

Create an NDVar by converting a data column to a dimension

table.difference(y, x, c1, c0, match[, sub, …])

Subtract data in one cell from another

table.frequencies(y[, x, of, sub, ds])

Calculate frequency of occurrence of the categories in y

table.melt(name, cells, cell_var_name, ds[, …])

Restructure a Dataset such that a measured variable is in a single column

table.melt_ndvar(ndvar[, dim, cells, ds, …])

Transform data to long format by converting an NDVar dimension into a variable

table.repmeas(y, x, match[, sub, ds])

Create a repeated-measures table

table.stats(y, row[, col, match, sub, fmt, …])

Make a table with statistics

Statistics

Univariate statistical tests:

test.Correlation(y, x[, sub, ds])

Pearson product moment correlation between y and x

test.TTestOneSample(y[, match, sub, ds, …])

One-sample t-test

test.TTestIndependent(y, x[, c1, c0, match, …])

Independent measures t-test

test.TTestRelated(y, x[, c1, c0, match, …])

Related-measures t-test

test.ANOVA(y, x[, sub, ds, title, caption])

Univariate ANOVA.

test.pairwise(y, x[, match, sub, ds, par, …])

Pairwise comparison table

test.ttest(y[, x, against, match, sub, …])

T-tests for one or more samples

test.correlations(y, x[, cat, sub, ds, asds])

Correlation with one or more predictors

test.pairwise_correlations(xs[, sub, ds, labels])

Pairwise correlation table

test.MannWhitneyU(y, x[, c1, c0, match, …])

Mann-Whitney U-test (non-parametric independent measures test)

test.WilcoxonSignedRank(y, x[, c1, c0, …])

Wilcoxon signed-rank test (non-parametric related measures test)

test.lilliefors(data[, formatted])

Lilliefors’ test for normal distribution

Mass-Univariate Statistics

testnd.TTestOneSample(y[, popmean, match, …])

Mass-univariate one sample t-test

testnd.TTestRelated(y, x[, c1, c0, match, …])

Mass-univariate related samples t-test

testnd.TTestIndependent(y, x[, c1, c0, …])

Mass-univariate independent samples t-test

testnd.TContrastRelated(y, x, contrast[, …])

Mass-univariate contrast based on t-values

testnd.ANOVA(y, x[, sub, ds, samples, pmin, …])

Mass-univariate ANOVA

testnd.Correlation(y, x[, norm, sub, ds, …])

Mass-univariate correlation

testnd.Vector(y[, match, sub, ds, samples, …])

Test a vector field for vectors with non-random direction

testnd.VectorDifferenceRelated(y, x[, c1, …])

Test difference between two vector fields for non-random direction

The tests in this module produce maps of statistical parameters, and implement different methods to compute corresponding maps of p-values that are corrected for multiple comparison:

Permutation for maximum statistic (samples=n)

Compute n parameter maps corresponding to n permutations of the data. In each permutation, store the maximum value of the test statistic across the map. This is the distribution of the maximum parameter in a map under the null hypothesis. Then, calculate a p-value for each data point in the original parameter map based on where it lies in this distribution.

Threshold-based clusters 1 (samples=n, pmin=p)

Find clusters of data points where the original statistic exceeds a value corresponding to an uncorrected p-value of p. For each cluster, calculate the sum of the statistic values that are part of the cluster (the cluster mass). Do the same in n permutations of the original data and retain for each permutation the value of the largest cluster. Evaluate all cluster values in the original data against the distributiom of maximum cluster values.

Threshold-free cluster enhancement 2 (samples=n, tfce=True)

Similar to permutation for maximum statstic, but each statistical parameter map is first processed with the cluster enhancement algorithm.

Two-stage tests

Two-stage tests proceed by estimating parameters for a fixed effects model for each subject, and then testing hypotheses on these parameter estimates on the group level. Two-stage tests are implemented by fitting an LM for each subject, and then combining them in a LMGroup to retrieve coefficients for group level statistics (see Two-stage test example).

testnd.LM(y, model[, ds, coding, subject, sub])

Fixed effects linear model

testnd.LMGroup(lms)

Group level analysis for linear model LM objects
References
1(1,2)

Maris, E., & Oostenveld, R. (2007). Nonparametric statistical testing of EEG- and MEG-data. Journal of Neuroscience Methods, 164(1), 177-190. 10.1016/j.jneumeth.2007.03.024

2

Smith, S. M., and Nichols, T. E. (2009). Threshold-Free Cluster Enhancement: Addressing Problems of Smoothing, Threshold Dependence and Localisation in Cluster Inference. NeuroImage, 44(1), 83-98. 10.1016/j.neuroimage.2008.03.061

Plotting

Plot univariate data (Var objects):

plot.Barplot(y[, x, match, sub, cells, …])

Barplot for a continuous variable

plot.BarplotHorizontal(y[, x, match, sub, …])

Horizontal barplot for a continuous variable

plot.Boxplot(y[, x, match, sub, cells, …])

Boxplot for a continuous variable

plot.PairwiseLegend([size, trend])

Legend for colors used in pairwise comparisons

plot.Scatter(y, x[, color, size, sub, ds, …])

Scatter-plot

plot.Histogram(y[, x, match, sub, ds, …])

Histogram plots with tests of normality

plot.Regression(y, x[, cat, match, sub, ds, …])

Plot the regression of y on x

plot.Timeplot(y, time[, categories, match, …])

Plot a variable over time

Color tools for plotting:

plot.colors_for_categorial(x[, hue_start, cmap])

Automatically select colors for a categorial model

plot.colors_for_oneway(cells[, hue_start, …])

Define colors for a single factor design

plot.colors_for_twoway(x1_cells, x2_cells[, …])

Define cell colors for a two-way design

plot.soft_threshold_colormap(cmap, …[, …])

Soft-threshold a colormap to make small values transparent

plot.two_step_colormap(left_max, left[, …])

Colormap using lightness to extend range

plot.Style([color, marker, hatch, …])

Control color/pattern by category

plot.ColorBar(cmap[, vmin, vmax, label, …])

A color-bar for a matplotlib color-map

plot.ColorGrid(row_cells, column_cells, colors)

Plot colors for a two-way design in a grid

plot.ColorList(colors[, cells, labels, size, h])

Plot colors with labels

See also

Example with Eelbrain Colormaps

Plot uniform time-series:

plot.LineStack(y[, x, sub, ds, offset, …])

Stack multiple lines vertically

plot.UTS(y[, xax, axtitle, ds, sub, xlabel, …])

Value by time plot for UTS data

plot.UTSClusters(res[, pmax, ptrend, …])

Plot permutation cluster test results

plot.UTSStat(y[, x, xax, match, sub, ds, …])

Plot statistics for a one-dimensional NDVar

Plot multidimensional uniform time series:

plot.Array(y[, xax, xlabel, ylabel, …])

Plot UTS data to a rectangular grid.

plot.Butterfly(y[, xax, sensors, axtitle, …])

Butterfly plot for NDVars

Plot topographic maps of sensor space data:

plot.TopoArray(y[, xax, ds, sub, vmax, …])

Channel by sample plots with topomaps for individual time points

plot.TopoButterfly(y[, xax, ds, sub, vmax, …])

Butterfly plot with corresponding topomaps

plot.Topomap(y[, xax, ds, sub, vmax, vmin, …])

Plot individual topogeraphies

plot.TopomapBins(y[, xax, ds, sub, …])

Topomaps in time-bins

Plot sensor layout maps:

plot.SensorMap(sensors[, labels, proj, …])

Plot sensor positions in 2 dimensions

plot.SensorMaps(sensors[, select, proj, …])

Multiple views on a sensor layout.

Method related plots:

plot.DataSplit(splits[, legend, colors])

xax parameter

Many plots have an xax parameter which is used to sort the data in y into different categories and plot them on separate axes. xax can be specified through categorial data, or as a dimension in y.

If a categorial data object is specified for xax, y is split into the categories in xax, and for every cell in xax a separate subplot is shown. For example, while

>>> plot.Butterfly('meg', ds=ds)

will create a single Butterfly plot of the average response,

>>> plot.Butterfly('meg', 'subject', ds=ds)

where 'subject' is the xax parameter, will create a separate subplot for every subject with its average response.

A dimension on y can be specified through a string starting with .. For example, to plot each case of meg separately, use:

>>> plot.Butterfly('meg', '.case', ds=ds)
General layout parameters

Most plots that also share certain layout keyword arguments. By default, all those parameters are determined automatically, but individual values can be specified manually by supplying them as keyword arguments.

h, wscalar

Height and width of the figure. Use a number ≤ 0 to defined the size relative to the screen (e.g., w=0 to use the full screen width).

axh, axwscalar

Height and width of the axes.

nrow, ncolNone | int

Limit number of rows/columns. If neither is specified, a square layout is produced

marginsdict

Absolute subplot parameters (in inches). Implies tight=False. If margins is specified, axw and axh are interpreted exclusive of the margins, i.e., axh=2, margins={'top': .5} for a plot with one axes will result in a total height of 2.5. Example:

margins={
    'top': 0.5,  # from top of figure to top axes
    'bottom': 1,  # from bottom of figure to bottom axes
    'hspace': 0.1,  # height of space between axes
    'left': 0.5,  #  from left of figure to left-most axes
    'right': 0.1,  # from right of figure to right-most axes
    'wspace': 0.1,  # width of space between axes
}
ax_aspectscalar

Width / height aspect of the axes.

framebool | str

How to frame the axes of the plot. Options:

  • True (default): normal matplotlib frame, spines around axes

  • False: omit top and right spines

  • 't': draw spines at x=0 and y=0, common for ERPs

  • 'none': no spines at all

namestr

Window title (not displayed on the figure itself).

titlestr

Figure title (displayed on the figure).

tightbool

Use matplotlib’s tight_layout to resize all axes to fill the figure.

right_ofeelbrain plot

Position the new figure to the right of this figure.

beloweelbrain plot

Position the new figure below this figure.

Plots that do take those parameters can be identified by the **layout in their function signature.

GUI Interaction

By default, new plots are automatically shown and, if the Python interpreter is in interactive mode the GUI main loop is started. This behavior can be controlled with 2 arguments when constructing a plot:

showbool

Show the figure in the GUI (default True). Use False for creating figures and saving them without displaying them on the screen.

runbool

Run the Eelbrain GUI app (default is True for interactive plotting and False in scripts).

The behavior can also be changed globally using configure().

By default, Eelbrain plots open in windows with enhance GUI features such as copying a figure to the OS clip-board. To plot figures in bare matplotlib windows, configure() Eelbrain with eelbrain.configure(frame=False).

Plotting Brains

The plot.brain module contains specialized functions to plot NDVar objects containing source space data. For this it uses a subclass of PySurfer’s surfer.Brain class. The functions below allow quick plotting. More specific control over the plots can be achieved through the Brain object that is returned.

plot.brain.brain(src[, cmap, vmin, vmax, …])

Create a Brain object with a data layer

plot.brain.butterfly(y[, cmap, vmin, vmax, …])

Shortcut for a Butterfly-plot with a time-linked brain plot

plot.brain.cluster(cluster[, vmax])

Plot a spatio-temporal cluster

plot.brain.dspm(src[, fmin, fmax, fmid])

Plot a source estimate with coloring for dSPM values (bipolar).

plot.brain.p_map(p_map[, param_map, p0, p1, …])

Plot a map of p-values in source space.

plot.brain.annot(annot[, subject, surf, …])

Plot the parcellation in an annotation file

plot.brain.annot_legend(lh, rh, *args, **kwargs)

Plot a legend for a freesurfer parcellation

Brain(subject, hemi[, surf, title, cortex, …])

PySurfer Brain subclass returned by plot.brain functions

plot.brain.SequencePlotter()

Grid of anatomical images in one figure

plot.GlassBrain(ndvar[, cmap, vmin, vmax, …])

Plot 2d projections of a brain volume

plot.GlassBrain.butterfly(y[, cmap, vmin, …])

Shortcut for a butterfly-plot with a time-linked glassbrain plot

In order to make custom plots, a Brain figure without any data added can be created with plot.brain.brain(ndvar.source, mask=False).

Surface options for plotting data on fsaverage:

white

surf_white

smoothwm

surf_smoothwm

inflated_pre

surf_inflated_pre

inflated

surf_inflated

inflated_avg

surf_inflated_avg

sphere

surf_sphere

GUIs

Tools with a graphical user interface (GUI):

gui.select_components(path, ds[, sysname, …])

GUI for selecting ICA-components

gui.select_epochs(ds[, data, accept, blink, …])

GUI for rejecting trials of MEG/EEG data

gui.load_stcs()

GUI for detecting and loading source estimates
Controlling the GUI Application

Eelbrain uses a wxPython based application to create GUIs. This GUI appears as a separate application with its own Dock icon. The way that control of this GUI is managed depends on the environment form which it is invoked.

When Eelbrain plots are created from within iPython, the GUI is managed in the background and control returns immediately to the terminal. There might be cases in which this is not desired, for example when running scripts. After execution of a script finishes, the interpreter is terminated and all associated plots are closed. To avoid this, the command gui.run(block=True) can be inserted at the end of the script, which will keep all gui elements open until the user quits the GUI application (see gui.run() below).

In interpreters other than iPython, input can not be processed from the GUI and the interpreter shell at the same time. In that case, the GUI application is activated by default whenever a GUI is created in interactive mode (this can be avoided by passing run=False to any plotting function). While the application is processing user input, the shell can not be used. In order to return to the shell, quit the application (the python/Quit Eelbrain menu command or Command-Q). In order to return to the terminal without closing all windows, use the alternative Go/Yield to Terminal command (Command-Alt-Q). To return to the application from the shell, run gui.run(). Beware that if you terminate the Python session from the terminal, the application is not given a chance to assure that information in open windows is saved.

gui.run([block])

Hand over command to the GUI (quit the GUI to return to the terminal)

Reports

The report submodule contains shortcuts for producing data summaries and visualizations. Meant to work in Jupyter notebooks.

report.scatter_table(xs[, color, sub, ds, …])

Table with pairwise scatter-plots (see plot.Scatter)
Formatted Text

The fmtxt submodule provides elements and tools for reports. Most eelbrain functions and methods that print tables in fact return fmtxt objects, which can be exported in different formats, for example:

>>> ds = datasets.get_uts()
>>> table.stats('Y', 'A', 'B', ds=ds)
           B
     ---------------
     b0       b1
--------------------
a0   0.8857   0.4525
a1   0.3425   1.377
>>> type(table.stats('Y', 'A', 'B', ds=ds))
eelbrain.fmtxt.Table

This means that the result can be exported as formatted text, for example:

>>> fmtxt.save_pdf(table.stats('Y', 'A', 'B', ds=ds))

See available export functions and the module documentation for details:

fmtxt

Document model for formatted text documents

fmtxt.copy_pdf(fmtext)

Copy an FMText object to the clipboard as PDF

fmtxt.copy_tex(fmtext)

Copy an FMText object to the clipboard as tex code

fmtxt.display(fmtext[, title])

Display the FMText rendered to HTML in a GUI window

fmtxt.save_html(fmtext[, path, …])

Save an FMText object in HTML format

fmtxt.save_pdf(fmtext[, path])

Save an FMText object as a PDF (requires LaTeX installation)

fmtxt.save_rtf(fmtext[, path])

Save an FMText object in Rich Text format

fmtxt.save_tex(fmtext[, path])

Save an FMText object as TeX code

Experiment Pipeline

The MneExperiment class provides a template for analyzing EEG and MEG data. The objects for specifying the analysis are all in the pipeline submodule.

See also

For the guide on working with the MneExperiment class see The MneExperiment Pipeline.

MneExperiment([root, find_subjects])

Analyze an MEG or EEG experiment

Result containers:

ROITestResult(subjects, samples, …)

Test results for temporal tests in one or more ROIs

ROI2StageResult(subjects, samples, …)

Test results for 2-stage tests in one or more ROIs

Participant groups:

Group(subjects)

Group defined as collection of subjects

SubGroup(base, exclude)

Group defined by removing subjects from a base group

Pre-processing:

RawSource([filename, reader, sysname, …])

Raw data source

RawFilter(source[, l_freq, h_freq, cache])

Filter raw pipe

RawICA(source, session[, method, …])

ICA raw pipe

RawMaxwell(source[, bad_condition, cache])

Maxwell filter raw pipe

RawReReference(source[, reference, add, …])

Re-reference EEG data

RawApplyICA(source, ica[, cache])

Apply ICA estimated in a RawICA pipe

Event variables:

EvalVar(code[, session])

Variable based on evaluating a statement

GroupVar(groups[, session])

Group membership for each subject

LabelVar(source, codes[, default, session])

Variable assigning labels to values

Epochs:

PrimaryEpoch(session[, sel])

Epoch based on selecting events from a raw file

SecondaryEpoch(base[, sel])

Epoch inheriting events from another epoch

SuperEpoch(sub_epochs, **kwargs)

Combine several other epochs

Tests:

ANOVA(x[, model, vars])

ANOVA test

TTestOneSample([tail])

One-sample t-test

TTestIndependent(model, c1, c0[, tail])

Independent measures t-test (comparing groups of subjects)

TTestRelated(model, c1, c0[, tail])

Related measures t-test

TContrastRelated(model, contrast[, tail])

Contrasts of T-maps (see eelbrain.testnd.TContrastRelated)

TwoStageTest(stage_1[, vars, model])

Two-stage test: T-test of regression coefficients

Brain parcellations:

SubParc(base, labels[, views])

A subset of labels in another parcellation

CombinationParc(base, labels[, views])

Recombine labels from an existing parcellation

FreeSurferParc([views])

Parcellation that is created outside Eelbrain for each subject

FSAverageParc([views])

Fsaverage parcellation that is morphed to individual subjects

SeededParc(seeds[, mask, surface, views])

Parcellation that is grown from seed coordinates

IndividualSeededParc(seeds[, mask, surface, …])

Seed parcellation with individual seeds for each subject

Datasets

Datasets for experimenting and testing:

datasets.get_loftus_masson_1994()

Dataset used for illustration purposes by Loftus and Masson (1994)

datasets.get_mne_sample([tmin, tmax, …])

Load events and epochs from the MNE sample data

datasets.get_uts([utsnd, seed, nrm, vector3d])

Create a sample Dataset with 60 cases and random data.

datasets.get_uv([seed, nrm, vector])

Dataset with random univariate data

datasets.simulate_erp([n_trials, seed])

Simulate event-related EEG data

Configuration

eelbrain.configure(n_workers=None, frame=None, autorun=None, show=None, format=None, figure_background=None, prompt_toolkit=None, animate=None, nice=None, tqdm=None)

Set basic configuration parameters for the current session

Parameters

  • n_workers (bool | int) – Number of worker processes to use in multiprocessing enabled computations. False to disable multiprocessing. True (default) to use as many processes as cores are available. Negative numbers to use all but n available CPUs.

  • frame (bool) – Open figures in the Eelbrain application. This provides additional functionality such as copying a figure to the clipboard. If False, open figures as normal matplotlib figures.

  • autorun (bool) – When a figure is created, automatically enter the GUI mainloop. By default, this is True when the figure is created in interactive mode but False when the figure is created in a script (in order to run the GUI at a specific point in a script, call eelbrain.gui.run()).

  • show (bool) – Show plots on the screen when they’re created (disable this to create plots and save them without showing them on the screen).

  • format (str) – Default format for plots (for example “png”, “svg”, …).

  • figure_background (bool | matplotlib color) – While matplotlib uses a gray figure background by default, Eelbrain uses white. Set this parameter to False to use the default from matplotlib.rcParams, or set it to a valid matplotblib color value to use an arbitrary color. True to revert to the default white.

  • prompt_toolkit (bool) – In IPython 5, prompt_toolkit allows running the GUI main loop in parallel to the Terminal, meaning that the IPython terminal and GUI windows can be used without explicitly switching between Terminal and GUI. This feature is enabled by default, but can be disabled by setting prompt_toolkit=False.

  • animate (bool) – Animate plot navigation (default True).

  • nice (int [-20, 19]) – Scheduling priority for muliprocessing (larger number yields more to other processes; negative numbers require root privileges).

  • tqdm (bool) – Enable or disable tqdm progress bars.

Examples

Note

Examples frequently use the print() function. In interactive work, this function can often be omitted. For example, to show a summary of a dataset, examples use >>> print(ds.summary()), but in interactive work simply entering >>> ds.summary() will print the summary.

Datasets

Examples of how to construct and manipulate Dataset objects.

Note

Click here to download the full example code

Introduction

Data are represented with there primary data-objects:

  • Factor for categorial variables

  • Var for scalar variables

  • NDVar for multidimensional data (e.g. a variable measured at different time points)

Multiple variables belonging to the same dataset can be grouped in a Dataset object.

Factor

A Factor is a container for one-dimensional, categorial data: Each case is described by a string label. The most obvious way to initialize a Factor is a list of strings:

# sphinx_gallery_thumbnail_number = 5
from eelbrain import *

a = Factor(['a', 'a', 'a', 'a', 'b', 'b', 'b', 'b'], name='A')
print(a)

Out:

Factor(['a', 'a', 'a', 'a', 'b', 'b', 'b', 'b'], name='A')

Since Factor initialization simply iterates over the given data, the same Factor could be initialized with:

a = Factor('aaaabbbb', name='A')
print(a)

Out:

Factor(['a', 'a', 'a', 'a', 'b', 'b', 'b', 'b'], name='A')

There are other shortcuts to initialize factors (see also the Factor class documentation):

a = Factor(['a', 'b', 'c'], repeat=4, name='A')
print(a)

Out:

Factor(['a', 'a', 'a', 'a', 'b', 'b', 'b', 'b', 'c', 'c', 'c', 'c'], name='A')

Indexing works like for arrays:

print(a[0])
print(a[0:6])

Out:

a
Factor(['a', 'a', 'a', 'a', 'b', 'b'], name='A')

All values present in a Factor are accessible in its Factor.cells attribute:

print(a.cells)

Out:

('a', 'b', 'c')

Based on the Factor’s cell values, boolean indexes can be generated:

print(a == 'a')
print(a.isany('a', 'b'))
print(a.isnot('a', 'b'))

Out:

[ True  True  True  True False False False False False False False False]
[ True  True  True  True  True  True  True  True False False False False]
[False False False False False False False False  True  True  True  True]

Interaction effects can be constructed from multiple factors with the % operator:

b = Factor(['d', 'e'], repeat=2, tile=3, name='B')
print(b)
i = a % b
print(i)

Out:

Factor(['d', 'd', 'e', 'e', 'd', 'd', 'e', 'e', 'd', 'd', 'e', 'e'], name='B')
A % B

Interaction effects are in many ways interchangeable with factors in places where a categorial model is required:

print(i.cells)
print(i == ('a', 'd'))

Out:

(('a', 'd'), ('a', 'e'), ('b', 'd'), ('b', 'e'), ('c', 'd'), ('c', 'e'))
[ True  True False False False False False False False False False False]
Var

The Var class is a container for one-dimensional numpy.ndarray:

y = Var([1, 2, 3, 4, 5, 6])
print(y)

Out:

Var([1, 2, 3, 4, 5, 6])

Indexing works as for factors

print(y[5])
print(y[2:])

Out:

6
Var([3, 4, 5, 6])

Many array operations can be performed on the object directly

print(y + 1)

Out:

Var([2, 3, 4, 5, 6, 7])

For any more complex operations the corresponding numpy.ndarray can be retrieved in the Var.x attribute:

print(y.x)

Out:

[1 2 3 4 5 6]

Note

The Var.x attribute is not intended to be replaced; rather, a new Var object should be created for a new array.

NDVar

NDVar objects are containers for multidimensional data, and manage the description of the dimensions along with the data. NDVar objects are often not constructed from scratch but imported from existing data. For example, mne source estimates can be imported with load.fiff.stc_ndvar(). As an example, consider data from a simulated EEG experiment:

ds = datasets.simulate_erp()
eeg = ds['eeg']
print(eeg)

Out:

<NDVar 'eeg': 80 case, 140 time, 65 sensor>

This representation shows that eeg contains 80 trials of data (cases), with 140 time points and 35 EEG sensors. Since eeg contains information on the dimensions like sensor locations, plotting functions can take advantage of that:

p = plot.TopoButterfly(eeg)
p.set_time(0.400)
intro

NDVar offer functionality similar to numpy.ndarray, but take into account the properties of the dimensions. For example, through the NDVar.sub() method, indexing can be done using meaningful descriptions, such as indexing a time slice in seconds:

eeg_400 = eeg.sub(time=0.400)
plot.Topomap(eeg_400)
intro

Out:

<Topomap: eeg>

Several methods allow aggregating data, for example an RMS over sensor:

eeg_rms = eeg.rms('sensor')
print(eeg_rms)
plot.UTSStat(eeg_rms)
intro

Out:

<NDVar 'eeg': 80 case, 140 time>

<UTSStat: eeg>

Or a mean in a time window:

eeg_400 = eeg.mean(time=(0.350, 0.450))
plot.Topomap(eeg_400)
intro

Out:

<Topomap: eeg>

As with a Var, the corresponding numpy.ndarray can always be accessed as array. The NDVar.get_data() method allows retrieving the data while being explicit about which axis represents which dimension:

array = eeg_400.get_data(('case', 'sensor'))
print(array.shape)

Out:

(80, 65)

NDVar objects can be constructed directly from an array and corresponding dimension objects, for example:

import numpy

frequency = Scalar('frequency', [1, 2, 3, 4])
time = UTS(0, 0.01, 50)
data = numpy.random.normal(0, 1, (4, 50))
ndvar = NDVar(data, (frequency, time))
print(ndvar)

Out:

<NDVar: 4 frequency, 50 time>

A case dimension can be added by including the bare Case class:

data = numpy.random.normal(0, 1, (10, 4, 50))
ndvar = NDVar(data, (Case, frequency, time))
print(ndvar)

Out:

<NDVar: 10 case, 4 frequency, 50 time>
Dataset

A Dataset is a container for multiple variables (Factor, Var and NDVar) that describe the same cases. It can be thought of as a data table with columns corresponding to different variables and rows to different cases. Variables can be assigned as to a dictionary:

ds = Dataset()
ds['x'] = Factor('aaabbb')
ds['y'] = Var([5, 4, 6, 2, 1, 3])
print(ds)

Out:

x   y
-----
a   5
a   4
a   6
b   2
b   1
b   3

A variable that’s equal in all cases can be assigned quickly:

ds[:, 'z'] = 0.

The string representation of a Dataset contains information on the variables stored in it:

# in an interactive shell this would be the output of just typing ``ds``
print(repr(ds))

Out:

<Dataset (6 cases) 'x':F, 'y':V, 'z':V>

n_cases=6 indicates that the Dataset contains 6 cases (rows). The subsequent dictionary-like representation shows the keys and the types of the corresponding values (F: Factor, V: Var, Vnd: NDVar).

A more extensive summary can be printed with the Dataset.summary() method:

print(ds.summary())

Out:

Key   Type     Values
-------------------------------
x     Factor   a:3, b:3
y     Var      1, 2, 3, 4, 5, 6
z     Var      0:6
-------------------------------
Dataset: 6 cases

Indexing a Dataset with strings returns the corresponding data-objects:

print(ds['x'])

Out:

Factor(['a', 'a', 'a', 'b', 'b', 'b'], name='x')

numpy.ndarray-like indexing on the Dataset can be used to access a subset of cases:

print(ds[2:])

Out:

x   y   z
---------
a   6   0
b   2   0
b   1   0
b   3   0

Row and column can be indexed simultaneously (in row, column order):

print(ds[2, 'x'])

Out:

a

Arry-based indexing also allows indexing based on the Dataset’s variables:

print(ds[ds['x'] == 'a'])

Out:

x   y   z
---------
a   5   0
a   4   0
a   6   0

Since the dataset acts as container for variable, there is a Dataset.eval() method for evaluatuing code strings in the namespace defined by the dataset, which means that dataset variables can be invoked with just their name:

print(ds.eval("x == 'a'"))

Out:

[ True  True  True False False False]

Many dataset methods allow using code strings as shortcuts for expressions involving dataset variables, for example indexing:

print(ds.sub("x == 'a'"))

Out:

x   y   z
---------
a   5   0
a   4   0
a   6   0
Example

Below is a simple example using data objects (for more, see the Examples):

y = numpy.empty(21)
y[:14] = numpy.random.normal(0, 1, 14)
y[14:] = numpy.random.normal(2, 1, 7)
ds = Dataset({
    'a': Factor('abc', 'A', repeat=7),
    'y': Var(y, 'Y'),
})
print(ds)

Out:

a   y
-------------
a   -1.0655
a   -1.4081
a   -0.11117
a   -1.1912
a   -0.2207
a   0.81642
a   0.37445
b   -0.036294
b   -1.1602
b   0.91944
b   -1.0213
b   0.57937
b   -0.42151
b   -1.8369
c   1.1327
c   2.4491
c   2.803
c   1.4325
c   2.8258
c   2.0877
c   3.6062

print(table.frequencies('a', ds=ds))

Out:

a   n
-----
a   7
b   7
c   7

print(test.ANOVA('y', 'a', ds=ds))

Out:

               SS   df      MS          F        p
--------------------------------------------------
a           35.22    2   17.61   21.79***   < .001
Residuals   14.55   18    0.81
--------------------------------------------------
Total       49.76   20

print(test.pairwise('y', 'a', ds=ds, corr='Hochberg'))

Out:

Pairwise T-Tests (independent samples)

    b                 c
--------------------------------------
a   t(12) = 0.05      t(12) = -6.01***
    p = .961          p < .001
    p(c) = .961       p(c) < .001
b                     t(12) = -5.58***
                      p < .001
                      p(c) < .001
--------------------------------------
(* Corrected after Hochberg, 1988)

t = test.pairwise('y', 'a', ds=ds, corr='Hochberg')
print(t.get_tex())

Out:

\begin{center}
\begin{tabular}{lll}
\toprule
 & b & c \\
\midrule
a & $t_{12} = 0.05^{   \ \ \ }$ & $t_{12} = -6.01^{***}$ \\
 & $p = .961$ & $p < .001$ \\
 & $p_{c} = .961$ & $p_{c} < .001$ \\
b &  & $t_{12} = -5.58^{***}$ \\
 &  & $p < .001$ \\
 &  & $p_{c} < .001$ \\
\bottomrule
\end{tabular}
\end{center}

plot.Boxplot('y', 'a', ds=ds, title="My Boxplot", ylabel="value", corr='Hochberg')
My Boxplot

Out:

<Boxplot: My Boxplot>

Gallery generated by Sphinx-Gallery

Note

Click here to download the full example code

Dataset basics
# Author: Christian Brodbeck <christianbrodbeck@nyu.edu>
from eelbrain import *
import numpy as np

A dataset can be constructed column by column, by adding one variable after another:

# initialize an empty Dataset:
ds = Dataset()
# numeric values are added as Var object:
ds['y'] = Var(np.random.normal(0, 1, 6))
# categorical data as represented in Factors:
ds['a'] = Factor(['a', 'b', 'c'], repeat=2)
# check the result:
print(ds)

Out:

y          a
------------
0.8465     a
-1.7586    a
-3.0142    b
-0.84743   b
0.46509    c
1.8053     c

For larger datasets it can be more convenient to print only the first few cases…

print(ds.head())

Out:

y          a
------------
0.8465     a
-1.7586    a
-3.0142    b
-0.84743   b
0.46509    c
1.8053     c

… or a summary of variables:

print(ds.summary())

Out:

Key   Type     Values
-------------------------------------------------------------------------
y     Var      -3.01424, -1.75855, -0.847432, 0.465085, 0.846497, 1.80532
a     Factor   a:2, b:2, c:2
-------------------------------------------------------------------------
Dataset: 6 cases

An alternative way of constructing a dataset is case by case (i.e., row by row):

rows = []
for i in range(6):
    subject = f'S{i}'
    y = np.random.normal(0, 1)
    a = 'abc'[i % 3]
    rows.append([subject, y, a])
ds = Dataset.from_caselist(['subject', 'y', 'a'], rows, random='subject')
print(ds)

Out:

subject   y          a
----------------------
S0        0.089567   a
S1        -0.64756   b
S2        1.3846     c
S3        -0.62415   a
S4        -1.7727    b
S5        0.85463    c

Gallery generated by Sphinx-Gallery

Note

Click here to download the full example code

Align datasets

Shows how to combine information from two datasets describing the same cases, but not necessarily in the same order.

# Author: Christian Brodbeck <christianbrodbeck@nyu.edu>
import random
import string

from eelbrain import *


# Generate a dataset with known sequence
ds = Dataset()
ds['ascii'] = Factor(string.ascii_lowercase)
# Add an index variable to the dataset to later identify the cases
ds.index()

# Generate two shuffled copies of the dataset (and print them to confirm that
# they are shuffled)
ds1 = ds[random.sample(range(ds.n_cases), 15)]
print(ds1.head())
ds2 = ds[random.sample(range(ds.n_cases), 16)]
print(ds2.head())

Out:

ascii   index
-------------
o       14
m       12
q       16
a       0
y       24
g       6
w       22
r       17
d       3
l       11
ascii   index
-------------
v       21
h       7
n       13
t       19
b       1
s       18
g       6
u       20
e       4
m       12
Align the datasets

Use the "index" variable added above to identify cases and align the two datasets

ds1_aligned, ds2_aligned = align(ds1, ds2, 'index')

# show the ascii sequences for the two datasets next to each other to
# demonstrate that they are aligned
ds1_aligned['ascii_ds2'] = ds2_aligned['ascii']
print(ds1_aligned)

Out:

ascii   index   ascii_ds2
-------------------------
o       14      o
m       12      m
a       0       a
g       6       g
d       3       d
x       23      x
e       4       e
s       18      s

Gallery generated by Sphinx-Gallery

Univariate Statistics

Note

Click here to download the full example code

Model coding

Illustrates how to inspect coding of regression models.

# Author: Christian Brodbeck <christianbrodbeck@nyu.edu>
from eelbrain import *

ds = Dataset()
ds['A'] = Factor(['a1', 'a0'], repeat=4)
ds['B'] = Factor(['b1', 'b0'], repeat=2, tile=2)

look at data

ds.head()

Out:

A    B
-------
a1   b1
a1   b1
a1   b0
a1   b0
a0   b1
a0   b1
a0   b0
a0   b0

create a fixed effects model

m = ds.eval('A * B')
print(repr(m))

Out:

A + B + A % B

show the model coding

print(m)

Out:

intercept   A   B   A x B
-------------------------
1           0   0   0
1           0   0   0
1           0   1   0
1           0   1   0
1           1   0   0
1           1   0   0
1           1   1   1
1           1   1   1

create random effects model

ds['subject'] = Factor(['s1', 's2'], tile=4, name='subject', random=True)
m = ds.eval('A * B * subject')
print(repr(m))

Out:

A + B + A % B + subject + A % subject + B % subject + A % B % subject

show the model coding

print(m)

Out:

intercept   A   B   A x B   subject   A x subject   B x subject   A x B x subject
---------------------------------------------------------------------------------
1           0   0   0       1         0             0             0
1           0   0   0       0         0             0             0
1           0   1   0       1         0             1             0
1           0   1   0       0         0             0             0
1           1   0   0       1         1             0             0
1           1   0   0       0         0             0             0
1           1   1   1       1         1             1             1
1           1   1   1       0         0             0             0

Gallery generated by Sphinx-Gallery

Note

Click here to download the full example code

Repeated measures ANOVA

Based on 1.

# Author: Christian Brodbeck <christianbrodbeck@nyu.edu>
from eelbrain import *

y = Var([7,  3,  6,  6,  5,  8,  6,  7,
         7, 11,  9, 11, 10, 10, 11, 11,
         8, 14, 10, 11, 12, 10, 11, 12],
        name='y')
a = Factor('abc', repeat=8, name='A')

Fixed effects ANOVA (independent measures, 1 p. 24):

print(test.ANOVA(y, a, title="Independent Measures"))

Out:

Independent Measures

                SS   df      MS          F        p
---------------------------------------------------
A           112.00    2   56.00   22.62***   < .001
Residuals    52.00   21    2.48
---------------------------------------------------
Total       164.00   23

Repeated measures ANOVA (1 p. 72): subject is defined as random effect and entered for model construction as completely crossed factor

subject = Factor(range(8), tile=3, name='subject', random=True)
print(test.ANOVA(y, a * subject, title="Repeated Measures"))

Out:

Repeated Measures

            SS   df      MS   MS(denom)   df(denom)          F        p
-----------------------------------------------------------------------
A       112.00    2   56.00        2.71          14   20.63***   < .001
-----------------------------------------------------------------------
Total   164.00   23
Two-way repeated measures ANOVA
y = Var([ 7,  3,  6,  6,  5,  8,  6,  7,
          7, 11,  9, 11, 10, 10, 11, 11,
          8, 14, 10, 11, 12, 10, 11, 12,
         16,  7, 11,  9, 10, 11,  8,  8,
         16, 10, 13, 10, 10, 14, 11, 12,
         24, 29, 10, 22, 25, 28, 22, 24])
a = Factor(['a0', 'a1'], repeat=3 * 8, name='A')
b = Factor(['b0', 'b1', 'b2'], tile=2, repeat=8, name='B')
subject = Factor(range(8), tile=6, name='subject', random=True)

print(test.ANOVA(y, a * b * subject, title="Repeated Measure:"))

Out:

Repeated Measure:

             SS   df       MS   MS(denom)   df(denom)          F        p
-------------------------------------------------------------------------
A        432.00    1   432.00       10.76           7   40.14***   < .001
B        672.00    2   336.00       11.50          14   29.22***   < .001
A x B    224.00    2   112.00        6.55          14   17.11***   < .001
-------------------------------------------------------------------------
Total   1708.00   47

Bar-plot with within-subject error bars and pairwise tests

plot.Barplot(y, a % b, match=subject)
RMANOVA

Out:

<Barplot: None ~ A x B>
References
1(1,2,3)

Rutherford, A. (2001). Introducing ANOVA and ANCOVA: A GLM Approach. Sage.

Gallery generated by Sphinx-Gallery

Note

Click here to download the full example code

ANCOVA
Example 1

Based on 1, Exercises (page 8).

# Author: Christian Brodbeck <christianbrodbeck@nyu.edu>
from eelbrain import *

y = Var([2, 3, 3, 4,
         3, 4, 5, 6,
         1, 2, 1, 2,
         1, 1, 2, 2,
         2, 2, 2, 2,
         1, 1, 2, 3], name="Growth Rate")

genotype = Factor(range(6), repeat=4, name="Genotype")

hours = Var([8, 12, 16, 24], tile=6, name="Hours")

Show the model

print(hours * genotype)

Out:

intercept   Hours   Genotype                 Hours x Genotype
-----------------------------------------------------------------------------
1           8       1    0    0    0    0    8      0      0      0      0
1           12      1    0    0    0    0    12     0      0      0      0
1           16      1    0    0    0    0    16     0      0      0      0
1           24      1    0    0    0    0    24     0      0      0      0
1           8       0    1    0    0    0    0      8      0      0      0
1           12      0    1    0    0    0    0      12     0      0      0
1           16      0    1    0    0    0    0      16     0      0      0
1           24      0    1    0    0    0    0      24     0      0      0
1           8       0    0    1    0    0    0      0      8      0      0
1           12      0    0    1    0    0    0      0      12     0      0
1           16      0    0    1    0    0    0      0      16     0      0
1           24      0    0    1    0    0    0      0      24     0      0
1           8       0    0    0    1    0    0      0      0      8      0
1           12      0    0    0    1    0    0      0      0      12     0
1           16      0    0    0    1    0    0      0      0      16     0
1           24      0    0    0    1    0    0      0      0      24     0
1           8       0    0    0    0    1    0      0      0      0      8
1           12      0    0    0    0    1    0      0      0      0      12
1           16      0    0    0    0    1    0      0      0      0      16
1           24      0    0    0    0    1    0      0      0      0      24
1           8       0    0    0    0    0    0      0      0      0      0
1           12      0    0    0    0    0    0      0      0      0      0
1           16      0    0    0    0    0    0      0      0      0      0
1           24      0    0    0    0    0    0      0      0      0      0

ANCOVA

print(test.ANOVA(y, hours * genotype, title="ANOVA"))

Out:

ANOVA

                      SS   df     MS          F        p
--------------------------------------------------------
Hours               7.06    1   7.06   54.90***   < .001
Genotype           27.88    5   5.58   43.36***   < .001
Hours x Genotype    3.15    5   0.63    4.90*       .011
Residuals           1.54   12   0.13
--------------------------------------------------------
Total              39.62   23

Plot the slopes

plot.Regression(y, hours, genotype)
ANCOVA

Out:

<Regression: Growth Rate ~ Hours | Genotype>
Example 2

Based on 2 (p. 118-20)

y = Var([16,  7, 11,  9, 10, 11,  8,  8,
         16, 10, 13, 10, 10, 14, 11, 12,
         24, 29, 10, 22, 25, 28, 22, 24])

cov = Var([9, 5, 6, 4, 6, 8, 3, 5,
           8, 5, 6, 5, 3, 6, 4, 6,
           5, 8, 3, 4, 6, 9, 4, 5], name='cov')

a = Factor([1, 2, 3], repeat=8, name='A')

Full model, with interaction

print(test.ANOVA(y, cov * a))
plot.Regression(y, cov, a)
ANCOVA

Out:

                 SS   df       MS          F        p
-----------------------------------------------------
cov          199.54    1   199.54   32.93***   < .001
A            807.82    2   403.91   66.66***   < .001
cov x A       19.39    2     9.70    1.60        .229
Residuals    109.07   18     6.06
-----------------------------------------------------
Total       1112.00   23

<Regression: None ~ cov | A>

Drop interaction term

print(test.ANOVA(y, a + cov))
plot.Regression(y, cov)
ANCOVA

Out:

                 SS   df       MS          F        p
-----------------------------------------------------
A            807.82    2   403.91   62.88***   < .001
cov          199.54    1   199.54   31.07***   < .001
Residuals    128.46   20     6.42
-----------------------------------------------------
Total       1112.00   23

<Regression: None ~ cov>
ANCOVA with multiple covariates

Based on 3, p. 139.

# Load data form a text file
ds = load.txt.tsv('Fox_Prestige_data.txt', delimiter=' ', skipinitialspace=True)
print(ds.head())

Out:

id                    education   income   women   prestige   census   type
---------------------------------------------------------------------------
GOV.ADMINISTRATORS    13.11       12351    11.16   68.8       1113     prof
GENERAL.MANAGERS      12.26       25879    4.02    69.1       1130     prof
ACCOUNTANTS           12.77       9271     15.7    63.4       1171     prof
PURCHASING.OFFICERS   11.42       8865     9.11    56.8       1175     prof
CHEMISTS              14.62       8403     11.68   73.5       2111     prof
PHYSICISTS            15.64       11030    5.13    77.6       2113     prof
BIOLOGISTS            15.09       8258     25.65   72.6       2133     prof
ARCHITECTS            15.44       14163    2.69    78.1       2141     prof
CIVIL.ENGINEERS       14.52       11377    1.03    73.1       2143     prof
MINING.ENGINEERS      14.64       11023    0.94    68.8       2153     prof

# Variable summary
print(ds.summary())

Out:

Key         Type     Values
-------------------------------------------------------------------------------------------------
id          Factor   ACCOUNTANTS, AIRCRAFT.REPAIRMEN, AIRCRAFT.WORKERS, ARCHITECTS... (102 cells)
education   Var      6.38 - 15.97
income      Var      611 - 25879
women       Var      0 - 97.51
prestige    Var      14.8 - 87.2
census      Var      1113 - 9517
type        Factor   NA:4, bc:44, prof:31, wc:23
-------------------------------------------------------------------------------------------------
Fox_Prestige_data.txt: 102 cases

# Exclude cases with missing type
ds2 = ds[ds['type'] != 'NA']

# ANOVA
print(test.ANOVA('prestige', '(income + education) * type', ds=ds2))

Out:

                         SS   df        MS          F        p
--------------------------------------------------------------
income              1131.90    1   1131.90   28.35***   < .001
education           1067.98    1   1067.98   26.75***   < .001
type                 591.16    2    295.58    7.40**      .001
income x type        951.77    2    475.89   11.92***   < .001
education x type     238.40    2    119.20    2.99        .056
Residuals           3552.86   89     39.92
--------------------------------------------------------------
Total              28346.88   97
References
1

Crawley, M. J. (2005). Statistics: an introduction using R. J Wiley.

2

Rutherford, A. (2001). Introducing ANOVA and ANCOVA: A GLM Approach. Sage.

3

Fox, J. (2008) Applied Regression Analysis and Generalized Linear Models, Second Edition. Sage.

Gallery generated by Sphinx-Gallery

Mass-univariate statistics

NDVar, n-dimensional variables, allow mass-univariate statistics analogous to univariate statistics.

Note

Click here to download the full example code

Generate NDVar (with artificial data)

Shows how to initialize an NDVar with the structure of EEG data from (randomly generate) data arrays. The data is intended for illustrating EEG analysis techniques and meant to vaguely resemble data from an N400 experiment, but it is not meant to be a physiologically realistic simulation.

# sphinx_gallery_thumbnail_number = 2
import numpy as np
import scipy.spatial
from eelbrain import *

Create a Sensor dimension from an actual montage

sensor = Sensor.from_montage('standard_alphabetic')
p = plot.SensorMap(sensor)
simulate N400

Generate N400-like topography

i_cz = sensor.names.index('Cz')
cz_loc = sensor.locs[i_cz]
dists = scipy.spatial.distance.cdist([cz_loc], sensor.locs)[0]
dists /= dists.max()
topo = -0.7 + dists
n400_topo = NDVar(topo, sensor)
p = plot.Topomap(n400_topo, clip='circle')
simulate N400

Generate N400-like timing

window = scipy.signal.windows.gaussian(200, 12)[:140]
time = UTS(-0.100, 0.005, 140)
n400_timecourse = NDVar(window, time)
p = plot.UTS(n400_timecourse)
simulate N400

Generate random values for the independent variable (call it “cloze probability”)

rng = np.random.RandomState(0)
n_trials = 100
cloze_x = np.concatenate([
    rng.uniform(0, 0.3, n_trials // 2),
    rng.uniform(0.8, 1.0, n_trials // 2),
])
rng.shuffle(cloze_x)
cloze = Var(cloze_x)
p = plot.Histogram(cloze)
simulate N400

Put all the dimensions together to simulate the EEG signal

signal = (1 - cloze) * n400_timecourse * n400_topo

# Add noise
noise = powerlaw_noise(signal, 1)
noise = noise.smooth('sensor', 0.02, 'gaussian')
signal += noise

# Apply the average mastoids reference
signal -= signal.mean(sensor=['M1', 'M2'])

# Store EEG data in a Dataset with trial information
ds = Dataset()
ds['eeg'] = signal
ds['cloze'] = Var(cloze_x)
ds['cloze_cat'] = Factor(cloze_x > 0.5, labels={True: 'high', False: 'low'})

Plot the average simulated response

p = plot.TopoButterfly('eeg', ds=ds, vmax=1.5, clip='circle', frame='t', axh=3)
p.set_time(0.400)
simulate N400

Plot averages separately for high and low cloze

p = plot.TopoButterfly('eeg', 'cloze_cat', ds=ds, vmax=1.5, clip='circle', frame='t', axh=3)
p.set_time(0.400)
simulate N400

Average over time in the N400 time window

p = plot.Topomap('eeg.mean(time=(0.300, 0.500))', 'cloze_cat', ds=ds, vmax=1, clip='circle')
high, low

Plot the first 20 trials, labeled with cloze propability

labels = [f'{i} ({c:.2f})' for i, c in enumerate(cloze[:20])]
p = plot.Butterfly('eeg[:20]', '.case', ds=ds, axtitle=labels)
0 (0.25), 1 (0.97), 2 (0.92), 3 (0.29), 4 (0.16), 5 (0.13), 6 (0.01), 7 (0.85), 8 (0.91), 9 (0.87), 10 (0.82), 11 (0.06), 12 (0.82), 13 (0.18), 14 (0.02), 15 (0.11), 16 (0.20), 17 (0.83), 18 (0.08), 19 (0.93)

Total running time of the script: ( 0 minutes 6.954 seconds)

Gallery generated by Sphinx-Gallery

Note

Click here to download the full example code

Cluster-based permutation t-test

This example show a cluster-based permutation test for a simple design (two conditions). The example uses simulated data meant to vaguely resemble data from an N400 experiment (not intended as a physiologically realistic simulation).

# sphinx_gallery_thumbnail_number = 3
from eelbrain import *
Simulated data

Each function call to datasets.simulate_erp() generates a dataset equivalent to an N400 experiment for one subject. The seed argument determines the random noise that is added to the data.

ds = datasets.simulate_erp(seed=0)
print(ds.summary())

Out:

Key         Type     Values
-------------------------------------------------------------------
eeg         NDVar    140 time, 65 sensor; -1.4998e-05 - 1.64215e-05
cloze       Var      0.00563694 - 0.997675
cloze_cat   Factor   high:40, low:40
n_chars     Var      0:2, 1:5, 2:6, 3:13, 4:23, 5:20, 6:7, 7:4
-------------------------------------------------------------------
Dataset: 80 cases

A singe trial of data:

p = plot.TopoButterfly('eeg[0]', ds=ds)
p.set_time(0.400)
sensor cluster based ttest

The Dataset.aggregate() method computes condition averages when sorting the data into conditions of interest. In our case, the cloze_cat variable specified conditions 'high' and 'low' cloze:

print(ds.aggregate('cloze_cat'))

Out:

n    cloze     cloze_cat   n_chars
----------------------------------
40   0.88051   high        4.125
40   0.17241   low         3.825
----------------------------------
NDVars: eeg
Group level data

This loop simulates a multi-subject experiment. It generates data and collects condition averages for 10 virtual subjects. For group level analysis, the collected data are combined in a Dataset:

dss = []
for subject in range(10):
    # generate data for one subject
    ds = datasets.simulate_erp(seed=subject)
    # average across trials to get condition means
    ds_agg = ds.aggregate('cloze_cat')
    # add the subject name as variable
    ds_agg[:, 'subject'] = f'S{subject:02}'
    dss.append(ds_agg)

ds = combine(dss)
# make subject a random factor (to treat it as random effect for ANOVA)
ds['subject'].random = True
print(ds.head())

Out:

n    cloze     cloze_cat   n_chars   subject
--------------------------------------------
40   0.88051   high        4.125     S00
40   0.17241   low         3.825     S00
40   0.89466   high        4         S01
40   0.13778   low         4.575     S01
40   0.90215   high        4.275     S02
40   0.12206   low         3.9       S02
40   0.88503   high        4         S03
40   0.14273   low         4.175     S03
40   0.90499   high        4.1       S04
40   0.15732   low         3.5       S04
--------------------------------------------
NDVars: eeg

Re-reference the EEG data (i.e., subtract the mean of the two mastoid channels):

ds['eeg'] -= ds['eeg'].mean(sensor=['M1', 'M2'])
Spatio-temporal cluster based test

Cluster-based tests are based on identifying clusters of meaningful effects, i.e., groups of adjacent sensors that show the same effect (see testnd for references). In order to find clusters, the algorithm needs to know which channels are neighbors. This information is refered to as the sensor connectivity (i.e., which sensors are connected). The connectivity graph can be visualized to confirm that it is set correctly.

p = plot.SensorMap(ds['eeg'], connectivity=True)
sensor cluster based ttest

With the correct connectivity, we can now compute a cluster-based permutation test for a related measures t-test:

res = testnd.TTestRelated(
    'eeg', 'cloze_cat', 'low', 'high', match='subject', ds=ds,
    pmin=0.05,  # Use uncorrected p = 0.05 as threshold for forming clusters
    tstart=0.100,  # Find clusters in the time window from 100 ...
    tstop=0.600,  # ... to 600 ms
)

Plot the test result. The top two rows show the two condition averages. The bottom butterfly plot shows the difference with only significant regions in color. The corresponding topomap shows the topography at the marked time point, with the significant region circled (in an interactive environment, the mouse can be used to update the time point shown).

p = plot.TopoButterfly(res, clip='circle')
p.set_time(0.400)
sensor cluster based ttest

Generate a table with all significant clusters:

clusters = res.find_clusters(0.05)
print(clusters)

Out:

id   tstart   tstop   duration   n_sensors   v         p   sig
--------------------------------------------------------------
51   0.35     0.46    0.11       59          -2657.7   0   ***

Retrieve the cluster map using its ID and visualize the spatio-temporal extent of the cluster:

cluster_id = clusters[0, 'id']
cluster = res.cluster(cluster_id)
p = plot.TopoArray(cluster, interpolation='nearest')
p.set_topo_ts(.350, 0.400, 0.450)
sensor cluster based ttest

Often it is desirable to summarize values in a cluster. This is especially useful in more complex designs. For example, after finding a signficant interaction effect in an ANOVA, one might want to follow up with a pairwise test of the value in the cluster. This can often be achieved using binary masks based on the cluster. Using the cluster identified above, generate a binary mask:

mask = cluster != 0
p = plot.TopoArray(mask, cmap='Wistia')
p.set_topo_ts(.350, 0.400, 0.450)
sensor cluster based ttest

Such a spatio-temporal boolean mask can be used to extract the value in the cluster for each condition/participant. Since mask contains both time and sensor dimensions, using it with the NDVar.mean() method collapses across these dimensions and returns a scalar for each case (i.e., for each condition/subject).

ds['cluster_mean'] = ds['eeg'].mean(mask)
p = plot.Barplot('cluster_mean', 'cloze_cat', match='subject', ds=ds, test=False)
sensor cluster based ttest

Similarly, a mask consisting of a cluster of sensors can be used to visualize the time course in that region of interest. A straight forward choice is to use all sensors that were part of the cluster (mask) at any point in time:

roi = mask.any('time')
p = plot.Topomap(roi, cmap='Wistia')
sensor cluster based ttest

When using a mask that ony contains a sensor dimension (roi), NDVar.mean() collapses across sensors and returns a value for each time point, i.e. the time course in sensors involved in the cluster:

ds['cluster_timecourse'] = ds['eeg'].mean(roi)
p = plot.UTSStat('cluster_timecourse', 'cloze_cat', match='subject', ds=ds, frame='t')
# mark the duration of the spatio-temporal cluster
p.set_clusters(clusters, y=0.25e-6)
sensor cluster based ttest

Now visualize the cluster topography, marking significant sensors:

time_window = (clusters[0, 'tstart'], clusters[0, 'tstop'])
c1_topo = res.c1_mean.mean(time=time_window)
c0_topo = res.c0_mean.mean(time=time_window)
diff_topo = res.difference.mean(time=time_window)
p = plot.Topomap([c1_topo, c0_topo, diff_topo], axtitle=['Low cloze', 'High cloze', 'Low - high'], ncol=3)
p.mark_sensors(roi, -1)
Low cloze, High cloze, Low - high
Temporal cluster based test

Alternatively, if a spatial region of interest exists, a univariate time course can be extracted and submitted to a temporal cluster based test. For example, the N400 is typically expected to be strong at sensor Cz:

ds['eeg_cz'] = ds['eeg'].sub(sensor='Cz')
res_timecoure = testnd.TTestRelated(
    'eeg_cz', 'cloze_cat', 'low', 'high', match='subject', ds=ds,
    pmin=0.05,  # Use uncorrected p = 0.05 as threshold for forming clusters
    tstart=0.100,  # Find clusters in the time window from 100 ...
    tstop=0.600,  # ... to 600 ms
)
clusters = res_timecoure.find_clusters(0.05)
print(clusters)

p = plot.UTSStat('eeg_cz', 'cloze_cat', match='subject', ds=ds, frame='t')
p.set_clusters(clusters, y=0.25e-6)
sensor cluster based ttest

Out:

id   tstart   tstop   duration   v         p            sig
-----------------------------------------------------------
5    0.35     0.455   0.105      -123.59   0.00097752   ***

Total running time of the script: ( 0 minutes 9.601 seconds)

Gallery generated by Sphinx-Gallery

Note

Click here to download the full example code

Two-stage test

When trials are associated with continuous predictor variables, averaging is often a poor solution that loses part of the data. In such cases, a two-stage design can be employed that allows using the continuous predictor variable to test hypotheses at the group level. A two-stage analysis involves:

  • Stage 1: fit a regression model to each individual subject’s data

  • Stage 2: test regression coefficients at the group level

The example uses simulated data meant to vaguely resemble data from an N400 experiment, but not intended as a physiologically realistic simulation.

# sphinx_gallery_thumbnail_number = 1
from eelbrain import *

Generate simulated data: each function call to datasets.simulate_erp() generates a dataset for one subject. For each subject, a regression model is fit using cloze probability as continuous predictor variables.

lms = []
for subject in range(10):
    # generate data for one subject
    ds = datasets.simulate_erp(seed=subject)
    # Re-reference EEG data
    ds['eeg'] -= ds['eeg'].mean(sensor=['M1', 'M2'])
    # Fit stage 1 model
    lm = testnd.LM('eeg', 'cloze', ds)
    lms.append(lm)

# Collect single-subject models for group analysis
stage2 = testnd.LMGroup(lms)

The testnd.LMGroup object allows quick access to t-tests of individual regression coefficients.

res = stage2.column_ttest('cloze', pmin=0.05, tstart=0.100, tstop=0.600)
print(res.find_clusters(0.05))
p = plot.TopoButterfly(res, frame='t')
p.set_time(0.400)
sensor two stage

Out:

id   tstart   tstop   duration   n_sensors   v        p   sig
-------------------------------------------------------------
1    0.34     0.46    0.12       59          2874.3   0   ***

Since this is a regression model, it also contains an intercept coefficient reflecting signal deflections shared in all trials.

res = stage2.column_ttest('intercept', pmin=0.05)
p = plot.TopoButterfly(res, frame='t')
p.set_time(0.120)
sensor two stage

The regression coefficients themselves can be retrieved as Dataset:

coeffs = stage2.coefficients_dataset()
print(coeffs.summary())

Out:

Key         Type     Values
----------------------------------------------------------------------------------------
intercept   NDVar    140 time, 65 sensor; -4.89413e-06 - 4.97167e-06
cloze       NDVar    140 time, 65 sensor; -4.91937e-06 - 6.47554e-06
subject     Factor   S000, S001, S002, S003, S004, S005, S006, S007, S008, S009 (random)
----------------------------------------------------------------------------------------
Dataset: 10 cases

Total running time of the script: ( 0 minutes 6.306 seconds)

Gallery generated by Sphinx-Gallery

Note

Click here to download the full example code

Compare topographies

This example shows how to compare EEG topographies, based on the method described by McCarthy & Wood 1.

# sphinx_gallery_thumbnail_number = 4
from eelbrain import *
Simulated data

Generate a simulated dataset (as in the Cluster-based permutation t-test example)

dss = []
for subject in range(10):
    # generate data for one subject
    ds = datasets.simulate_erp(seed=subject)
    # average across trials to get condition means
    ds_agg = ds.aggregate('cloze_cat')
    # add the subject name as variable
    ds_agg[:, 'subject'] = f'S{subject:02}'
    dss.append(ds_agg)

ds = combine(dss)
# make subject a random factor (to treat it as random effect for ANOVA)
ds['subject'].random = True
# Re-reference the EEG data (i.e., subtract the mean of the two mastoid channels):
ds['eeg'] -= ds['eeg'].mean(sensor=['M1', 'M2'])
print(ds.head())

Out:

n    cloze     cloze_cat   n_chars   subject
--------------------------------------------
40   0.88051   high        4.125     S00
40   0.17241   low         3.825     S00
40   0.89466   high        4         S01
40   0.13778   low         4.575     S01
40   0.90215   high        4.275     S02
40   0.12206   low         3.9       S02
40   0.88503   high        4         S03
40   0.14273   low         4.175     S03
40   0.90499   high        4.1       S04
40   0.15732   low         3.5       S04
--------------------------------------------
NDVars: eeg

The simulated data in the two conditions:

p = plot.TopoArray('eeg', 'cloze_cat', ds=ds, axh=2, axw=10, t=[0.120, 0.200, 0.280])
high, low
Test between conditions

Test whether the 120 ms topography differs between the two cloze conditions. The dataset already includes one row per cell (i.e., per cloze condition and subject). Consequently, we can just index the topography at the desired time point:

topography = ds['eeg'].sub(time=0.120)
# normalize the data in accordance with McCarth & Wood (1985)
topography = normalize_in_cells(topography, 'sensor', 'cloze_cat', ds)
# "melt" the topography NDVar to turn the sensor dimension into a Factor
ds_topography = table.melt_ndvar(topography, 'sensor', ds=ds)
# Note EEG is a single column, and the last column indicates the sensor
ds_topography.head()

Out:

n    cloze     cloze_cat   n_chars   subject   eeg        sensor
----------------------------------------------------------------
40   0.88051   high        4.125     S00       -0.55174   Fp1
40   0.17241   low         3.825     S00       -0.94975   Fp1
40   0.89466   high        4         S01       -1.0303    Fp1
40   0.13778   low         4.575     S01       -1.5701    Fp1
40   0.90215   high        4.275     S02       -1.3857    Fp1
40   0.12206   low         3.9       S02       -1.5753    Fp1
40   0.88503   high        4         S03       -1.1537    Fp1
40   0.14273   low         4.175     S03       -1.6285    Fp1
40   0.90499   high        4.1       S04       -1.9244    Fp1
40   0.15732   low         3.5       S04       -1.5299    Fp1

ANOVA to test whether the effect of cloze_cat differs between sensors:

test.ANOVA('eeg', 'cloze_cat * sensor * subject', ds=ds_topography)

Out:

<ANOVA: eeg ~ cloze_cat + sensor + cloze_cat x sensor + subject + cloze_cat x subject + sensor x subject + cloze_cat x sensor x subject
                            SS     df      MS   MS(denom)   df(denom)           F        p
  ----------------------------------------------------------------------------------------
  cloze_cat              -0.00      1   -0.00        3.64           9    -0.00       1.000
  sensor               1295.07     64   20.24        0.08         576   255.94***   < .001
  cloze_cat x sensor      4.93     64    0.08        0.07         576     1.06        .366
  ----------------------------------------------------------------------------------------
  Total                1468.53   1299
>

The non-significant interaction suggests that the effect of cloze_cat does not differ between sensors, i.e., the topographies do not differ, which is consistent with being generated by the same underlying neural sources.

Test two time points

Since we’re not interested in condition here, we first average across conditions, i.e., with the goal of having one row per subject:

ds_average = ds.aggregate('subject', drop_bad=True)
print(ds_average)

Out:

n    cloze     n_chars   subject
--------------------------------
40   0.52646   3.975     S00
40   0.51622   4.2875    S01
40   0.51211   4.0875    S02
40   0.51388   4.0875    S03
40   0.53115   3.8       S04
40   0.52163   4.0375    S05
40   0.53789   4.2       S06
40   0.52491   4.1125    S07
40   0.52464   4.15      S08
40   0.52559   3.8875    S09
--------------------------------
NDVars: eeg

In order to compare two time points, we need to construct a new dataset with time point as Factor:

dss = []
for time in [0.120, 0.280]:
    ds_time = ds_average['subject',]  # A new dataset with the 'subject' variable only
    ds_time['eeg'] = ds_average['eeg'].sub(time=time)
    ds_time[:, 'time'] = f'{time*1000:.0f} ms'
    dss.append(ds_time)
ds_times = combine(dss)
ds_times.summary()

Out:

Key       Type     Values
------------------------------------------------------------------------------------------------
subject   Factor   S00:2, S01:2, S02:2, S03:2, S04:2, S05:2, S06:2, S07:2, S08:2, S09:2 (random)
eeg       NDVar    65 sensor; -3.62446e-06 - 3.90635e-06
time      Factor   120 ms:10, 280 ms:10
------------------------------------------------------------------------------------------------
None: 20 cases

Then, normalize the data in accordance with McCarth & Wood (1985)

topography = normalize_in_cells('eeg', 'sensor', 'time', ds=ds_times)
# "melt" the topography NDVar to turn the sensor dimension into a Factor
ds_topography = table.melt_ndvar(topography, 'sensor', ds=ds_times)
# Note EEG is a single column, and the last column indicates the sensor
ds_topography.head()

Out:

subject   time     eeg        sensor
------------------------------------
S00       120 ms   -0.74756   Fp1
S01       120 ms   -1.2964    Fp1
S02       120 ms   -1.4811    Fp1
S03       120 ms   -1.3882    Fp1
S04       120 ms   -1.735     Fp1
S05       120 ms   -1.3842    Fp1
S06       120 ms   -1.2591    Fp1
S07       120 ms   -1.3146    Fp1
S08       120 ms   -1.4126    Fp1
S09       120 ms   -1.0279    Fp1

Plot the topographies before and after normalization:

p = plot.Topomap('eeg', 'time', ds=ds_times, ncol=2, title="Original")
p = plot.Topomap(topography, 'time', ds=ds_times, ncol=2, title="Normalized")
  • Original, 120 ms, 280 ms
  • Normalized, 120 ms, 280 ms

Compre the topographies with the ANOVA – test whether the effect of time differs between sensors:

test.ANOVA('eeg', 'time * sensor * subject', ds=ds_topography)

Out:

<ANOVA: eeg ~ time + sensor + time x sensor + subject + time x subject + sensor x subject + time x sensor x subject
                       SS     df      MS   MS(denom)   df(denom)           F        p
  -----------------------------------------------------------------------------------
  time               0.00      1    0.00       12.55           9
  sensor          1269.88     64   19.84        0.19         576   105.37***   < .001
  time x sensor     30.12     64    0.47        0.18         576     2.57***   < .001
  -----------------------------------------------------------------------------------
  Total           1754.93   1299
>

Visualize the difference

res = testnd.TTestRelated(topography, 'time', match='subject', ds=ds_times)
p = plot.Topomap(res, ncol=3, title="Normalized topography differences")
Normalized topography differences, 120 ms, 280 ms, (120 ms) - (280 ms)
References
1

McCarthy, G., & Wood, C. C. (1985). Scalp Distributions of Event-Related Potentials—An Ambiguity Associated with Analysis of Variance Models. Electroencephalography and Clinical Neurophysiology, 61, S226–S227. 10.1016/0013-4694(85)90858-2

Total running time of the script: ( 0 minutes 12.612 seconds)

Gallery generated by Sphinx-Gallery

Plotting

Examples illustrating the use of the plot functionality.

Note

Click here to download the full example code

Boxplot

Much of the functionality is shared with plot.Barplot.

# %matplotlib inline
from eelbrain import *
Basic boxplot
ds = datasets.get_uv()
ds.summary()
p = plot.Boxplot('fltvar', 'A % B', match='rm', ds=ds)
boxplot
Test against population mean

Instead of pairwise tests, test the values in each category against a population mean. Use the test parameter to set the value to test against. In the example, test=0 will test whether each sample is significantly different form zero. Use the tail parameter to determine the tailedness of the test (tail=1 to test for a value greater than 0).

p = plot.Boxplot('fltvar', 'A % B', match='rm', ds=ds, test=0, tail=1)
boxplot
Appearance: colored boxes and category labels
colors = {
    ('a1', 'b1'): (1.0, 0.0, 0.0),
    ('a1', 'b2'): (1.0, 0.7, 0.0),
    ('a2', 'b1'): (0.0, 0.0, 1.0),
    ('a2', 'b2'): (0.0, 0.8, 1.0),
}
labels = {
    ('a1', 'b1'): 'A-1 (B-1)',
    ('a1', 'b2'): 'A-1 (B-2)',
    ('a2', 'b1'): 'A-2 (B-1)',
    ('a2', 'b2'): 'A-2 (B-2)',
}
p = plot.Boxplot('fltvar', 'A % B', match='rm', ds=ds, colors=colors, xticks=labels)

colors = {
    ('a1', 'b1'): plot.Style((1, .5, .5), hatch=''),
    ('a1', 'b2'): plot.Style((1, .5, .5), hatch='///'),
    ('a2', 'b1'): plot.Style((.5, .5, 1), hatch=''),
    ('a2', 'b2'): plot.Style((.5, .5, 1), hatch='///'),
}
p = plot.Boxplot('fltvar', 'A % B', match='rm', ds=ds, colors=colors)

colors = {
    ('a1', 'b1'): plot.Style('w', hatch='///'),
    ('a1', 'b2'): plot.Style('w', hatch='/'),
    ('a2', 'b1'): plot.Style('w', hatch='O'),
    ('a2', 'b2'): plot.Style('w', hatch='.'),
}
p = plot.Boxplot('fltvar', 'A % B', match='rm', ds=ds, colors=colors)
  • boxplot
  • boxplot
  • boxplot
Label fliers

The label_fliers=True option is used to identify outlier points (labels are based on the match argument).

p = plot.Boxplot('fltvar', 'A % B', match='rm', ds=ds, label_fliers=True)
boxplot

Gallery generated by Sphinx-Gallery

Note

Click here to download the full example code

Colors

In general, colors for categorial variables are set as {cell: color} dictionaries. For factors, cells are strings, and for interactions they are tuples of strings. Colors can be set to any color understood by matplotlib.

# sphinx_gallery_thumbnail_number = 3
from eelbrain import *
import matplotlib.style

ds = datasets.get_uv()


colors = {'a1': (1, 0, 0), 'a2': (0, 0, 1)}
p = plot.Barplot('fltvar', 'A', ds=ds, w=2, colors=colors)

colors = {
    ('a1', 'b1'): 'red',
    ('a1', 'b2'): (1, 1, 0),
    ('a2', 'b1'): (0, 0, 1),
    ('a2', 'b2'): (0, 1, 1),
}
p = plot.Barplot('fltvar - 0.5', 'A % B', ds=ds, w=2, colors=colors)
  • colors
  • colors
Unambiguous colors

Eelbrain also comes with utilities to quickly generate such color dictionaries. Unambiguous colors are colors for up to 8 categories that are universally distinguishable.

# add a Factor with more categories
ds['AB'] = ds.eval('A%B').as_factor()

colors = plot.colors_for_oneway(ds['AB'].cells, unambiguous=True)
p = plot.Barplot('fltvar - 0.5', 'AB', ds=ds, w=2, xticks=False, colors=colors)
p_colors = plot.ColorList(colors)

colors = plot.colors_for_oneway(ds['AB'].cells, unambiguous=[2, 5, 3, 6])
p = plot.Barplot('fltvar - 0.5', 'AB', ds=ds, w=2, xticks=False, colors=colors)
p_colors = plot.ColorList(colors)
  • colors
  • colors
  • colors
  • colors
Hue-based colors

Colors can also be generated based on a color wheel.

cells = 'abc'
colors = plot.colors_for_oneway(cells)
p = plot.ColorList(colors)

cells = 'abc'
colors = plot.colors_for_oneway(cells, hue_start=0.3)
p = plot.ColorList(colors)

cells = 'abcdefgh'
colors = plot.colors_for_oneway(cells)
p = plot.ColorList(colors)
  • colors
  • colors
  • colors

Hue and lightness can be used to distinguish two factors in an interaction.

cells_a = ('a', 'b')
cells_b = ('1', '2')
colors = plot.colors_for_twoway(cells_a, cells_b)
p = plot.ColorGrid(cells_a, cells_b, colors, h=1, w=1)

cells_a = ('a', 'b')
cells_b = ('1', '2', '3')
colors = plot.colors_for_twoway(cells_a, cells_b)
p = plot.ColorGrid(cells_a, cells_b, colors, h=1, w=1)

cells_a = ('a', 'b', 'c')
cells_b = ('1', '2')
colors = plot.colors_for_twoway(cells_a, cells_b, hue_start=0.3)
p = plot.ColorGrid(cells_a, cells_b, colors, h=1, w=1)
  • colors
  • colors
  • colors

Gallery generated by Sphinx-Gallery

Note

Click here to download the full example code

Eelbrain Colormaps

Eelbrain adds some colormaps to matplotlib’s standard colormaps. In general, the colormaps come with an *-a (“alpha”) version which uses transparency for low values instead of a solid color. For example, to use the polar map with transparency, use polar-a.

Eelbrain colormaps 
import numpy as np
import matplotlib.pyplot as plt


gradient = np.linspace(0, 1, 256)
gradient = np.vstack((gradient, gradient))

cmap_list = [
    'polar',
    'xpolar',
    'lpolar',
    'lux',
    'polar-lux',
    'lux-purple',
    'polar-lux-purple',
]

nrows = len(cmap_list)
figh = 0.35 + 0.15 + (nrows + (nrows-1)*0.1)*0.22
fig, axs = plt.subplots(nrows=nrows, figsize=(6.4, figh))
fig.subplots_adjust(top=1-.35/figh, bottom=.15/figh, left=0.2, right=0.99)

axs[0].set_title('Eelbrain colormaps', fontsize=14)

for ax, name in zip(axs, cmap_list):
    ax.imshow(gradient, aspect='auto', cmap=plt.get_cmap(name))
    ax.text(-.01, .5, name, va='center', ha='right', fontsize=10,
            transform=ax.transAxes)

for ax in axs:
    ax.set_axis_off()

Gallery generated by Sphinx-Gallery

Note

Click here to download the full example code

Customizing plots

With the exception of plot.brain plots, Eelbrain’s plots are all based on matplotlib. A lot of fine control over the plots can be achieved through two means:

  • Customizing Matplotlib globally, before calling Eelbrain plotting functions, through styles or ``rcParams` <https://matplotlib.org/tutorials/ introductory/customizing.html#customizing-matplotlib-with-style-sheets-and-rcparams>`_

  • Accessing and modifying components of the plots after calling Eelbrain plotting functions

# sphinx_gallery_thumbnail_number = 3
from eelbrain import *
import matplotlib.style

ds = datasets.get_uv()
Styles

Matplotlib offers several styles

p = plot.Boxplot('fltvar', 'A % B', match='rm', ds=ds, w=2)

# Apply a style
matplotlib.style.use('ggplot')
p = plot.Boxplot('fltvar', 'A % B', match='rm', ds=ds, w=2)

matplotlib.style.use('bmh')
p = plot.Boxplot('fltvar', 'A % B', match='rm', ds=ds, w=2)

matplotlib.style.use('seaborn')
p = plot.Boxplot('fltvar', 'A % B', match='rm', ds=ds, w=2)
  • customizing
  • customizing
  • customizing
  • customizing
rcParams

Individual styles parameters can be modified directly in``rcParams``

matplotlib.rcParams['font.family'] = 'serif'
matplotlib.rcParams['font.size'] = 8
p = plot.Boxplot('fltvar', 'A % B', match='rm', ds=ds, w=2)

# revert back to the default style
matplotlib.style.use('default')
customizing
Modifying components
p = plot.Boxplot('fltvar', 'A % B', match='rm', ds=ds, w=2)

p = plot.Boxplot('fltvar', 'A % B', match='rm', ds=ds, w=2, xlabel=False)
ax = p.figure.axes[0]
ax.set_xticklabels(['A long label', 'An even longer label', 'Another label', 'And yet another one'], rotation=45, ha='right')
ax.grid(axis='y')
ax.set_yticks([-2, 0, 2])
ax.tick_params('y', left=False)
for spine in ax.spines.values():
    spine.set_visible(False)
  • customizing
  • customizing

Gallery generated by Sphinx-Gallery

Reverese Correlation

Note

Click here to download the full example code

Introduction to Temporal Response Functions (TRFs)

A temporal response functions (TRFs) is a linear model of the stimulus-response depency. The response is predicted by a linear convolution of the stimulus with the TRF. This means that for every non-zero element in the stimulus, there will be a response in the shape of the TRF:

from eelbrain import *
import numpy as np

# Construct a 10 s long stimulus
time = UTS(0, 0.01, 1000)
x = NDVar(np.zeros(len(time)), time)
# add a few impulses
x[1] = 1
x[3] = 1
x[5] = 1

# Construct a TRF of length 500 ms
trf_time = UTS(0, 0.01, 50)
trf = gaussian(0.200, 0.050, trf_time) - gaussian(0.300, 0.050, trf_time)

# The response is the convolution of the stimulus with the TRF
y = convolve(trf, x)

plot_args = dict(ncol=1, axh=2, w=10, frame='t', legend=False, colors='r')
plot.UTS([x, trf, y], ylabel=['Stimulus (x)', 'TRF', 'Response (y)'], **plot_args)
trf intro

Out:

<UTS: None>

Since The convolution is linear,

  • scaled stimuli cause scaled responses

  • overlapping responses add up

x[2] = 0.66
x[5.2] = 1
x[7] = 0.8
x[7.2] = 0.3
x[8.8] = 0.2
x[9] = 0.5
x[9.2] = 0.7

y = convolve(trf, x)
plot.UTS([x, trf, y], ylabel=['Stimulus (x)', 'TRF', 'Response (y)'], **plot_args)
trf intro

Out:

<UTS: None>

When the stimulus contains only non-zero elements this works just the same, but the TRF might not be apparent in the response anymore:

x += np.random.normal(0, 0.1, x.shape)
filter_data(x, 1, 40)
y = convolve(trf, x)
plot.UTS([x, trf, y], ylabel=['Stimulus (x)', 'TRF', 'Response (y)'], **plot_args)
trf intro

Out:

<UTS: None>

Given a stimulus and a response, there are different methods to reconstruct the TRF. Eelbrain comes with an implementation of the boosting() coordinate descent algorithm:

fit = boosting(y, x, 0.000, 0.500, basis=0.050, partitions=2)
plot_args = {**plot_args, 'ncol': 3}
plot.UTS([trf, fit.h, None], axtitle=["Model TRF", "Reconstructed TRF"], **plot_args)
Model TRF, Reconstructed TRF

Out:

Fitting models:   0%|          | 0/2 [00:00<?, ?it/s]

<UTS: None>

Gallery generated by Sphinx-Gallery

Note

Click here to download the full example code

EEG speech envelope TRF

Analyze continuous speech data from the mTRF dataset 1: use the boosting algorithm for estimating temporal response functions (TRFs) to the acoustic envelope.

# Author: Christian Brodbeck <christianbrodbeck@nyu.edu>
# sphinx_gallery_thumbnail_number = 4
import os

from scipy.io import loadmat
import mne
from eelbrain import *

# Load the mTRF speech dataset and convert data to NDVars
root = mne.datasets.mtrf.data_path()
speech_path = os.path.join(root, 'speech_data.mat')
mdata = loadmat(speech_path)

# Time axis
tstep = 1. / mdata['Fs'][0, 0]
n_times = mdata['envelope'].shape[0]
time = UTS(0, tstep, n_times)
# Load the EEG sensor coordinates (drop fiducials coordinates, which are stored
# after sensor 128)
sensor = Sensor.from_montage('biosemi128')[:128]
# Frequency dimension for the spectrogram
band = Scalar('frequency', range(16))
# Create variables
envelope = NDVar(mdata['envelope'][:, 0], (time,), name='envelope')
eeg = NDVar(mdata['EEG'], (time, sensor), name='EEG', info={'unit': 'µV'})
spectrogram = NDVar(mdata['spectrogram'], (time, band), name='spectrogram')
# Exclude a bad channel
eeg = eeg[sensor.index(exclude='A13')]
Data

Plot the spectrogram of the speech stimulus:

plot.Array(spectrogram, xlim=5, w=6, h=2)
mtrf

Out:

<Array: spectrogram>

Plot the envelope used as stimulus representation for reverse correlation:

plot.UTS(envelope, xlim=5, w=6, h=2)
mtrf

Out:

<UTS: envelope>

Plot the corresponding EEG data:

p = plot.TopoButterfly(eeg, xlim=5, w=7, h=2)
p.set_time(1.200)
mtrf
Reverse correlation

TRF for the envelope using boosting:

  • TRF from -100 to 400 ms

  • Basis of 100 ms Hamming windows

  • Use 4 partitionings of the data for cross-validation based early stopping

res = boosting(eeg, envelope, -0.100, 0.400, basis=0.100, partitions=4)
p = plot.TopoButterfly(res.h_scaled, w=6, h=2)
p.set_time(.180)
mtrf

Out:

Fitting models:   0%|          | 0/508 [00:00<?, ?it/s]
Fitting models:   0%|          | 1/508 [00:00<02:14,  3.76it/s]
Fitting models:   2%|1         | 9/508 [00:00<00:16, 30.04it/s]
Fitting models:   3%|3         | 17/508 [00:00<00:11, 44.04it/s]
Fitting models:   5%|4         | 25/508 [00:00<00:09, 52.64it/s]
Fitting models:   6%|6         | 33/508 [00:00<00:08, 56.48it/s]
Fitting models:   8%|7         | 40/508 [00:00<00:08, 57.68it/s]
Fitting models:   9%|9         | 48/508 [00:00<00:07, 62.12it/s]
Fitting models:  11%|#         | 55/508 [00:01<00:07, 63.68it/s]
Fitting models:  12%|#2        | 62/508 [00:01<00:07, 60.29it/s]
Fitting models:  14%|#3        | 69/508 [00:01<00:07, 59.07it/s]
Fitting models:  15%|#5        | 77/508 [00:01<00:07, 59.83it/s]
Fitting models:  17%|#6        | 84/508 [00:01<00:06, 61.35it/s]
Fitting models:  18%|#8        | 92/508 [00:01<00:06, 63.71it/s]
Fitting models:  20%|#9        | 100/508 [00:01<00:06, 66.67it/s]
Fitting models:  21%|##1       | 107/508 [00:01<00:06, 66.01it/s]
Fitting models:  22%|##2       | 114/508 [00:02<00:06, 60.86it/s]
Fitting models:  24%|##3       | 121/508 [00:02<00:06, 56.22it/s]
Fitting models:  25%|##5       | 129/508 [00:02<00:06, 59.31it/s]
Fitting models:  27%|##6       | 137/508 [00:02<00:06, 61.19it/s]
Fitting models:  28%|##8       | 144/508 [00:02<00:06, 60.22it/s]
Fitting models:  30%|##9       | 151/508 [00:02<00:06, 58.78it/s]
Fitting models:  31%|###1      | 159/508 [00:02<00:05, 62.97it/s]
Fitting models:  33%|###2      | 167/508 [00:02<00:05, 64.67it/s]
Fitting models:  34%|###4      | 175/508 [00:02<00:05, 65.83it/s]
Fitting models:  36%|###6      | 183/508 [00:03<00:04, 68.05it/s]
Fitting models:  38%|###7      | 191/508 [00:03<00:04, 69.70it/s]
Fitting models:  39%|###9      | 199/508 [00:03<00:04, 69.49it/s]
Fitting models:  41%|####      | 207/508 [00:03<00:04, 69.89it/s]
Fitting models:  42%|####2     | 215/508 [00:03<00:04, 69.88it/s]
Fitting models:  44%|####3     | 223/508 [00:03<00:04, 71.21it/s]
Fitting models:  45%|####5     | 231/508 [00:03<00:03, 70.24it/s]
Fitting models:  47%|####7     | 239/508 [00:03<00:03, 68.37it/s]
Fitting models:  48%|####8     | 246/508 [00:04<00:03, 65.56it/s]
Fitting models:  50%|#####     | 254/508 [00:04<00:03, 67.66it/s]
Fitting models:  52%|#####1    | 262/508 [00:04<00:03, 68.71it/s]
Fitting models:  53%|#####2    | 269/508 [00:04<00:03, 68.13it/s]
Fitting models:  54%|#####4    | 276/508 [00:04<00:03, 64.12it/s]
Fitting models:  56%|#####5    | 283/508 [00:04<00:03, 64.07it/s]
Fitting models:  57%|#####7    | 291/508 [00:04<00:03, 65.07it/s]
Fitting models:  59%|#####8    | 299/508 [00:04<00:03, 68.29it/s]
Fitting models:  60%|######    | 306/508 [00:04<00:03, 66.94it/s]
Fitting models:  62%|######1   | 313/508 [00:05<00:03, 64.56it/s]
Fitting models:  63%|######2   | 320/508 [00:05<00:03, 55.11it/s]
Fitting models:  64%|######4   | 326/508 [00:05<00:03, 51.21it/s]
Fitting models:  65%|######5   | 332/508 [00:05<00:03, 49.35it/s]
Fitting models:  67%|######6   | 339/508 [00:05<00:03, 54.23it/s]
Fitting models:  68%|######8   | 346/508 [00:05<00:02, 57.06it/s]
Fitting models:  69%|######9   | 353/508 [00:05<00:02, 59.23it/s]
Fitting models:  71%|#######   | 360/508 [00:05<00:02, 55.32it/s]
Fitting models:  72%|#######2  | 366/508 [00:06<00:02, 55.16it/s]
Fitting models:  73%|#######3  | 373/508 [00:06<00:02, 56.83it/s]
Fitting models:  75%|#######4  | 379/508 [00:06<00:02, 51.54it/s]
Fitting models:  76%|#######5  | 386/508 [00:06<00:02, 54.10it/s]
Fitting models:  78%|#######7  | 394/508 [00:06<00:01, 57.69it/s]
Fitting models:  79%|#######8  | 400/508 [00:06<00:01, 54.05it/s]
Fitting models:  80%|#######9  | 406/508 [00:06<00:01, 54.78it/s]
Fitting models:  81%|########1 | 412/508 [00:06<00:01, 48.17it/s]
Fitting models:  82%|########2 | 417/508 [00:07<00:01, 47.74it/s]
Fitting models:  83%|########3 | 422/508 [00:07<00:01, 44.97it/s]
Fitting models:  84%|########4 | 428/508 [00:07<00:01, 47.28it/s]
Fitting models:  85%|########5 | 434/508 [00:07<00:01, 50.58it/s]
Fitting models:  87%|########6 | 440/508 [00:07<00:01, 51.74it/s]
Fitting models:  88%|########7 | 447/508 [00:07<00:01, 55.25it/s]
Fitting models:  89%|########9 | 454/508 [00:07<00:00, 55.08it/s]
Fitting models:  91%|######### | 460/508 [00:07<00:00, 51.28it/s]
Fitting models:  92%|#########1| 466/508 [00:08<00:00, 50.42it/s]
Fitting models:  93%|#########2| 472/508 [00:08<00:00, 49.11it/s]
Fitting models:  94%|#########4| 479/508 [00:08<00:00, 53.07it/s]
Fitting models:  95%|#########5| 485/508 [00:08<00:00, 51.51it/s]
Fitting models:  97%|#########6| 491/508 [00:08<00:00, 50.19it/s]
Fitting models:  98%|#########8| 499/508 [00:08<00:00, 56.13it/s]
Fitting models:  99%|#########9| 505/508 [00:08<00:00, 52.95it/s]
Multiple predictors

Multiple predictors additively explain the signal:

# Derive acoustic onsets from the envelope
onset = envelope.diff('time', name='onset').clip(0)
onset *= envelope.max() / onset.max()
plot.UTS([[envelope, onset]], xlim=5, w=6, h=2)
mtrf

Out:

<UTS: envelope, onset>

res_onset = boosting(eeg, [onset, envelope], -0.100, 0.400, basis=0.100, partitions=4)
p = plot.TopoButterfly(res_onset.h_scaled, w=6, h=3)
p.set_time(.150)
mtrf

Out:

Fitting models:   0%|          | 0/508 [00:00<?, ?it/s]
Fitting models:   0%|          | 1/508 [00:00<02:25,  3.48it/s]
Fitting models:   1%|1         | 7/508 [00:00<00:25, 19.58it/s]
Fitting models:   2%|2         | 11/508 [00:00<00:20, 24.69it/s]
Fitting models:   3%|2         | 15/508 [00:00<00:17, 28.90it/s]
Fitting models:   4%|3         | 19/508 [00:00<00:15, 30.68it/s]
Fitting models:   5%|4         | 23/508 [00:00<00:15, 30.70it/s]
Fitting models:   5%|5         | 27/508 [00:01<00:15, 30.64it/s]
Fitting models:   6%|6         | 31/508 [00:01<00:14, 32.11it/s]
Fitting models:   7%|6         | 35/508 [00:01<00:14, 32.56it/s]
Fitting models:   8%|7         | 39/508 [00:01<00:14, 33.31it/s]
Fitting models:   8%|8         | 43/508 [00:01<00:13, 34.37it/s]
Fitting models:   9%|9         | 47/508 [00:01<00:12, 35.67it/s]
Fitting models:  10%|#         | 51/508 [00:01<00:12, 36.63it/s]
Fitting models:  11%|#         | 55/508 [00:01<00:12, 36.36it/s]
Fitting models:  12%|#1        | 59/508 [00:01<00:13, 33.26it/s]
Fitting models:  12%|#2        | 63/508 [00:02<00:13, 33.16it/s]
Fitting models:  13%|#3        | 67/508 [00:02<00:13, 33.50it/s]
Fitting models:  14%|#3        | 71/508 [00:02<00:13, 33.00it/s]
Fitting models:  15%|#4        | 75/508 [00:02<00:12, 33.73it/s]
Fitting models:  16%|#5        | 79/508 [00:02<00:13, 32.47it/s]
Fitting models:  16%|#6        | 83/508 [00:02<00:12, 34.14it/s]
Fitting models:  17%|#7        | 87/508 [00:02<00:12, 34.44it/s]
Fitting models:  18%|#7        | 91/508 [00:02<00:11, 35.63it/s]
Fitting models:  19%|#8        | 95/508 [00:02<00:11, 35.71it/s]
Fitting models:  19%|#9        | 99/508 [00:03<00:11, 35.81it/s]
Fitting models:  20%|##        | 103/508 [00:03<00:11, 34.98it/s]
Fitting models:  21%|##1       | 107/508 [00:03<00:11, 33.56it/s]
Fitting models:  22%|##1       | 111/508 [00:03<00:12, 32.77it/s]
Fitting models:  23%|##2       | 115/508 [00:03<00:11, 32.86it/s]
Fitting models:  23%|##3       | 119/508 [00:03<00:11, 33.91it/s]
Fitting models:  24%|##4       | 123/508 [00:03<00:11, 34.61it/s]
Fitting models:  25%|##5       | 127/508 [00:03<00:10, 35.11it/s]
Fitting models:  26%|##5       | 131/508 [00:04<00:10, 35.84it/s]
Fitting models:  27%|##6       | 135/508 [00:04<00:10, 34.31it/s]
Fitting models:  27%|##7       | 139/508 [00:04<00:10, 34.10it/s]
Fitting models:  28%|##8       | 143/508 [00:04<00:11, 30.81it/s]
Fitting models:  29%|##8       | 147/508 [00:04<00:11, 31.23it/s]
Fitting models:  30%|##9       | 151/508 [00:04<00:11, 30.97it/s]
Fitting models:  31%|###       | 155/508 [00:04<00:10, 32.94it/s]
Fitting models:  31%|###1      | 159/508 [00:04<00:10, 33.97it/s]
Fitting models:  32%|###2      | 163/508 [00:05<00:09, 34.78it/s]
Fitting models:  33%|###2      | 167/508 [00:05<00:09, 35.06it/s]
Fitting models:  34%|###3      | 171/508 [00:05<00:09, 35.67it/s]
Fitting models:  34%|###4      | 175/508 [00:05<00:09, 35.15it/s]
Fitting models:  35%|###5      | 179/508 [00:05<00:09, 35.27it/s]
Fitting models:  36%|###6      | 183/508 [00:05<00:09, 35.45it/s]
Fitting models:  37%|###6      | 187/508 [00:05<00:08, 36.11it/s]
Fitting models:  38%|###7      | 191/508 [00:05<00:08, 35.96it/s]
Fitting models:  38%|###8      | 195/508 [00:05<00:08, 35.74it/s]
Fitting models:  39%|###9      | 199/508 [00:06<00:09, 33.13it/s]
Fitting models:  40%|###9      | 203/508 [00:06<00:09, 32.39it/s]
Fitting models:  41%|####      | 207/508 [00:06<00:09, 32.98it/s]
Fitting models:  42%|####1     | 211/508 [00:06<00:08, 34.17it/s]
Fitting models:  42%|####2     | 215/508 [00:06<00:08, 35.03it/s]
Fitting models:  43%|####3     | 219/508 [00:06<00:07, 36.13it/s]
Fitting models:  44%|####3     | 223/508 [00:06<00:07, 36.50it/s]
Fitting models:  45%|####4     | 227/508 [00:06<00:07, 37.12it/s]
Fitting models:  45%|####5     | 231/508 [00:06<00:07, 36.89it/s]
Fitting models:  46%|####6     | 235/508 [00:07<00:07, 34.34it/s]
Fitting models:  47%|####7     | 239/508 [00:07<00:07, 35.70it/s]
Fitting models:  48%|####7     | 243/508 [00:07<00:07, 35.31it/s]
Fitting models:  49%|####8     | 247/508 [00:07<00:07, 35.73it/s]
Fitting models:  49%|####9     | 251/508 [00:07<00:07, 34.37it/s]
Fitting models:  50%|#####     | 255/508 [00:07<00:07, 34.80it/s]
Fitting models:  51%|#####     | 259/508 [00:07<00:07, 34.99it/s]
Fitting models:  52%|#####1    | 263/508 [00:07<00:06, 35.22it/s]
Fitting models:  53%|#####2    | 267/508 [00:07<00:06, 36.08it/s]
Fitting models:  53%|#####3    | 271/508 [00:08<00:06, 34.89it/s]
Fitting models:  54%|#####4    | 275/508 [00:08<00:07, 33.28it/s]
Fitting models:  55%|#####4    | 279/508 [00:08<00:07, 32.67it/s]
Fitting models:  56%|#####5    | 283/508 [00:08<00:06, 33.11it/s]
Fitting models:  56%|#####6    | 287/508 [00:08<00:06, 33.11it/s]
Fitting models:  57%|#####7    | 291/508 [00:08<00:06, 34.18it/s]
Fitting models:  58%|#####8    | 295/508 [00:08<00:06, 34.44it/s]
Fitting models:  59%|#####8    | 299/508 [00:08<00:05, 35.06it/s]
Fitting models:  60%|#####9    | 303/508 [00:09<00:06, 32.98it/s]
Fitting models:  60%|######    | 307/508 [00:09<00:06, 33.08it/s]
Fitting models:  61%|######1   | 311/508 [00:09<00:05, 34.26it/s]
Fitting models:  62%|######2   | 315/508 [00:09<00:05, 33.82it/s]
Fitting models:  63%|######2   | 319/508 [00:09<00:05, 33.15it/s]
Fitting models:  64%|######3   | 323/508 [00:09<00:05, 31.38it/s]
Fitting models:  64%|######4   | 327/508 [00:09<00:05, 30.76it/s]
Fitting models:  65%|######5   | 331/508 [00:09<00:05, 30.61it/s]
Fitting models:  66%|######5   | 335/508 [00:10<00:05, 32.40it/s]
Fitting models:  67%|######6   | 339/508 [00:10<00:05, 33.11it/s]
Fitting models:  68%|######7   | 343/508 [00:10<00:05, 31.23it/s]
Fitting models:  68%|######8   | 347/508 [00:10<00:05, 29.60it/s]
Fitting models:  69%|######9   | 351/508 [00:10<00:04, 31.87it/s]
Fitting models:  70%|######9   | 355/508 [00:10<00:04, 31.26it/s]
Fitting models:  71%|#######   | 359/508 [00:10<00:04, 32.42it/s]
Fitting models:  71%|#######1  | 363/508 [00:10<00:04, 32.56it/s]
Fitting models:  72%|#######2  | 367/508 [00:11<00:04, 28.47it/s]
Fitting models:  73%|#######3  | 371/508 [00:11<00:04, 29.98it/s]
Fitting models:  74%|#######3  | 375/508 [00:11<00:04, 29.47it/s]
Fitting models:  75%|#######4  | 379/508 [00:11<00:04, 27.28it/s]
Fitting models:  75%|#######5  | 383/508 [00:11<00:04, 29.13it/s]
Fitting models:  76%|#######6  | 387/508 [00:11<00:04, 28.41it/s]
Fitting models:  77%|#######6  | 391/508 [00:11<00:04, 26.05it/s]
Fitting models:  78%|#######7  | 395/508 [00:12<00:04, 25.04it/s]
Fitting models:  78%|#######8  | 398/508 [00:12<00:04, 25.67it/s]
Fitting models:  79%|#######8  | 401/508 [00:12<00:04, 23.86it/s]
Fitting models:  80%|#######9  | 404/508 [00:12<00:04, 21.32it/s]
Fitting models:  80%|########  | 407/508 [00:12<00:04, 22.53it/s]
Fitting models:  81%|########  | 410/508 [00:12<00:04, 24.06it/s]
Fitting models:  81%|########1 | 413/508 [00:13<00:04, 20.78it/s]
Fitting models:  82%|########1 | 416/508 [00:13<00:04, 20.64it/s]
Fitting models:  82%|########2 | 419/508 [00:13<00:04, 21.11it/s]
Fitting models:  83%|########3 | 422/508 [00:13<00:03, 22.35it/s]
Fitting models:  84%|########3 | 425/508 [00:13<00:03, 22.79it/s]
Fitting models:  84%|########4 | 428/508 [00:13<00:03, 20.02it/s]
Fitting models:  85%|########4 | 431/508 [00:13<00:03, 22.14it/s]
Fitting models:  86%|########5 | 435/508 [00:13<00:02, 24.87it/s]
Fitting models:  86%|########6 | 439/508 [00:14<00:02, 27.70it/s]
Fitting models:  87%|########7 | 443/508 [00:14<00:02, 29.71it/s]
Fitting models:  88%|########7 | 447/508 [00:14<00:01, 30.84it/s]
Fitting models:  89%|########8 | 451/508 [00:14<00:01, 30.13it/s]
Fitting models:  90%|########9 | 455/508 [00:14<00:01, 30.89it/s]
Fitting models:  90%|######### | 459/508 [00:14<00:01, 27.88it/s]
Fitting models:  91%|#########1| 463/508 [00:14<00:01, 29.52it/s]
Fitting models:  92%|#########1| 467/508 [00:15<00:01, 28.42it/s]
Fitting models:  93%|#########2| 470/508 [00:15<00:01, 26.00it/s]
Fitting models:  93%|#########3| 474/508 [00:15<00:01, 27.93it/s]
Fitting models:  94%|#########3| 477/508 [00:15<00:01, 27.86it/s]
Fitting models:  94%|#########4| 480/508 [00:15<00:00, 28.22it/s]
Fitting models:  95%|#########5| 483/508 [00:15<00:00, 28.18it/s]
Fitting models:  96%|#########5| 486/508 [00:15<00:00, 28.66it/s]
Fitting models:  96%|#########6| 489/508 [00:15<00:00, 26.38it/s]
Fitting models:  97%|#########7| 493/508 [00:15<00:00, 27.09it/s]
Fitting models:  98%|#########7| 497/508 [00:16<00:00, 29.67it/s]
Fitting models:  99%|#########8| 501/508 [00:16<00:00, 31.36it/s]
Fitting models:  99%|#########9| 505/508 [00:16<00:00, 28.85it/s]
Compare models

Compare model quality through the correlation between measured and predicted responses:

plot.Topomap([res.r, res_onset.r], w=4, h=2, ncol=2, axtitle=['envelope', 'envelope + onset'])
envelope, envelope + onset

Out:

<Topomap: Correlation>
References
1

Crosse, M. J., Liberto, D., M, G., Bednar, A., & Lalor, E. C. (2016). The Multivariate Temporal Response Function (mTRF) Toolbox: A MATLAB Toolbox for Relating Neural Signals to Continuous Stimuli. Frontiers in Human Neuroscience, 10. https://doi.org/10.3389/fnhum.2016.00604

Total running time of the script: ( 0 minutes 30.834 seconds)

Gallery generated by Sphinx-Gallery

Note

Click here to download the full example code

Impulse predictors for epochs

epoch_impulse_predictor() generates predictor variables for reverse correlation in trial-based experiments with discrete events. The function generates one impulse per trial, and these impulsescan be of varaiable magnitude and have variable latency.

The example uses simulated data meant to vaguely resemble data from an N400 experiment, but not intended as a physiologically realistic simulation.

# sphinx_gallery_thumbnail_number = 2
from eelbrain import *

ds = datasets.simulate_erp()
print(ds.summary())

Out:

Key         Type     Values
-------------------------------------------------------------------
eeg         NDVar    140 time, 65 sensor; -1.4998e-05 - 1.64215e-05
cloze       Var      0.00563694 - 0.997675
cloze_cat   Factor   high:40, low:40
n_chars     Var      0:2, 1:5, 2:6, 3:13, 4:23, 5:20, 6:7, 7:4
-------------------------------------------------------------------
Dataset: 80 cases
Discrete events

Computing a TRF for an impulse at trial onset is very similar to averaging:

any_trial = epoch_impulse_predictor('eeg', 1, ds=ds)
fit = boosting('eeg', any_trial, -0.100, 0.600, basis=0.050, ds=ds, partitions=2, delta=0.05)
average = ds['eeg'].mean('case')
trf = fit.h.sub(time=(average.time.tmin, average.time.tstop))
p = plot.TopoButterfly([trf, average], axtitle=['Impulse response', 'Average'])
p.set_time(0.400)
epoch impulse

Out:

Fitting models:   0%|          | 0/130 [00:00<?, ?it/s]
Fitting models:   1%|          | 1/130 [00:00<00:24,  5.35it/s]
Fitting models:   7%|6         | 9/130 [00:00<00:03, 35.70it/s]
Fitting models:  13%|#3        | 17/130 [00:00<00:02, 48.53it/s]
Fitting models:  19%|#9        | 25/130 [00:00<00:01, 54.84it/s]
Fitting models:  25%|##4       | 32/130 [00:00<00:01, 59.36it/s]
Fitting models:  31%|###       | 40/130 [00:00<00:01, 64.37it/s]
Fitting models:  37%|###6      | 48/130 [00:00<00:01, 66.97it/s]
Fitting models:  43%|####3     | 56/130 [00:00<00:01, 70.65it/s]
Fitting models:  49%|####9     | 64/130 [00:01<00:00, 71.70it/s]
Fitting models:  55%|#####5    | 72/130 [00:01<00:00, 73.52it/s]
Fitting models:  62%|######1   | 80/130 [00:01<00:00, 72.62it/s]
Fitting models:  68%|######8   | 89/130 [00:01<00:00, 74.28it/s]
Fitting models:  75%|#######4  | 97/130 [00:01<00:00, 74.05it/s]
Fitting models:  82%|########1 | 106/130 [00:01<00:00, 76.18it/s]
Fitting models:  88%|########7 | 114/130 [00:01<00:00, 73.18it/s]
Fitting models:  94%|#########3| 122/130 [00:01<00:00, 71.88it/s]
Categorial coding

Impulse predictors can be used like dummy codes in a regression model. Use one impulse to code for occurrence of any word (any_word), and a second impulse to code for unpredictable words only (cloze):

any_word = epoch_impulse_predictor('eeg', 1, ds=ds, name='any_word')
# effect code for cloze (1 for low cloze, -1 for high cloze)
cloze_code = Var.from_dict(ds['cloze_cat'], {'high': 0, 'low': 1})
low_cloze = epoch_impulse_predictor('eeg', cloze_code, ds=ds, name='low_cloze')

# plot the predictors for each trial
plot.UTS([any_word, low_cloze], '.case')
0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79

Out:

<UTS: any_word, low_cloze ~ .case>

Estimate response functions for these two predictors. Based on the coding, any_word reflects the response to predictable words, and low_cloze reflects how unpredictable words differ from predictable words:

fit = boosting('eeg', [any_word, low_cloze], 0, 0.5, basis=0.050, model='cloze_cat', ds=ds, partitions=2, delta=0.05)
p = plot.TopoButterfly(fit.h)
p.set_time(0.400)
epoch impulse

Out:

Fitting models:   0%|          | 0/130 [00:00<?, ?it/s]
Fitting models:   1%|          | 1/130 [00:00<00:20,  6.40it/s]
Fitting models:   5%|5         | 7/130 [00:00<00:04, 28.79it/s]
Fitting models:  10%|#         | 13/130 [00:00<00:03, 37.44it/s]
Fitting models:  15%|#4        | 19/130 [00:00<00:02, 42.69it/s]
Fitting models:  18%|#8        | 24/130 [00:00<00:02, 43.10it/s]
Fitting models:  23%|##3       | 30/130 [00:00<00:02, 44.42it/s]
Fitting models:  28%|##7       | 36/130 [00:00<00:02, 45.61it/s]
Fitting models:  32%|###2      | 42/130 [00:01<00:01, 46.59it/s]
Fitting models:  37%|###6      | 48/130 [00:01<00:01, 47.88it/s]
Fitting models:  42%|####1     | 54/130 [00:01<00:01, 47.43it/s]
Fitting models:  46%|####6     | 60/130 [00:01<00:01, 47.29it/s]
Fitting models:  50%|#####     | 65/130 [00:01<00:01, 47.97it/s]
Fitting models:  55%|#####4    | 71/130 [00:01<00:01, 48.99it/s]
Fitting models:  59%|#####9    | 77/130 [00:01<00:01, 49.37it/s]
Fitting models:  64%|######3   | 83/130 [00:01<00:00, 49.06it/s]
Fitting models:  68%|######8   | 89/130 [00:01<00:00, 48.08it/s]
Fitting models:  73%|#######3  | 95/130 [00:02<00:00, 48.79it/s]
Fitting models:  78%|#######7  | 101/130 [00:02<00:00, 49.57it/s]
Fitting models:  82%|########1 | 106/130 [00:02<00:00, 48.91it/s]
Fitting models:  86%|########6 | 112/130 [00:02<00:00, 49.49it/s]
Fitting models:  91%|######### | 118/130 [00:02<00:00, 49.97it/s]
Fitting models:  95%|#########5| 124/130 [00:02<00:00, 50.20it/s]
Fitting models: 100%|##########| 130/130 [00:02<00:00, 50.36it/s]

Total running time of the script: ( 0 minutes 13.319 seconds)

Gallery generated by Sphinx-Gallery

Gallery generated by Sphinx-Gallery

Recipes

Plots for Publication

See also

In order to produce multiple plots with consistent style it is useful to set some matplotlib options globally One way this can be done is by updating rcParams inside a script, e.g.:

import matplotlib as mpl

mpl.rcParams['font.family'] = 'sans-serif'
mpl.rcParams['font.size'] = 8
# Change all line-widths
for key in mpl.rcParams:
    if 'linewidth' in key:
        mpl.rcParams[key] *= 0.5
# The size of saved figures depends on the figure size (w, h) and DPI
mpl.rcParams['figure.dpi'] = 300
mpl.rcParams['savefig.dpi'] = 300

Matplotlib’s tight_layout() functionality provides an easy way for plots to use the available space, and most Eelbrain plots use it by default. However, when trying to produce multiple plots with identical scaling it can lead to unwanted discrepancies. Subplot placement can be specified in absolute scale through the margins parameter (see General layout parameters).

If a script produces several plots and there is no need to show them on the screen, they can be kept hidden through the show=False argument:

p = plot.UTSStat('uts', 'A', ds=ds, w=5, tight=False, show=False)

When working on more complex figures, it is often desirable to save the legend separately and combine it in a layout application:

p = plot.UTSStat('uts', 'A', ds=ds, w=5, tight=False, show=False, legend=False)
p.save('plot.pdf', transparent=True)
p.save_legend('legend.pdf', transparent=True)

The MneExperiment Pipeline

See also

Introduction

The MneExperiment class is a template for an MEG/EEG analysis pipeline. The pipeline is adapted to a specific experiment by creating a subclass, and specifying properties of the experiment as attributes.

Once set up, an MneExperiment subclass instance provides access into the pipeline at different stages of analysis through its methods:

  • .load_... methods are for loading data.

  • .make_... methods are for generating various intermediate results. Most of these methods don’t have to be called by the user, but they are used internally when needed. The exception are those that require user input, like ICA component selection, which are mentioned below.

  • .show_... methods are for retrieving and displaying information at different stages.

  • .plot_... methods are for generating plots of the data.

An MneExperiment instance has a state, which determines what data and settings it is currently using. Not all settings are always relevant. For example, subject is relevenat for steps applied separately to each subject, like make_ica_selection(), whereas group defines the group of subjects in group level analysis, such as in load_test(). For more information, see State Parameters.

Step by step

Setting up the file structure
MneExperiment.sessions

The first step is to define an MneExperiment subclass with the name of the experiment:

from eelbrain import *

class WordExperiment(MneExperiment):

    sessions = 'words'

Where sessions is the name which you included in your raw data files after the subject identifier.

The pipeline expects input files in a strictly determined folder/file structure. In the schema below, curly brackets indicate slots to be replaced with specific names, for example '{subject}' should be replaced with each specific subject’s label:

root
mri-sdir                                /mri
mri-dir                                    /{mrisubject}
meg-sdir                                /meg
meg-dir                                    /{subject}
raw-dir
trans-file                                       /{mrisubject}-trans.fif
raw-file                                         /{subject}_{session}-raw.fif

This schema shows path templates according to which the input files should be organized. Assuming that root="/files", for a subject called “R0001” this includes:

  • MRI-directory at /files/mri/R0001

  • the raw data file at /files/meg/R0001/R0001_words-raw.fif (the session is called “words” which is specified in WordExperiment.sessions)

  • the trans-file from the coregistration at /files/meg/R0001/R0001-trans.fif

Once the required files are placed in this structure, the experiment class can be initialized with the proper root parameter, pointing to where the files are located:

>>> e = WordExperiment("/files")

The setup can be tested using MneExperiment.show_subjects(), which shows a list of the subjects that were discovered and the MRIs used:

>>> e.show_subjects()
#    subject   mri
-----------------------------------------
0    R0026     R0026
1    R0040     fsaverage * 0.92
2    R0176     fsaverage * 0.954746600461
...
MneExperiment.visits

Note

If participants come back for the experiment on multiple occasions, a visits attribute might also be needed. For details see the corresponding wiki page.

Pre-processing

Make sure an appropriate pre-processing pipeline is defined as MneExperiment.raw.

To inspect raw data for a given pre-processing stage use:

>>> e.set(raw='1-40')
>>> y = e.load_raw(ndvar=True)
>>> p = plot.TopoButterfly(y, xlim=10, w=0)

Which will plot a 10 s excerpt and allow scrolling through the data.

Events

If needed, set MneExperiment.merge_triggers to handle spurious events. Then, add event labels. Initially, events are only labeled with the trigger ID. Use the MneExperiment.variables settings to add labels. Events are represented as Dataset objects and can be inspected with corresponding methods and functions, for example:

>>> e = WordExperiment("/files")
>>> ds = e.load_events()
>>> ds.head()
>>> print(table.frequencies('trigger', ds=ds))

For more complex designs and variables, you can override methods that provide complete control over the events. These are the transformations applied to the triggers extracted from raw files (in this order):

Defining data epochs

Once events are properly labeled, define MneExperiment.epochs.

There is one special epoch to define, which is called 'cov'. This is the data epoch that will be used to estimate the sensor noise covariance matrix for source estimation.

In order to find the right sel epoch parameter, it can be useful to actually load the events with MneExperiment.load_events() and test different selection strings. The epoch selection is determined by selection = event_ds.eval(epoch['sel']). Thus, a specific setting could be tested with:

>>> ds = e.load_events()
>>> print(ds.sub("event == 'value'"))
Bad channels

Flat channels are automatically excluded from the analysis.

An initial check for noisy channels can be done by looking at the raw data (see Pre-processing above). If this inspection reveals bad channels, they can be excluded using MneExperiment.make_bad_channels().

Another good check for bad channels is plotting the average evoked response, and looking for channels which are uncorrelated with neighboring channels. To plot the average before trial rejection, use:

>>> ds = e.load_epochs(epoch='epoch', reject=False)
>>> plot.TopoButterfly('meg', ds=ds)

The neighbor correlation can also be quantified, using:

>>> nc = neighbor_correlation(concatenate(ds['meg']))
>>> nc.sensor.names[nc < 0.3]
Datalist(['MEG 099'])

A simple way to cycle through subjects when performing a given pre-processing step is MneExperiment.next(). If a general threshold is adequate, the selection of bad channels based on neighbor-correlation can be automated using the MneExperiment.make_bad_channels_neighbor_correlation() method:

>>> for subject in e:
...     e.make_bad_channels_neighbor_correlation()
ICA

If preprocessing includes ICA, select which ICA components should be removed. To open the ICA selection GUI, The experiment raw state needs to be set to the ICA stage of the pipeline:

>>> e.set(raw='ica')
>>> e.make_ica_selection()

See MneExperiment.make_ica_selection() for more information on display options and on how to precompute ICA decomposition for all subjects.

When selecting ICA components for multiple subject, a simple way to cycle through subjects is MneExperiment.next(), like:

>>> e.make_ica_selection(epoch='epoch', decim=10)
>>> e.next()
subject: 'R1801' -> 'R2079'
>>> e.make_ica_selection(epoch='epoch', decim=10)
>>> e.next()
subject: 'R2079' -> 'R2085'
...
Trial selection

For each primary epoch that is defined, bad trials can be rejected using MneExperiment.make_epoch_selection(). Rejections are specific to a given raw state:

>>> e.set(raw='ica1-40')
>>> e.make_epoch_selection()
>>> e.next()
subject: 'R1801' -> 'R2079'
>>> e.make_epoch_selection()
...

To reject trials based on a pre-determined threshold, a loop can be used:

>>> for subject in e:
...     e.make_epoch_selection(auto=1e-12)
...
Analysis

With preprocessing completed, there are different options for analyzing the data.

The most flexible option is loading data from the desired processing stage using one of the many .load_... methods of the MneExperiment. For example, load a Dataset with source-localized condition averages using MneExperiment.load_evoked_stc(), then test a hypothesis using one of the mass-univariate test from the testnd module. To make this kind of analysis replicable, it is probably useful to write the complete analysis as a separate script that imports the experiment (see the example experiment folder).

Many statistical comparisons can also be specified in the MneExperiment.tests attribute, and then loaded directly using the MneExperiment.load_test() method. This has the advantage that the tests will be cached automatically and, once computed, can be loaded very quickly. However, these definitions are not quite as flexible as writing a custom script.

Finally, for tests defined in MneExperiment.tests, the MneExperiment can generate HTML report files. These are generated with the MneExperiment.make_report() and MneExperiment.make_report_rois() methods.

Warning

If source files are changed (raw files, epoch rejection or bad channel files, …) reports are not updated automatically unless the corresponding MneExperiment.make_report() function is called again. For this reason it is useful to have a script to generate all desired reports. Running the script ensures that all reports are up-to-date, and will only take seconds if nothing has to be recomputed (for an example see make-reports.py in the example experiment folder).

Example

The following is a complete example for an experiment class definition file (the source file can be found in the Eelbrain examples folder at examples/mouse/mouse.py):

# skip test: data unavailable
from eelbrain.pipeline import *


class Mouse(MneExperiment):

    # Name of the experimental session(s), used to locate *-raw.fif files
    sessions = 'CAT'

    # Pre-processing pipeline: each entry in `raw` specifies one processing step. The first parameter
    # of each entry specifies the source (another processing step or 'raw' for raw input data).
    raw = {
        # Maxwell filter as first step (taking input from raw data, 'raw')
        'tsss': RawMaxwell('raw', st_duration=10., ignore_ref=True, st_correlation=0.9, st_only=True),
        # Band-pass filter data between 1 and 40 Hz (taking Maxwell-filtered data as input, 'tsss)
        '1-40': RawFilter('tsss', 1, 40),
        # Perform ICA on filtered data
        'ica': RawICA('1-40', 'CAT', n_components=0.99),
    }

    # Variables determine how event triggeres are mapped to meaningful labels. Events are represented
    # as data-table in which each row corresponds to one event (i.e., one trigger). Each variable
    # defined here adds one column in that data-table, assigning a label or value to each event.
    variables = {
        # The first parameter specifies the source variable (here the trigger values),
        # the second parameter a mapping from source to target labels/values
        'stimulus': LabelVar('trigger', {(162, 163): 'target', (166, 167): 'prime'}),
        'prediction': LabelVar('trigger', {(162, 166): 'expected', (163, 167): 'unexpected'}),
    }

    # Epochs specify how to extract time-locked data segments ("epochs") from the continuous data.
    epochs = {
        # A PrimaryEpoch definition extracts epochs directly from continuous data. The first argument
        # specifies the recording session from which to extract the data (here: 'CAT'). The second
        # argument specifies which events to extract the data from (here: all events at which the
        # 'stimulus' variable, defined above, has a value of either 'prime' or 'target').
        'word': PrimaryEpoch('CAT', "stimulus.isin(('prime', 'target'))", samplingrate=200),
        # A secondary epoch inherits its properties from the base epoch ("word") unless they are
        # explicitly modified (here, selecting a subset of events)
        'prime': SecondaryEpoch('word', "stimulus == 'prime'"),
        'target': SecondaryEpoch('word', "stimulus == 'target'"),
        # The 'cov' epoch defines the data segments used to compute the noise covariance matrix for
        # source localization
        'cov': SecondaryEpoch('prime', tmax=0),
    }

    tests = {
        '=0': TTestOneSample(),
        'surprise': TTestRelated('prediction', 'unexpected', 'expected'),
        'anova': ANOVA('prediction * subject'),
    }

    parcs = {
        'frontotemporal-lh': CombinationParc('aparc', {
            'frontal-lh': 'parsorbitalis + parstriangularis + parsopercularis',
            'temporal-lh': 'transversetemporal + superiortemporal + '
                           'middletemporal + inferiortemporal + bankssts',
            }, views='lateral'),
    }


root = '~/Data/Mouse'
e = Mouse(root)

The event structure is illustrated by looking at the first few events:

>>> from mouse import *
>>> ds = e.load_events()
>>> ds.head()
trigger   i_start   T        SOA     subject   stimulus   prediction
--------------------------------------------------------------------
182       104273    104.27   12.04   S0001
182       116313    116.31   1.313   S0001
166       117626    117.63   0.598   S0001     prime      expected
162       118224    118.22   2.197   S0001     target     expected
166       120421    120.42   0.595   S0001     prime      expected
162       121016    121.02   2.195   S0001     target     expected
167       123211    123.21   0.596   S0001     prime      unexpected
163       123807    123.81   2.194   S0001     target     unexpected
167       126001    126      0.598   S0001     prime      unexpected
163       126599    126.6    2.195   S0001     target     unexpected

Experiment Definition

Basic setup
MneExperiment.owner

Set MneExperiment.owner to your email address if you want to be able to receive notifications. Whenever you run a sequence of commands with mne_experiment.notification: you will get an email once the respective code has finished executing or run into an error, for example:

>>> e = MyExperiment()
>>> with e.notification:
...     e.make_report('mytest', tstart=0.1, tstop=0.3)
...

will send you an email as soon as the report is finished (or the program encountered an error)

MneExperiment.auto_delete_results

Whenever a MneExperiment instance is initialized with a valid root path, it checks whether changes in the class definition invalidate previously computed results. By default, the user is prompted to confirm the deletion of invalidated results. Set auto_delete_results to True to delete them automatically without interrupting initialization.

MneExperiment.auto_delete_cache

MneExperiment caches various intermediate results. By default, if a change in the experiment definition would make cache files invalid, the outdated files are automatically deleted. Set auto_delete_cache to 'ask' to ask for confirmation before deleting files. This can be useful to prevent accidentally deleting files that take long to compute when editing the pipeline definition. When using this option, set MneExperiment.screen_log_level to 'debug' to learn about what change caused the cache to be invalid.

MneExperiment.screen_log_level

Determines the amount of information displayed on the screen while using an MneExperiment (see logging).

MneExperiment.meg_system

Starting with mne 0.13, fiff files converted from KIT files store information about the system they were collected with. For files converted earlier, the MneExperiment.meg_system attribute needs to specify the system the data were collected with. For data from NYU New York, the correct value is meg_system="KIT-157".

MneExperiment.merge_triggers

Use a non-default merge parameter for load.fiff.events().

MneExperiment.trigger_shift

Set this attribute to shift all trigger times by a constant (in seconds). For example, with trigger_shift = 0.03 a trigger that originally occurred 35.10 seconds into the recording will be shifted to 35.13. If the trigger delay differs between subjects, this attribute can also be a dictionary mapping subject names to shift values, e.g. trigger_shift = {'R0001': 0.02, 'R0002': 0.05, ...}.

Subjects
MneExperiment.subject_re

Subjects are identified on initialization by looking for folders in the data directory (meg by default) whose name matches the MneExperiment.subject_re regular expression. By default, this is '(R|A|Y|AD|QP)(\d{3,})$', which matches R-numbers like R1234, but also numbers prefixed by A, Y, AD or QP (for information about how to define a different regular expression, see re).

Defaults
MneExperiment.defaults

The defaults dictionary can contain default settings for experiment analysis parameters (see State Parameters), e.g.:

defaults = {'epoch': 'my_epoch',
            'cov': 'noreg',
            'raw': '1-40'}
Pre-processing (raw)
MneExperiment.raw

Define a pre-processing pipeline as a series of linked processing steps (mne refers to data that is not time-locked to specific events as Raw, with filenames matching *-raw.fif):

RawFilter(source[, l_freq, h_freq, cache])

Filter raw pipe

RawICA(source, session[, method, …])

ICA raw pipe

RawApplyICA(source, ica[, cache])

Apply ICA estimated in a RawICA pipe

RawMaxwell(source[, bad_condition, cache])

Maxwell filter raw pipe

RawSource([filename, reader, sysname, …])

Raw data source

RawReReference(source[, reference, add, …])

Re-reference EEG data

The raw data that constitutes the input to the pipeline can be accessed in a pipe named "raw" (the input data can be customized by adding a RawSource pipe). Each subsequent preprocessing step is defined with its input as first argument (source).

For example, the following definition sets up a pipeline using TSSS, a band-pass filter and ICA:

class Experiment(MneExperiment):

    sessions = 'session'

    raw = {
        'tsss': RawMaxwell('raw', st_duration=10., ignore_ref=True, st_correlation=0.9, st_only=True),
        '1-40': RawFilter('tsss', 1, 40),
        'ica': RawICA('1-40', 'session', 'extended-infomax', n_components=0.99),
    }

To use the raw --> TSSS --> 1-40 Hz band-pass pipeline, use e.set(raw="1-40"). To use raw --> TSSS --> 1-40 Hz band-pass --> ICA, select e.set(raw="ica").

Note

Continuous files take up a lot of hard drive space. By default, files for most pre-processing steps are cached This can be controlled with the cache parameter. To delete files correspoding to a specific step (e.g., raw='1-40'), use the MneExperiment.rm() method:

>>> e.rm('cached-raw-file', True, raw='1-40')
Event variables
MneExperiment.variables

Event variables add labels and variables to the events:

LabelVar(source, codes[, default, session])

Variable assigning labels to values

EvalVar(code[, session])

Variable based on evaluating a statement

GroupVar(groups[, session])

Group membership for each subject

Most of the time, the main purpose of this attribute is to turn trigger values into meaningful labels:

class Mouse(MneExperiment):

    variables = {
        'stimulus': LabelVar('trigger', {(162, 163): 'target', (166, 167): 'prime'}),
        'prediction': LabelVar('trigger', {162: 'expected', 163: 'unexpected'}),
    }

This defines a variable called “stimulus”, and on this variable all events that have triggers 162 and 163 have the value "target", and events with trigger 166 and 167 have the value "prime". The “prediction” variable only labels triggers 162 and 163. Unmentioned trigger values are assigned the empty string ('').

Epochs
MneExperiment.epochs

Epochs are specified as a {name: epoch_definition} dictionary. Names are str, and epoch_definition are instances of the classes described below:

PrimaryEpoch(session[, sel])

Epoch based on selecting events from a raw file

SecondaryEpoch(base[, sel])

Epoch inheriting events from another epoch

SuperEpoch(sub_epochs, **kwargs)

Combine several other epochs

Examples:

epochs = {
    # some primary epochs:
    'picture': PrimaryEpoch('words', "stimulus == 'picture'"),
    'word': PrimaryEpoch('words', "stimulus == 'word'"),
    # use the picture baseline for the sensor covariance estimate
    'cov': SecondaryEpoch('picture', tmax=0),
    # another secondary epoch:
    'animal_words': SecondaryEpoch('noun', sel="word_type == 'animal'"),
    # a superset-epoch:
    'all_stimuli': SuperEpoch(('picture', 'word')),
}
Tests
MneExperiment.tests

Statistical tests are defined as {name: test_definition} dictionary. Test- definitions are defined from the following:

TTestOneSample([tail])

One-sample t-test

TTestRelated(model, c1, c0[, tail])

Related measures t-test

TTestIndependent(model, c1, c0[, tail])

Independent measures t-test (comparing groups of subjects)

ANOVA(x[, model, vars])

ANOVA test

TContrastRelated(model, contrast[, tail])

Contrasts of T-maps (see eelbrain.testnd.TContrastRelated)

TwoStageTest(stage_1[, vars, model])

Two-stage test: T-test of regression coefficients

Example:

tests = {
    'my_anova': ANOVA('noise * word_type * subject'),
    'my_ttest': TTestRelated('noise', 'a_lot_of_noise', 'no_noise'),
}
Subject groups
MneExperiment.groups

A subject group called 'all' containing all subjects is always implicitly defined. Additional subject groups can be defined in MneExperiment.groups with {name: group_definition} entries:

Group(subjects)

Group defined as collection of subjects

SubGroup(base, exclude)

Group defined by removing subjects from a base group

Example:

groups = {
    'good': SubGroup('all', ['R0013', 'R0666']),
    'bad': Group(['R0013', 'R0666']),
}
Parcellations (parcs)
MneExperiment.parcs

The parcellation determines how the brain surface is divided into regions. A number of standard parcellations are automatically defined (see parc/mask (parcellations) below). Additional parcellations can be defined in the MneExperiment.parcs dictionary with {name: parc_definition} entries. There are a couple of different ways in which parcellations can be defined, described below.

SubParc(base, labels[, views])

A subset of labels in another parcellation

CombinationParc(base, labels[, views])

Recombine labels from an existing parcellation

SeededParc(seeds[, mask, surface, views])

Parcellation that is grown from seed coordinates

IndividualSeededParc(seeds[, mask, surface, …])

Seed parcellation with individual seeds for each subject

FreeSurferParc([views])

Parcellation that is created outside Eelbrain for each subject

FSAverageParc([views])

Fsaverage parcellation that is morphed to individual subjects
Visualization defaults
MneExperiment.brain_plot_defaults

The MneExperiment.brain_plot_defaults dictionary can contain options that changes defaults for brain plots (for reports and movies). The following options are available:

surf‘inflated’ | ‘pial’ | ‘smoothwm’ | ‘sphere’ | ‘white’

Freesurfer surface to use as brain geometry.

viewsstr | iterator of str

View or views to show in the figure. Can also be set for each parcellation, see MneExperiment.parc.

foregroundmayavi color

Figure foreground color (i.e., the text color).

backgroundmayavi color

Figure background color.

smoothing_stepsNone | int

Number of smoothing steps to display data.

State Parameters

These are parameters that can be set after an MneExperiment has been initialized to affect the analysis, for example:

>>> my_experiment = MneExperiment()
>>> my_experiment.set(raw='1-40', cov='noreg')

sets up my_experiment to use a 1-40 Hz band-pass filter as preprocessing, and to use sensor covariance matrices without regularization. Most methods also accept state parameters, so MneExperiment.set() does not have to be used separately.

session

Which raw session to work with (one of MneExperiment.sessions; usually set automatically when epoch is set)

visit

Which visit to work with (one of MneExperiment.visits)

raw

Select the preprocessing pipeline applied to the continuous data. Options are all the processing steps defined in MneExperiment.raw, as well as "raw" for using unprocessed raw data.

subject

Any subject in the experiment (subjects are identified based on MneExperiment.subject_re).

group

Any group defined in MneExperiment.groups. Will restrict the analysis to that group of subjects.

epoch

Any epoch defined in MneExperiment.epochs. Specify the epoch on which the analysis should be conducted.

rej (trial rejection)

Trial rejection can be turned off e.set(rej=''), meaning that no trials are rejected, and back on, meaning that the corresponding rejection files are used e.set(rej='man').

model

While the epoch state parameter determines which events are included when loading data, the model parameter determines how these events are split into different condition cells. The parameter should be set to the name of a categorial event variable which defines the desired cells. In the Example, e.load_evoked(epoch='target', model='prediction') would load responses to the target, averaged for expected and unexpected trials.

Cells can also be defined based on crossing two variables using the % sign. In the Example, to load corresponding primes together with the targets, you would use e.load_evoked(epoch='word', model='stimulus % prediction').

equalize_evoked_count

By default, the analysis uses all epochs marked as good during rejection. Set equalize_evoked_count='eq' to discard trials to make sure the same number of epochs goes into each cell of the model (see equal_count parameter to Dataset.aggregate()).

‘’ (default)

Use all epochs.

‘eq’

Make sure the same number of epochs n is used in each cell by discarding epochs. The first n epochs are used for each condition (assuming that habituation increases by condition).

cov

The method for correcting the sensor covariance.

‘noreg’

Use raw covariance as estimated from the data (do not regularize).

‘bestreg’ (default)

Find the regularization parameter that leads to optimal whitening of the baseline.

‘reg’

Use the default regularization parameter (0.1).

‘auto’

Use automatic selection of the optimal regularization method.

src

The source space to use.

  • ico-x: Surface source space based on icosahedral subdivision of the white matter surface x steps (e.g., ico-4, the default).

  • vol-x: Volume source space based on a volume grid with x mm resolution (x is the distance between sources, e.g. vol-10 for a 10 mm grid).

inv

What inverse solution to use for source localization. This parameter can also be set with MneExperiment.set_inv(), which has a more detailed description of the options. The inverse solution can be set directly using the appropriate string as in e.set(inv='fixed-1-MNE').

parc/mask (parcellations)

The parcellation determines how the brain surface is divided into regions. There are a number of built-in parcellations:

  • FreeSurfer Parcellations: aparc.a2005s, aparc.a2009s, aparc, aparc.DKTatlas, PALS_B12_Brodmann, PALS_B12_Lobes, PALS_B12_OrbitoFrontal, PALS_B12_Visuotopic.

  • cortex: All sources in cortex, based on the FreeSurfer “cortex” label.

  • lobes: Modified version of PALS_B12_Lobes in which the limbic lobe is merged into the other 4 lobes.

  • lobes-op: One large region encompassing occipital and parietal lobe in each hemisphere.

  • lobes-ot: One large region encompassing occipital and temporal lobe in each hemisphere.

Additional parcellation can be defined in the MneExperiment.parcs attribute. Parcellations are used in different contexts:

  • When loading source space data, the current parc state determines the parcellation of the souce space (change the state parameter with e.set(parc='aparc')).

  • When loading tests, setting the parc parameter treats each label as a separate ROI. For spatial cluster-based tests that means that no clusters can cross the boundary between two labels. On the other hand, using the mask parameter treats all named labels as connected surface, but discards any sources labeled as "unknown". For example, loading a test with mask='lobes' will perform a whole-brain test on the cortex, while discarding subcortical sources.

Parcellations are set with their name, with the expception of SeededParc: for those, the name is followed by the radious in mm, for example, to use seeds defined in a parcellation named 'myparc' with a radius of 25 mm around the seed, use e.set(parc='myparc-25').

connectivity

Possible values: '', 'link-midline'

Connectivity refers to the edges connecting data channels (sensors for sensor space data and sources for source space data). These edges are used to find clusters in cluster-based permutation tests. For source spaces, the default is to use FreeSurfer surfaces in which the two hemispheres are unconnected. By setting connectivity='link-midline', this default connectivity can be modified so that the midline gyri of the two hemispheres get linked at sources that are at most 15 mm apart. This parameter currently does not affect sensor space connectivity.

select_clusters (cluster selection criteria)

In thresholded cluster test, clusters are initially filtered with a minimum size criterion. This can be changed with the select_clusters analysis parameter with the following options:

Name

Min time

Min sources

Min sensors

"all"

"10ms"

10 ms

10

4

"" (default)

25 ms

10

4

"large"

25 ms

20

8

To change the cluster selection criterion use for example:

>>> e.set(select_clusters='all')

See also

Indices and tables

Eelbrain relies on NumPy, SciPy, Matplotlib, MNE-Python, PySurfer, WxPython, Cython and incorporates icons from the Tango Desktop Project.

Current funding: National Institutes of Health (NIH) grant R01-DC-014085 (since 2016). Past funding: NYU Abu Dhabi Institute grant G1001 (2011-2016).