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 channejupyls 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)

../../../_images/1a33359f7e36c6f0e3021b861ab9fb0e1e342ad040500e178b6b3ba20b50b85a.png

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

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 ]]])

Coordinates: (3)

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])

Attributes: (36)
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()

../../../_images/89ea6341d363741a43ff5d5130b80147da7d4212fd0c518a0780f4db141f0279.png

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()

../../../_images/c63be37ef82f5b9f4cd2e93d413b230d2bf2c9ec4389afa420748b35e0113054.png

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()

../../../_images/7dfe4668d2a17721cbe28581541766ea75e71b1d0bb840c1fdb5af174dedd237.png

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()

../../../_images/e1446b2497ff3fb1e0183cfa2eb972550f7c9809279dfed6ded7a4c150334caa.png

Cleaning up our HFA data#

Next let’s clean up this high-frequency activity (HFA) data.

Z-scoring our array#

Instead of simple de-meaning the data like before, we’ll re-scale our data using the same baseline timepoints. What we’d like to do is the following:

Calculate the mean and standard deviation across trials of all pre-baseline data values, per channel
Z-score each channel using this mean and standard deviation

Once again we’ll use the groupby / apply combination to apply our function to subsets of the data.

# For each channel, apply a z-score that uses the mean/std of pre-baseline activity for all trials
def z_score(activity):
    """Take a DataArray and apply a z-score using the baseline"""
    baseline = activity.where(activity.time < -.1 )
    return (activity - np.nanmean(baseline)) / np.nanstd(baseline)

da_hf_zscored = da_hf_power.groupby('channel').apply(z_score)

Taking a look at the result, we can see a much cleaner separation of activity for some of the channels after time==0.

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

../../../_images/2482864996fbf2da7b325ee48bbabf7880924fa1e2ca95fd9f5f268a9a8e80c2.png

Smoothing our HFA data#

Finally, let’s smooth this HFA so it has less jitter to it, and pick a smaller window that removes some of the filtering artifacts at the edges.

We’ll use the same filter_data function as before, but this time applied with the .groupby and .apply combination to show two ways of accomplishing the same thing. We’ll also use .sel to pick a subset of time for visualization

da_hf_zscored_lowpass = da_hf_zscored.groupby('epoch').apply(
    mne.filter.filter_data,
    sfreq=da.sfreq,
    l_freq=None,
    h_freq=10,
    verbose=False
)

Note that quickly selecting a subset of timepoints if we used numpy is much more verbose. Here’s a quick comparison:

# Numpy alone
mask_time = (times > -.8) * (times < 2.8)
epoch_dim = 0
da_hf_zscored_lowpass[..., mask_time].mean(epoch_dim)

# xarray
da_hf_zscored_lowpass.sel(time=slice(-.8, 2.8)).mean(dim='epoch')

fig, ax = plt.subplots(figsize=(15, 5))
da_hf_zscored_lowpass.mean(dim='epoch').sel(time=slice(-.8, 2.8)).plot.line(x='time', hue='channel')
ax.get_legend().remove()

../../../_images/0b44d99cd8f37bb91c76781c339d9db6a0f706f9eb0ed61b70b6351322a5bfed.png

Now we can see there are clearly some channels that become active just after t==0. We can reduce our dataarray to a single dimension of “mean post-baseline activity in each channel” and convert it to a DataFrame for further processing:

# Find the channel with the most activity by first converting to a dataframe
total_activity = da_hf_zscored_lowpass.sel(time=slice(0, 2)).mean(dim=['epoch', 'time'])
total_activity = total_activity.to_dataframe()
total_activity.head()

	Sample dataset
channel
ch_0	0.011101
ch_1	0.012152
ch_2	0.051911
ch_3	0.176188
ch_4	0.246021

Let’s grab the channel with maximal activation to look into a bit further.

max_chan = total_activity.squeeze().sort_values(ascending=False).index[0]

Time frequency analysis#

As a final step, let’s expand our DataArray and add another dimension. In the above steps we specifically focused on high-frequency activity. A more common approach is to first create a spectrogram of your data to see activity across many frequencies.

To do this, we’ll use another MNE function for creating a Time-Frequency Representation or TFR. We’ll define a range of frequencies, and apply MNE’s function directly on our DataArray. This will return a NumPy array with the filtered values.

frequencies = [2**ii for ii in np.arange(2, 9, .5)]
tfr = mne.time_frequency.tfr_array_morlet(
    da,
    sfreq=da.sfreq,
    freqs=frequencies,
    n_cycles=4,
)

# Take the absolute value to throw out the non-real parts of the numbers
tfr = np.abs(tfr)
tfr[:2, :2, :2, :2]

