Hybrid scattering-LSTM networks for automated detection of sleep arousals

Historically, most published methods for detecting transients in biomedical signals are composed of two stages: domain-specific feature extraction and general-purpose machine learning. In the context of polysomnography, features are engineered either in the time domain, the frequency domain, or the time-scale domains. Furthermore, the most widespread machine learning methods are decision trees (De Carli et al 1999, Agarwal 2006, Shmiel et al 2009), multi-layer perceptrons (MLP) (Huupponen et al 1996), and support vector machines (SVM) (Cho et al 2006).

The design of handcrafted features and finding their optimal combination for improving classifier performance can be difficult and time-consuming; to overcome this issue in recent years, the feature-engineering step is often automated using deep neural networks (DNN) such as CNNs (Tsinalis et al 2016, Aggarwal et al 2018, Zhang and Wu 2018).

De Carli et al (1999) recorded EEG (F4-C4 and C4-O2 channels) and a chin EMG to develop their arousal detector, using the American Academy of Sleep Medicine criteria defining all types of arousal (AASM 1992). Their dataset was relatively small, consisting of eleven overnight recordings of patients with various pathological conditions. They established their ground truth from an aggregation of two human experts. Then, they compared the arousal detection of their linear-discriminant classifier to the markings of the same two human experts on independent test data. They found that the sensitivity was higher (88.1% versus 72.4% and 78.4%) while the precision was lower (74.5% versus 83.0% and 82.0%). Since ground truth was derived from the marking of these same humans, the reported human scores are possibly too optimistic to reflect a real-world use case. Despite this caveat and the small dataset size, the results of De Carli et al (1999) offer a useful insight that automating the detection of sleep arousals has great application potential, yet requires the development of advanced techniques in signal processing and machine learning.

Although numerous publications have proposed to apply deep learning to sleep stage classification (wakefulness, stage 1, stage 2, stage 3 and REM) (Tsinalis et al 2016, Aggarwal et al 2018, Zhang and Wu 2018), fewer works have applied them to sleep arousal data. Recently, Saeed et al (2017) explored the potential of CNN for classifying sleep arousals into three categories: under-aroused, normal, and over-aroused. They measured raw physiological signals collected from wrist-wearable devices, as worn by eleven subjects. They found that their CNN outperformed a non-convolutional baseline, with respective F-scores of 0.82 and 0.75.

Our approach submitted to the 2018 Challenge was rated sixth out of 19 teams during the official phase. These top-six approaches all used DNNs, and four of these, including ours, also used an ensemble method. Therefore, this section is divided in two subsections: the first one concerns approaches that used a single deep model, while the second focuses on those approaches that used an ensemble method. All performance measures refer to the area under the precision-recall curve (AUPRC) of the algorithm on the hidden test set. Table 1 shows the various studies conducted on the automated detection of arousal regions in PSG signals, along with the results yielded by the present study which are discussed in section 4.

2.2.1. Deep architecture

2.2.2. Deep ensemble classifier

The development and evaluation of our proposed approach was based on the 2018 Challenge dataset. It consists of 1985 subjects monitored at the Massachusetts General Hospital (MGH) sleep laboratory for the diagnosis of sleep disorders. The data were partitioned into training (n  =  994) and hidden test sets (n  =  989). For each subject, 13 different physiological signals were collected during PSG sleep studies. The PSG recordings include six channels of EEG (F3-M2, F4-M1, C3-M2, C4-M1, O1-M2 and O2-M1), left eye EOG, an EMG lead placed at the chin (Chin1-2), respiratory movements placed on the chest and abdomen (CHEST and ABD), a single lead of ECG, SaO₂ and AIRFLOW. The sampling rate for all signals and labels was 200 Hz, sufficient to capture most of the PSG signal frequency bands mentioned in the introduction.

