Wavelet Analysis of Overnight Airflow to Detect Obstructive Sleep Apnea in Children

This study focused on the automatic analysis of the airflow signal (AF) to aid in the diagnosis of pediatric obstructive sleep apnea (OSA). Thus, our aims were: (i) to characterize the overnight AF characteristics using discrete wavelet transform (DWT) approach, (ii) to evaluate its diagnostic utility, and (iii) to assess its complementarity with the 3% oxygen desaturation index (ODI3). In order to reach these goals, we analyzed 946 overnight pediatric AF recordings in three stages: (i) DWT-derived feature extraction, (ii) feature selection, and (iii) pattern recognition. AF recordings from OSA patients showed both lower detail coefficients and decreased activity associated with the normal breathing band. Wavelet analysis also revealed that OSA disturbed the frequency and energy distribution of the AF signal, increasing its irregularity. Moreover, the information obtained from the wavelet analysis was complementary to ODI3. In this regard, the combination of both wavelet information and ODI3 achieved high diagnostic accuracy using the common OSA-positive cutoffs: 77.97%, 81.91%, and 90.99% (AdaBoost.M2), and 81.96%, 82.14%, and 90.69% (Bayesian multi-layer perceptron) for 1, 5, and 10 apneic events/hour, respectively. Hence, these findings suggest that DWT properly characterizes OSA-related severity as embedded in nocturnal AF, and could simplify the diagnosis of pediatric OSA.


Introduction
Pediatric obstructive sleep apnea (OSA) is a major sleep-related breathing disorder, affecting a large number of children (5%) and increasing the risk of negative health consequences among those affected [1,2]. Several studies have reported that the characteristic nocturnal respiratory disruptions manifesting as either cessation or reductions in airflow lead to inadequate gas exchange and fragmented sleep that affects not only other physiological processes, but also cognitive development [2][3][4]. Consequently, the morbidities can worsen and become irreversible if OSA is not timely treated [5].
Despite the potential seriousness of OSA-related complications, it remains an underdiagnosed disorder due to relative unawareness of both parents and primary care physicians and the inherent difficulties in accessing the diagnostic test [6,7]. In this regard, the gold standard approach for diagnosing OSA is overnight polysomnography (PSG), which is technically complex, labor intensive, expensive, potentially distressing to the child and uncomfortable to the parent, and relatively unavailable [6,7]. These drawbacks AF and SpO 2 recordings were obtained from the PSG recordings by means of a thermistor and a pulse oximeter, respectively. According to AASM, the AF and SpO 2 recordings used in our study were resampled at the recommended frequency: 100 Hz for AF and 25 Hz for SpO 2 [15]. Both signals were subjected to a preprocessing stage in order to automatically remove possible artifacts. This stage was conducted following the artifact removal methods proposed in previous studies [12,17]. Signals whose duration was less than 3 h after artifact removal were excluded from our study [17,30]. Moreover, AF signals were normalized to minimize the inter-individual differences related to particular physiological characteristics other than OSA [31]. Figure 1 shows the scheme of the stages followed in this study: (i) feature extraction to characterize AF by means of WT, as well as to compute ODI3, (ii) feature selection using the fast correlation-based filter (FCBF) to obtain an optimal feature subset, and (iii) application of machine-learning approaches to carry out multiclass classification and regression. The multiclass classification was performed to determine the OSA severity degree through AdaBoost.M2 with decision trees as base classifier. Regarding to the regression process, it was carried out to estimate the AHI of each child by means of multi-layer perceptron neural network with a Bayesian approach (BY-MLP). In this sense, recent studies have shown the usefulness of these selection and machine-learning methods in the context of pediatric OSA diagnosis [12,17].

