Wavelet analysis of oximetry recordings to assist in the automated detection of moderate-to-severe pediatric sleep apnea-hypopnea syndrome

Abstract

Background

The gold standard for pediatric sleep apnea hypopnea syndrome (SAHS) is overnight polysomnography, which has several limitations. Thus, simplified diagnosis techniques become necessary.

Objective

The aim of this study is twofold: (i) to analyze the blood oxygen saturation (SpO2) signal from nocturnal oximetry by means of features from the wavelet transform in order to characterize pediatric SAHS; (ii) to evaluate the usefulness of the extracted features to assist in the detection of pediatric SAHS.

Methods

981 SpO2 signals from children ranging 2–13 years of age were used. Discrete wavelet transform (DWT) was employed due to its suitability to deal with non-stationary signals as well as the ability to analyze the SAHS-related low frequency components of the SpO2 signal with high resolution. In addition, 3% oxygen desaturation index (ODI3), statistical moments and power spectral density (PSD) features were computed. Fast correlation-based filter was applied to select a feature subset. This subset fed three classifiers (logistic regression, support vector machines (SVM), and multilayer perceptron) trained to determine the presence of moderate-to-severe pediatric SAHS (apnea-hypopnea index cutoff ≥ 5 events per hour).

Results

The wavelet entropy and features computed in the D9 detail level of the DWT reached significant differences associated with the presence of SAHS. All the proposed classifiers fed with a selected feature subset composed of ODI3, statistical moments, PSD, and DWT features outperformed every single feature. SVM reached the highest performance. It achieved 84.0% accuracy (71.9% sensitivity, 91.1% specificity), outperforming state-of-the-art studies in the detection of moderate-to-severe SAHS using the SpO2 signal alone.

Conclusion

Wavelet analysis could be a reliable tool to analyze the oximetry signal in order to assist in the automated detection of moderate-to-severe pediatric SAHS. Hence, pediatric subjects suffering from moderate-to-severe SAHS could benefit from an accurate simplified screening test only using the SpO2 signal.

About the Authors:

Fernando Vaquerizo-Villar

Roles Investigation, Methodology, Software, Validation, Writing – original draft

* E-mail: [email protected]

Affiliation: Biomedical Engineering Group, Universidad de Valladolid, Valladolid, Spain

ORCID logo http://orcid.org/0000-0002-5898-2006

Daniel Álvarez

Roles Investigation, Methodology, Software, Supervision, Writing – original draft

Affiliations Biomedical Engineering Group, Universidad de Valladolid, Valladolid, Spain, Pneumology Service, Hospital Universitario Río Hortega, Valladolid, Spain

ORCID logo http://orcid.org/0000-0003-1027-2395

Leila Kheirandish-Gozal

Roles Investigation, Methodology, Validation, Writing – original draft

Affiliation: Department of Child Health, The University of Missouri School of Medicine, Columbia, Missouri, United States of America

Gonzalo C. Gutiérrez-Tobal

Roles Investigation, Methodology, Software, Validation, Writing – original draft

Affiliation: Biomedical Engineering Group, Universidad de Valladolid, Valladolid, Spain

Verónica Barroso-García

Roles Investigation, Methodology, Software, Validation, Writing – original draft

Affiliation: Biomedical Engineering Group, Universidad de Valladolid, Valladolid, Spain

Andrea Crespo

Roles Investigation, Methodology, Validation, Writing – original draft

Affiliations Biomedical Engineering Group, Universidad de Valladolid, Valladolid, Spain, Pneumology Service, Hospital Universitario Río Hortega, Valladolid, Spain

Félix del Campo

Roles Funding acquisition, Investigation, Methodology, Software, Validation, Writing – original draft

Affiliations Biomedical Engineering Group, Universidad de Valladolid, Valladolid, Spain, Pneumology Service, Hospital Universitario Río Hortega, Valladolid, Spain

ORCID logo http://orcid.org/0000-0002-4554-2167

David Gozal

Roles Investigation, Methodology, Project administration, Validation, Writing – original draft

Affiliation: Department of Child Health, The University of Missouri School of Medicine, Columbia, Missouri, United States of America

Roberto Hornero

Roles Funding acquisition, Investigation, Methodology, Project administration, Supervision, Writing – original draft

Affiliations Biomedical Engineering Group, Universidad de Valladolid, Valladolid, Spain, IMUVA, Instituto de Investigación en Matemáticas, Universidad de Valladolid, Valladolid, Spain

Introduction

The American Academy of Pediatrics (AAP) defines pediatric sleep apnea-hypopnea syndrome (SAHS) as a breathing disorder characterized by recurrent episodes of complete cessation (apnea) and/or significant reduction (hypopnea) of airflow during sleep [1]. SAHS is a highly prevalent condition among children (in the range of 1% to 5%) that may lead to many adverse consequences on the overall health and quality of life, such as cognitive deficits, behavioral abnormalities, sleepiness, systemic inflammation, and cardiac and metabolic derangements [1].

The gold standard technique for pediatric SAHS diagnosis is overnight polysomnography (PSG). It involves recording a wide range of biomedical signals in a specialized sleep laboratory [2,3]. These recordings are used to score apneas and hypopneas in order to compute the apnea-hypopnea index (AHI), defined as the number of apneas and hypopneas per hour (e/h) of sleep. AHI is the clinical variable used to establish SAHS. The diagnosis of moderate-to-severe pediatric SAHS is confirmed when they present an AHI ≥5 e/h, irrespective of other co-morbidities [1]. These children are at increased risk of suffering from the major negative consequences of the disease [3–5]. Thus, to expedite the diagnosis and treatment is essential in these patients. In this sense, surgical treatment with adenotonsillectomy is consistently recommended for children suffering from SAHS with an AHI ≥5 e/h [6]. This treatment leads to an improvement in the condition in the majority of pediatric patients who suffer from moderate-to-severe childhood SAHS [1]. However, in spite of the PSG serving as the current recommended diagnostic gold standard, it is costly and complex due to the necessary equipment and trained staff, as well as highly intrusive due to the use of multiple sensors. In addition, it is a time-demanding method that shows limited availability and absent scalability, thereby delaying the diagnosis and treatment of SAHS patients [7,8].