array([[[[160.04103045, 160.38413909],
         [171.73704543, 175.09249553]],

        [[283.28699104, 285.21855726],
         [241.65630528, 245.77098295]]],


       [[[ 93.99546124,  94.4406537 ],
         [ 78.02050045,  79.15324341]],

        [[ 47.85148993,  49.90845151],
         [ 53.97221461,  54.55793674]]]])

Convert this data into a DataArray with `.expand_dims`#

Next, we’ll convert this into a DataArray by using the metadata from our original DataArray. We can use the expand_dims method to create a new dimension for our DataArray. We’ll use this to store frequency information.

We’ll then reshape our new DataArray so that it matches the output of the MNE function, and use the copy method to create a new DataArray. By supplying the data= argument to copy, we directly insert the new data inside the generated DataArray.

da_tfr = (da
    .expand_dims(frequency=frequencies)
    .transpose('epoch', 'channel', 'frequency', 'time')
    .copy(data=np.log(tfr))
)

We can now visualize this time-frequency representation over time

fig, ax = plt.subplots(figsize=(15, 5))
(da_tfr
    .sel({'frequency': slice(None, 180), 'channel': max_chan})
    .mean('epoch')
    .plot.imshow(x='time', y='frequency')
)

<matplotlib.image.AxesImage at 0x7feb35626c88>

../../../_images/026ae0809f92b4be072bd8ba84caa22e3fa8b40ba05550e201724d674d0cabb4.png

Similar to our one-dimensional visualizations above, it can be hard to visualize relative changes in activity over a baseline (particularly because the amplitude scales inversely with the frequency).

Let’s apply a re-scaling function to our data so that we can see things more clearly. This time we’ll use MNE’s rescale function, which acts similarly to our zscore function above.

da_tfr_baselined = xr.apply_ufunc(
    mne.baseline.rescale,
    da_tfr,
    kwargs={'times': da_tfr.coords['time'], 'baseline': (None, -.1), "mode": 'zscore'}
)

Applying baseline correction (mode: zscore)

again, the result should be a DataArray, so we can directly visualize it:

(da_tfr_baselined
    .sel({'frequency': slice(None, 180), 'channel': max_chan, 'time': slice(-.8, 2.5)})
    .mean('epoch')
    .plot.imshow(x='time', y='frequency')
)

<matplotlib.image.AxesImage at 0x7feb06faaa58>

../../../_images/c6642bf543d37385bee49d402fd762abe7397f03891588d36c60ef95544dc3a8.png

Now we can see a clear increase in activity in the higher frequencies at t==0.

Wrapping up#

In all, I was pretty happy with what you can do using xarray’s DataArray structure. It’s pretty nice to be able to refer to axes by their names, and to make more intelligent selection / slicing operations using their coordinate values. Moreover, this post is just scratching the surface for how to use this information in a way that speeds up the exploration and analysis post.

For example, we might have sped-up some feature extraction steps by using a distributed processing framework like Dask in the operations above. Dask integrates nicely with xarray, and offers a lot of interesting opportunities to parallelize interactive computation. I’ll explore that in another blog post.

Finally - the goal of this post has largely been to learn a bit more about xarray. This means I might be totally mis-using functionality, or missing something that would have made the above process much easier. If anybody has tips or thoughts on the code above, please do reach out!

Draft

Recent Posts

Archives

Analyzing intracranial electrophysiology data with xarray#

A quick background on our data#

Converting to xarray#

Subsetting out data with `da.sel`#

De-meaning the data with `da.where`#

Extracting a more useful feature with `xr.apply_ufunc`#

High-pass filtering the signal#

Calculate the envelope of this signal with `da.groupby`#

Cleaning up our HFA data#

Z-scoring our array#

Smoothing our HFA data#

Time frequency analysis#

Convert this data into a DataArray with `.expand_dims`#

Combining the two with `xr.merge`#

Wrapping up#

Draft

Recent Posts

Archives

Analyzing intracranial electrophysiology data with xarray#

A quick background on our data#

Converting to xarray#

Subsetting out data with da.sel#

De-meaning the data with da.where#

Extracting a more useful feature with xr.apply_ufunc#

High-pass filtering the signal#

Calculate the envelope of this signal with da.groupby#

Cleaning up our HFA data#

Z-scoring our array#

Smoothing our HFA data#

Time frequency analysis#

Convert this data into a DataArray with .expand_dims#

Combining the two with xr.merge#

Wrapping up#

Subsetting out data with `da.sel`#

De-meaning the data with `da.where`#

Extracting a more useful feature with `xr.apply_ufunc`#

Calculate the envelope of this signal with `da.groupby`#

Convert this data into a DataArray with `.expand_dims`#

Combining the two with `xr.merge`#