Fuzzy Approximate Entropy of Extrema Based on Multiple Moving Averages as a Novel Approach in Obstructive Sleep Apnea Screening

Objective: Obstructive sleep apnea (OSA) is a respiratory disease associated with autonomic nervous system dysfunction. As a novel method for analyzing OSA depending on heart rate variability, fuzzy approximate entropy of extrema based on multiple moving averages (Emma-fApEn) can effectively assess the sympathetic tension limits, thereby realizing a good performance in the disease severity screening. Method: Sixty 6-h electrocardiogram recordings (20 healthy, 16 mild/moderate OSA and 34 severe OSA) from the PhysioNet database were used in this study. The performances of minima of Emma-fApEn (fApEn-minima), maxima of Emma-fApEn (fApEn-maxima) and classic time-frequency domain indices for each recording were assessed by significance analysis, correlation analysis, parameter optimization and OSA screening. Results: fApEn-minima and fApEn-maxima had significant differences between the severe OSA group and the other two groups, while the mean value (Mean) and the ratio of low-frequency power and high-frequency power (LH) could significantly differentiate OSA recordings from healthy recordings. The correlation coefficient between fApEn-minima and apnea-hypopnea index was the highest (|R| = 0.705). Machine learning methods were used to evaluate the performances of the above four indices. Random forest (RF) achieved the highest accuracy of 96.67% in OSA screening and 91.67% in severe OSA screening, with a good balance in both. Conclusion: Emma-fApEn may be used as a simple preliminary detection tool to assess the severity of OSA prior to polysomnography analysis.

of the upper respiratory tract during sleep, which will lead 23 to poor quality of sleep and thus daytime sleepiness [1]. 24 The occurrence mechanism of OSA includes upper airway  27 [5], etc. It's reported that OSA is closely associated with 28 diabetes [6], atherosclerosis [7], stroke [8], hypertension [9] 29 and a variety of cardiovascular diseases [10]. The severity 30 of OSA can also be used as a predictor of mortality [11]. 31 At present, there are at least about 1 billion OSA patients 32 in the world [12], and the prevalence of OSA is showing an 33 upward trend [13]. Moreover, many OSA patients are still 34 undiagnosed [14]. Therefore, there is an urgent need for a 35 simple and effective method of detecting OSA. 36 Polysomnography is considered as a gold standard for 37 OSA detection [15]. However, it has some disadvantages 38 including being expensive, causing patients' discomfort, 39 requiring expert judgment and having obscure results [16]. 40 In recent years, many researchers have proposed more effi-41 cient methods for OSA detection. OSA can affect a patient's 42 breathing state, nerve activity, muscle activity, heart func-43 tion and so on. Consequently, some studies developed OSA 44 detection algorithms based on oxygen saturation (SpO 2 ) [17], 45 VOLUME 10, 2022 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ [18], [19], [20], electroencephalogram (EEG) [21], [22], [23], 46 [24], electromyography (EMG) [25], [26], [27], electrocardiogram (ECG) [28], [29], [30], [31], etc. The results of these The performance of OSA classification is not only affected 91 by features, but also related to the quality of the classifier. 92 Janbakhshi et al. [46] proposed an ECG method to identify 93 OSA, and combined it with support vector machine (SVM), to classify apnea events. Alireza et al. [48] proposed a 102 random forest classifier to classify features after dimension-103 ality reduction by principal component analysis and linear 104 discriminant analysis, achieving an accuracy of 95.01%. 105 In above studies, most of them were analyzed by compar-106 ing the indices of the overall data of OSA groups and control 107 groups, and less attention was paid to the information of 108 some special points in the data. In HRV analysis, the general 109 trend of the RR interval sequences (RRs) reflects the general 110 changes in the ANS tension, and the extrema of the RRs are 111 related to the ANS tension limits.

112
Therefore, in order to analyze the changes of the ANS 113 tension limits in OSA patients through the extrema in the 114 general trend of the RRs, this study proposed the fuzzy 115 approximate entropy of extrema based on multiple moving 116 averages (Emma-fApEn). Due to the small partial fluctua-117 tions of the RRs, the multiple moving average method was 118 introduced into the experiment, which reduces the fluctuation 119 effect of the sequence by obtaining segments' averages of 120 the data several times. Sample entropy has the advantages 121 of anti-noise and anti-interference, but it depends on small 122 tolerances and forward matching of data length [49]. Per-123 mutation entropy can reflect the regularity of time series, 124 but it is not sensitive to the outliers in internal series [50]. 125 Fuzzy approximate entropy (fApEn) shows better monotonic-126 ity, consistency and robustness when describing signals of 127 different complexities [51]. So fApEn was applied to repre-128 sent the complexity of extrema. Moreover, the Emma-fApEn 129 method was combined with random forest (RF) classifier, 130 to achieve OSA screening. Therefore, the multiple moving 131 average method combined with the standardization method 132 was used to alleviate the small partial self-fluctuation of 133 RRs. Then the extremum sequences reflecting ANS tension 134 limits were extracted. By applying fApEn, Emma-fApEn was 135 obtained from extremum sequences. To evaluate the complex-136 ity of fluctuations of the ANS tension limits, we conducted 137 the significance analysis, the correlation analysis, the param-138 eter optimization, and the OSA screening for Emma-fApEn. 139

