DWFS: A Wrapper Feature Selection Tool Based on a Parallel Genetic Algorithm

Many scientific problems can be formulated as classification tasks. Data that harbor relevant information are usually described by a large number of features. Frequently, many of these features are irrelevant for the class prediction. The efficient implementation of classification models requires identification of suitable combinations of features. The smaller number of features reduces the problem’s dimensionality and may result in higher classification performance. We developed DWFS, a web-based tool that allows for efficient selection of features for a variety of problems. DWFS follows the wrapper paradigm and applies a search strategy based on Genetic Algorithms (GAs). A parallel GA implementation examines and evaluates simultaneously large number of candidate collections of features. DWFS also integrates various filtering methods that may be applied as a pre-processing step in the feature selection process. Furthermore, weights and parameters in the fitness function of GA can be adjusted according to the application requirements. Experiments using heterogeneous datasets from different biomedical applications demonstrate that DWFS is fast and leads to a significant reduction of the number of features without sacrificing performance as compared to several widely used existing methods. DWFS can be accessed online at www.cbrc.kaust.edu.sa/dwfs.


Introduction
In the last decade, the leading high-throughput experimental techniques in biology, such as next generation sequencing, mass spectrometry, array-based methods and others, let to the massively increasing volume and dimensionality of the produced data in this field. This expansion in volume and dimensionality of data is also occurring in many other domains such as, for example, web content, social networks, graphical information systems, etc. From a Data Mining perspective, the need for faster, more reliable and more cost-effective classification models based on such data requires the extraction of a smaller and optimized set of features that can be obtained by removing largely redundant, irrelevant and unnecessary features for the class prediction [1]. A reduced number of features may also in some cases improve the classification performance [1]. In addition, the reduced number of features may lead to a better description of the underlying process from which data is generated and, thus, may contribute to better interpretation of the results [2]. Consequently, the feature selection (FS) problem is a From the category of filters, mRMR [14] belongs to multivariate filtering methods and a web tool is available. However, the web-tool has some limitations since only input files with size less than 2 MB are supported. FSelector [15] presents an extensive collection of existing filtering algorithms implemented in a Ruby package whereas, the CBFS algorithm [16] introduces the clearness-based scoring scheme for describing features that maximize separability between the classes.
There are also many dataset-specific FS methods, e.g. ArrayMining [17] that targets selection of genes from microarray datasets, Multi-Relief that [18] ranks crucial residues for protein functionality or MetaboloAnalyst [19] that selects important metabolites from metabolic network data. Table 1 summarizes the advantages and disadvantages of the above-mentioned FS techniques.