Wavelet Features
WT allows to conduct a multiresolution analysis of non-stationary signals, i.e., it allows to analyze a time series in different frequency ranges with variable resolution [24]. In our study, the analysis is performed using the discrete wavelet transform (DWT) to deal with the redundancy and computational complexity issues of the continuous approach [23,24]. Moreover, DWT has been successfully used in previous studies aimed at OSA detection [25,26,28,29]. As can be seen in Figure 2, multiresolution analysis with DWT consists of performing N = log2(M) decomposition steps, where M is the size of the time series x(n) [23,25]. At each step, the decomposition is conducted by using a scaling ϕj,n and a wavelet function ψj,n, which are generated by scaling and translation of basis functions and where j denotes the decomposition level [32]: where s = 2 and τ = 1 (dyadic sampling) are the scale and translation parameters, respectively, j ∈ Z is the decomposition level, k ∈ Z is the coefficient position within each subband of the decomposition, ϕ is the basis scaling function, and ψ is the basis wavelet function or mother wavelet. The scaling and wavelet functions allow to characterize the approximation and detail spaces at different resolutions, respectively. It cannot be affirmed that a generic optimal wavelet function exists, since in each particular case there will be a mother wavelet that better adapts to the signal under study. 3.1. Feature Extraction 3.1.1. Wavelet Features WT allows to conduct a multiresolution analysis of non-stationary signals, i.e., it allows to analyze a time series in different frequency ranges with variable resolution [24]. In our study, the analysis is performed using the discrete wavelet transform (DWT) to deal with the redundancy and computational complexity issues of the continuous approach [23,24]. Moreover, DWT has been successfully used in previous studies aimed at OSA detection [25,26,28,29]. As can be seen in Figure 2, multiresolution analysis with DWT consists of performing N = log 2 (M) decomposition steps, where M is the size of the time series x(n) [23,25]. At each step, the decomposition is conducted by using a scaling φ j,n and a wavelet function ψ j,n , which are generated by scaling and translation of basis functions and where j denotes the decomposition level [32]: where s = 2 and τ = 1 (dyadic sampling) are the scale and translation parameters, respectively, j ∈ Z is the decomposition level, k ∈ Z is the coefficient position within each sub-band of the decomposition, φ is the basis scaling function, and ψ is the basis wavelet function or mother wavelet. The scaling and wavelet functions allow to characterize the approximation and detail spaces at different resolutions, respectively. It cannot be affirmed that a generic optimal wavelet function exists, since in each particular case there will be a mother wavelet that better adapts to the signal under study. In our work, Haar and Daubechies-5 wavelets have been evaluated due to their previous suitability to the AF signal [33,34]. Figure 3 shows the mother wavelets Haar and Daubechies-5, as well as a 10-min segment of AF signal. Haar is the simpler orthonormal wavelet, which does not cause edge effect due to it uses a single vanishing moment and a support width = 1 (Figure 3) [35,36]. It would avoid ignoring the relevant information contained in these regions. Moreover, the stepped shape of the Haar wavelet would allow to detect abrupt reductions of AF caused by apneic events. Regarding Daubechies-5, it uses 5 vanishing moments and a support width = 9 ( Figure 3). When increasing the vanishing moment and the support width, the time-frequency localization improves [37,38], but the edge effect also increases [36]. Hence, the Daubechies-5 wavelet would provide a better localization while causing less edge effect than other higher order variants, such as Daubechies-10 or Daubechies-20 [36]. In addition, its shape resembles an AF signal, which would allow a better setting.  In practice, the scaling (φ) and wavelet (ψ) functions are considered as half-band low-pass h(n) and half-band high-pass g(n) filters, respectively, such as [24,32]: Sensors 2021, 21, 1491 6 of 19 frequency band can be split again, thus generating new levels of decomposition. At each decomposition level j, after the corresponding filter and subsampling, the approximation Aj (low frequency) and the detail Dj (high frequency) signals are obtained [24,25,29]: where A0 is the time series x(n). This decomposition process finishes when level j = N is reached.  Hence, as shown in Figure 2a, DWT can be implemented as a cascade of recursive filters, each of them followed by a subsampling process by 2 (dyadic sampling) to reduce the sampling frequency and increase the spectral resolution [24,29,32]. Thereby, the low frequency band can be split again, thus generating new levels of decomposition. At each decomposition level j, after the corresponding filter and subsampling, the approximation A j (low frequency) and the detail D j (high frequency) signals are obtained [24,25,29]: where A 0 is the time series x(n). This decomposition process finishes when level j = N is reached.
The AF signal was segmented into epochs of length 2 16 samples (≈10-min), as conducted in previous studies [12,39]. Thus, DWT decomposition for each segment was carried out at 16 levels (N = log 2 (2 16 )), using Haar and Daubechies-5 as mother wavelets [33,34]. An example of decomposition of AF segment by means of Daubechies-5 can be observed in Figure 2b.
Our DWT analysis mainly focused on D 8 (0.1953-0.3906 Hz), which corresponds to the normal breathing frequency band in sleeping children [16,40]. This choice let us characterize the alterations caused by OSA in normal nocturnal respiration. Thereby, after obtaining the detail coefficients of D 8 from each of the 946 AF signals, the following features were extracted to quantify the information contained in them [25]: Energy (E D8 ). This feature measures the quadratic amplitude of the detail signal D 8 , providing information about the activity produced in the resolution level associated to the representative frequency band of the normal breathing [25,41]. It is computed as the sum of the modulus of the detail coefficients squared [22,23]: In addition to the D 8 -derived features, the wavelet entropy from all detail levels has been obtained to also characterize the OSA global effects on the complete AF signal:

•
Wavelet entropy (WE). It is an extension of the well-known Shannon's entropy. Therefore, this feature allows quantifying the energy distribution changes generated in the decomposition process, offering information about the underlying dynamical behavior and the irregularity of the signal [22,25,41]: where p j is the normalized energy distribution at the decomposition level j: Due to the amplitude reductions of the AF signal produced by apneic events [15], as well as the ability of the DWT to assign low coefficients to the flatter signal parts and high coefficients to the steeper [42], it is expected to find higher coefficients in the signal D 8 of the subjects without OSA. Consequently, and according to previous spectral analyses of AF [16,43], higher values of M 1D8 , M 2D8 , Max D8 , and Min D8 , as well as lower values of M 3D8 and M 4D8 , are expected in these subjects. In addition, it is also expected that there will be greater activity in the normal breathing band in the absence of apneas and hypopneas, i.e., that subjects without OSA have a higher E D8 . Regarding WE, it has been observed that apneic events introduce changes in time and frequency domains of AF signal of children [11,12,16]. Thus, AF signals are expected to be more irregular (higher values of WE) as OSA severity increases. Thereby, since DWT allows to obtain higher frequency resolution at low frequencies, it is expected to provide a more detailed analysis than classical spectral analysis.

Oximetry Index
ODI3 was obtained from the SpO 2 signal in order to compare the AF information with a common clinical index that usually acts as a surrogate of the full PSG, when the latter is not available. The guideline followed to calculate ODI3 was established according to the definition provided by Taha et al. [44]. Thus, the SpO 2 reductions greater than or equal to 3% and lasting at least 10 s were considered desaturation events. Finally, ODI3 was computed dividing the total number of desaturations presented in the SpO 2 signal by the recording time expressed in hours. Due to the blood oxygen desaturations are closely related to the occurrence of apneic events [15], ODI3 is expected to be higher as OSA severity increases.

Feature Selection
The fast correlation-based filter (FCBF) was used to obtain an optimal subset of relevant and non-redundant features [45]. Some motivations for using this algorithm are that it does not depend on subsequent analyzes and it reduces the complexity and dimensionality of the predictive models [45,46]. Moreover, this method has been successfully applied in the pediatric OSA context [12,13,17,25,47].
A bootstrapping procedure was conducted (1000 replicates) to obtain a stable and optimal subset [46]. The average significance was used as selection threshold T s [30]. It was computed as the average number of times that all features were selected [30]. Hence, the features selected a number of times ≥T s constituted the optimal set of features that fed the predictive models proposed in our study.

