Sleep stage classification with ECG and respiratory effort

The data set was the same as used in earlier work (Long et al 2014d, 2014) and comprised full single-night polysomnographic (PSG) recordings of 48 subjects (27 females) acquired in the SIESTA project (Klosch et al 2001). All subjects were healthy sleepers with a Pittsburgh Sleep Quality Index (Buysse et al 1989) of less than 6 and had no regular sleep complaints nor earlier diagnosis of sleep disorders. The subjects had an average age of $41.3\left(\pm 16.1\right)$ years at the time of the recording. Full subject demographics can be found in our earlier work (Long et al 2014d). Sleep stages were scored by trained sleep technicians in six classes according to the R and K rules (Rechtschaffen and Kales 1968). In the scope of this study, S1 and S2 were merged in a single L (light sleep) class and S3 and S4 were merged in a single D (deep sleep) class. Each PSG recording comprised, besides the standard signals required for sleep scoring, modified lead II ECG, and (thoracic) respiratory effort recorded with respiratory inductance plethysmography (RIP). QRS complexes were detected and localized from ECG signals using a combination of a Hamilton–Tompkins detector (Hamilton and Tompkins 1986, Hamilton 2002) and a post-processing localization algorithm (Fonseca et al 2014). Prior to feature extraction, RIP signals were filtered with a 10th order Butterworth low-pass filter with a cut-off frequency of 0.6 Hz, after which baseline was removed by subtracting the median peak-to-through amplitude (Long et al 2014d).

We extracted a set of 142 features from cardiac and respiratory activity, and from cardiorespiratory interaction (CRI) using a sliding window centered on each 30 s epoch, guaranteeing sufficient data to capture the changes in autonomic activity (Malik et al 1996). Since some features are computed based on windows which exceed the epoch length, epochs at the start and end of each recording required a special handling: for each such feature, all epochs for which the window crosses the boundaries of the recording were marked as invalid; the feature values for these epochs were interpolated during post-processing using spline fitting (section 2.2.4).

2.2.1. Cardiac features.

2.2.2. Respiratory features.

2.2.3. Cardiorespiratory interaction features.

2.2.4. Feature post-processing.

In order to reduce the impact of physiological differences and equipment-related variations from subject to subject, the features of each subject were first Z-score normalized by subtracting their mean and dividing by their SD. Further, it is known that the sleep pattern of healthy adults progresses with several cycles throughout the night (Carskadon and Dement 2011). For example, REM and NREM sleep alternate with 4–6 cycles of about 90–110 min with deep sleep usually dominating the NREM periods during the first half of the night. This suggests that the autonomic physiological response with its associated sleep stage is time-variant across the night for each subject. For this reason, we were motivated to smoothen each feature for each subject by means of a cubic spline fitting method (De Boor 2001). This is also expected to help reduce signal measurement noise and variability within subjects for each sleep stage conveyed by the feature values. The latter can be caused by body movements, conscious breathing control, internal physiological variations, or other external factors such as changes in environmental noise and temperature during bedtime sleep. Instead of other simpler low-pass filters, spline fitting was chosen since it can interpolate feature values which could not be computed, for example due to motion artifacts (about 10% observed in our data set) or at the start and end of each recording for features computed with windows exceeding the epoch duration. This procedure allows all epochs in each recording to be classified.

Let t represent a sequence of feature values v   =   {v₁, v₂, ..., v_n} at their corresponding time (or epoch) indices t   =   {t₁,t₂, ...,t_n} (in 30 s), then a relation between them can be modeled by

$\begin{eqnarray}&&{{v}_{i}}=h\left({{t}_{i}}\right)+{{\varepsilon}_{i}}\ \left(i=1,2,...,n\right),\end{eqnarray} \tag{ 1 }$

where h is a smoothing (spline) function, ${{\varepsilon}_{i}}$ are independent and identically distributed residuals. The smoothing function can be estimated by minimizing the objective function (i.e. penalized sum of square) such that