The training data was provided with five-class sleep-state labels and thirteen-class arousal-type labels. The sleep stages were annotated by clinical staff at the MGH according to the American Academy of Sleep Medicine manual for the scoring of sleep. The following six sleep stages were annotated in 30 s contiguous intervals: wakefulness, stage 1, stage 2, stage 3, rapid eye movement (REM), and undefined. Certified sleep technologists at the MGH annotated waveforms for the presence of arousals that interrupted the sleep of the subjects. The annotated arousals were classified as either: spontaneous arousals, respiratory effort related arousals (RERA), bruxisms, hypoventilations, hypopneas, apneas (central, obstructive and mixed), vocalizations, snores, periodic leg movements, Cheyne–Stokes breathing or partial airway obstructions.

Target arousals were then defined as intervals from 2 s before a RERA arousal begins, up to 10 s after it ends or from 2 s before a non-RERA, non-apnea arousal begins, up to 2 s after it ends. Accordingly, the records for the training data were provided with three-class labels: arousal (A, +1), non-arousal (NA, 0) and unscored (NS, −1) regions.

The main architecture of the base deep classifier shown in figure 1 consists of two main components; representation learning and sequencing learning, which are described below.

Zoom In Zoom Out Reset image size

Our solution is designed in the spirit of full representation learning, i.e. learning directly from the raw data without any preprocessing. But because the training data corpus was 134 GB in size and we wanted a working solution on 'modest' computing resources (e.g. a system with 32 GB RAM and a single GPU with 12 GB of memory), we needed to make some tradeoffs. In the extreme, full representation learning from 200 Hz data is conceptually possible, but because the recordings are long (approximately 5M samples), training time would be onerous. Therefore our original motivation to use ST, a fixed nonlinear operation, was to provide a data reduction step before learning. LSTM training time is also a function of sequence length, so this first decimation step by a factor of 512 to a sampling rate of 0.5 Hz for the detection function was critical to a working solution.

The scattering transform is a multidimensional, nonlinear function of real-valued signals (Bruna and Mallat 2013). In the past, this function has found practical applications in two different domains of signal processing: time series forecasting (Andreux 2018) and signal classification (Andén et al 2015). Interestingly, the task of arousal detection in PSG recordings lies at the interaction between these two perspectives. Indeed, sleep arousal can either be regarded, first, as a disruption in the predictability of the time series; or, alternatively, as a nonstationary pattern emerging from a stationary background. From this observation, we give hereafter two intuitive justifications of the interest behind adopting the scattering transform as a frontend to our deep learning system for sleep arousal detection. We refer to Lostanlen (2017) for further details in the argumentation.

First, in the context of time series forecasting, the scattering transform aims to encode, at every time step, past scalar values of a given sequence, and embed them into a feature vector, that is, another time series of fixed dimensionality. The purpose of this feature vector is to extract trends in the past values so that the prediction of the next few future values is relatively simple. In this respect, there is empirical evidence that the scattering transform reaches a favorable tradeoff between simplicity and accuracy: on a large class of real-world time series, the regression of scattering coefficients has a relatively high forecast accuracy, yet a relatively low intrinsic dimensionality (Andreux 2018).

In the specific case of the scattering representation, the feature space aggregates two subspaces, hereafter called 'orders'. The first order, which is comparable to a discrete wavelet transform (DWT), decomposes past trends into 11 typical time scales of periodicity which grow by factors of 2 from 10 ms to 10 s. As a result, if the stationary process underlying the observed time series were to consist of two independent linear trends of comparable periodicity, these would interfere within the same first-order subband. Resolving this interference is precisely the purpose of the second layer (Chudáček et al 2014). More precisely, second-order scattering decomposes, within each first-order subband, the temporal alternation between in-phase, constructive interference and out-of-phase, destructive interference, again into multiple time scales of periodicity. These second-order time scales are systematically greater than the 'baseband', first-order time scales. Therefore, the number of pairs of time scales that lead to a non-negligible outcome grows logarithmically with the second-order time scale. In other words, in consideration of the dimensionality of its feature vector, the scattering representation encodes its recent history with finer details, and its distant history with coarser details. This procedure of adaptation temporal resolution as a function of time lag is known in wavelet theory as foveation (Chang et al 2000).