140
A. DATA

141
The database used in this paper was downloaded from 142 http://www.physionet.org/physiobank/database/apnea-ecg/, 143 which was applied in the Computers in Cardiology Challenge 144 2000 [52]. The database includes 70 single-channel nocturnal 145 ECG recordings from 32 OSA patients and normal subjects. 146 In this database, the sampling frequency of each recording is 147 100Hz, and the recording time is 401∼587min. Each minute 148 of each recording has a corresponding reference note, which 149 identifies whether apnea occurs in that minute.

150
The number of apnea and hypopnea events per hour during 151 sleep is defined as the apnea-hypopnea index (AHI), which 152 is an indicator of the severity of sleep apnea. Recordings 153 that contained at least 1 hour with AHI of 10 or more, and 154 at least 100 minutes labeled apnea were classified as the 155 apnea group. Recordings that contained at least 1 hour with 156 AHI of 5 or more, and between 5 and 99 minutes labeled 157 apnea were classified as the borderline group. Recordings that 158 contained at least 1 hour with AHI of less than 5, and fewer 159 than 5 minutes labeled apnea were classified as the control 160 group. In this study, the apnea and control groups were used 161 for analysis. According to AHI, the selected recordings were     the mean of the sum of the squares of differences between 189 adjacent the normal-to-normal intervals (RMSSD) of con-190 secutive 5-minute RR intervals. For N points RRs consisting 191 of consecutive 5-minute RR intervals, they are defined as 192 follows: The frequency domain analysis is to analyze the law of heart 198 rate change based on the power spectral density of RRs The extremum sequences can reflect the tension limits of the 211 ANS, but it cannot be extracted directly from the RRs due to 212 the small partial fluctuation of the RRs. Therefore, the RRs 213 were filtered by the multiple moving average method first, 214 and then the extremum sequences were extracted from them. 215 Next, fApEn of the minima (fApEn-minima) and fApEn 216 of the maxima (fApEn-maxima) were obtained by applying 217 fApEn to calculate the complexity of the minimum and max-218 imum sequences. The calculation method of the indices is 219 shown in Fig. 2. First of all, the gross errors of the RRs 220 were removed by the Pauta criterion [55]. And the Z-score 221 method was used to standardize the RRs. Secondly, in order 222 to extract the extremum sequences accurately, the multi-223 ple moving average method was used to reduce the local 224 slight fluctuation of the RRs. Finally, fApEn-minima and 225 fApEn-maxima were calculated to measure the complexity 226 of extremum sequences.

227
The raw RRs were divided into 1-minute non-repeated 228 sequence subsets. First, since gross errors of the RRs can 229 greatly affect the performance of the extremum sequences 230 reflecting ANS tension limits, it is necessary to further 231 VOLUME 10, 2022

Algorithm 1 (Error Reduction and Standardization)
Input: Initialize the following parameters.
remove them by Pauta criterion [55]. Then, the RRs were 2) the repeat counter K ; 3) the interval length M ; Output:the maximum sequence U and the minimum sequence V calculate the average length of Y : Algorithm 3 (Fuzzy Approximate Entropy, fApEn) Input: Initialize the following parameters. 1) the non-stationary time sequence X {x 1 , x 2 , . . . , x n }; 2) the parameters m, n, r; Output:the index fApEn for c = m;c ≤ m + 1;c + + do for i = 1;i ≤ n − c + 1;i + + do for j = 1;j ≤ c;j + + do The self-matching was not eliminated in this study, which 248 means d c ij = 0(i = j) in algorithm 3. In the calculation 249   The design of decision tree needs proper attribute selection 271 measure and pruning method. Most methods choose attribute 272 assignment quality measures to construct decision trees. Gini 273 index, the most commonly used attribute selection measure in 274 decision tree induction, is used to measure the total variance 275 across two OSA classes. It is defined as: In (7), p N (m) represents the proportion of observation 281 of the m th node in the tree belonging to the normal and 282 mild/moderate OSA groups; in (8), p M (m) represents the pro-283 portion of observation of the m th node in the tree belonging 284 to the mild/moderate OSA and severe OSA groups.   In the calculation of Emma-fApEn, the multiple moving 323 average method was used to reduce the interference of fine 324 fluctuations in the RRs to the general trend. The moving 325 average method reduced the interference while losing some 326 information contained in the RRs, so the repeat counter K and 327 the interval length M had an important impact on the results. 328 The relationships between the absolute value of the cor-329 relation coefficient |R| and K about fApEn-minima and 330 fApEn-maxima when M =3 are shown in Fig. 7(a). It was 331 The above analysis and the results in Table 1 show that 341 Mean and LH have significant differences between the nor-  Table 2. 358 As shown in Fig 8, the corresponding ROC curve of RF 359 classifier achieves the best performance, and the AUC of RF 360 classifier is obviously higher than others. The classification 361 performance of RF classifier reached 96.67% accuracy to 362 distinguish the normal group from the other two groups and 363 91.67% accuracy to distinguish the severe OSA group from 364 the other two groups. They provided the highest accuracy and 365 a good balance between sensitivity and specificity. Moreover, 366 RF achieved the best AUC in OSA and severe OSA detection. 367 The above results indicate that using RF classifier in the 368 detection of OSA and severe OSA is significantly better than 369 others. Therefore, RF classifier is chosen to achieve the OSA 370 screening.