Multiclass Classification
In order to carry out a classification into 4-classes (no-OSA, mild, moderate, and severe OSA), we used the adaptive boosting (AdaBoost.M2) ensemble-learning algorithm due to previous success in classifying AF signals in the OSA context [12,48]. AdaBoost.M2 is based on iterative training of multiple 'weak classifiers' (also, base classifiers) of the same type, so that each new one focuses on the data misclassified by its predecessors [49,50]. Thereby, this method calls the base classification algorithm L times, giving it each time a different weight distribution for the training. The same weight is initially assigned to all instances and then it is updated in each iteration: more weight to wrongly classified instances and less weight to those correctly classified [49,50]. This allows the algorithm to adapt to the data, minimizing the expected error and focusing on correctly classifying the instances with more weight [50]. The weight update also involves the use of the learning rate α, a regularization parameter to deal with overfitting [48,50]. Finally, the weighted vote of all the previously trained classifiers is computed, obtaining a more robust prediction of each class [49,50].
In our study, decision trees were used as base classifiers, which is a common choice when using ensemble-learning methods [49]. The optimal number of weak classifiers L, as well as the learning rate α, require to be tuned. In this regard, both hyperparameters were optimized by applying bootstrapping (1000 replicates) in the training group and using 0.632 bootstrap to estimate the Cohen's kappa (kappa) for each L/α pair [12,48].

Regression
A regression process by means of multi-layer perceptron neural network with a Bayesian approach (BY-MLP) was conducted to estimate AHI. BY-MLP has already shown its utility in OSA diagnosis [17,51,52]. This method is based on a set of perceptrons organized in layers, so that each perceptron is connected to all those of the next layer with a certain weight [53]. In order to obtain a universal approximation, BY-MLP is typically formed by 3 layers: (i) an input layer with N I perceptrons that receive the input patterns and propagate them to all the perceptrons of the next layer, (ii) a hidden layer with N H perceptrons that perform a non-linear processing of the received patterns and propagate it to the next layer, and (iii) an output layer with N O perceptrons that process the information received from the hidden layer and provide the response of the neural network [53]. The learning process of the network involves adjusting the weights associated to the connections between perceptrons. The technique applied to carry out this adjustment is the Bayesian inference, which allows finding the optimal weights that minimize the error function [52].
In our study, N I is equal to the number of features selected in the previous stage, N H is a parameter to be tuned, and N O is equal to one perceptron since the network purpose is to estimate the AHI. As in Adaboost.M2, the value of N H was optimized by applying bootstrapping (1000 replicates) in the training group and using 0.632 bootstrap to estimate the kappa [51].

Statistical Analysis
The data distribution of each biomedical index was evaluated by means of the Lilliefors test. The results showed that wavelet features did not follow a normal distribution. Consequently, the existence of statistically significant differences (p-value < 0.01) among OSA severity groups was assessed by means of the non-parametric Kruskal-Wallis test. Moreover, Spearman's correlation was used to evaluate the relationship between AHI and the features under study. Cohen's kappa of two-class (kappa 2 ) and four-class (kappa 4 ), as well as the four-class accuracy (Acc 4 ), assessed the agreement between predicted and actual diagnosis [54]. In addition, the metrics used to evaluate the diagnostic performance of the machine-learning approaches for the common AHI thresholds 1 e/h, 5 e/h, and 10 e/h were sensitivity (Se), specificity (Sp), accuracy (Acc), positive and negative predictive values (PPV and NPV, respectively), and positive and negative likelihood ratios (LR+ and LR-, respectively). The statistically significant differences (p-value < 0.001) between diagnostic metrics of the models were evaluated using the Mann-Whitney U test for pairwise comparison with Bonferroni correction.

Extracted Features
The coefficients of D 8 with and without sign computed by means of Haar obtained features with an average Spearman's correlation of 0.2075 and 0.3190, respectively. Regarding Daubechies-5, the features achieved an average correlation of 0.2379 and 0.3447 for coefficients with and without sign, respectively. Consequently, the features extracted from D 8 in absolute value and with Daubechies-5 were used in the following stages of our study. Figure 4 displays the averaged detail signal D 8 for each severity group in the training dataset. In this figure, higher amplitude values of D 8 can be observed in subjects without OSA. In addition, the frequency distribution of the values of the D 8 coefficients as OSA severity increases can be visualized in Figure 5. According to this figure, coefficients close to 0 are more frequent in children with OSA, increasing the asymmetry and the sharpness of the peak of the distribution as the AHI increases.          Table 2 shows the median and interquartile range values by OSA severity group of each feature, as well as its Spearman's correlation coefficient with the AHI and the p-value obtained by means of the Kruskal-Wallis test in the training group. M 1D8 , M 2D8 , Max D8 , Min D8 , and E D8 showed a decreasing trend as OSA severity increases. In contrast, the tendency of M 3D8 , M 4D8 , and ODI3 was towards higher values. Regarding WE, it showed less separability between groups, although with a notable increase in the most severe subjects. In addition, all the extracted features showed significant differences among OSA severity groups (p-value < 0.01 after Bonferroni correction).