Secondly, in the context of signal classification, the scattering transform performs a mapping of raw signal samples which reduces geometric, intra-class variability but retains informative, inter-class variability (Mallat 2016). In the case of PSG, the presence of small temporal lags and phase offsets between modalities is an example of such uninformative intra-class variability. Indeed, although PSG samples all modalities at a common rate of 200 Hz, spurious factors of experimental acquisition may cause a slight asynchrony between any two of the 13 channels. Likewise, swapping a pair of EEG electrodes leads to an inversion of electrical polarity, and thus a half-wave phase shift. None of these random modifications of data acquisition across different trials and subjects affect the task of interest. Therefore, it so appears that engineering a mid-level representation which, by design, remains invariant to the action of channel-dependent phase shifts is beneficial to statistical robustness and generalization.

Given a time series of low-level, short-term features, the simplest way to ensure local invariance to time shifts is to apply a moving average. However, this local averaging procedure comes at the expense of losing fine details in highly oscillatory data. The rationale behind the scattering transform is, instead of building a linear invariant to translation in a single step, to perform multiple nonlinear operations of continuous wavelet transform and pointwise complex modulus, so as to demodulate the input signal into progressively slower modulation subbands, hereafter called scattering 'paths' (Mallat 2011). In the case of PSG, we apply two of such nonlinear wavelet modulus operators in a cascade before eventually convolving each path with a Gaussian window function. For all modalities, we set the time constant of this Gaussian window equal to 1 s. We used the open-source MATLAB library scattering.m (Lostanlen 2019) to compute the scattering coefficients.

In recent years, various authors in have proposed to train CNN models on scattering transform features for audio signal processing, thus forming hybrid ST-CNN pipelines (Peddinti et al 2014, Fousek et al 2015, Zeghidour et al 2016, Andreux and Mallat 2018). Conversely, there is a growing amount of literature on biomedical applications of the scattering transform, in conjunction with shallow learning methods (Chudáček et al 2014, Leonarduzzi et al 2018). Yet, prior to Warrick and Homsi (2018b), there had not been an empirical study on the applicability of deep learning from scattering coefficients of biomedical signals. Furthermore, there had not been any deep learning architecture, in any application domain, which trained long short-term memory networks (LSTM) on scattering coefficients. We refer to Oyallon et al (2017) for a review of the state of the art in deep hybrid networks featuring both wavelet convolutions and trainable convolutions.

Another dimension in which our work is novel is the fact that we apply the scattering transform on multimodal data. Although one prior study has applied scattering transforms to stereophonic audio (Lostanlen and Andén 2016), the question of training a machine learning model on the scattering representations of various sources of physiological data in simultaneity has remained considerably under-studied. In particular, the resort to depthwise-separable convolutions as an interface between multiple single-channel scattering coefficients with one multichannel LSTM network is a novel architectural contribution of ours.

Because the generated ST coefficients had long exponential right tails, we performed a normalization step to transform them to quasi-normal distributions. This was done in two steps using the statistics of the training data for each ST path p  of each channel e. For each coefficient $x[i,j,e,p]$ where i is the recording and j  is the time sample, scaling by the median (calculated over all i and j ) was followed by a log-like transformation $\mathcal{K}$ as in (1):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eqn:Norm} {\rm Norm}(x[i,j,e,p])=\mathcal{K}\Big\{\frac{\mu \cdot x[i,j,e,p]}{{\rm median}(x[:,:,e,p])}\Big\}. \nonumber \end{align} \tag{ 1 }$$

In our 2018 submission we had used $\mathcal{K}(x)=\log(x+1)$. However, the Gaussian filter used to obtain the low frequency ST paths sometimes had transient negative values. We therefore used a modified transform $\mathcal{K}(x)=\sinh^{-1}(x)= \log[x + (x^2+1)]$ which admitted negative values. The constant $\mu$ was selected by searching for minimal skewness (i.e. quasi-normality). We first performed a search with coarse-grained resolution and then repeated the search with a finer resolution. Finally, we confirmed by visual inspection that the distributions were close to normal.

The arousal detection task from the 13 channel PSG recordings was carried out through two stages: building the base deep classifier and building the ensemble classifiers.

3.5.1. Depthwise-separable convolution

