Exploring xarray - Analyzing iEEG data

40 minute read

Published:

Analyzing intracranial electrophysiology data with xarray

Over the last few years, it has been exciting to see the xarray project evolve, add new functionality, and mature. This post is an attempt at giving xarray another visit to see how it could integrate into electrophysiology workflows.

A quick background on our data

It is common in neuroscience to ask individuals to perform a task over and over again. You record the activity in the brain each time they perform the task (called an "epoch" or a "trial"). Time is recorded relative to some onset when the task begins. That is t==0. The result is usually a matrix of epochs x channels x time. You can do a lot of stuff with this data, but our task in this paper is to detect changes in neural activity at trial onset (t==0).

In our case, we've got a small dataset from an old paper of mine. The repository contains several tutorial notebooks and sample data to describe predictive modeling in cognitive neuroscience. You can find the repository here. The task that individuals were performing was passively listening to spoken sentences through a speaker. While they did this, we recorded electrical activity at the surface of their brain (these were surgical patients, and had implanted electrodes under their scalp).

In the Feature Extraction notebook, I covered how to do some simple data manipulation and feature extraction with timeseries analysis. Let's try to re-create some of the main steps in that tutorial, but now using xarray as an in-memory structure for our data.

Note: The goal here is to learn a bit about xarray moreso than to discuss ecog modeling, so I'll spend more time talking about my thoughts on the various functions/methods/etc in Xarray than talking about neuroscience.

In this post, we'll perform a few common processing and extraction steps. The goal is to do a few munging operations that require manipulating data and visualizing simple statistics.

# Imports we'll use later
import mne
import numpy as np
import matplotlib.pyplot as plt
from download import download
import os
from sklearn.preprocessing import scale
import xarray as xr
xr.set_options(display_style="html")

import warnings
warnings.simplefilter('ignore')
%matplotlib inline

We'll load the data from my GitHub repository (probably not the most efficient way to store or retrieve the data, but hey, this was 3 years ago :-) ).

url_epochs = "https://github.com/choldgraf/paper-encoding_decoding_electrophysiology/blob/master/raw_data/ecog-epo.fif?raw=true"

path_data = download(url_epochs, './ecog-epo.fif', replace=True)
ecog = mne.read_epochs(path_data, preload=True)
os.remove(path_data)

Here's what the raw data looks like - each horizontal line is electrical activity in a channel over time. The faint vertical green lines show the onset of each trial (they are concatenated together, but in reality there's a bit of time between trials). This will be one of the last times we use MNE hopefully.

_ = ecog.plot(scalings='auto', n_epochs=5, n_channels=10)

Converting to xarray

First off, we'll define a helper function that converts the MNE Epochs object into an xarray DataArray object. DataArrays provide an N-Dimensional representation of data, but with the option to include a lot of extra metadata.

DataArrays are useful because you can include information about each dimension of the data. For example, we can tell our DataArray the name, values, and units of each dimension. In this case, in our case one dimension is "time" so we can label it as such.

def epochs_to_dataarray(epochs):
    """A simple function to convert an Epochs object to DataArray"""
    da = xr.DataArray(
    epochs._data,
    dims=['epoch', 'channel', 'time'],
    coords={
        'time': ecog.times,
        'channel': ecog.ch_names,
        'epoch': range(ecog._data.shape[0])
    },
    name='Sample dataset',
    attrs=dict(ecog.info)
    )
    return da

Just look at all the metadata that we were able to pack into the DataArray. Almost all of MNE's metadata fit nicely into .attrs.

