1. Introduction

Diabetes mellitus and diabetic complications are predicted to increase world-wide. Cardiovascular complications of diabetes account for 65% of diabetic deaths and cardiac autonomic neuropathy (CAN) belongs to the major contributors to sudden cardiac deaths (SCD). The outcomes of the Action to Control Cardiovascular Risk in Diabetes (ACCORD) and National Diabetes Strategy (NDS) trial recommended regular screening of people with diabetes, especially for comorbidities including autonomic nervous system (ANS) dysfunction [1,2,3]. The increased risk of arrhythmias and SCD associated with CAN makes screening of people with diabetes an important health issue as early detection and intervention improve treatment effectiveness [4,5].

Mortality associated with CAN has been estimated to be approximately five times higher in diabetes with autonomic neuropathy (29%) compared to diabetes without autonomic neuropathy [5,6]. Prevalence of CAN lies between 22% and 60% in people with diabetes depending on duration of diabetes and treatment efficacy. CAN, even at the asymptomatic stage can lead to significant changes in the heart rate variability (HRV) metrics reflecting impaired heart rate and blood pressure regulation. Early detection of CAN and timely treatment are important for risk stratification and reducing the effects of diabetic neuropathy, which affects the function of all major body organs, leading to kidney failure, urinary dysfunction and cardiac arrhythmias.

In the present paper we propose a new machine learning algorithm for the diagnosis of CAN progression based on HRV attributes, the Multi-Layer Attribute Selection and Classification (MLASC) algorithm. It begins with attribute extraction and incorporates our new automated attribute selection procedure, Double Wrapper Subset Evaluator with Particle Swarm Optimization (DWSE-PSO) combined with a multi-layer classifier. We present the results of a comprehensive collection of experiments, which compare MLASC with other simpler versions and counterpart methods.

It is a new task to detect CAN automatically on the basis of HRV attributes, which are often available in relation to other health investigations. We conducted comprehensive experiments comparing several automated feature selection and optimization methods that have not been considered before for the task of classification of CAN progression using HRV attributes. In addition, we introduce a new feature selection method, and apply it for solving this task. Our experiments compare the effectiveness of these automated feature selection methods in conjunction with the application of hierarchical classifiers for the classification of CAN progression.

Our previous article [7] studied multi-level classifiers for neurological diagnostics using different attributes including traditional tests for identification of CAN and available in the Diabetes Screening Complications Research Initiative (DiScRi) database. The Ewing battery is commonly used for detecting CAN but is often not conclusive and therefore more sensitive and accurate tests are required. ECGs are routinely assessed in clinical practice and although they do not directly indicate CAN, HRV can be determined from the interbeat interval tachogram or from a continuous heart rate recording [8]. HRV as a clinical tool using ECG recordings has been shown to be a sensitive marker for risk of future arrhythmias or CAN and easier to use clinically compared to the Ewing battery.

Here we concentrate on developing optimal combinations of options to be employed at the different layers of MLASC for the diagnostics of CAN progression based on HRV attributes in people with diabetes. Background information on HRV attributes, CAN and the DiScRi database are included in Section 3. HRV information can include as many as 20-30 measures sensitive to different characteristics of the ECG time series that can be divided into time, frequency, and nonlinear measures. The use of HRV in diagnostics has been verified in identification of neonatal sepsis [9], arrhythmia risk [10], and in the risk of future adverse cardiac events. Which of the 20-30 measures are most sensitive for the diagnostics of CAN in diabetes progression has not been addressed in [8,11] and will be the subject of subsequent studies.

The present work was motivated by other recent results indicating the utility of incorporating several HRV measures in the classification of CAN [12,13]. Several other studies handled other different multi-stage systems [14,15,16,17]. However, this is the first article where multi-layer attribute selection is included in one of the layers of the main algorithm. The ability to use only HRV for accurate identification of CAN and CAN progression reduces the need for invasive testing such as obtaining cholesterol, BGL and HbA1c results.

2. Heart Rate Variability Attributes as Analysis Tool

HRV attributes are very important and are often collected for patients with various pathologies such as sepsis, heart failure and kidney failure. In order to investigate the role of HRV attributes and the capability of multi-level classifiers in improving diagnostic accuracy, we used a large database of health-related parameters and tests amalgamated in DiScRi [18] organized by Charles Sturt University.

Let us refer to [19,20,21] and [7] for preliminaries on the classification of disease progression associated with CAN. In particular, this article uses the classical set of five Ewing attributes specified in the following table (Table 1).

Table 1. Tests in the Ewing battery.