The ST output is hierarchically structured by channel $e\in\{1,...,E\}$ and ST path $p\in\{1,...,P\}$. In our 2018 Challenge entry we used J  =  8 octave ST giving eight first-order and ${8 \choose 2}=26$ second-order coefficients, such that P  =  36. Our 2018 Challenge entry flattened the ST output as input to the subsequent LSTM layer, ignoring this structure.

Lower-frequency information, down to 0.1 Hz, is referred to in the EEG literature as the 'slow-frequency' band (Hiltunen et al 2014). In this new work we wished to capture this information which required J  =  11. This increased P to 11  +  ${11 \choose 2}=66$. This near two-fold increase in P would significantly augment the DNN model weight count. Furthermore, we wished to experiment with third-order ST which would add another ${11 \choose 3}=165$ coefficients.

Therefore, we considered strategies to overcome the impact of the ST complexity on the model parameterization. We included this ST structure by factorizing the weights of the first LSTM layer into E electrodes $\times$ ST paths. Such an approach should also reduce rank, speed up training, and offer statistical regularization.

We chose a depthwise-separable convolution (DSC) (Chollet 2017) to address these issues. A DSC separates a standard convolution into two steps as shown in figure 2 for our configuration. We first perform the depthwise convolution X to mix the electrodes for each ST path with the $P \times E$ filter map F:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eqn:dsc} X[\,p] = \sum_{e} S[e,p]F[\,p,e]. \nonumber \end{align} \tag{ 2 }$$

Zoom In Zoom Out Reset image size

Then, the pointwise convolution mixes these transformed paths N times with $P \times N$ feature map G. N can be chosen arbitrarily and we choose N  =  P to allow at least one output for each path p . Output Y results from applying the bias B and activation function $\rho$:

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \label{eqn:dsc2} Y[n] = \rho\big(B[n] + \sum_{p} X[\,p]G[\,p,n]\big). \nonumber \end{align} \tag{ 3 }$$

The total number of convolution coefficients including the bias weights is therefore $P\times E + (P+1) \times N$. For N  =  P  =  66 ST paths per channel, the DSC weight count is 5280. The main impact in overall weight count should occur in the subsequent BiLSTM layer. LSTM weight counts can be implementation-specific, so we tested with our Keras-Tensorflow subsystem for a precise comparison. We compared the weight counts of direct ST-BiLSTM layers to ST-DSC-BiLSTM layers using N  =  P  =  66 and 100 cells per BiLSTM in each case. The former had 456 000 weights in the BiLSTM layer while the latter had 134 400. Therefore, inserting the DSC layer between the ST and BiLSTM layers had the effect of reducing the weight count by 69% for this architecture.

3.5.2. Sequence learning

3.5.3. Class imbalance

3.5.4. Auxiliary targets

3.5.5. Ensemble classifier

The proposed classifier was evaluated with two metrics that give complementary insights into the classification performance: AUPRC and area under the receiver operating characteristic (AUROC). We relied on the Python code provided by the 2018 Challenge to compute these metrics, which calculates histograms to estimate scaled versions of the two conditional prediction probabilities $Prob[arousal|(truth={{\rm arousal}})]$ and $Prob[arousal|(truth={{\rm non\mbox{-}arousal}}))$. These histograms are used to determine the false-positive rates, recall and precision values for each probability threshold, taken from corresponding histogram bins; we used the default setting of 1000 histogram bins. These allow us to generate the ROC and PR curves and estimate their areas.

For performance reasons, we compared targets and predictions at the 0.5 Hz decision sampling rate in all our experiments. For final evaluation on the Challenge server, we up-sampled the predictions by sample and hold to the expected 200 Hz rate.

As a baseline comparison of ST processing with conventional spectral methods, we also performed an experiment using power spectral density (PSD) analysis. This was done by first detrending the 200 Hz signals with a FIR high-pass filter having a cutoff frequency of 1 Hz. Then we performed sliding-window analysis with a window size of 512 samples that was advanced by the same amount, thus conforming to the decision sampling rate of the ST analysis. For each window we created a model using the matlab function $\textsf{pcov}$ using an FFT length of 64, generating approximately the same number of coefficients (33) as the ST36 processing, A fixed AR model order of 20 was chosen empirically. The frequency resolution provided by the PSD coefficients was therefore 200 Hz / 64  =  3.1 Hz. For the very slow-changing SaO₂ signal, the above processing was inappropriate, and the coefficient was set to the window average. We performed the normalization as before using the $\log(x+1)$ transform. Finally, we trained a detector in cross-validation using a three-layer BiLSTM for fair comparison.