DWFS web-tool
Search strategy description. DWFS takes as input text files with tab-delimited format where rows represent data samples and columns represent features. The last column represents class labels. The multiclass classification tasks (classes> = 2) are supported. The wrapper FS search strategy is based on GA that is a well-known heuristic optimization technique inspired by the principles of natural selection [6]. Initially, a starting population Θ 0 = {θ 01 , θ 02 ,. . ., θ 0r } (the collection of initial solutions) of individuals θ 0i ,i = 1,2,. . .,r, where r is the number of individuals, is formed. The individuals represent chromosomes. Note that in our case the chromosomes are always described by vectors of length n whose components correspond to genes (features); n is the number of all features. The values of the components are binary 1s and 0s, where '1' encodes the existence and '0' the absence of a gene. In the initial population chromosomes are initialized by a random selection of 1s and 0s. In GA, the population changes in each evolutionary cycle forming different generations. An evolutionary cycle k, is characterized by the population Θ k = {θ k1 , θ k2 ,. . ., θ kr }; θ ki represents the i th chromosome in population Θ k . In every evolutionary cycle, a fitness value is computed based on a predefined fitness (objective) function. DWFS uses the fitness function described by Equation (1), where Performance(θ ki ) is the classification performance achieved by using chromosome θ ki , size(θ ki ) is the sum of binary 1s in chromosome θ ki , and α is a weighting parameter (default 0.15).
Note that in our case the actual Performance function is determined by the selected criterion from the web-tool. Then, GA applies two operations, crossover and mutation, and the best solutions are selected to survive to the next generation. The optimization process is run for a specified  [12] mRMR no yes filter [14] Arraymining no yes filter-wrapper-embedded [17] WEKA no no filter-wrapper [13] Multi-Relief no yes filter [18] MetaboloAnalyst no yes filter [19] FSelector no no filter-wrapper [15] CBFS no no filter [16] FeaLect no no wrapper [35] Proposed method DWFS yes yes hybrid-filter-wrapper number of evolutionary cycles, during which time GA attempts to increase the value of the fitness function. The chromosome that produces the largest value of the fitness function during this specified number of evolutionary cycles is considered an optimized solution and it may or may not be the global (optimal) solution. It is difficult to know if the obtained optimized solution is the global one. The procedure terminates whenever termination criteria are satisfied [6]. Implementation. DWFS search strategy is based on the PGAPack software [20], which is a collection of libraries that execute all GA steps. PGAPack deploys a parallel master/slave single population GA using the message-passing interface (MPI). Master node stores the initial population and applies the GA operations. Slave nodes are responsible for evaluating the fitness function for every chromosome. The master/slave paradigm is very efficient in terms of execution time because the most costly part is the fitness function evaluation, whereas the communication overhead is minimal [21]. Fig. 1 demonstrates the GA algorithm under the master/slave framework of PGAPack. In principle, increasing the population size does not necessarily improve the classification performance as there is no guarantee that the algorithm will be able to find any better solution than it does with the original population size. We use a population size of 100 individuals, which is a good compromise considering the trade-off between the performance and over-fitting [22]. Note, that DWFS master/slave model was deployed on a cluster with nodes having 64-cores. In our experimentation in one experiment 63 cores are used for replacing the old solutions for every new generation and the 64 th core (i.e. the master) was responsible only for executing the GA operators. In the first population we initialize genes to 0 and 1 randomly with equal probability and we ensure that a chromosome containing all genes set to 1 is present (i.e., this represents the solution that contains all features). In our implementation we adopt constant rates for crossover and mutation. The crossover operator produces new offspring by combining genes of two chromosomes (they are considered parents in the GA terminology). We use two-point crossover technique that swaps all genes between two crossover points for both chromosomes. The following example illustrates the two-point crossover operation over a pair of chromosomes with 9 genes. The underlined features in the chosen parent solutions indicate the points within which crossover will take place (i.e. genes between these points will be crossed-over to generate child solutions). A pre-defined crossover rate will determine the number of genes where crossover will take place.
Mutation modifies specific genes of a given chromosome. We use the simple bit-string mutation that flips specific genes at random positions. Both GA operations simulate parts of the natural evolutionary process and are used to control transfer of information between different generations of population. The evolution of the population follows a binary tournament selection, where the best two chromosomes in the population are selected to serve as parents of new offspring. The GA operators including selection, mutation and crossover are then carried over in an iterative process through a predefined number of generations.
Tuning search strategy parameters. Depending on the application requirements users can specify mutation and crossover rates, as well as the number of generations for the optimization procedure. After experimentation with different mutation and crossover rates we observed that values close to 0.8 for crossover and 0.01 for mutation are effective in various classification tasks (see also [5]). These values we selected as default for DWFS. The effects of extensive experimentation with various mutation and crossover rates for all datasets are captured in Fig. 2 and Fig. 3. Both Fig. 2 and Fig. 3, indicate a considerable sensitivity of GA towards mutation and crossover rates. Overall, increasing either one of these rates leads to not exploiting surrounding regions of good solutions and mostly exploring the search space in a random manner. This is apparent in Fig. 2 and Fig. 3 for most of the datasets. More results about the effects of mutation and crossover rates over stability and reduction in size of selected features are provided in S1 Table and S2 Table. The default number of evolutionary cycles in DWFS is 100. Also, in order to reduce the processing time we included additional stopping criteria in which we consider that the population has reached the steady-state (i.e. it has converged) if there is no difference in the fitness function value for 50 consecutive generations.
Modifying the fitness function. DWFS supports different validation methods and users can choose between automatic split of data to training and testing sets, or the custom option in which they can upload explicitly the training and testing datasets. DWFS offers different performance metrics for the fitness function from Equation (1) including accuracy, geometric mean of sensitivity and specificity (GMean), F1-score and Mathews Correlation Coefficient (MCC). Possibility to select different performance metrics for the optimization process increases the utility of DWFS. For instance, GMean is preferable in cases where the data classes are highly imbalanced [23]. Three different classification models can be used: a/ Naïve Bayes Classifier (NBC-in-house implementation), b/ K-Nearest-Neighbors (KNN) available from the AlgLib library (http://www. alglib.net/) [24], and c/ and the combination of a/ and b/. In principle, selecting more than one classifier tends to minimize potential bias towards a particular type of classifier. The idea behind this is relatively simple: features selected using a single classifier in the wrapper FS usually result in higher performance for that particular classification model and weaker performance for other classifiers. Including a combination of classifiers in the wrapper FS setting, may lead to the selection of a feature subset that is less biased, but this is not guaranteed.
Hybrid combination of filtering and wrapping. When the number of features is very large, the GA becomes very computationally demanding and its run time can be prohibitive. To reduce the time required for training, DWFS supports a two-phase hybrid combination of filtering and wrapping. In the first phase, filtering is applied as a pre-processing step and the top-ranked features are selected based on a pre-defined but tunable threshold. In the second phase, the GAbased wrapper FS is applied to the selected features, which in practice means that the FS process uses a much smaller search space instead of the original large search space. This heuristic optimization technique may or may not result in combinations of features that have high discriminative power. However, this cannot be known in advance. In practice, however, the experiments show that this heuristics performs well and that it is fast and memory efficient for datasets described by large number of features and having large number of samples. Of course, the user may choose not to use filtering at all. Current implementation of DWFS incorporates the following well-known and widely-used filtering techniques: minimum redundancy maximum relevance (mRMR) [14], joint mutual information (JMI) [25], conditional mutual information maximization (CMIM) [26] and interaction capping (ICAP) [27]. mRMR selects features based on the minimum redundancymaximum relevancy criterion. Both JMI and CMIM are based on different notions of mutual information (joint vs conditional), whereas ICAP selects features based on the concept of "interaction" of the class label with different sets of features. Recently, a study has proposed a unifying framework for information theoretic feature selection that derive the previous criteria from conditional likelihood of the training labels [28]. In particular, the scoring criteria including mRMR, JMI, CMIM and ICAP can be derived from the following function F described by Equation (2).
The first term is a likelihood ratio between the true distribution p(y|x θ of a class given some particular selected features x θ and the predicted class distribution q(y|x θ , τ) given the selected features and a model parameter τ, averaged over the input space. Given the selected features X θ , IðX ỹ; YjX y Þ defines conditional mutual information between the class label Y and the unselected features X ỹ. The final term H(Y|X) is the conditional entropy of labels Y given all original set of features X. In [28] a rigorous derivation of many information theoretic measures for feature selection from Equation (2) is provided to eventually reveal either a direct or conditionally dependent relationship between the target label and the selected features of dataset.

Experimental setup
Datasets. In our experiments we use nine datasets related to different biomedical problems. Table 2 summarizes the number of samples, number of features, ratio between positive and negative samples and corresponding data sources. We chose on purpose datasets from a broad spectrum of different problems to demonstrate that DWFS is a general tool for FS and that it performs well over a variety of different types of data.
The methods used in comparison. To estimate the effectiveness of DWFS and quantify the contribution of the filtering (pre-processing) step, we select features using wrapper FS only (denoted as DWFS), wrapper FS combined with mRMR (denoted as mRMR+DWFS) and wrapper FS combined with JMI (denoted as JMI+DWFS). To provide better cross-benchmarking results, we compare the classification performance of DWFS and its variants with the most effective wrappers and the filtering approaches presented in Table 1. Specifically, we include in the comparison WEKA (version 3.6.6) wrapper that uses an integrated method of forward and backward search (Bi-directional BestFirst), FST3 that implements a powerful wrapper-based method called sequential forward floating search (SFFS), as well as mRMR and JMI filters. In our comparison analysis, we view WEKA and FST3 wrappers as representative implementations for other tools that use similar sophisticated approaches for FS. Moreover, mRMR and JMI filters are also considered representatives of a wider class of advanced filtering approaches. In particular, JMI stands as a representative for single-objective filtering methods that examine joint relevance between features with the target variable, while mRMR represent methods that consider multi-objective formulation for filtering features with high relevance to the target and minimum redundancy between each other. In other words, comparison with WEKA and FST3 wrappers and mRMR and JMI filters, is expected to reveal how our approach may compare to other filters or wrappers. We used two baselines, one where the classification performance is obtained utilizing all features (the initial/original feature vector), and the other that uses top 10% of features with the highest class correlation (we call it Correlation-baseline). We did not include in the comparison the datasetspecific FS tools such as Arraymining, Multi-Relief and MetaboloAnalyst, or tools that at the time of this study were not fully functional (e.g. that crushed during experiments).
Performance assessment. Classification performance is assessed using nested 5-fold crossvalidation. This technique splits data to five folds and for every fold FS is applied. To evaluate the performance of the model inside a single fold we split again this particular fold to five folds, assign the performances and average the classification performance of the folds. Such cross-validation setup has been proposed in [29] for avoiding potential over-fitting effects. Given the large size of our experimental datasets, 5-fold setup seems to be a proper choice for computing a representative (i.e. not biased) estimate while being computationally faster compared to a setup with a larger number of folds [30,31].
Classification performance is measured using five classification performance metrics, namely: Sensitivity that measures the true positive rate; Specificity that measures true negative rate; GMean that measures the combined effect of Sensitivity and Specificity; Positive Predictive Value (PPV) that measures proportion of correctly predicted positives out of all predicted positives; and F1-measure that estimates the combined effect of PPV and Sensitivity. These performance metrics are frequently used performance indicators for classification systems [23,32]. Let TP be the number of true positives, FP the number of false positives, TN the number of true negatives and FN the number of false negatives, Equations 3 to 7 show formulae of these performance metrics: In addition to classification performance indicators, we also include a measure of robustness called stability [33]. The stability S(θ i , θ j ) between two selected feature subsets θ i , θ j over two different folds of the data is defined by Equation (8). The sizes of these two feature subsets are represented by k i , k j , respectively. The size of common features (i.e. the number of elements in the intersection between the two sets) is c: We compute the pair-wise scores in Equation (8) for different folds (i.e. 5 folds in our setup) and report their average. The stability measure in Equation (8) has a range of (−1,1] were a value of 0 indicates an output from a random FS. A positive range states that the FS method is more stable than a random one and a negative value indicate the opposite [33].

Comparison analysis
Classification performance and ranking. We do not know in advance the specifics of the classifier implementation that the different tools have. For this reason, we select features with each of the wrapper-based tools included in the comparison, using their in-house implementations. However, we evaluate all programs using the same MATLAB R2012b implementations as these were not part of any of the wrapper-based FS tools used in the comparison analysis. This ensures fairness in the comparison. Finally, for JMI, mRMR and Correlation-baseline, we select the 10% of the top-ranked features.
Tables 3-20 present the actual average and standard deviation of classification performance for all studied methods using NBC and KNN classifiers. We observe that, using different performance indicators the studied programs show different advantages and disadvantages and perform differently on different datasets. In terms of GMean, DWFS is comparable in many cases to the other methods. For the KNN classifier, in 6 out of 9 datasets, DWFS hybrid with a filtering method enhanced the filtering performance. This, nevertheless, is not the case for the NBC classifier where such enhancement happens in 4 out of 9 datasets. In 3 out of 18 cases (9 datasets with 2 classifiers), WEKA achieved the best results with the NBC and KNN classifiers. The same was also the case for FST3. Interestingly, the Correlation-baseline achieved the best result for the Microarray dataset with the KNN classifier (Table 16). These experimental results suggest that DWFS or one of its variants are better suited for datasets with large number of samples (>200).
To draw conclusions about their relative performance, we rank the FS methods following the method presented in [34] as follows: for every dataset and for every performance indicator, the FS method that achieves the best result is ranked first and obtains 1 point, the second best obtains 2 points and so on, up to 9 points. The overall ranking for each FS method is reflected in the sum of obtained points across all datasets and all performance metrics using KNN and NBC, which we call the final ranking score. At the end, the FS method with the smallest sum of ranking scores is ranked at position 1 and so on. One can also divide the final ranking score with the total number of tests and obtain the average rank position. This however does not change the final ranking order. The ranking that we present and conclusions that we derive from it are relative: they are based on the collection of methods that are compared, the collection of criteria used, and the collection of datasets used. Based on our experiments, we observe that the simple DWFS wrapper FS has the best (smallest) average rank across all tests. Yet, for a specific performance metric or a particular dataset, this need not be the case. WEKA is ranked second followed by JMI combined with DWFS. Table 21 presents these ranked positions. Note that in the ranking process all the tested cases equally contribute to the ranking score.
Effectiveness of selected features and feature subset size reduction. Next, to provide more insights about the performance of the studied methods we compute the ratio of classification performance over the number of selected features. This ratio measures the effectiveness of the FS process and depicts the maximum classification performance that can be achieved using the reduced number of features selected by an FS method. As performance indicator, we choose GMean since it combines two other important performance metrics (note that other selections were also possible). The results for the NBC and KNN classifiers are presented in Fig. 4. DWFS and its variants achieve the best results in 10 out of 18 tested cases. For the remaining 8 cases Run time. In Table 24 and Table 25 we report the run time required for selecting features and evaluating classification performance. From this comparison we exclude mRMR and JMI since they are filtering methods and their run time is much shorter than for any of the wrapper models. We observe that on average and across different datasets, DWFS is the slowest. However, the combination of DWFS wrapper with mRMR reduces significantly the run time and achieves the fastest wrapper setting compared to WEKA and FST3. Stability results. Stability quantifies the degree of presence of selected features over different splits of the training data. Features that appear more frequently in selected subsets may represent more significant ones for the class prediction, but this factor is not explicitly correlated to classification performance. Note, that most stability metrics are defined for feature subsets of equal size. However, this is not the case for wrappers, where the final selected subsets have different size for every split of data. In our case, we use a definition proposed by Lustgarten et al. [33] (see Equation (8)), which meets our need for measuring stability between feature subsets of different sizes. It should be noted that higher stability may not necessarily reflect a better classification performance. Fig. 5 depicts the results of the comparison analysis using the NBC classifier for all the methods used and across all studied datasets. We observe that filtering approaches are more stable than wrapper-based ones in FS. WEKA with forward and backward search achieves more stable results than the other wrappers. For specific datasets (i.e. WDBC, TIS and Pre-miRNAs), DWFS combined with any of the filters achieves the most stable results.

Conclusion
In this study we present DWFS, a web-based FS tool based on the wrapper model combined with a randomized search strategy. DWFS offers a variety of options to enhance the optimized FS including incorporation of filtering methods as pre-processing step and many other options for tailoring the objective function according to application requirements. From our extensive experimentation with several benchmark datasets from different biological and biomedical problems we observe that, DWFS integrated in a hybrid setup with a filtering approach for FS is capable of reducing significantly the number of features without sacrificing classification performance. Moreover, using different performance indicators we show that DWFS performs well relative to the existing FS methods. We hope that DWFS will find good use in the analysis of many complex biological and biomedical data complementing other available tools for FS that biologists may use.
Supporting Information S1 Table. Effects of changing crossover rate on stability of NBC and reduction in the size of selected features. (XLS) S2 Table. Effect of changing mutation rate on stability of NBC and reduction in the size of selected features. (XLS)