$\begin{eqnarray}&&\hat{h}=\text{arg}\underset{h}{\mathop{\min}}\,\left[\underset{i=1}{\overset{n}{\mathop \sum}}\,{{\left[{{v}_{i}}-h\left({{t}_{i}}\right)\right]}^{2}}+\lambda \int_{{{t}_{1}}}^{{{t}_{n}}}{{h}^{\prime \prime}}{{(t)}^{2}}\text{d}t\right],\end{eqnarray} \tag{ 2 }$

where $\lambda$ is a smoothing parameter that controls the trade-off between residual and local variation. The smoothing function can be expressed by cubic B-splines as basis functions and determined via least squares approximation (Unser 1999, De Boor 2001).

For a specific overnight recording with a total of m epochs, it is divided in s continuous segments ($s=\lceil m/n\rceil$), designated as smoothing splines. Each segment can then be modeled by the spline function, yielding a general spline fitting for the epochs over the entire recording. n represents the smoothing window size where a larger n translates to a smoother fitting curve. In this work, a window size of nine epochs for modeling splines was experimentally found to be appropriate for the task of sleep stage classification.

This work used a multi-class Bayesian linear discriminant with time-varying prior probabilities (Redmond et al 2007), similar to that used in previous work (Long et al 2014d). For each epoch, the selected class (D, L, R, or W) is the class ${{\omega}_{i}}$ that maximizes the posterior probability given an feature vector x (Duda et al 2000),

$\begin{eqnarray}&&{{\omega}_{i}}\left(\mathbf{x}\right)=\underset{i}{\mathop{\text{arg}\,\text{max}}}\,\left[{{\text{g}}_{i}}\left(\mathbf{x}\right)\right]\end{eqnarray} \tag{ 3 }$

with the the discriminant function g_i for each class given by

$\begin{eqnarray}&&{{\text{g}}_{i}}\left(\mathbf{x}\right)=-\frac{1}{2}{{\left(\mathbf{x}-{{\mu}_{i}}\right)}^{t}}{{\Sigma}^{-1}}\left(\mathbf{x}-{{\mu}_{i}}\right)+\ln \text{P}\left({{\omega}_{i}},t\right)\end{eqnarray} \tag{ 4 }$

where ${{\mu}_{i}}$ is the average feature vector for class i, $\Sigma$ is the pooled covariance matrix for all classes, and $\text{P}\left({{\omega}_{i}},t\right)$ is the prior probability for class i at time (since lights off) t. All parameters were estimated during training.

To select the final list of features we used a wrapper feature selection method based on sequential forward selection (SFS) (Whitney 1971) using as criterion the Cohen's kappa coefficient of agreement $\kappa$ (Cohen 1960) on the training set. This measure of agreement between the classification predictions and the ground-truth annotations is more adequate than traditional measures of accuracy for this problem since there is a strong imbalance between classes (L epochs, for example, account for more than 50% of all epochs in the data set) and this coefficient factors out chance agreement, compensating for class imbalance.

In many machine learning studies supervised feature selection is often applied on the entire data set, even if the training and validation are kept separate (for example using cross-validation). This common pitfall is known to introduce a bias in the evaluation of a classifier's performance, which will often be overestimated (Smialowski et al 2010). Although keeping a hold-out set for validation would solve this problem, the limited size of the data set would either mean that the model learning would be based on potentially insufficient examples, or that the classifier would be evaluated on a very small sample, potentially unrepresentative of the problem at hand. Instead, the feature selection procedure was executed by strictly separating, on an iterative procedure akin to cross-validation, the training and validation sets. For each iteration, unbound SFS was applied using as criteria the classification performance obtained in the training set of each iteration. The final number of features was chosen as the smallest number S that yield a certain percentage of the maximum training kappa obtained across all iterations. The final list of selected features was chosen as the S features most often selected during the process.

The discriminative power of selected features was evaluated with the absolute standardized mean distance (ASMD) between the feature values of two classes, computed as