These drawbacks have led to extensive exploration of the use of simplified diagnostic techniques [9,10]. One common approach is the analysis of a reduced set of cardiorespiratory signals involved in PSG. In this regard, overnight oximetry is a common alternative due to its reliability, simplicity, and suitability for children [7,11]. Nocturnal oximetry records the blood oxygen saturation (SpO2) signal, which provides a numerical measure of the oxygen content in hemoglobin [12]. Apneic events result in decreases in blood oxygen levels and such events are termed oxyhemoglobin desaturations [12]. Hence, the SpO2 signal contains useful information to detect pediatric SAHS. Previous studies have shown the usefulness of automated analysis of the SpO2 signal from nocturnal oximetry to assist in the screening of moderate-to-severe pediatric SAHS [13–19]. However, the results obtained in these studies indicate that an accurate diagnosis of pediatric SAHS is difficult, and in fact, substantially more difficult than in adults, particularly because the frequency of apneic events and reductions in SpO2 is markedly lower in children. Thus, further scientific evidence is still necessary before the diagnostic ability of the SpO2 signal can be widely implemented as a pragmatic tool to assist in an automated detection of childhood SAHS.

Different signal processing techniques have already been applied to characterize the changes produced in the SpO2 signal as elicited by apneic events. Conventional oximetry indices, statistical measures, nonlinear parameters, and spectral analysis from the SpO2 recordings have all been evaluated [13–19]. Among these approaches, the use of spectral analysis is a common choice due to the recurrence of apneic events. In this sense, previous studies have assessed features extracted from the power spectral density (PSD) and bispectrum [17–19]. However, these methods are based on the Short-Time Fourier Transform (STFT), thus having a fixed time-frequency resolution [20]. In contrast, wavelet transform (WT) offers high frequency resolution at low frequencies as well as high time resolution at high frequencies [20,21]. This property makes WT a potentially more suitable technique to accurately detect low frequency components, such as those associated with the duration of SpO2 desaturations. Additionally, WT is also suitable to analyze non-stationarities like those occurring in the SpO2 signal by apnea-hypopnea events. In this sense, wavelet analysis has proven its usefulness to detect changes produced in biomedical signals by apneic events among adult SAHS patients [22–27]. Nevertheless, only two single preliminary studies by our group evaluated the usefulness of the wavelet analysis in the detection of pediatric SAHS using the SpO2 signal [28,29]. Therefore, additional research is clearly needed to further corroborate previous findings in a small cohort and to assess the usefulness of wavelet analysis of SpO2 in the diagnosis of pediatric SAHS. Thus, we propose to develop a more exhaustive wavelet analysis with a larger database of 981 overnight SpO2 recordings.

We hypothesized that the multiresolution analysis afforded by the WT could provide a set of useful features to precisely characterize changes occurring in the SpO2 signal associated with pediatric SAHS. Consequently, the aim of this study was twofold: (i) to analyze oximetry dynamics by means of WT-derived features in order to characterize differences associated with the presence of SAHS; and (ii) to assess the usefulness of these features to assist in an automated detection of moderate-to-severe pediatric SAHS.

Materials and methods

Subjects and signals under study

