Epilepsy Detection Using EEG Data
IPython notebook: download
In this example we’ll use the cesium library to compare
various techniques for epilepsy detection using a classic EEG time series dataset from
Andrzejak et al..
The raw data are separated into five classes: Z, O, N, F, and S; we will consider a
three-class classification problem of distinguishing normal (Z, O), interictal (N, F), and
ictal (S) signals.
The overall workflow consists of three steps: first, we “featurize” the time series by selecting some set of mathematical functions to apply to each; next, we build some classification models which use these features to distinguish between classes; finally, we validate our models by generating predictions for some unseen holdout set and comparing them to the true class labels.
First, we’ll load the data and inspect a representative time series from each class:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import seaborn; seaborn.set()
from cesium import datasets
eeg = datasets.fetch_andrzejak()
# Group together classes (Z, O), (N, F), (S) as normal, interictal, ictal
eeg["classes"] = eeg["classes"].astype('U16') # allocate memory for longer class names
eeg["classes"][np.logical_or(eeg["classes"]=="Z", eeg["classes"]=="O")] = "Normal"
eeg["classes"][np.logical_or(eeg["classes"]=="N", eeg["classes"]=="F")] = "Interictal"
eeg["classes"][eeg["classes"]=="S"] = "Ictal"
fig, ax = plt.subplots(1, len(np.unique(eeg["classes"])), sharey=True)
for label, subplot in zip(np.unique(eeg["classes"]), ax):
i = np.where(eeg["classes"] == label)[0][0]
subplot.plot(eeg["times"][i], eeg["measurements"][i])
subplot.set(xlabel="time (s)", ylabel="signal", title=label)
Featurization
Once the data is loaded, we can generate features for each time series using the
cesium.featurize
module. The featurize module includes many built-in choices of features which can be applied
for any type of time series data; here we’ve chosen a few generic features that do not have
any special biological significance.
If Celery is running, the time series will automatically be split among the available workers
and featurized in parallel; setting use_celery=False will cause the time series to be
featurized serially.
from cesium import featurize
features_to_use = ['amplitude',
'percent_beyond_1_std',
'maximum',
'max_slope',
'median',
'median_absolute_deviation',
'percent_close_to_median',
'minimum',
'skew',
'std',
'weighted_average']
fset_cesium = featurize.featurize_time_series(times=eeg["times"],
values=eeg["measurements"],
errors=None,
features_to_use=features_to_use,
targets=eeg["classes"], use_celery=True)
print(fset_cesium)
<xarray.Dataset>
Dimensions: (channel: 1, name: 500)
Coordinates:
* channel (channel) int64 0
* name (name) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 ...
target (name) object 'Normal' 'Normal' 'Normal' ...
Data variables:
minimum (name, channel) float64 -146.0 -254.0 -146.0 ...
amplitude (name, channel) float64 143.5 211.5 165.0 ...
median_absolute_deviation (name, channel) float64 28.0 32.0 31.0 31.0 ...
percent_beyond_1_std (name, channel) float64 0.1626 0.1455 0.1523 ...
maximum (name, channel) float64 141.0 169.0 184.0 ...
median (name, channel) float64 -4.0 -51.0 13.0 -4.0 ...
percent_close_to_median (name, channel) float64 0.505 0.6405 0.516 ...
max_slope (name, channel) float64 1.111e+04 2.065e+04 ...
skew (name, channel) float64 0.0328 -0.09271 ...
weighted_average (name, channel) float64 -4.132 -52.44 12.71 ...
std (name, channel) float64 40.41 48.81 47.14 ...
The output of
featurize_time_series
is an xarray.Dataset which contains all the feature information needed to train a machine
learning model: feature values are stored as data variables, and the time series index/class
label are stored as coordinates (a channel coordinate will also be used later for
multi-channel data).
Custom feature functions
Custom feature functions not built into cesium may be passed in using the
custom_functions keyword, either as a dictionary {feature_name: function}, or as a
dask graph. Functions should take
three arrays times, measurements, errors as inputs; details can be found in the
cesium.featurize
documentation.
Here we’ll compute five standard features for EEG analysis provided by Guo et al.
(2012):
import numpy as np
import scipy.stats
def mean_signal(t, m, e):
return np.mean(m)
def std_signal(t, m, e):
return np.std(m)
def mean_square_signal(t, m, e):
return np.mean(m ** 2)
def abs_diffs_signal(t, m, e):
return np.sum(np.abs(np.diff(m)))
def skew_signal(t, m, e):
return scipy.stats.skew(m)
Now we’ll pass the desired feature functions as a dictionary via the custom_functions
keyword argument.
guo_features = {
'mean': mean_signal,
'std': std_signal,
'mean2': mean_square_signal,
'abs_diffs': abs_diffs_signal,
'skew': skew_signal
}
fset_guo = featurize.featurize_time_series(times=eeg["times"], values=eeg["measurements"],
errors=None, targets=eeg["classes"],
features_to_use=list(guo_features.keys()),
custom_functions=guo_features,
use_celery=True)
print(fset_guo)
<xarray.Dataset>
Dimensions: (channel: 1, name: 500)
Coordinates:
* channel (channel) int64 0
* name (name) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 ...
target (name) object 'Normal' 'Normal' 'Normal' 'Normal' 'Normal' ...
Data variables:
abs_diffs (name, channel) float64 4.695e+04 6.112e+04 5.127e+04 ...
mean (name, channel) float64 -4.132 -52.44 12.71 -3.992 -18.0 ...
mean2 (name, channel) float64 1.65e+03 5.133e+03 2.384e+03 ...
skew (name, channel) float64 0.0328 -0.09271 -0.0041 0.06368 ...
std (name, channel) float64 40.41 48.81 47.14 47.07 44.91 45.02 ...
Multi-channel time series
The EEG time series considered here consist of univariate signal measurements along a
uniform time grid. But
featurize_time_series
also accepts multi-channel data; to demonstrate this, we will decompose each signal into
five frequency bands using a discrete wavelet transform as suggested by Subasi
(2005), and then
featurize each band separately using the five functions from above.
import pywt
n_channels = 5
eeg["dwts"] = [pywt.wavedec(m, pywt.Wavelet('db1'), level=n_channels-1)
for m in eeg["measurements"]]
fset_dwt = featurize.featurize_time_series(times=None, values=eeg["dwts"], errors=None,
features_to_use=list(guo_features.keys()),
targets=eeg["classes"],
custom_functions=guo_features)
print(fset_dwt)
<xarray.Dataset>
Dimensions: (channel: 5, name: 500)
Coordinates:
* channel (channel) int64 0 1 2 3 4
* name (name) int64 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 ...
target (name) object 'Normal' 'Normal' 'Normal' 'Normal' 'Normal' ...
Data variables:
abs_diffs (name, channel) float64 2.513e+04 1.806e+04 3.241e+04 ...
skew (name, channel) float64 -0.0433 0.06578 0.2999 0.1239 0.1179 ...
mean2 (name, channel) float64 1.294e+04 5.362e+03 2.321e+03 664.4 ...
mean (name, channel) float64 -17.08 -6.067 -0.9793 0.1546 0.03555 ...
std (name, channel) float64 112.5 72.97 48.17 25.77 10.15 119.8 ...
The output featureset has the same form as before, except now the channel coordinate is
used to index the features by the corresponding frequency band. The functions in
cesium.build_model
and cesium.predict
all accept featuresets from single- or multi-channel data, so no additional steps are
required to train models or make predictions for multichannel featuresets using the
cesium library.
Model Building
Model building in cesium is handled by the
build_model_from_featureset
function in the cesium.build_model submodule. The featureset output by
featurize_time_series
contains both the feature and target information needed to train a
model; build_model_from_featureset is simply a wrapper that calls the fit method of a
given scikit-learn model with the appropriate inputs. In the case of multichannel
features, it also handles reshaping the featureset into a (rectangular) form that is
compatible with scikit-learn.
For this example, we’ll test a random forest classifier for the built-in cesium features,
and a 3-nearest neighbors classifier for the others, as suggested by Guo et al.
(2012).
from cesium.build_model import build_model_from_featureset
from sklearn.ensemble import RandomForestClassifier
from sklearn.neighbors import KNeighborsClassifier
from sklearn.cross_validation import train_test_split
train, test = train_test_split(np.arange(len(eeg["classes"])), random_state=0)
rfc_param_grid = {'n_estimators': [8, 16, 32, 64, 128, 256, 512, 1024]}
model_cesium = build_model_from_featureset(fset_cesium.isel(name=train),
RandomForestClassifier(max_features='auto',
random_state=0),
params_to_optimize=rfc_param_grid)
knn_param_grid = {'n_neighbors': [1, 2, 3, 4]}
model_guo = build_model_from_featureset(fset_guo.isel(name=train),
KNeighborsClassifier(),
params_to_optimize=knn_param_grid)
model_dwt = build_model_from_featureset(fset_dwt.isel(name=train),
KNeighborsClassifier(),
params_to_optimize=knn_param_grid)
Prediction
Making predictions for new time series based on these models follows the same pattern:
first the time series are featurized using
featurize_timeseries,
and then predictions are made based on these features using
predict.model_predictions,
from sklearn.metrics import accuracy_score
from cesium.predict import model_predictions
preds_cesium = model_predictions(fset_cesium, model_cesium, return_probs=False)
preds_guo = model_predictions(fset_guo, model_guo, return_probs=False)
preds_dwt = model_predictions(fset_dwt, model_dwt, return_probs=False)
print("Built-in cesium features: training accuracy={:.2%}, test accuracy={:.2%}".format(
accuracy_score(preds_cesium[train], eeg["classes"][train]),
accuracy_score(preds_cesium[test], eeg["classes"][test])))
print("Guo et al. features: training accuracy={:.2%}, test accuracy={:.2%}".format(
accuracy_score(preds_guo[train], eeg["classes"][train]),
accuracy_score(preds_guo[test], eeg["classes"][test])))
print("Wavelet transform features: training accuracy={:.2%}, test accuracy={:.2%}".format(
accuracy_score(preds_dwt[train], eeg["classes"][train]),
accuracy_score(preds_dwt[test], eeg["classes"][test])))
Built-in cesium features: training accuracy=100.00%, test accuracy=83.20%
Guo et al. features: training accuracy=90.93%, test accuracy=84.80%
Wavelet transform features: training accuracy=100.00%, test accuracy=95.20%
The workflow presented here is intentionally simplistic and omits many important steps
such as feature selection, model parameter selection, etc., which may all be
incorporated just as they would for any other scikit-learn analysis.
But with essentially three function calls (featurize_time_series,
build_model_from_featureset, and model_predictions), we are able to build a
model from a set of time series and make predictions on new, unlabeled data. In
upcoming posts we’ll introduce the web frontend for cesium and describe how
the same analysis can be performed in a browser with no setup or coding required.