# There's quite a lot of output, so keep scrolling down!
da = epochs_to_dataarray(ecog)
da
Show/Hide data repr Show/Hide attributes
xarray.DataArray
'Sample dataset'
  • epoch: 29
  • channel: 32
  • time: 2250
  • -118.1 -121.3 -117.6 -111.4 -107.3 ... -108.2 -105.9 -105.4 -102.5
    array([[[-118.10999298, -121.30046844, -117.64116669, ...,
               20.78869057,   11.38018703,    8.12319088],
            [-126.95540619, -120.84265137, -105.75999451, ...,
               25.57149315,   28.6223793 ,   24.82316971],
            [ -28.1061306 ,  -16.39686775,  -21.64873505, ...,
              -31.78049469,  -24.28160286,  -25.21776772],
            ...,
            [  18.55451965,   22.87419701,   25.66038895, ...,
               25.84903336,   19.71798515,   14.92882633],
            [ -44.79590607,  -49.5848465 ,  -52.08906555, ...,
              -48.92469025,  -54.11177826,  -60.23026276],
            [-105.09076691, -108.88887024, -113.57571411, ...,
               -3.28264022,    5.5660615 ,   13.18798733]],
    
           [[  20.84839821,   18.30262566,   20.12442017, ...,
              -19.40711594,  -18.51371574,  -20.21596718],
            [  18.0386982 ,    5.35804272,   -6.45382166, ...,
              -42.1801033 ,  -44.59093857,  -47.80654144],
            [  76.62892914,   73.66176605,   64.21563721, ...,
              -12.31007576,  -21.26686478,  -26.83379555],
            ...,
            [  58.06867599,   53.49727631,   51.89107895, ...,
               42.84626007,   52.36317062,   55.75439453],
            [ -24.79452324,  -18.58691406,   -8.12459183, ...,
              -45.95447159,  -42.66623306,  -41.06335449],
            [-110.1750946 , -101.40830994,  -88.24209595, ...,
              -26.7741375 ,  -32.31249619,  -37.07465744]],
    
           [[  44.85819626,   57.10462952,   68.87900543, ...,
               13.73396015,   28.33650208,   40.76254654],
            [  10.10721684,   19.94246101,   31.81379509, ...,
              120.01114655,  119.50100708,  128.0358429 ],
            [ -96.35434723,  -89.9883194 ,  -72.41400909, ...,
              140.28143311,  131.52635193,  133.30158997],
            ...,
            [   1.81244755,   -1.55214143,   -2.54425907, ...,
              -30.9526062 ,  -31.50832367,  -32.25312042],
            [  44.95462799,   38.7992363 ,   33.4192009 , ...,
               64.5262146 ,   62.74245834,   61.11261368],
            [ -37.66891861,  -37.9514122 ,  -40.97248077, ...,
               61.204319  ,   58.57411194,   56.27007675]],
    
           ...,
    
           [[  88.75709534,   83.42398071,   76.10377502, ...,
              -23.7650528 ,  -23.43914986,  -19.55894852],
            [  63.89505005,   63.79450226,   61.43893433, ...,
              -54.35200119,  -59.907547  ,  -69.36830139],
            [  39.31270599,   46.15828705,   51.12403107, ...,
              -30.19447899,  -27.1382637 ,  -19.75760651],
            ...,
            [  -2.09485602,  -13.05238819,  -19.73264503, ...,
              150.31195068,  129.28895569,  114.28322601],
            [-142.40054321, -142.72599792, -140.20773315, ...,
               90.81705475,   82.73035431,   77.06429291],
            [-161.34713745, -160.04364014, -151.092453  , ...,
              238.87109375,  230.41633606,  222.64045715]],
    
           [[  67.85353088,   72.97834778,   70.99682617, ...,
               99.13801575,   95.14775085,   86.59616089],
            [ 169.86994934,  172.12055969,  173.22143555, ...,
              132.95147705,  135.66346741,  141.64862061],
            [  80.4973526 ,   77.63601685,   73.43898773, ...,
               91.37036896,   84.76886749,   82.48825836],
            ...,
            [ -13.24606895,   -8.65564728,    6.53395557, ...,
             -105.62326813, -104.58720398, -107.99065399],
            [  22.20972443,   49.33536148,   85.349823  , ...,
             -103.54721069, -102.84375   , -104.26371002],
            [ -33.03855133,   -9.5574255 ,   10.97984695, ...,
             -148.20103455, -148.95999146, -154.4241333 ]],
    
           [[ -98.52385712,  -99.22042847,  -99.56749725, ...,
               49.72173309,   41.6667099 ,   31.01161194],
            [ -70.90259552,  -66.83739471,  -68.66448975, ...,
               35.73180008,   33.55708694,   28.1260643 ],
            [ -29.78354645,  -34.5358963 ,  -40.23687363, ...,
              -39.39603043,  -42.75005341,  -50.70994186],
            ...,
            [ 110.57292175,  103.35498047,   95.09171295, ...,
             -162.01602173, -164.45585632, -168.76620483],
            [  38.4355278 ,   39.19929504,   43.87767792, ...,
             -108.41282654, -106.09828949, -103.17415619],
            [  98.9158783 ,  105.16780853,  116.0725174 , ...,
             -105.856987  , -105.39962769, -102.4755249 ]]])
    • time
      (time)
      float64
      -1.5 -1.497 -1.493 ... 5.993 5.997
      array([-1.5     , -1.496667, -1.493333, ...,  5.99    ,  5.993333,  5.996667])
    • channel
      (channel)
      'ch_0' 'ch_1' ... 'ch_30' 'ch_31'
      array(['ch_0', 'ch_1', 'ch_2', 'ch_3', 'ch_4', 'ch_5', 'ch_6', 'ch_7', 'ch_8',
             'ch_9', 'ch_10', 'ch_11', 'ch_12', 'ch_13', 'ch_14', 'ch_15', 'ch_16',
             'ch_17', 'ch_18', 'ch_19', 'ch_20', 'ch_21', 'ch_22', 'ch_23', 'ch_24',
             'ch_25', 'ch_26', 'ch_27', 'ch_28', 'ch_29', 'ch_30', 'ch_31'],
            dtype='<U5')
    • epoch
      (epoch)
      int64
      0 1 2 3 4 5 6 ... 23 24 25 26 27 28
      array([ 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])
  • file_id :
    {'version': 65539, 'machid': array([808661043, 808661043], dtype=int32), 'secs': 1471644431, 'usecs': 91610}
    events :
    []
    hpi_results :
    []
    hpi_meas :
    []
    subject_info :
    None
    device_info :
    None
    helium_info :
    None
    hpi_subsystem :
    None
    proc_history :
    []
    meas_id :
    {'version': 65539, 'machid': array([808661043, 808661043], dtype=int32), 'secs': 1471644431, 'usecs': 91754}
    experimenter :
    None
    description :
    None
    proj_id :
    None
    proj_name :
    None
    meas_date :
    (0, 0)
    utc_offset :
    None
    sfreq :
    300.0
    highpass :
    0.0
    lowpass :
    500.0
    line_freq :
    None
    gantry_angle :
    None
    chs :
    [{'scanno': 1, 'logno': 1, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_0', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 2, 'logno': 2, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_1', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 3, 'logno': 3, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_2', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 4, 'logno': 4, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_3', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 5, 'logno': 5, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_4', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 6, 'logno': 6, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_5', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 7, 'logno': 7, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_6', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 8, 'logno': 8, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_7', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 9, 'logno': 9, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_8', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 10, 'logno': 10, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_9', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 11, 'logno': 11, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_10', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 12, 'logno': 12, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_11', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 13, 'logno': 13, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_12', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 14, 'logno': 14, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_13', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 15, 'logno': 15, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_14', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 16, 'logno': 16, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_15', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 17, 'logno': 17, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_16', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 18, 'logno': 18, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_17', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 19, 'logno': 19, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_18', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 20, 'logno': 20, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_19', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 21, 'logno': 21, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_20', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 22, 'logno': 22, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_21', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 23, 'logno': 23, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_22', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 24, 'logno': 24, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_23', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 25, 'logno': 25, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_24', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 26, 'logno': 26, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_25', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 27, 'logno': 27, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_26', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 28, 'logno': 28, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_27', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 29, 'logno': 29, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_28', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 30, 'logno': 30, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_29', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 31, 'logno': 31, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_30', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}, {'scanno': 32, 'logno': 32, 'kind': 902, 'range': 1.0, 'cal': 1.0, 'coil_type': 1, 'loc': array([0., 0., 0., 1., 0., 0., 0., 1., 0., 0., 0., 1.]), 'unit': 107, 'unit_mul': 0, 'ch_name': 'ch_31', 'coord_frame': 0 (FIFFV_COORD_UNKNOWN)}]
    dev_head_t :
    <Transform | MEG device->head> [[1. 0. 0. 0.] [0. 1. 0. 0.] [0. 0. 1. 0.] [0. 0. 0. 1.]]
    ctf_head_t :
    None
    dev_ctf_t :
    None
    dig :
    []
    bads :
    []
    ch_names :
    ['ch_0', 'ch_1', 'ch_2', 'ch_3', 'ch_4', 'ch_5', 'ch_6', 'ch_7', 'ch_8', 'ch_9', 'ch_10', 'ch_11', 'ch_12', 'ch_13', 'ch_14', 'ch_15', 'ch_16', 'ch_17', 'ch_18', 'ch_19', 'ch_20', 'ch_21', 'ch_22', 'ch_23', 'ch_24', 'ch_25', 'ch_26', 'ch_27', 'ch_28', 'ch_29', 'ch_30', 'ch_31']
    nchan :
    32
    projs :
    []
    comps :
    []
    acq_pars :
    None
    acq_stim :
    None
    custom_ref_applied :
    False
    xplotter_layout :
    None
    kit_system_id :
    None

The data consists of many trials, channels, and timepoints. Let's start by selecting a time region within each trial that we can visualize more cleanly.

Subsetting out data with da.sel

In xarray, we select items with the sel and isel method. This behaves kind of like the pandas loc and iloc methods, however because we have named dimensions, we can directly specify them in our call.

# We'll drop a subset of timepoints for visualization
da = da.sel(time=slice(-1, 3))

Now let's calculate the average across all epochs for each electrode/time point. This is a reduction of our data array, in that it reduces the number of dimensions. Xarray has many of the same statistical methods that NumPy does. An interesting twist is that you can specify named dimensions instead of simply an axis=<integer> argument. In addition, we'll choose the colors that we'll use for cycling through our channels - because we can quickly reference the channels axis by name, we don't need to remember which axis corresponds to channels.

fig, ax = plt.subplots(figsize=(15, 5))
n_channels = da['channel'].shape[0]
ax.set_prop_cycle(color=plt.cm.viridis(np.linspace(0, 1, n_channels)))
da.mean(dim='epoch').plot.line(x='time', hue='channel')
ax.get_legend().remove()

It doesn't look like much is going on...let's see if we can clean it up a bit.

De-meaning the data with da.where

First off - we'll subtract the "pre-baseline mean" from each trial. This makes it easier to visualize how each channel's activity changed at time == 0.

To accomplish this we'll use da.where. This takes some kind of boolean-style mask, does a bunch of clever projections according to the names of coordinates, and returns the dataarray masked values removed (as NaNs) and other values unchanged. We can use this to calculate the mean of each channel / epoch only for the pre-baseline timepoints.

# This returns a version of the data array with NaNs where the query is False
# The dimensions will intelligently broadcast 
prebaseline_mean = da.where(da.time < 0).mean(dim='time')
da_demeaned = da - prebaseline_mean

Now we can visualize the de-baseline-meaned data

fig, ax = plt.subplots(figsize=(15, 5))
ax.set_prop_cycle(color=plt.cm.viridis(np.linspace(0, 1, da['channel'].shape[0])))
da_demeaned.mean(dim='epoch').plot.line(x='time', hue='channel')
ax.get_legend().remove()

Hmmm, there still doesn't seem to be much going on (that channel down at the bottom looks noisy to me, rather than having a meaningful signal) so let's transform this signal into something with a bit more SNR to it.

Extracting a more useful feature with xr.apply_ufunc

Without going into too much details on the neuroscience, iEEG data is particularly useful because there is information about neuronal activity in the higher frequency parts of the signal (AKA, parts of the electrical signal that change very quickly, but have very low amplitude). To pull that out, we'll do the following:

  • High-pass filter the signal, which will remove all the slow-moving components
  • Calculate the envelope of the signal, which will tell us the power of high-frequency activity over time.

High-pass filtering the signal

MNE has a lot of nice functions for filtering a timeseries. Most of these operate on numpy arrays instead of MNE objects. We'll use xarray's apply_ufunc function to simply map that function onto our dataarray. xarray should keep track of the metadata (e.g. coordinates etc) and output a new DataArray with updated values.

flow = 80
fhigh = 140
da_lowpass = xr.apply_ufunc(
    mne.filter.filter_data, da,
   kwargs=dict(
       sfreq=da.sfreq,
       l_freq=flow,
       h_freq=fhigh,
   )
)
Setting up band-pass filter from 80 - 1.4e+02 Hz

FIR filter parameters
---------------------
Designing a one-pass, zero-phase, non-causal bandpass filter:
- Windowed time-domain design (firwin) method
- Hamming window with 0.0194 passband ripple and 53 dB stopband attenuation
- Lower passband edge: 80.00
- Lower transition bandwidth: 20.00 Hz (-6 dB cutoff frequency: 70.00 Hz)
- Upper passband edge: 140.00 Hz
- Upper transition bandwidth: 10.00 Hz (-6 dB cutoff frequency: 145.00 Hz)
- Filter length: 99 samples (0.330 sec)

Visualizing our data, we can see all the slower fluctuations (e.g. long arcs over time) are gone.

fig, ax = plt.subplots(figsize=(15, 5))
da_lowpass.mean(dim='epoch').plot.line(x='time')
ax.get_legend().remove()

Calculate the envelope of this signal with da.groupby

Next, we'll calculate the envelope of the high-pass-filtered data. This is roughly the power that is present in these high frequencies over time. We do so by using something called a hilbert transform.

MNE also has a function for applying Hilbert transforms to data, but it has a weird quirk that expects the data to be of a particular shape. We can work around this by using our DataArray's groupby method. This works similar to DataFrame.groupby - we'll iterate through each channel, which will return a DataArray with shape epochs x timepoints. We can then calculate the Hilbert transform in each and re-combine into the original shape.

Note: This can be an expensive operation depending on the number of channels/epochs and the length of each trial. This might be a good place to insert paralellization via Dask.

def hilbert_2d(array):
    """Perform a Hilbert transforms on an (n_channels, n_times) array."""
    for ii, channel in enumerate(array):
        array[ii] = mne.filter._my_hilbert(channel, envelope=True)
    return array

da_hf_power = da_lowpass.groupby(da.coords['epoch']).apply(hilbert_2d)

The output dataarray should be the exact same shape, because we haven't done any dimensional reductions. If we take a look at the resulting data, we can see what seems to be more structure in there:

fig, ax = plt.subplots(figsize=(15, 5))
da_hf_power.mean(dim='epoch').plot.line(x='time', hue='channel')
ax.get_legend().remove()