Differences between experiment results were often subtle (within a few percent mean AUPRC and AUROC), so we needed a principled way to compare two experiments. We did this by performing a 2-sided t-test of test AUPRC over all ten folds to test the null hypothesis that the mean AUPRCs were equal.

Figure 3(a) shows the normalized ST coefficients for all channels for a typical recording. Figure 3(b) is a magnified view that includes a region of arousals. It is apparent from these overall and magnified views that the onset and termination of the arousals appear as higher ST values (i.e. towards red on the color maps). It is also clear that there are numerous elevated values where there are no target arousals. These clearly appear itermittently in non-scored (noisy) regions (e.g. between sample 3.4 $\times$ 10⁶ and 4.4 $\times$ 10⁶) but they may also indicate sensitivity to non-target arousals. Finally it is apparent from figure 3(b) that the time support of the second-order coefficients (near the bottom of each channel map) is wider than those of the first-order coefficients (at the top of the map), and especially with decreasing frequency band (lower in the map).

Zoom In Zoom Out Reset image size

Figures 4 and 5 show the DSC depthwise and pointwise filter-weight maps, respectively, for an example fold of the cross-validation. They are presented as 11-by-11 lower-triangular images with first-order paths along the diagonal and second-order paths beneath; table 2 gives the Morlet filter octave(s) associated with each path in the matrix.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For the depthwise convolution, all channels have non-zero weights for at least some paths; however, the four channels ABD, CHEST, AIRFLOW and SaO₂ seem to dominate over the seven EEG channels, Chin1-2 and ECG: their weights tends to more towards the red (positive) and blue (negative) extremes of the colormap. While these four channels also had non-zero weights in the higher octaves (towards the upper left), their lower octave weights (toward the lower right) were especially strong.

AIRFLOW, for example had strong first-order 0.097 Hz (intense red) and 0.39 Hz (intense blue) weights as well as strong second-order weights for the 12.5/25 Hz (red) and 1.56/25 Hz (red) and 0.781/3.125 Hz (blue) octave pairs. SaO₂ had higher-octave contributions in addition to its strong lower octave contributions, possibly reflecting sensitivity to sensor disturbance.

Weights were diversely spread across the seven EEG channels. In particular, lower octaves had significant weights over several electrodes: the first-order 0.197 Hz octave was present in E1-M2, and the second order 0.097/0.197 octave pair was present for F3-M2 and F4-M1 and to a lesser degree, C3-M2. ECG was especially strong for the 0.39/6.25 Hz (red) and 0.781/25 Hz (blue) octave pairs. For the Chin1-2 EMG signal, the first-order 25 Hz (red) and second order 0.39/50 Hz octave pair (blue) stood out.

For brevity we show only 12 of the P  =  66 mixtures of the pointwise filter in figure 5. These maps are more difficult to interpret, with strong peaks (red) and valleys (blue) occuring in many of the first and second order paths, and with much diversity between mixes.

Table 4 shows the results of our classifier training experimentation. The models were compared incrementally to assess the impact of the introduction of each change. We used a slightly modified version of our 2018 Challenge entry for ST36-UniLSTM (Exp. 1). It consisted of second-order ST cofficients spanning J  =  8 octaves (P  =  36 ST paths per channel), three unidirectional LSTM layers, each followed by a BN layer, and a dense-softmax layer. Unlike the Challenge entry, however, we used all three arousal targets {A, NA, NS} rather than just two {A, NA} to allow more noise-related expressiveness in the model. The mean number of training epochs at convergence was 140, with training time on the order of 16 hrs per fold. The mean $\pm$ standard deviation of the fold AUPRC was 0.2913 $\pm$ 0.0398.