Feature Selection
We carried out 2 selection trials, one only with wavelet features from AF and another that also included the ODI3. FCBF was applied to 1000 bootstrap replicates obtained from the training group. In the trial with wavelet features from AF, only M 3D8 was selected more than T s times (T s = 125.25). As can be seen in Figure 6, M 3D8 , Min D8 , and ODI3 exceeded this threshold (T s = 205.33) when ODI3 was included in the selection process. We carried out 2 selection trials, one only with wavelet features from AF and another that also included the ODI3. FCBF was applied to 1000 bootstrap replicates obtained from the training group. In the trial with wavelet features from AF, only M3D8 was selected more than Ts times (Ts = 125.25). As can be seen in Figure 6, M3D8, MinD8, and ODI3 exceeded this threshold (Ts = 205.33) when ODI3 was included in the selection process.

Test Group
In order to improve the generalization of our results, the trained models were assessed in 1000 bootstrap replicates from the test group. The diagnostic performance metrics were obtained by means of the bootstrap 0.632 procedure and the statistically significant differences (p-value < 0.001) between models were evaluated using the Mann-Whitney U test for pairwise comparison with the Bonferroni correction. Tables 3 and 4 show the diagnostic performance (median [95% confidence interval]) achieved by each of the machine-learning models proposed in our study, as well as the ODI3. ODI3 showed a severity underestimation in 1 and 5 e/h. Regarding AB AF and BY-MLP AF models, these obtained an unbalanced Se-Sp pair, with a severity overestimation in 1 and 5 e/h and an underestimation in 10 e/h. When combining wavelet features from AF with ODI3 (AB AF,ODI3 and BY-MLP AF,ODI3 ) these negative effects of overestimation and underestimation were reduced. BY-MLP AF,ODI3 obtained highest Acc for 1 and 5 e/h and AB AF,ODI3 for 10 e/h, significantly outperforming (p-value < 0.001) the individual approaches. Moreover, these models also achieved higher performance than AB AF , BY-MLP AF , and ODI3 in terms of kappa2, kappa4, and Acc4.