Test	Normal	Borderline	Abnormal
(i) Heart rate response to standing (ratio of HR at 15 and 30 sec)	≥ 1.04	1.01-1.03	≤ 1.00
(ii) Blood pressure response to standing (mmHg)	≤ 10	11-29	≥ 30
(iii) Heart rate response to deep breathing (BPM)	≥ 15	11-14	≤ 10
(iv) Valsalva maneuver and heart rate response (BPM)	≥ 1.21	1.11-1.20	≤ 1.10
(v) Blood pressure response to sustained handgrip (mmHg)	≥ 16	11-15	≤ 10

 | Show Table

DownLoad: CSV

Following [7], we investigated three classifications of CAN progression introduced in [19,20]. They are the classifications into 2,3 and 4 classes: the presence or absence of CAN (2 classes), absence of CAN, early and definite CAN (3 classes) or adding the severe class.

The motivation to use HRV data is that HRV data are more often available and are easier to obtain in clinical practice than the Ewing battery features. HRV measures also provide many more variables compared to the five attributes in the Ewing battery.

Testing for CAN in people with diabetes has progressed from the simple Ewing tests to HRV analysis and analysis of more complex ECG features [22,23,24,25]. These new analysis techniques determine the changes in heart rate or ECG characteristics over time and have resulted in CAN being diagnosed before clinical signs are present with greater accuracy and reproducibility (see [26,27,28] for more details and further references). HRV analysis involves determining the interbeat intervals (RR intervals) between successive QRS complexes on an ECG. The next table (Table 2) summarizes various HRV analysis methods, see [29,30,31].

Table 2. HRV analysis methods.

HRV method	Measure	Description
Time Domain	SDNN	The standard deviation of the NN intervals
	RMSSD	The square root of the mean squared difference of the NN intervals
Frequency Domain	Total Power	Variance of N-N intervals over the temporal segment (freq ＜ 0.4)
	VLF	Power in very low frequency range (freq ＜ 0.04)
	LF	Power in low frequency range (freq 0.04 to 0.15)
	HF	Power in high frequency range (freq 0.15 to 0.4)
Nonlinear	SD1, SD2	The standard deviations perpendicular to and along the line-of-identity of the Poincaré plot
	ApEn	Approximate entropy
	SampEn	Sample entropy
	D2	Correlation dimension
	DFA	Detrended fluctuation analysis:
		α1 - Short-term fluctuation slope; α2- Long-term fluctuation slope
	ShanEn	Shannon entropy
	Rényi En	Rényi entropy

 | Show Table

DownLoad: CSV

Nonlinear HRV measures have become popular in recent times as they are more robust against the nonstationarity and nonlinearity characteristics of the RR tachogram and are able to detect how aging and pathological conditions affect the nonstationarity of the signal [32,30]. Non-linear HRV features such as Detrended Fluctuation Analysis (DFA), estimate complexity inherent in the signal. The correlation dimension (D₂) can also be applied [33]. Several entropy measures have been proposed such as approximate entropy, sample entropy, and tone-entropy [34,35]. These measures have subsequently led to the multiscale entropy measures including the multi-scale Rényi entropy, which is a generalization of the Shannon entropy [13].

3. Methods

To prepare data for experiments in this paper we kept the Ewing attributes and the HRV attributes in the data files and deleted all other attributes including biomarkers such as cholesterol profile, co-morbidities, HbA1c and blood sugar levels, which are associated with CAN in [8,36].

3.1. Classifiers

We used the following classifiers implemented in WEKA.

· Bayes Net (BN) provides network structure, conditional probability distributions and other data structures and facilities for Bayes Network learning. We used it with the K2 search algorithm.

· J48 classifier is one of the most efficient decision tree classifiers in WEKA employing C4.5 tree as described in [37].

· Random Forest (RF) applies a collection of random trees [38]. RF is also implemented in WEKA and Rattle [39].

· SMO trains a support vector classifier by applying sequential minimal optimization [40].

· Decorate ensemble based on Random Forest (DRF). This classifier outperformed several other ensemble classifiers with two tiers for the classification of CAN in [7]. Decorate constructs special artificial training examples to build diverse ensembles of classifiers.

Let us refer to [41,42,43,44] for preliminaries on these methods, and to [39,45] for more information on Random Forest.

3.2. Attribute Selection Methods

Wrapper Subset Evaluator (WSE) in WEKA searches for a set of attributes tailored for the use with a particular classifier. It evaluates attribute sets by using a classifier indicated as an input parameter for the Wrapper as explained in [46]. WSE employs a number of search algorithms based on different strategies. The results produced by WSE were then compared with three search methods:

· Wrapper Subset Evaluator with Best First Search (WSE-BFS). This method produced good results in [46].

· Wrapper Subset Evaluator with Genetic Search (WSE-GS) using genetic algorithm described in [47] and applied previously by our research group [48].

· Wrapper Subset Evaluator with Particle Swarm Optimization search algorithm (WSE-PSO) explores the attribute space using the Particle Swarm Optimization (PSO) algorithm as explained in [49].