The database is composed of 981 pediatric subjects (602 males and 379 females) ranging from 2 to 10 years of age. All children were referred to the Pediatric Sleep Unit at the University of Chicago Medicine-Comer Children’s Hospital (Chicago, IL, USA) in the context of clinical suspicion of SAHS. All legal caretakers of the children gave their informed consent as a prerequisite to be part of the study and the Ethics Committee of the hospital approved the protocols (#11-0268-AM017, # 09-115-B-AM031, and # IRB14-1241).

Children’s sleep was monitored using a digital polysomnography system (Nihon Kohden America Inc., CA, USA). SpO2 recordings were acquired during overnight polysomnography at sampling rates of 25, 200, or 500 Hz. In a preprocessing stage, artifacts were removed by discarding those SpO2 values below 50% and those intervals with a slope higher than 4%/s [30]. Then, SpO2 recordings were resampled to a common rate of 25 Hz, as recommended by the American Academy of Sleep Medicine (AASM) [12], and were rounded to the second decimal place in order to have the same resolution [31]. The guidelines of the AASM were used by a certified pediatric sleep specialist to quantify sleep and cardiorespiratory events. The AHI was subsequently derived in order to diagnose pediatric SAHS. An AHI of 5 e/h was the threshold used to establish moderate-to-severe SAHS because of the enhanced risk of morbidity and thus the importance of an early detection and treatment in these cases. According to this AHI-based cutoff, 405 children were in the group AHI ≥5 e/h, whereas 576 children were in the group AHI <5 e/h.

The dataset was randomly divided into an optimization set (60%) and a cross-validation set (40%) [19]. Table 1 shows demographic and clinical data of the population under study (median [interquartile range] or n (%)). No statistically significant differences (p-value < .01) emerged in either age or body mass index (BMI) between optimization and cross-validation groups.

[Figure omitted. See PDF.]

Table 1. Demographic and clinical characteristics of the patient groups under study.

https://doi.org/10.1371/journal.pone.0208502.t001

Methods

Our methodology is divided into three steps: feature extraction, selection, and classification. In the first step, the wavelet transform was applied to analyze each SpO2 signal. A set of features was computed using the discrete wavelet transform (DWT) to characterize the changes produced in SpO2 recordings due to SAHS. In addition, 3% oxygen desaturation index (ODI3), statistical moments in the time domain and PSD features, which are common features from the SpO2 signal [17,19], were obtained to compose a wide initial feature set with relevant as well as complementary information. In the second step, a feature subset was selected using the fast correlation-based filter (FCBF) method [32]. Finally, binary logistic regression (LR) [33], support vector machines (SVM) [34] and multi-layer perceptron (MLP) neural network [35] classifiers were trained using this selected feature subset in order to detect moderate-to-severe pediatric SAHS.

Fig 1 shows the validation approach employed in each methodological step. The first set (optimization set) was employed to perform descriptive analysis of the extracted features, select a subset of features with FCBF, and select the optimal design parameters of the SVM and MLP classifiers. Bootstrapping has been employed in the feature selection stage, in order to avoid overfitting [36]. In the same way, 10-fold stratified cross validation has been applied to optimize the design parameters of SVM and MLP. The second set (cross-validation set) was used to evaluate the diagnostic performance of the single features and classifiers. Stratified K-fold cross validation (K = 5) was applied for this purpose [37].

[Figure omitted. See PDF.]

Fig 1. Validation approach employed in each methodological step of the study.

https://doi.org/10.1371/journal.pone.0208502.g001

Feature extraction.

Discrete Wavelet Transform

WT can be seen as the decomposition of a signal x(t) onto a set of basis functions, called wavelets [20]. Wavelets are obtained by time translations and scaling of a unique function called the mother wavelet. WT can be seen as an extension of the Fourier transform where, instead of analyzing a single scale, a multiscale analysis is performed. This multiscale property of the WT allows decomposing a signal into a set of scales, where each scale analyzes a different frequency range of the signal. WT can be continuous (Continuous Wavelet Transform, CWT) or discrete (DWT), depending on the scale and translation values [20]. CWT computes WT for each scale, whereas DWT only computes WT for dyadic (power of 2) scales, thus presenting lower complexity and higher computational efficiency than CWT [38]. Consequently, DWT was chosen in this study. In addition, it has previously shown its usefulness to detect different frequency components in physiological signals associated to SAHS events in adult patients [22–27].

Fig 2 shows how DWT is computed. In Fig 2A, the decomposition process of a SpO2 signal x[n] using DWT, the so-called subband coding scheme, is illustrated. It is a filter-bank tree where each stage consists of a high pass-filter g[n] (the mother wavelet) and a low pass filter h[n] (the mirror version of the mother wavelet), followed by a subsampling process of factor two [20]. The relationship between these two filters is as follows [20]:(1)where L, an even number, is the length of the filter. First, x[n] is decomposed in an approximation signal (lowpass version), A1, and a detail signal (highpass version), D1. Then, A1 is further decomposed into another approximation signal, A2, and another detail signal, D2. Each iteration increases the frequency resolution of the approximation and the detail version by two, as well as decreases the number of samples of both approximation and detail signals. This process continues until the maximum detail level of the signal, N = log2(M) is reached, being M the length of x[n] [39]. At each level (i = 1, 2, …, N), the approximation signal, Ai, and the detail signal, Di, can be computed as follows:(2)(3)where Ai-1 is the approximation signal in the level i-1. In the level 1, A0 is the original signal x[n]. Fig 2B shows an example of SpO2 signal, x[n], the detail signal Di obtained at each level i of the DWT decomposition, and the approximation signal AN obtained at the level N of the DWT decomposition.

[Figure omitted. See PDF.]

Fig 2. DWT computation.

(A) Decomposition process of a signal using DWT. (B) Original SpO2 signal, detail signals at each decomposition level and approximation signal at the maximum level of the decomposition.

https://doi.org/10.1371/journal.pone.0208502.g002

DWT was applied to the upper power of 2 for 5 minute segments (M = 213 samples (5.5 minutes)) and, consequently, N = 13 [23]. In this study, the Haar wavelet was chosen as mother wavelet. The reason is twofold [27]: (i) its suitability for picking up abrupt changes, which is appropriate to detect the changes produced in the SpO2 values due to apneic events; and (ii) its smoothing feature, which does not distort the original form of the SpO2 signal. At each level of the decomposition, detail coefficients contain information about a different frequency band, as stated in Fig 2B. We focused on the detail coefficients of the level 9 (D9, i.e., 0.0244–0.0488 Hz), since it is the level which is contained in the band of interest previously related to the recurrence of apneic events [18]. SpO2 signal presents both drops and rises associated to apneic events, which result in decreased and increased values in D9 coefficients, respectively. The information contained in the D9 coefficients may be canceled due to the presence of both positive and negative values, such as mean or skewness. To avoid this, the absolute values of the DWT coefficients were used. The following seven features were extracted from the DWT coefficients:

* Statistical moments of the D9 coefficients (Mean (M1D9), variance (M2D9), skewness (M3D9) and kurtosis (M4D9)). M1D9–M4D9 measure the central tendency, dispersion, asymmetry and peakedness of the data, respectively.

* Maximum amplitude of the D9 coefficients (MaxD9). It quantifies the highest amplitude in this frequency band.

* Energy of the D9 coefficients (EnD9). It measures the averaged quadratic amplitude of the signal in D9. It is computed as follows:

(4)

* Wavelet Entropy (WE), which measures the irregularity introduced in the DWT. It was extracted in order to obtain information about the changes produced in the energy distribution of the different detail levels of the DWT of the SpO2 signal by apneic events [39]. It is computed as follows:

(5)where pi is the relative wavelet energy at the detail level Di:(6)

Where EnDi is the wavelet energy at the detail level Di:(7)

Conventional features from the SpO2 signal

In order to enhance the diagnostic ability of our proposal, the following features, that are common parameters of the oximetry signal [17,19], were computed:

* ODI3. It was estimated as the number of desaturations of at least 3% from preceding baseline per hour of recording [40]. This parameter has shown its usefulness in clinical studies, even though it underestimates AHI [13–15].

* Statistical moments. First-to-fourth order statistical moments were computed from the SpO2 signal in the time domain (M1T -M4T): mean (M1T), variance (M2T), skewness (M3T), and kurtosis (M4T) [17,19]. These features measure the central tendency, dispersion, asymmetry, and peakedness of the data, respectively.

* PSD features. PSD was estimated using the Welch’s method (213-sample Hamming window, 50% overlap and 214-points DFT) [41]. The following features were obtained: first-to-fourth order statistical moments (M1PSD- M4PSD) and maximum amplitude (MaxPSD) from the band of interest determined in [18] (0.018–0.050 Hz) and spectral entropy (SEPSD) in the full spectrum. These features provide information about the recurrence and duration of apneic events.

Feature selection: Fast Correlation-Based Filter (FCBF).

The FCBF method was applied to select a subset of relevant and non-redundant features [32]. FCBF is a feature selection algorithm that has previously shown its usefulness in the context of pediatric SAHS [18,19]. First, FCBF computes the symmetrical uncertainty (SU) between each feature (xi) and the AHI (y). SU is a normalization of the information gain between two variables. SU is computed as follows [32]:(8)where IG(xi|y) = H(xi)—H(xi|y), N is the total number of features extracted and H refers to Shannon’s entropy [32]. According to their SU value (between 0 and 1), features are ranked from the most relevant (highest SU with the AHI) to the least relevant one (lowest SU with the AHI). Then, a redundancy analysis is performed. SU between each pair of features (xj, xi) is computed. Features xj sharing more information with a more relevant one than with the AHI (SU (xj| xi) ≥ SU(xj|y)) were discarded. Finally, an optimum subset composed of the features not discarded in this process is obtained.

A bootstrap approach was employed in order to obtain a subset of features independent of a particular dataset. In this regard, FCBF was applied to 1000 bootstrap replicates built with a sample with replacement procedure from the optimization set [42,43]. Those variables that were selected with FCBF more than 500 times (50% of runs) formed the feature subset [18,19].

Feature classification.

In this study, we employed LR, SVM, and MLP, which are well-known algorithms in the context of binary classification. Particularly, these algorithms were applied to assign each subject to the groups AHI <5 e/h and AHI ≥5 e/h [33–35].

Logistic regression

LR is a standard machine learning approach for binary classification. Given a set of input features, LR estimates the posterior probability of a given instance (subject) belonging to one of two mutually exclusive groups (AHI <5 e/h and AHI ≥5 e/h) by the use of the logistic function [33]:(9)where Cl represents the two groups (AHI <5 e/h and AHI ≥5 e/h), β = β0, β1, …, βN are the coefficients of the model for each input feature, xk = x1,k, …, xN,k, is the input pattern for the instance k, and N is the number of features. A Bernoulli distribution is used to model the probability density function and β coefficients are optimized using the maximum likelihood ratio [33].

Support vector machines

A SVM is a binary classifier that searches for the best hyperplane that separates instances from the classes under study [34]. The hyperplane has the following expression [34]:(10)where x ϵ ℝN is the input pattern of dimension N (number of features), ɸ(x) ϵ ℝP transforms the data into a high-dimensional space P>N, and w is the weight vector. The weight vector w is optimized in order to maximize the margin of separation between the two groups [34]. A regularization parameter C was applied to control the trade-off between maximizing the margin of separation between groups and obtaining a good generalization ability in an independent set [34]. The optimization problem of SVM is formulated using Lagrange multipliers:(11)where S is a subset of the indices {1, …, L} corresponding to the non-zero Lagrange multipliers (support vectors) ηi, L is the number of observations in the training set, ti are the output labels (±1 for the AHI ≥5 and AHI <5 e/h groups), and K(∙,∙) is the kernel function in the transformed space. In this study, a linear kernel was used, which has previously shown its usefulness in the context of adult SAHS [44]. The value of C was optimized by means of 10-fold cross-validation using the optimization set.

Multi-layer perceptron neural network

A MLP is an artificial neural network arranged in several fully connected layers: input, hidden, and output layers [35]. These layers are composed of computing units called perceptrons or neurons. Each neuron consists of an activation function gk{.} and adaptive weights wkj that interconnect the neuron with neurons from the subsequent layer [35]. The input layer was composed of one neuron for each input feature. Additionally, a configuration with one single hidden layer with a hyperbolic tangent activation function was applied since it provides a fast convergence for the training algorithm [35]. This configuration can provide universal approximation to any continuous function with the only condition that there are enough hidden units [35,45]. Finally, two neurons composed the output layer, since our problem is a binary classification task. A logistic sigmoid activation function has been used in the output layer, because it allows the output neurons to be interpreted probabilistically [35]:(12)where yk are the outputs neurons, wkj are the weights connecting the hidden layer to the output layer, wji are the weights connecting the input layer to the hidden layer, bj and bk are the bias associated to the hidden and the output units, respectively, xi is the feature i, gk{.} and gj{.} are the activation functions of the output and hidden layer, respectively, NH is the number of neurons in the hidden layer, and N is the number of input features [35]. Random initialization was performed for the weights of the network. Then, the scaled conjugate gradient algorithm with weight-decay regularization was used to optimize the weights [35]. NH and the regularization parameter (α) were optimized by means of 10-fold cross-validation using the optimization set.

Statistical analysis.

The software tools Matlab version R2017a was used for performing signal processing and statistical analyses. Normality and homoscedasticity tests showed that extracted parameters were not normality distributed and had different variances. Consequently, the Mann-Whitney U test was applied to search for statistical significant differences in the extracted features (p-value <0.01) between groups. Diagnostic performance was assessed by means of sensitivity (Se, percentage of patients with an AHI ≥5 e/h correctly classified), specificity (Sp, percentage of children with an AHI <5 e/h correctly classified), positive predictive value (PPV, proportion of subjects classified as positive that are true positives), negative predictive value (NPV, proportion of subjects classified as negative that are true negatives), positive likelihood ratio (LR+, likelihood ratio for subjects classified as positive), negative likelihood ratio (LR-, likelihood ratio for subjects classified as negative), and accuracy (Acc, percentage of subjects correctly classified).

K-fold stratified cross validation (K = 5) was applied to assess the performance of the extracted features and the binary classifiers [37]. The cross-validation set was randomly divided into K subsets, preserving the proportion of subjects belonging to the groups AHI <5 e/h and AHI ≥5 e/h. K-1 folds formed the training folds (80% of the cross-validation set), whereas the remaining one formed the test fold (20% of the cross-validation set). Accordingly, Receiver Operating Characteristics (ROC) curves were used to obtain optimum classification cut-off points for the single features using the K-1 training folds. Similarly, the classification algorithms were trained using the training folds. Then, the diagnostic performance of the single features and the LR, SVM, and MLP classifiers was measured using the test fold. This process was repeated K times, so each fold was considered once as the test fold. Finally, all the metrics are averaged across the K = 5 iterations.

Results

Feature separability

A total of seven DWT-derived features were obtained for each SpO2 recording (S1 Table). Fig 3 shows the histogram of the D9 coefficients in the optimization set for the groups AHI <5 e/h and AHI ≥5 e/h. It can be observed that D9 coefficients are more concentrated near zero in the AHI <5 e/h group, whereas in the group AHI ≥5 e/h these coefficients are more disperse. Table 2 shows the median and interquartile range of all these extracted features in the optimization set for both groups. All features showed significant statistical differences (p-value <0.01) between groups. M1D9, M2D9, MaxD9, EnD9, and WE showed higher values in the AHI ≥5 e/h group, whereas M3D9 and M4D9 showed higher values in the AHI <5 e/h group. ODI3, statistical moments and PSD features were also computed for each SpO2 recording (S1 Table).ODI3, 3 out of 4 statistical moments (M1T, M2T, and M3T) and 3 out of 6 spectral features (M1PSD, M2PSD, and MaxPSD) also showed significant statistical differences (p-value <0.01), which agrees with previous studies [17,18].

[Figure omitted. See PDF.]

Fig 3. Histogram of the D9 coefficients for each group in the optimization set.

https://doi.org/10.1371/journal.pone.0208502.g003

[Figure omitted. See PDF.]

Table 2. DWT-derived features for each group in the optimization set.

https://doi.org/10.1371/journal.pone.0208502.t002

Optimum feature subset

FCBF was applied to each bootstrap replicate from the optimization set, each one composed of all the extracted features (ODI3, statistical moments, PSD, and DWT features). ODI3, 1 statistical moment (M2T), 3 features from PSD (M2PSD, M3PSD, and MaxPSD), and 3 DWT-derived features (M3D9, EnD9 and WE) were selected more than 50% of times (500) (S2 Table). Thus, these features formed the selected feature subset [18,19]. Notice that features from all the different methodological approaches were selected.

Classification models optimization

LR, SVM, and MLP classifiers were designed using the selected feature subset obtained with FCBF (ODI3, M2T, M2PSD, M3PSD, and MaxPSD, M3D9, EnD9, and WE). Optimum values for the design parameters of the SVM (regularization parameter: C) and MLP classifiers (number of neurons in the hidden layer: NH; regularization parameter: α) were obtained as those for which the Acc of the classifiers was the highest in the optimization set. Concerning SVM, the following values of C were assessed: 10−5, 10−4, 10−3, …, 104, 105. The optimum value of the input parameter C was 103, which maximizes Acc. Regarding MLP, NH was varied from 2 up to 50 and α was varied from 0 up to 10. Since the network depends on the initial random values of the weights, the accuracy was computed and averaged for a total of 10 runs for each pair NH-α. Finally, user-dependent network parameters NH = 5 and α = 1 were chosen since this pair reached the highest accuracy.

Diagnostic performance

The value of all the extracted features (ODI3, statistical moments, PSD, and DWT features) and the classification score of the LR, SVM, and MLP classifiers were obtained for each subject in the cross-validation set (S3 Table). Table 3 shows the diagnostic ability of each single feature in the cross-validation set obtained using optimum cut-off point obtained from the ROC curve. Most of the DWT-derived features (5 out of 7) showed accuracies near 80%. In this regard, MaxD9 achieved the highest performance (81.7±5.6% Acc, with 75.4±7.1% Se and 85.4±6.8% Sp), outperforming statistical moments and PSD features. Only ODI3 achieved slightly higher Acc than MaxD9, reaching 81.9±7.2% Acc (78.1±7.3% Se and 84.2±8.1% Sp). Table 4 shows the diagnostic performance of LR, SVM, and MLP classifiers, which were trained using the selected feature subset (ODI3, M2T, M2PSD, M3PSD, and MaxPSD, M3D9, EnD9, and WE) obtained with FCBF, in the cross validation set. These classifiers showed high diagnostic performance, outperforming all the extracted features in terms of Sp, PPV, LR+, and Acc. SVM achieved the highest accuracy (84.0±5.2% Acc, with 71.9±4.4% Se and 91.1±7.2% Sp) for the cutoff of 5 e/h.

[Figure omitted. See PDF.]

Table 3. Diagnostic ability of the proposed features (ODI3, statistical moments, PSD, and DWT) in the cross-validation set.

https://doi.org/10.1371/journal.pone.0208502.t003

[Figure omitted. See PDF.]

Table 4. Diagnostic ability of the LR, SVM, and MLP models in the cross-validation set.

https://doi.org/10.1371/journal.pone.0208502.t004

Discussion

In the present study, we examined the usefulness of wavelet analysis to identify features that characterize oximetry dynamics in order to expedite detection of moderate-to-severe pediatric SAHS. WE and features from the coefficients in D9 (M1D9-M4D9, MaxD9, and EnD9) were obtained from the DWT of each SpO2 recording. D9 (0.0244–0.0488 Hz) was chosen according to a previous study in the context of pediatric SAHS [18], and is related to the duration and frequency of the SpO2 desaturations associated with apneic events [40]. Statistically significant differences (p-value < 0.01) emerged in all DWT-derived features between the groups AHI <5 e/h and AHI ≥5 e/h in the optimization set (Table 2). The higher values showed by M1D9, MaxD9, and EnD9 in the AHI ≥5 e/h group agree with a higher amplitude of the histogram for high values of the D9 coefficients in this group. In addition, the SpO2 drops and rises caused by apneic events are reflected in a higher dispersion in the histogram of D9 coefficients, as reported by the higher values of M2D9 in the AHI ≥5 e/h group. In contrast, the lower values that M3D9 and M4D9 as reflected in the AHI ≥5 e/h group indicate that the variations produced in the SpO2 signal due to apneic events result in values less proximal to zero in the histogram of the D9 coefficients. Finally, the higher irregularity reported by WE in the SAHS positive group suggests that apnea-hypopnea events alter the energy distribution of the whole DWT profile of the SpO2 signal.

Regarding the diagnostic performance of the proposed features, ODI3 and MaxD9 reached similar Acc in the cross-validation set, higher than the remaining features. In addition, higher accuracies were generally obtained with the DWT-derived features with respect to statistical moments and features from PSD. This suggests that DWT is a useful approach to analyze the changes produced in the SpO2 signal associated to SAHS. In the feature selection stage, a feature subset composed of ODI3 (conventional oximetric index); M2T (time); M2PSD, M3PSD, and MaxPSD (PSD), and M3D9, EnD9 y WE (DWT) was obtained with FCBF. LR, SVM, and MLP models built with this subset obtained high diagnostic performance for the detection of moderate-to-severe SAHS (AHI ≥5 e/h), improving the diagnostic ability of the single features (Table 3) in terms of Sp, PPV, LR+, and Acc. It is worthy to note that the SVM model achieved the highest average Acc (84.0%), Sp (91.1%), PPV (83.8%), and LR+ (14.6) among the single features and binary classifiers. In addition, SVM reached similar NPV and LR- to LR, MLP, ODI3 and the remaining features. A high LR+ is especially important for screening tests [17,46]. In this sense, a LR+ greater than 10 is considered to provide strong evidence to confirm diagnoses [46]. Thus, our method is especially useful to confirm the presence of pediatric SAHS.

Three DWT features were involved in the feature subset obtained with FCBF: M3D9, EnD9 and WE. As aforementioned, these features provide information about the concentration of the D9 coefficients near zero (M3D9), the amplitude of the D9 coefficients (EnD9), and the irregularity of the distribution of the whole DWT profile of the SpO2 signal (WE). According to our results, M3D9, EnD9 and WE provide both relevant and complementary (non-redundant) information on the changes occurring in the SpO2 signal due to SAHS. This is consistent with the different properties of the SpO2 signal these DWT-derived features quantify. The fact that a high performance was reached with the three classification algorithms reinforces the notion that DWT is a useful method to analyze the SpO2 signal in the context of pediatric SAHS.

To the best of our knowledge, this is the first study assessing wavelet analysis of SpO2 recordings in the context of pediatric SAHS. Our results suggest that DWT is an appropriate tool to analyze the low frequency components of the SpO2 signals related to the duration of the desaturations caused by apnea-hypopnea events since it provides high resolution at low frequencies of the power spectrum [20,21]. This assumption is further supported by previous studies, whereby DWT was also applied to quantify the frequency components of different biomedical signals associated to respiratory events in the context of adult SAHS [23,24]. Additionally, the favorable performance of our approach may be due to the suitability of the WT to analyze non stationary properties of a signal [20,21], which is appropriate to events such as the non-stationary changes of the SpO2 signal associated with apneic events. The high resolution afforded by WT at low frequencies, as well as its suitability to analyze non-stationary signals clearly support the contention that DWT is more appropriate than conventional spectral analysis techniques to analyze the SpO2 signal [20,21].

Table 5 shows the performance of previous studies focused on the automated analysis of SpO2 as an alternative to PSG in the screening of moderate-to-severe pediatric SAHS [13–19]. Oxygen desaturation index and clusters of desaturations have been employed for this task [13–16]. Kirk et al. [13] applied ODI3, reaching 67.0% Se, 60.0% Sp, and 64.0% Acc. Tsai et al. [14] obtained 83.8% Se, 86.5% Sp, and 85.1% Acc using 4% ODI (ODI4). However, ODI4 cutoff values were optimized and validated using the same population, such that no true post-hoc verification was achieved. Chang et al. [15] combined ODI3 with common symptoms to assess a discriminative score, reaching 60% Se, 86% Sp, and 72% Acc. Pia-Villa et al. [16] reported 69.4% Acc (40.6% Se and 97.9% Sp) combining clusters of desaturations and clinical history in a discriminative score. Our approach achieved a high the diagnostic performance while also strengthening its validity since the methods were derived using not only a much larger sample size, but also applying a cross validation approach to validate the results.

[Figure omitted. See PDF.]

Table 5. Summary of the state-of-the-art studies in the context of detection of moderate-to-severe pediatric SAHS using SpO2 recordings.

https://doi.org/10.1371/journal.pone.0208502.t005

In order to increase the diagnostic ability of the SpO2 signal, conventional oximetric indices have been combined with features from other signal processing approaches in studies developed by our group [17–19]. Álvarez et al. [17] assessed LR models fed with conventional oximetric indices, statistical parameters, PSD, and nonlinear features. These models were validated using a bootstrap procedure, reaching 82.8% Acc (82.2% Se and 83.6% Sp). Vaquerizo-Villar et al. [18] assessed the usefulness of oximetry bispectrum. A multiclass multi-layer perceptron (MLP) model fed with ODI3, anthropometrical variables, PSD, and bispectral features reached 61.8% Se, 97.6% Sp, and 81.3% Acc in an independent test set, outperforming a MLP classifier built without bispectral features. Finally, Hornero et al. [19] analyzed 4,191 SpO2 recordings obtained from 13 sleep laboratories in a multicenter international study. A MLP regression model with ODI3 and the skewness of the PSD reached 68.2% Se, 87.2% Sp, and 81.7% Acc. In contrast with the findings of these studies, our current results achieved improved diagnostic ability for the screening of moderate-to-severe SAHS with the use of DWT-derived features. This suggests that wavelet analysis could enhance the detection of this clinically important and vulnerable group of SAHS severity from single-channel oximetry recordings. In these patients, it is essential to early detect this condition, since they are more likely to suffer from morbidities such as decreases in cognitive performance [3,4], as well as an increased C-reactive protein level due to systemic inflammation [5]. Moreover, an AHI ≥5 e/h is also associated with increased systemic blood pressure measurements and an increased risk for cardiac strain [3]. All these important negative consequences highlight the necessity of an early detection of moderate-to-severe pediatric SAHS [3].

Notwithstanding the highly promising results of our current approach, several limitations must be considered. First, the exclusive use of the SpO2 signal to detect SAHS may restrict the spectrum of physiological perturbations being detected by the oximetry signal, such as electroencephalographic arousals or reductions in airflow and increased intrathoracic pressure swings [1]. In this regard, the combination of SpO2 with other physiological signals from PSG could potentially enhance the performance of our proposed method but at the cost of adding significant complexity to the test. In addition, future research efforts may prospectively focus on identifying a specific mother wavelet for this task. However, our proposed approach achieved high performance with the Haar’s mother wavelet. Of note, the lack of universally accepted AHI severity cutoffs is another limitation that affects our study. Nevertheless, we have assessed the diagnostic ability of our proposal using an AHI cutoff of 5 e/h, a widely used criterion in the clinical decision making leading to the recommendation of surgical treatment [3,10]. Finally, it would be an interesting future goal to further validate our methodology in a larger sample of unattended oximetry recordings obtained at patients’ homes.

Conclusions

The application of WT has enabled the identification of features with the ability to characterize the effects of SAHS in the overnight oximetry profile of children. Features computed in the D9 detail level of the DWT as well as WE reached significant differences associated with the presence of SAHS. DWT has been found to provide complementary information to conventional approaches. Additionally, high diagnostic performance was reached using different reference binary classifiers, which emphasizes the usefulness of the DWT to provide discriminant information from oximetry signals. These results suggest that wavelet analysis could be useful to further characterize the oximetry signal and improve the diagnostic performance and implementation of abbreviated screening test for pediatric SAHS.

Supporting information

[Figure omitted. See PDF.]

S1 Table. Actual AHI from PSG and values of all the extracted features from the SpO2 signal (ODI3, statistical moments, PSD, and DWT features) of each subject in the optimization set.

https://doi.org/10.1371/journal.pone.0208502.s001

(XLSX)

S2 Table. The number of times each feature was selected with FCBF in the optimization set.

https://doi.org/10.1371/journal.pone.0208502.s002

(XLSX)

S3 Table. Actual AHI from PSG, values of all the extracted features from the SpO2 signal (ODI3, statistical moments, PSD, and DWT features), and the classification scores of LR, SVM, and MLP of each subject in the cross-validation set.

https://doi.org/10.1371/journal.pone.0208502.s003

(XLSX)

Citation: Vaquerizo-Villar F, Álvarez D, Kheirandish-Gozal L, Gutiérrez-Tobal GC, Barroso-García V, Crespo A, et al. (2018) Wavelet analysis of oximetry recordings to assist in the automated detection of moderate-to-severe pediatric sleep apnea-hypopnea syndrome. PLoS ONE 13(12): e0208502. https://doi.org/10.1371/journal.pone.0208502

References

1. Marcus CL, Brooks LJ, Ward SD, Draper KA, Gozal D, Halbower AC, et al. Diagnosis and Management of Childhood Obstructive Sleep Apnea Syndrome. Pediatrics. 2012;130: e714–e755. pmid:22926176

2. Alonso-Álvarez ML, Canet T, Cubell-Alarco M, Estivill E, Fernández E, Gozal D, et al. Documento de consenso del síndrome de apneas-hipopneas durante el sueño en niños. Arch Bronconeumol. 2011;47: 2–18.

3. Kaditis A, Kheirandish-Gozal L, Gozal D. Pediatric OSAS: Oximetry can provide answers when polysomnography is not available. Sleep Med Rev. Elsevier Ltd; 2016;27: 96–105. pmid:26146027

4. Hunter SJ, Gozal D, Smith DL, Philby MF, Kaylegian J, Kheirandish-Gozal L. Effect of sleep-disordered breathing severity on cognitive performance measures in a large community cohort of young school-aged children. Am J Respir Crit Care Med. 2016;194: 739–747. pmid:26930303

5. Church GD. The Role of Polysomnography in Diagnosing and Treating Obstructive Sleep Apnea in Pediatric Patients. Curr Probl Pediatr Adolesc Health Care. 2012;42: 22–25. pmid:22221590

6. Tan H-L, Gozal D, Ramirez HM, Bandla HPR, Kheirandish-Gozal L. Overnight polysomnography versus respiratory polygraphy in the diagnosis of pediatric obstructive sleep apnea. Sleep. 2014;37: 255–260. pmid:24497654

7. Nixon GM, Kermack AS, Davis GM, Manoukian JJ, Brown A, Brouillette RT. Planning adenotonsillectomy in children with obstructive sleep apnea: the role of overnight oximetry. Pediatrics. 2004;113: e19–e25. pmid:14702490

8. Katz ES, Mitchell RB, Ambrosio CMD. Obstructive Sleep Apnea in Infants. Am J Respir Crit Care Med. 2012;185: 805–816. pmid:22135346

9. Brockmann PE, Schaefer C, Poets A, Poets CF, Urschitz MS. Diagnosis of obstructive sleep apnea in children: A systematic review. Sleep Med Rev. Elsevier Ltd; 2013;17: 331–340. pmid:23375659

10. Kaditis AG, Alvarez MLA, Boudewyns A, Alexopoulos EI, Ersu R, Joosten K, et al. Obstructive sleep disordered breathing in 2- to 18-year-old children: Diagnosis and management. Eur Respir J. 2016;47: 69–94. pmid:26541535

11. Garde A, Dehkordi P, Karlen W, Wensley D, Ansermino JM, Dumont GA. Development of a screening tool for sleep disordered breathing in children using the phone oximeterTM. PLoS One. 2014;9. pmid:25401696

12. Berry RB, Budhiraja R, Gottlieb DJ, Gozal D, Iber C, Kapur VK, et al. Rules for scoring respiratory events in sleep: update of the 2007 AASM manual for the scoring of sleep and associated events: deliberations of the sleep apnea definitions task force of the American Academy of Sleep Medicine. J Clin sleep Med. American Academy of Sleep Medicine; 2012;8: 597. pmid:23066376

13. Kirk VG, Bohn SG, Flemons WW, Remmers JE. Comparison of Home Oximetry Monitoring With Laboratory Polysomnography in Children. CHEST J. 2003;124: 1702–1708.

14. Tsai C, Kang C, Su M, Lin H, Huang E, Chen C, et al. Usefulness of desaturation index for the assessment of obstructive sleep apnea syndrome in children. Int J Pediatr Otorhinolaryngol. Elsevier Ireland Ltd; 2013;77: 1286–1290. pmid:23732021

15. Chang L, Wu J, Cao L. Combination of symptoms and oxygen desaturation index in predicting childhood obstructive sleep apnea. Int J Pediatr Otorhinolaryngol. 2013;77: 365–371. pmid:23246417

16. Villa MP, Pietropaoli N, Supino MC, Vitelli O, Rabasco J, Evangelisti M, et al. Diagnosis of Pediatric Obstructive Sleep Apnea Syndrome in Settings With Limited Resources. JAMA Otolaryngol Neck Surg. 2015;141: 990–996. pmid:26540025

17. Álvarez D, Alonso-Álvarez ML, Gutiérrez-Tobal GC, Crespo A, Kheirandish-Gozal L, Gozal D, et al. Automated Screening of Children With Obstructive Sleep Apnea Using Nocturnal Oximetry: An Alternative to Respiratory Polygraphy in Unattended Settings. J Clin Sleep Med. 2017;13: 7–11.

18. Vaquerizo-Villar F, Álvarez D, Kheirandish-Gozal L, Gutiérrez-Tobal GC, Barroso-García V, Crespo A, et al. Utility of bispectrum in the screening of pediatric sleep apnea-hypopnea syndrome using oximetry recordings. Comput Methods Programs Biomed. Elsevier B.V.; 2018;156: 141–149. pmid:29428066

19. Hornero R, Kheirandish-Gozal L, Gutiérrez-Tobal GC, Philby MF, Alonso-Álvarez ML, Álvarez D, et al. Nocturnal Oximetry-based Evaluation of Habitually Snoring Children. Am J Respir Crit Care Med. 2017;196: 1591–1598. pmid:28759260

20. Rioul O, Vetterli M. Wavelets and signal processing. IEEE signal processing magazine. 19918: 14–38.

21. Daubechies I. The Wavelet Transform, Time-Frequency Localization and Signal Analysis. IEEE Trans Inf theory. 1990;36: 961–1005.

22. Fontenla-Romero O, Guijarro-Berdin B, Alonso-Betanzos A, Moret-Bonillo V. A new method for sleep apnea classification using wavelets and feedforward neural networks. Artif Intell Med. 2005;34: 65–76. pmid:15885567

23. Emin Tagluk M, Akin M, Sezgin N. Classification of sleep apnea by using wavelet transform and artificial neural networks. Expert Syst Appl. Elsevier Ltd; 2010;37: 1600–1607.

24. Khandoker AH, Palaniswami M, Karmakar CK. Support Vector Machines for Automated Recognition of Obstructive Sleep Apnoea Syndrome from Electrocardiogram Recordings. IEEE Trans Inf Technol Biomed. 2008;13: 1–32.

25. Jiang J-A, Lin R, Lee R-G, Tseng C-L, Zhou H-K, Chao C-F. A New Approach for Identifying Sleep Apnea Syndrome Using Wavelet Transform and Neural Networks. Biomed Eng Appl Basis Comm. 2006;18: 138–143.

26. Mendez MO, Corthout J, Van Huffel S, Matteucci M, Penzel T, Cerutti S, et al. Automatic screening of obstructive sleep apnea from the ECG based on empirical mode decomposition and wavelet analysis. Physiol Meas. 2010;31: 273–289. pmid:20086277

27. Lee YK, Bister M, Blanchfield P, Salleh YM. Automated detection of obstructive apnea and hypopnea events from oxygen saturation signal. Conf Proc IEEE Eng Med Biol Soc. 2004;1: 321–324. pmid:17271675

28. Vaquerizo-Villar F, Alvarez D, Gutierrez-Tobal GC, Barroso-Garcia V, Kheirandish-Gozal L, Crespo A, et al. Usefulness of discrete wavelet transform in the analysis of oximetry signals to assist in childhood sleep apnea-hypopnea syndrome diagnosis. Proceedings of the Annual International Conference of the IEEE Engineering in Medicine and Biology Society, EMBS. 2017. pp. 3753–3756. pmid:29060714

29. Sedano AC, Vaquerizo-Villar F, Álvarez D, Gutiérrez-Tobal GC, Barroso-Garc\’\ia V, Cerezo A, et al. Automated detection of childhood sleep apnea using discrete wavelet transform of nocturnal oximetry and anthropometric variables. Eur Respiratory Soc; 2017.

30. Magalang UJ, Dmochowski J, Veeramachaneni S, Draw A, Mador J, El-Solh A, et al. Prediction of the Apnea-Hypopnea Index From Overnight Pulse Oximetry. CHEST J. 2003;124: 1694–1701.

31. Garde A, Karlen W, Dehkordi P, Ansermino JM, Dumont GA. Oxygen Saturation Resolution Influences Regularity Measurements. 36th Annual International Conference of the IEEE In Engineering in Medicine and Biology Society (EMBC). 2014. pp. 2257–2260.

32. Yu L, Liu H. Efficient Feature Selection via Analysis of Relevance and Redundancy. J Mach Learn Res. 2004;5: 1205–1224.

33. Hosmer D, Lemeshow S. Applied logistic regression. John Wiley & Sons; 2004.

34. Bishop CM. Pattern Recognition and Machine Learning. 2006.

35. Bishop CM. Neural Networks for Pattern Recognition. 1995.

36. Witten I, Frank E, Hall M. Data Mining. 2011.

37. Steyerberg EW, Vergouwe Y. Towards better clinical prediction models: Seven steps for development and an ABCD for validation. Eur Heart J. 2014;35: 1925–1931. pmid:24898551

38. Cvetkovic D, Übelyi ED, Cosic I. Wavelet transform feature extraction from human PPG, ECG, and EEG signal responses to ELF PEMF exposures: A pilot study. Digit Signal Process. 2008;18: 861–874.

39. Rosso OA, Blanco S, Yordanova J, Kolev V, Figliola A, Schürmann M, et al. Wavelet entropy: a new tool for analysis of short duration brain electrical signals. J Neurosci Methods. 2001;105: 65–75. pmid:11166367

40. Taha BH, Dempsey JA, Weber SM, Badr MS, Skatrud JB, Young TB., et al. Automated detection and classification of sleep-disordered breathing from conventional polysomnography data. Sleep. 1997;20: 991–1001. pmid:9456464

41. Welch PD. The Use of Fast Fourier Transform fot he Estimation of Power Spectra: A Method Based on Time Averaging Over Short, Modified Periodograms. IEEE Trans audio Electroacoust. 1967;15: 70–73.

42. Efron B, Tibshirani RJ. Introduction to the Bootstrap. CRC press. 1994.

43. Guyon I. An Introduction to Variable and Feature Selection. J Mach Learn Res. 2003;3: 1157–1182.

44. Álvarez D, Hornero R, Marcos JVí, Wessel N, Penzel T, Glos M, et al. Assessment of Feature Selection and Classification Approaches To Enhance Information From Overnight Oximetry in the Context of Apnea Diagnosis. Int J Neural Syst. 2013;23: 1350020. pmid:23924411

45. Hornik K. Approximation capabilities of multilayer feedforward networks. Neural Networks. 1991;4: 251–257.

46. Deeks JJ. Diagnostic tests 4: likelihood ratios. Bmj. 2004;329: 168–169. pmid:15258077