$\begin{eqnarray}&&\text{ASMD}=\left| \frac{{{{\bar{x}}}_{1}}-{{{\bar{x}}}_{2}}}{\sigma}\right| \end{eqnarray} \tag{ 5 }$

where ${{\bar{x}}_{1}}$ and ${{\bar{x}}_{2}}$ are the sample means for class 1 and 2 and $\sigma$ is the pooled sample SD.

After feature selection is performed and the set of features is chosen, the classification results per subject were evaluated using a ten-fold cross-validation procedure. The kappa coefficient for all subjects in the data set as well as the average and pooled performance were then calculated. In addition to the kappa coefficient, we also computed the traditional metric accuracy, i.e. the percentage of correctly identified epochs. For kappa and accuracy, the results were computed both after pooling the predictions over all epochs of all subjects and after averaging the performance for each subject.

Figure 1 indicates the kappa coefficient obtained in the training set of each iteration of the feature selection procedure, for a varying number of features. As illustrated, the maximum average training performance is obtained for 105 features, with an average kappa of 0.58. Also clear in the figure, is a plateau in performance between 70 and 100 features. This suggests that the number of features can be greatly decreased without affecting the training performance. A small feature set is often desirable to prevent over-fitting to the training data, as long as it is not so small that the model cannot learn the characteristics of the problem.

Zoom In Zoom Out Reset image size

By choosing different operating points in figure 1, we can choose a smaller feature set at the expense of a reduction in training performance. By allowing a reduction of 1% in training performance (from the maximum kappa of 0.58–0.57), we can reduce the number of features from 150 to 80 features (a reduction of 23.8%). Below this number there is a statistically significant (after a Wilcoxon signed-rank test, with p  <  0.05) decrease in training performance and further decreasing the number of features will likely lead to a decrease in classification performance after cross-validation. Using as criteria the smallest number of features which does not cause a statistically significant decrease below the maximum training performance, we chose a total of S   =   80 features.

After ranking all features by the number of times they were selected during the iterative feature selection procedure and selecting the 80 features with the highest count, we found that all features in the final feature set were selected in at least 5 of the 10 iterations (with a mean count of 7.67) with 14 features having been selected in all 10 iterations. This illustrates the robustness of the modified SFS method: despite their simplicity and computational efficiency, sequential selection algorithms are known to suffer from a so-called 'nesting effect', potentially leading to sub-optimal feature sets (Pudil et al 1994). By iteratively performing several unbound SFS searches on different training sets and keeping only the features that are selected most often, this effect was reduced, as attested by the large number of iterations each feature in the final set was selected.

For brevity only the 14 features selected in all iterations will be discussed further. Table 1 indicates the discriminative power of each feature using the pooled ASMD. It was computed for each pair of classes after aggregating the feature values for all subjects and also the 90th percentile of the ASMD (in parenthesis) obtained for each feature, for all individual subjects. Pooled ASMD values below 0.5 were omitted and 90th percentile ASMD values below 1 were omitted.

The top features are clearly discriminative for different pairs of classes which helps explain the relatively large number of features selected. Additionally, it is interesting to observe that there is one feature (median likelihood ratio) which does not have a pooled ASMD above 0.5 for any class pair. However, its 90th percentile ASMD value is larger than 1 for the pairs D/W and L/W. This is a good example of a feature which is discriminative for only a subset of the subjects (at least 10%) but not for all subjects. The fact that it was selected in every single iteration using the wrapper method described in section 2.4 suggests that it is complementary to other features for some subjects, helping raise the overall training performance.