Test Group
In order to improve the generalization of our results, the trained models were assessed in 1000 bootstrap replicates from the test group. The diagnostic performance metrics were obtained by means of the bootstrap 0.632 procedure and the statistically significant differences (p-value < 0.001) between models were evaluated using the Mann-Whitney U test for pairwise comparison with the Bonferroni correction. Tables 3 and 4 show the diagnostic performance (median [95% confidence interval]) achieved by each of the machine-learning models proposed in our study, as well as the ODI3. ODI3 showed a severity underestimation in 1 and 5 e/h. Regarding AB AF and BY-MLP AF models, these obtained an unbalanced Se-Sp pair, with a severity overestimation in 1 and 5 e/h and an underestimation in 10 e/h. When combining wavelet features from AF with ODI3 (AB AF,ODI3 and BY-MLP AF,ODI3 ) these negative effects of overestimation and underestimation were reduced. BY-MLP AF,ODI3 obtained highest Acc for 1 and 5 e/h and AB AF,ODI3 for 10 e/h, significantly outperforming (p-value < 0.001) the individual approaches. Moreover, these models also achieved higher performance than AB AF , BY-MLP AF , and ODI3 in terms of kappa 2 , kappa 4 , and Acc 4 . Se = sensitivity; Sp = specificity; Acc = accuracy; PPV = positive predictive value; NPV = negative predictive value; LR+ = positive likelihood ratio; LR-= negative likelihood ratio; kappa 2 = Cohen's kappa of two-class; 95%CI = 95% confidence interval; ND = Non defined; ODI3 = 3% oxygen desaturation index; AB AF = Adaboost.M2 fed with optimal wavelet features from AF; AB AF,ODI3 = Adaboost.M2 fed with optimal wavelet features from AF and ODI3; BY-MLP AF = Bayesian multi-layer perceptron fed with optimal wavelet features from AF; BY-MLP AF,ODI3 = Bayesian multi-layer perceptron fed with optimal wavelet features from AF and ODI3, a Significant differences (p-value < 0.001) between AB AF and AB AF,ODI3 ; b Significant differences (p-value < 0.001) between AB AF and BY-MLP AF ; c Significant differences (p-value < 0.001) between AB AF and BY-MLP AF,ODI3 ; d Significant differences (p-value < 0.001) between AB AF and ODI3; e Significant differences (p-value < 0.001) between AB AF,ODI3 and BY-MLP AF ; f Significant differences (p-value < 0.001) between AB AF,ODI3 and BY-MLP AF,ODI3 ; g Significant differences (p-value < 0.001) between AB AF,ODI3 and ODI3; h Significant differences (p-value < 0.001) between BY-MLP AF and BY-MLP AF,ODI3 ; i Significant differences (p-value < 0.001) between BY-MLP AF and ODI3; j Significant differences (p-value < 0.001) between BY-MLP AF,ODI3 and ODI3. Adaboost.M2 fed with optimal wavelet features from AF and ODI3; BY-MLP AF = Bayesian multi-layer perceptron fed with optimal wavelet features from AF; BY-MLP AF,ODI3 = Bayesian multi-layer perceptron fed with optimal wavelet features from AF and ODI3; a Significant differences (p-value < 0.001) between AB AF and AB AF,ODI3 ; b Significant differences (p-value < 0.001) between AB AF and BY-MLP AF ; c Significant differences (p-value < 0.001) between AB AF and BY-MLP AF,ODI3 ; d Significant differences (p-value < 0.001) between AB AF and ODI3; e Significant differences (p-value < 0.001) between AB AF,ODI3 and BY-MLP AF ; f Significant differences (p-value < 0.001) between AB AF,ODI3 and BY-MLP AF,ODI3 ; g Significant differences (p-value < 0.001) between AB AF,ODI3 and ODI3; h Significant differences (p-value < 0.001) between BY-MLP AF and BY-MLP AF,ODI3 ; i Significant differences (p-value < 0.001) between BY-MLP AF and ODI3; j Significant differences (p-value < 0.001) between BY-MLP AF,ODI3 and ODI3.

Discussion
In the present work, we characterized overnight pediatric AF by means of wavelet features obtained from the normal respiratory band. In addition, we showed the complementarity between the wavelet analysis conducted on AF and ODI3, obtaining two machine-learning models with high performance to diagnose pediatric OSA. The interpretation of our findings is detailed below.

Training Group
Our results revealed that the mother wavelet Daubechies-5 is better suited to the nocturnal AF behavior of children than the Haar wavelet. This may be because Daubechies-5 uses more vanishing moments and a larger support width than Haar (Figure 3), which allowed to obtain a more precise phase space and time-frequency location, as well as to focus and detect the singularities of the AF signal [37,38].
According to Figures 4 and 5, the subjects without OSA presented higher values of D 8 , as well as a less asymmetric and less peaked distribution than the subjects with OSA. These amplitude and distribution differences agree with the information provided by the statistical analysis carried out in the training group. Thereby, M 1D8 and M 2D8 showed a decreasing tendency as the OSA severity increased. This fact is consistent with a less activity in the normal breathing band as AHI increased, which causes a notable reduction of the coefficients of signal D 8 from AF and a narrower dispersion range. Regarding M 3D8 and M 4D8 , these experienced an increasing tendency, i.e., greater positive skewness and greater sharpness of the distribution peak in lower values of the coefficients of D 8 as the AHI is higher. This indicates that apneas and hypopneas change the frequency distribution of AF signal and reduce its frequency components in the normal breathing band, which leads to fewer high coefficients and more coefficients close to zero in this band. According to this reduction of coefficients, the maximum and minimum values of D 8 (Max D8 and Min D8 ), as well as the energy of this level (E D8 ), were lower as the severity of OSA increased. This decrease revealed that apneic events reduce the detail signal amplitude and the activity produced in the resolution level associated to the normal breathing band, which agrees with a lower occurrence of normal breathing patterns. Moreover, an increase of WE could also be observed in the severely affected children. This fact suggests that severe OSA disturbs the energy distribution of AF signal at the different levels of decomposition. In this way, the wavelet energy is redistributed in other frequency bands associated with apneic events instead of concentrating in the normal breathing band. Consequently, the AF signal becomes more irregular in these cases. Thus, we showed that DWT provided useful information to characterize the nocturnal AF of children, as well as its behavior in presence of apneic events.