373
The performance of the time-frequency domain indices was 374 also verified by previous studies. The decreases of SDNN 375 and RMSSD in OSA patients reported by Roche et al. [37] 376 were consistent with our results. Lack of monotonicity with 377 the degree of disease may be one of the reasons for their 378 significance. Gula et al. [38] presented that LH was the 379 most useful indicator for OSA detection, whose screening 380 performance has been verified in our work. Although LH 381 was a reliable index for screening OSA patients, there was 382 no significant difference in LH between the mild/moderate 383 VOLUME 10, 2022    Table 3, fApEn-427 maxima and fApEn-minima have a higher correlation (|R| ≥ 428 0.2) than other indices, whose violin plots are shown in Fig. 9. 429 fApEn-minima was statistically significant (R = −0.538), 430 demonstrating the relationship between the nonlinear indices 431 and AHI is more significant.  each decision tree is built by random attribute selections and 489 random samples from the original training set [63]. Outliers 490 caused by noise in the samples will reduce the accuracy of the 491 corresponding decision tree. The introduction of two types 492 of randomnesses allows outliers to affect some decision trees 493 instead of all the decision trees in a random forest. The forest 494 is a vote of all decision trees, which is less affected by the 495 decision trees with low accuracy. So RF does not overfit and 496 has strong anti-noise performance. Unbiased estimation of 497 generalization error can be realized while balancing errors, 498 which enables the method to obtain the strong generalization 499 ability [64]. The imbalance between the sympathetic and vagus nerve 503 system can lead to ANS dysfunction [36]. HRV is proved 504 to be an effective indicator to evaluate ANS function, and 505 its reduction is one of the features of OSA occurrence [34]. 506 As verified in this study, multiple classical time domain 507 indices can be used to represent HRV changes [37]. Low 508 frequency power reflects sympathetic tension, and high fre-509 quency power reflects parasympathetic tension. As their ratio, 510 LH can reflect the balance between parasympathetic tone and 511 sympathetic tone.

512
However, time-frequency domain analysis is linear and 513 more suitable for stationary signal, while the process of 514 regulating heart rate through ANS is nonlinear and non-515 stationary. In contrast, non-linear indices of HRV are bet-516 ter to analyze the complex dynamics of the heartbeat [41]. 517 Therefore, nonlinear methods are usually used to analyze 518 HRV [40]. In this study, multiple moving average method 519 was used to weaken the small partial fluctuations of the RRs 520 to obtain the general trend. The general trend of the RRs 521 reflected the changes of physiological sympathetic tension 522 during sleep. The extremum sequences in the general trend 523 can reflect the changes of the relative tension limits of the 524 sympathetic nerve. As shown in Fig. 3, the periodicity of the 525 general trend of the RRs in OSA patients was obvious, and 526 the fluctuation ranges of maximum and minimum sequences 527 were significantly slighter than that of normal subjects. These 528 results indicate that sympathetic tension in OSA patients 529 enters a nearly periodic change when sleep apnea occurs. The 530 spectrum of HRV is narrowed and the high frequency com-531 ponent decreases, which reduces the variation in the maxima 532 and minima. Moreover, it reduces the ability of OSA patients 533 to adapt to external changes [65]. These may be related to 534 decreased tension and activity of ANS in OSA patients.  Denoising method is important when obtaining extremum 560 sequences. The data must conform to independent identically 561 normal distribution [67], but the actual conditions do not gen-562 erally meet this requirement. Therefore, the algorithm used in 563 this study may not remove gross errors well. Second, there are 564 differences in age and BMI among different groups in the data 565 set we used. Individual differences are not taken into account 566 in the method, which may affect the experimental results.

567
Finally, potential cardiovascular diseases are not mentioned 568 in this study, which may have influenced the HRV analysis.

569
In the future studies, we will adopt more effective methods to