In addition, we applied the following attribute ranking methods implemented in WEKA.

· Classifier Attribute Evaluator (CAE) evaluates attributes applying a classifier and ranks all attributes according to their significance for the particular classifier.

· Information Gain Attribute Evaluator (IGAE) assesses attributes by calculating the information gain.

Each of the ranking methods orders all attributes according to its own algorithm for determining significance of the attributes. For each classifier, we used the resulting ranking to decide which attributes to select. To this end, we started with the most significant attribute and determined the effectiveness of the classifier based on this attribute alone. Then we incrementally added the next most significant attribute on the ranking list, and continued this process until the effectiveness of the classifier stopped improving.

4. Multi-layer Attribute Selection and Classification Algorithm

MLASC algorithm has three layers explained in this section. Note that this is the first article where multi-layer attribute selection method is introduced and included in one of the layers of the main algorithm.

Layer 1 of MLASC begins by extracting a large collection of HRV attributes. A C++ program was written by the author to extract attributes. In particular, Rényi entropy was calculated with α = −5, −4, …, 4,5, and estimated probabilities p_i of sequences of RR intervals of length π = 1,2,4,8,16.

Layer 2 of MLASC uses a hierarchical ensemble, Bagging of Decorate ensemble based on Random Forest (BDRF) and applies a new attribute selection method introduced in this paper.

The BDRF classifier was introduces in our previous article [7], where it produced the best outcomes for a different task. It combines Bagging with Decorate and is based on Random Forest. Bagging (bootstrap aggregating) generates a collection of new sets by resampling the given training set at random and with replacement. New classifiers are then trained, one for each of these new training sets. They are amalgamated via a majority vote.

Our new attribute selection method is a Double Wrapper Subset Evaluator with Particle Swarm Optimization (DWSE-PSO). The first sublayer of this method applies WSE-PSO to select the first subset of attributes for BDRF. All instances are then classified, using the selected first subset of features. The misclassified instances, where incorrect outcomes had been produced, are collected in a new dataset. After that in the second sublayer WSE-PSO is applied to the new dataset of misclassified instances again using BDRF, which results in a second subset of attributes. Two subsets are then combined together to compose the final set of chosen attributes.

Layer 3 of MLASC extracts from the whole training set all attributes chosen in Layer 2 and trains the BDRF classifier using all these attributes. The BDRF is then applied to validate the set using all attributes chosen in Layer 2.

5. Simulation Details

To prepare the dataset for experiments, we selected all patients in the DiScRi database with available HRV attributes and recorded all DiScRi parameters for these patients in a csv file. We have kept only the HRV attributes and the Ewing battery results, deleting all other features, and created 3 copies of the data file to explore the progression of CAN indicated in the DiScRi database as “no CAN”, “early CAN” and “definite CAN”. In the first copy, a two-class paradigm was investigated with the last column containing the class value for classification including “definite CAN” and “no CAN”. In the second copy, we added three CAN classes as “no CAN”, “early CAN” or “definite CAN”. In the third copy, as the last column we added one of four CAN classes including “no CAN”, “early CAN”, “definite CAN” or “severe CAN”. In order to enable all classifiers available in WEKA Explorer to process these three files, the files were reformatted into ARFF format, which is the standard format used by all classifiers in WEKA. These three files were used in all experiments presented in the paper. For each classifier in every test, the data file was divided in two parts of equal size. The first part was used to apply feature selection strategy for the classifier. After the optimal set of features had been selected using the first part, we selected the attributes of the chosen optimal set of features in the second part of the file. This second part with chosen features was then used to determine the effectiveness of the classifier.

6. Results

We carried out experiments to determine the optimum classifier algorithm and an optimum attribute selection algorithm using WEKA [42,43] and R [39,45] in all experiments.

To compare the effectiveness of classifiers we used the Area Under Curve (AUC), as the measure of performance. The value of AUC is included in the output of all classifiers in WEKA. Since WEKA is open source software, its source code has also always been available for inspection by the researchers. We refer the readers to [42] and [43] for more details.

We carried out experiments to compare the effectiveness of MLASC algorithm with other counterpart methods, and also with several simplified versions of MLASC, where the main ingredients of MLASC are replaced by simpler alternatives. In the figures with the results we include the average values of AUC obtained in the 5 × 2 cross-validation (5 × 2 CV), because it is shown in Table 1 of [52] that the average value is a standard measure of performance which is one of the most frequently used indications of effectiveness of classifiers used in the literature. To prevent over-fitting we used 5 × 2 cross-validation introduced and recommended in [51] for comparison of classifiers. This method carries out five iterations of twofold cross-validation. The results of cross-validation implemented in WEKA are included in the output of all classifiers automatically.