Next, in Exp. 2 we replaced the LSTM units with BiLSTM. This time the AUPRC was 0.3753 $\pm$ 0.0398, an increase of 28.8% over baseline (p   =  0.0002). Training time was similar to baseline, but the mean epochs at minimum loss reduced to 112.6.

The baseline PSD-BiLSTM detector of Exp. 3 was inferior to the equivalent ST36-BiLSTM detector of Exp. 2, with AUPRC decreasing from 0.3753 to 0.2839 (p   <  0.0001).

We then augmented the ST analysis to capture the 'slow-frequency' band (0.1 Hz), increasing P to 66 (Exp. 4). The AUPRC was 0.4074 $\pm$ 0.0502, an increase of 8.5% over the Exp. 2 (p   =  0.0141). The rate of training convergence was similar.

Then a DSC layer with N  =  P  =  66 and linear activation was inserted between the ST and first BiLSTM layers. The change in AUPRC was insignificant (Exp. 5, p   =  0.9078); however, when the DSC activation was changed to a ReLU unit and a BN layer was added to the input, the AUPRC increased to 0.4344 $\pm$ 0.0369, an increase of 6.6% (Exp. 6, p   =  0.0395). As well, addition of the DSC reduced training time considerably in both cases; in the latter case, it reduced to a mean of 24 epochs at minimum loss and 10 hrs per fold. The rationale for adding the initial BN layer was that while normalization successfully transformed the ST coefficients to a quasi-normal distribution, the mean was considerably offset (approximately in the range of 6–10). Incorporating mean removal into the normalization step would have been the ideal approach, but to avoid renormalizing the data before training, we used the BN layer to accomplish the same task during training.

Removing the increased loss weighting on arousal regions by setting this value to 1 increased the AUPRC to 0.4547 $\pm$ 0.0308, an increase of 4.7% (Exp. 7, p   =  0.0391). Convergence was slightly longer, however, with mean epochs at minimum loss of 31.

As a further study in model ablation, we reduced the number of BiLSTM layers from three to two. The AUPRC did not change significantly (Exp. 8, p   =  0.8121), but convergence was slower, with the mean epoch at minimum loss increasing slightly to 42. We retained this simpler model in subsequent experiments.

Finally, we added the auxiliary target of arousal type to the network architecture, giving it equal weight in the loss function to the primary arousal target. This increased AUPRC to 0.4675 $\pm$ 0.0444, an increase of 2.4% (Exp. 9, p   =  0.0741). The mean epoch at minimum loss increased to 50, likely due mostly to the additional calculation of gradients for this 15-class, one-hot auxiliary target. Training and validation AUPRC were 0.6092 $\pm$ 0.0246 and 0.4608 $\pm$ 0.0175, respectively. The test, training and validation AUROCs were 0.9180 $\pm$ 0.0066, 0.9500 $\pm$ 0.0047 and 0.9230 $\pm$ 0.0050, respectively.

We considered this to be our best architecture. We used the models generated from 10-fold cross-validation in our proposed ensemble. In a 2018 Challenge follow-up entry, this ensemble scored 0.50 on the hidden test set.

Other experiments that did not significantly change performance included adding third-order ST coefficients, inserting 20% dropout layers before or after the DSC layer and using the sleep stages as auxiliary targets.

Figure 6 shows sample receiver operating characteristic (ROC) curves and precision-recall (PR) curves using a single fold of our final architecture. Since our final classifier model had relatively small variance (standard deviation of 0.0444) across the ten folds of cross-validation, as shown in table 4, we considered the ROC and PR curves of one of these folds (fold 1) to be representative; we chose a decision probability threshold (p   =  0.17) based on an example criteria of maximum F1 (0.48) using the validation data set. F1 is the harmonic mean of precision and recall, defined as $2\cdot\frac{precision \cdot recall}{precision + recall}$. Applying this threshold to the test data resulted in a false-positive rate of 0.06, a recall of 0.61 and a precision of 0.44. Table 3 tabulates these performance measures on the validation data as a function of probability threshold at 0.05 recall intervals.

Zoom In Zoom Out Reset image size