Feature Selection and Diagnostic Performance
Despite the clear tendencies shown by the 8 wavelet features extracted from AF, only the asymmetry of the distribution of the coefficients of D 8 (M 3D8 ) was relevant and not redundant with respect to the rest. This is coherent with the statistical analysis carried out in the training group, where M 3D8 not only showed significant differences among severity groups, but was also the wavelet feature that obtained highest Spearman's correlation with the AHI. Thus, the selection of M 3D8 suggests that its individual utility to characterize the pediatric OSA is greater than individual and joint usefulness of the other wavelet features.
However, FCBF revealed that there is complementarity between ODI3 and the wavelet information obtained from both M 3D8 and Min D8 . In this regard, M 3D8 , Min D8 , and ODI3 showed a clear separability among OSA severity groups (lower p-value), as well as a higher Spearman's correlation with the AHI. Notice that FCBF did not select Min D8 in the trial with wavelet features from AF, but it was selected when ODI3 was included in the selection process. Due to the final subset is obtained as the union of the feature subsets selected in ≥T s bootstrap replicates [46], this fact suggests that Min D8 contributes with additional information to ODI 3 and different from that provided by M 3D8 , which highlights the joint usefulness of these features. Therefore, the information about the occurrence of apneic events provided by DWT through the distribution asymmetry (M 3D8 ) and the minimum amplitude of D 8 (Min D8 ) from AF is complementary to the information provided by ODI3 about the occurrence of desaturations.
This complementarity was also reflected in the machine-learning models used to classify children in 4-classes of OSA severity and estimate their AHI. AB AF , BY-MLP AF , and ODI3 achieved a moderate diagnostic performance in the testing group. These approaches were significantly outperformed (p-value < 0.001) by BY-MLP AF,ODI3 in kappa 2 for the 3 AHI cut-offs, kappa 4 , and Acc 4 , as well as by AB AF,ODI3 in kappa 2 for 1 and 10 e/h and kappa 4 . This fact revealed that the agreement between the predicted and actual OSA severity improves when wavelet information and ODI3 are combined, which confirmed the complementarity of both approaches. Moreover, the diagnostic accuracies reached by BY-MLP AF,ODI3 for 1, 5, and 10 e/h and AB AF,ODI3 for 1 and 10 e/h were also significantly higher (p-value < 0.001) than those obtained by the individual approaches. In this regard, BY-MLP AF,ODI3 achieved a more balanced Se-Sp pair in 5 e/h, as well as significantly higher Se and Acc (p-value < 0.001) than AB AF,ODI3 in 1 and 5 e/h. Hence, the AHI estimated by BY-MLP AF,ODI3 could be used by medical personnel to discriminate between mildly to moderately affected subjects. In addition, AB AF,ODI3 obtained a statistically significant higher Acc (p-value < 0.001) for 10 e/h. Furthermore, as LR+ > 10 is a robust indicator to determine the presence of a disease [55], the significantly higher LR+ (p-value < 0.001) obtained with AB AF,ODI3 for 10 e/h (LR+ = 18.99 [14.60, 51.76]) indicated that this model provides greater evidence than BY-MLP AF,ODI3 to detect severe OSA. Consequently, this multiclass classifier could be used as a potential automatic tool to this purpose. Therefore, the methodology proposed in the present study would be a useful alternative to nocturnal polysomnography, since it would help to simplify and streamline the pediatric OSA diagnosis in a timely fashion that potentially prevents the generation and worsening of deleterious consequences.