Figure 1 presents the results of our experiments comparing the performance of MLASC with other algorithms. In this figure BDRF, BN, DRF, J48, RF, SMO are applied to process the data directly and not as a part of larger multi-layer process. The 5 × 2 CV test recommended in [52] established that MLASC produced better outcomes and these differences are significant at the 95% confidence level.

Figure 1. Comparison of MLASC and other classifiers.

DownLoad: Full-Size ImgPowerPoint

Figure 2 presents the results of experiments comparing the performance of MLASC with other simplified versions, where both BDRF and DWSE-PSO, or at least one of these main ingredients of MLASC, are replaced by simpler alternatives.

Let us introduce notation to refer to such simpler versions in the diagram. If A is an attribute selection method and C is a classifier, then we use the symbol SV (A, C) to denote a simplified version of MLASC algorithm, where BDRF is replaced by C and the application of DWSE-PSO is replaced by A. Our experiments compare MLASC with the following simplified versions: SV(WSE-BFS, BN), SV(WSE-BFS, J48), SV(WSE-BFS, RF), SV(WSE-BFS, SMO), SV(WSE-BFS, DRF), SV(WSE-BFS, DRF), SV(WSE-GS, BN), SV(WSE-GS, J48), SV(WSE-GS, RF), SV(WSE-GS, SMO), SV(WSE-GS, DRF), SV(WSE-GS, BDRF), SV(WSE-PSO, BN), SV(WSE-PSO, J48), SV(WSE-PSO, RF), SV(WSE-PSO, SMO), SV(WSE-PSO, DRF), SV(WSE-PSO, BDRF), SV(CAE, BN), SV(CAE, J48), SV(CAE, RF), SV(CAE, SMO), SV(CAE, DRF), SV(CAE, BDRF), SV(IGAE, BN), SV(IGAE, J48), SV(IGAE, RF), SV(IGAE, SMO), SV(IGAE, DRF), SV(IGAE, BDRF), SV(DWSE-PSO, BN), SV(DWSE-PSO, J48), SV(DWSE-PSO, RF), SV(DWSE-PSO, SMO), SV(DWSE-PSO, DRF). This list contains all simplified versions of MLASC, where at least one of the main ingredients is replaced by a simpler version. Figure 2 presents the outcomes of these classifiers and MLASC grouped according to the attribute selection applied at Layer 2: a) applying WSA with best first search, b) genetic search, c) PSO, d) DWSA with PSO, e) CAE ranking, and f) IGAE ranking.

Figure 2. AUC achieved by MLASC classifier and simplified versions grouped according to the attribute selection employed at Layer 2: (a) WSA with BFS, (b) WSA with genetic search, (c) WSA with PSO, (d) DWSA with PSO, (e) CAE ranking, (f) IGAE ranking.

DownLoad: Full-Size ImgPowerPoint

In general, it can be observed that selecting features using PSO (Figure 1c) provided the highest values of AUC and that the BDRF classifier provided the highest accuracy, regardless of which method was used for feature selection. The best outcome was obtained by the MLASC algorithm with an AUC equal to 0.96. The 5 × 2 CV test recommended in [52] established that these results are significant at the 95% confidence level.

7. Conclusion

In this article, we introduced Multi-Layer Attribute Selection and Classification (MLASC) algorithm for the diagnosis of CAN progression based on HRV attributes. We carried out a comprehensive set of experiments comparing the effectiveness of MLASC with other simplified versions and counterpart methods. The results obtained demonstrated that the MLASC algorithm produced best results. Diagnosis by this method achieved an AUC level of 0.96 when applying DWSA with PSO attribute selection. Classification of presence of CAN and CAN progression remained above 0.90 in the 3 and 4 class paradigms, which is high enough for clinical diagnostics.

The variation of success using different classifier algorithms when using a separate subset for each classifier supports the findings of Kohavi [46], by showing that the attributes chosen should be regarded as part of the classifier algorithm. The current application involving the Rényi entropy addresses methodological issues inherent in previous multiscale applications [52].

The results obtained by applying MLASC algorithm show that recognizing CAN progression from HRV data can provide an effective test which would benefit diabetic patients by early diagnosis and better treatment, and consequently a reduction in hospitalization and length of stay [53].

Acknowledgments

The authors wish to thank the three reviewers for comments that have helped to improve this article. We acknowledge the many students that have contributed to the data analyzed and reported in this paper. Part of the research presented was funded by the Diabetes Australia Research Trust, Albury and Tallangatta Council Funding, Deakin-Ballarat Collaboration Funding and Charles Sturt University Research Compacts grants. Roche Australia Pty is gratefully acknowledged for provided glucose meters and test strips.

Conflict of Interest

All authors declare no conflicts of interest in this paper.