A note should be made regarding long-term cardiac features such as ${{\alpha}_{2}}$ or larger-scale MSE features. None of these features were part of the final set of features selected with our method. One possible explanation for this is related to the length of the time series used to compute them. The choice of window sizes (for MSE, 17 epochs, i.e. 8.5 min) represents a compromise between having as much data as possible to accurately calculate the features, while at the same time not exceeding the average length of a given sleep stage (in our data set, the average length of deep and REM sleep periods was found to be 5.1 and 8.7 min, respectively). Although theoretically the window sizes we used are sufficient to calculate these features (a window of 17 epochs corresponds, at an average heart rate of 60 bpm, to 510 samples), it has been suggested that the estimation of sample entropy is low in series shorter than 10^m (where m is the pattern length, in samples) (Richman and Moorman 2000, Yentes et al 2013). This means that for a pattern length of m   =   2 and scales higher than 5, the coarse-grained time series used to calculate the sample entropy will have less data points than the suggested limit and the features might not be accurate, and therefore, not representative of the autonomic characteristics of different sleep stages.

Table 2 indicates the overall classification performance obtained after 10-fold cross-validation using the selected set of 80 features. In addition, it indicates the performance per class, obtained by considering each class as the positive class and merging the remaining in a single negative class. The highest performance is obtained for R detection, followed by W. The lowest performance is obtained for L. This is further confirmed by the confusion matrix of table 3 which shows that the largest proportion of errors occurs when trying to distinguish L from the other classes. For all other classes, the percentage of misclassified epochs (relative to the total number of epochs) is below 1% except for L.

Table 2. Cross-validation performance for 3 and 4 classes.

	Pooled kappa	Pooled acc.	Mean kappa	Mean acc.
WRLD	0.49	0.69	$0.49\pm 0.13$	$0.69\pm 0.08$
D	0.51	0.89	$0.50\pm 0.17$	$0.89\pm 0.04$
L	0.40	0.71	$0.41\pm 0.14$	$0.71\pm 0.07$
R	0.57	0.87	$0.58\pm 0.19$	$0.87\pm 0.08$
W	0.54	0.91	$0.51\pm 0.18$	$0.91\pm 0.04$
WRN	0.56	0.80	$0.56\pm 0.15$	$0.80\pm 0.08$

Note: the pooled performance was computed after aggregating all epochs of all subjects. The mean and SD were calculated based on the performance for each individual subject

Table 3. Confusion matrix after cross-validation.

Pred.$\downarrow$ ref.$\to$	D		L		R		W
D	3431	(7.6%)	1949	(4.3%)	5	(0.0%)	97	(0.2%)
L	2969	(6.6%)	19165	(42.6%)	2947	(6.5%)	2302	(5.1%)
R	86	(0.2%)	2071	(4.6%)	5383	(12.0%)	404	(0.9%)
W	31	(0.1%)	952	(2.1%)	243	(0.5%)	2996	(6.7%)

In order to evaluate the performance of the classifier in a three-class task (WRN), classes D and L were merged in a single N (non-REM) class. Table 2 indicates the resulting performance. Analyzing the performance of the classifier we see that the classification performance rises substantially, to a kappa of 0.56 and an accuracy of 80%. This was expected since a large number of classification errors occurred between D and L, and in a WNR task these two classes no longer need to be distinguished.

Figure 2 illustrates examples of predicted hypnograms, as compared with the reference, for three subjects in the data set: the subject with the worst performance (with $\kappa =0.17$), with the median performance (with $\kappa =0.50$) and with the best performance (with $\kappa =0.69$). A possible explanation for the poor performance obtained for the worst subject is that the model trained with the characteristics of the general sample population does not fully capture this subject's cardiac and respiratory expression of different sleep stages. However, despite the low kappa coefficient, the predicted hypnogram still exhibits some correct features, namely, most REM intervals were detected, albeit with the incorrect length, and the two deep sleep periods were also detected. As the performance improves, we see that the predicted hypnograms match better the characteristics of the reference hypnogram, and in the best case the most obvious mistakes are in the missed detection of brief periods of wake during the night while the rest of the sleep stages are correctly predicted. This is likely caused by the use of spline smoothing during feature post-processing, which is adequate to capture the slow-changing characteristics of most sleep stages, but penalizes short, abrupt changes such as brief periods of awakening.

Zoom In Zoom Out Reset image size