Comparison with Other Studies
As shown in Table 5, several studies focused on the simplification of pediatric OSA have applied automatic analysis techniques to polysomnographic signals such as ECG, PPG, SpO 2 , and AF [8][9][10][11][12][13][14]17,25,47]. In this regard, Shouldice et al. [10] applied temporal and spectral analysis methods to 50 ECG signals, obtaining 84.00% Acc for 1 e/h. Other studies such as those carried out by Gil et al. [8] and Dehkordi et al. [9] analyzed 21 and 146 PPG signals, respectively. They evaluated their proposal using 5 e/h as AHI cut-off for positive OSA and reached 80.00% Acc and 71.00% Acc, respectively. However, despite the high diagnostic performance achieved by these studies, the small sample size makes their results difficult to generalize.
The work of Vaquerizo-Villar et al. [25], a previous study from our research group, was based on wavelet analysis combined with other features extracted from 981 SpO 2 signals. This study used only the AHI cut-off 5 e/h to determine the presence of OSA, achieving a high Acc. However, our current proposal was evaluated according to the different severity categories, which allowed us to also identify children without OSA as well as the severely affected ones. In this regard, our methodology showed great usefulness for detecting severe cases (Acc = 90.99% [89. 29,92.61] and LR+ = 18.99 [14.60,51.76] with AB AF,ODI3 for 10 e/h). As these pediatric subjects have an increased risk of developing comorbidities and neurocognitive deficits [1,5], our tool would allow to diagnose and treat these cases early on, before the consequences are irreversible. Moreover, automatic detection of no-OSA cases would reduce the waiting lists and streamline the diagnosis of children affected by this disease. In addition, it is worth noting the achieved improvement with respect to our previous study based on recurrence plots from AF [17]. Although the Acc in 1 e/h was lower, a more balanced Se-Sp pair should be noted, as well as the higher diagnostic performance obtained in 5 e/h. Moreover, we outperformed the Acc and LR+ obtained in 10 e/h, thus our current proposal is more robust to detect severe OSA. Regarding our latest work focused on the bispectral analysis of AF [56], it was improved in Sp, PPV, and LR+ for 1 and 5 e/h, as well as in all diagnostic metrics for 10 e/h by the current proposal. Hence, the DWT would be a more potentially useful tool than bispectrum to diagnose severe OSA cases in a timely fashion. Furthermore, DWT revealed changes of energy and frequency components of AF, that could not have been detected with bispectrum or recurrence plots.

Limitations and Future Work
One of the limitations of this study is the size of the population under study, since it would have been desirable to analyze a larger set of AF signals originating from multiple sources to ensure more generalizable results. Although we have shown that Daubechies-5 better reflects the behavior of AF than Haar, other mother wavelets that have not yet been applied in this context could be analyzed in future studies, and the results obtained with each of them could be compared. Moreover, it would also be interesting to search for mother wavelets that more closely resemble the shape of AF signal. The use of DWT and the analysis of the detail level D 8 have also shown its usefulness in the diagnosis of the disease. However, wavelet packet decomposition could be applied in future research to obtain detail levels more fitting to the frequency features of the AF signal [57]. In addition, another future goal is to perform multiclass classification and/or AHI estimation through other advanced techniques of pattern recognition, such as deep-learning [58], and pediatric AF at-home recordings. Finally, the analysis proposed in this study could be used along with other PSG-derived signals, such as thoracic or abdominal effort, to distinguish among obstructive, central, and mixed apneic events in future research [15].

Conclusions
In this study, nocturnal AF of children has been characterized by means of features from wavelet analysis. We found that apneic events decrease the detail coefficients of the DWT associated to the normal breathing band (lower values of M 1D8 , Max D8 , and Min D8 ), as well as the activity produced in the frequency range 0.1953-0.3906 Hz of AF (lower values of E D8 ). Wavelet analysis also revealed that OSA changes the frequency and energy distribution of AF signal, reducing its frequency components in the normal breathing band (higher values of M 3D8 and M 4D8 ) and generating more global irregularity in the signal (higher values of WE). Our study also found complementarity between DWT features from AF and ODI3. These findings allowed us to obtain a BY-MLP model with high diagnostic accuracy to estimate the AHI of mildly to moderately affected children, as well as an AdaBoost.M2 model particularly useful for classifying severe pediatric subjects. Therefore, we conclude that DWT can characterize the nocturnal AF, and that such approach could be jointly used with ODI3 to simplify the diagnosis of OSA in children.