Table 4 compares the results of our work with other studies reported in literature. As indicated, only a few studies focused on WRLD classification based on cardiac and/or respiratory signals and our results are amongst the best performing. The first observation is that the results of our previous work (Long et al 2014), which used only respiratory features, are worse than those produced in the present work, indicating that combining cardiac and respiratory activity can lead to an improved classification performance. The study of Hedner et al (2011) achieved similar results but they used more signal modalities including peripheral arterial tone, actigraphy, and pulse oximetry. The recent study by Willemen et al (2014) also achieved a good performance, although it was validated with a younger sample population, excluded 12% of the epochs from validation and used a basis of 60 s epochs instead of the standard scoring basis of 30 s which makes the results incomparable.

Table 4. Performance comparison with state-of-the-art.

Reference	Modalitiesa	N	Age (year)	Average $\kappa$	Average accuracy
WRLD
Our work	RIP, ECG	$48$	$\mathbf{41}.\mathbf{3}\pm \mathbf{16}.\mathbf{1}$	$\mathbf{0}.\mathbf{49}$	$\mathbf{0}.\mathbf{69}$
Our previous workb	RIP	48	$41.3\pm 16.1$	0.41	0.65
Isa et al (2011)	ECG	16	$\left[32-56\right]$	0.26	0.60
Hedner et al (2011)	PAT, PO, ACT	227 OSA	$49\pm 14$	0.48	0.66
Willemen et al (2014)c	RIP, ECG, ACT	85 (36 subj.)	$22.1\pm 3.2$	0.56	0.69
WNR
Our work	RIP, ECG	$48$	$\mathbf{41}.\mathbf{3}\pm \mathbf{16}.\mathbf{1}$	$\mathbf{0}.\mathbf{56}$	$\mathbf{0}.\mathbf{80}$
Our previous workb	RIP	48	$41.3\pm 16.1$	0.48	0.77
Redmond and Heneghan (2006)	RIP, ECG	37 OSA	$46.7\pm 10.4$	0.32	0.67
Redmond et al (2007)	RIP, ECG	31	$42.0\pm 7.4$	0.45	0.76
Mendez et al (2010)	BCG	22 (11 subj.)	n.a.	0.42	0.72
Kortelainen et al (2010)	BCG	18 (9 subj.)	$\left[20-54\right]$	0.44	0.79
Migliorini et al (2010)	BCG	22 (11 subj.)	n.a.	0.55	0.77
Kurihara and Watanabe (2012)	BCG	20	22.2	0.48	0.78
Xiao et al (2013)	ECG	45	$\left[16-61\right]$	0.46	0.73
Domingues et al (2014)	RIP, ECG, ACT	20	$42.1\pm 9$	0.58	0.78
Willemen et al (2014)c	RIP, ECG, ACT	85 (36 subj.)	$22.1\pm 3.2$	0.62	0.81

^aECG: electrocardiography, RIP: respiratory inductance plethysmography, PAT: peripheral arterial tone, PO: pulse oximetry, ACT: actigraphy, BCG: ballistocardiography. ^bLong et al (2014). ^c60 s epochs; 12% of all epochs excluded from validation. Note: all data sets comprise healthy subjects unless indicated. There are other studies in literature which present subject-dependent classification results. In this comparison only results obtained with subject-independent schemes were considered.

For WRN classification we see that, to the best of our knowledge, our results also outperform those reported in almost all of the previous studies. In comparison with one of the best performing studies (Domingues et al 2014), we obtain a higher accuracy (albeit a slightly smaller kappa) but require one less modality (actigraphy). Regarding the work of Willemen et al (2014) it is again important to note that the results in that study were obtained on basis of 60 s epochs.

These results suggest that our choice for a Bayesian linear discriminant was appropriate for this task. Besides its simplicity, it offers the benefit of a probabilistic framework which allows, for instance, the direct use of time-varying prior probabilities to improve classification. In comparison with increasingly popular black-box approaches, this classifier has the additional advantage that it does not require the tuning of critical parameters such as kernels for support vector machines, number of nodes in classification trees or number of hidden layers in neural networks.