Abstract
Background
Classification and variable selection play an important role in knowledge discovery in highdimensional data. Although Support Vector Machine (SVM) algorithms are among the most powerful classification and prediction methods with a wide range of scientific applications, the SVM does not include automatic feature selection and therefore a number of feature selection procedures have been developed. Regularisation approaches extend SVM to a feature selection method in a flexible way using penalty functions like LASSO, SCAD and Elastic Net.
We propose a novel penalty function for SVM classification tasks, Elastic SCAD, a combination of SCAD and ridge penalties which overcomes the limitations of each penalty alone.
Since SVM models are extremely sensitive to the choice of tuning parameters, we adopted an interval search algorithm, which in comparison to a fixed grid search finds rapidly and more precisely a global optimal solution.
Results
Feature selection methods with combined penalties (Elastic Net and Elastic SCAD SVMs) are more robust to a change of the model complexity than methods using single penalties. Our simulation study showed that Elastic SCAD SVM outperformed LASSO (L_{1}) and SCAD SVMs. Moreover, Elastic SCAD SVM provided sparser classifiers in terms of median number of features selected than Elastic Net SVM and often better predicted than Elastic Net in terms of misclassification error.
Finally, we applied the penalization methods described above on four publicly available breast cancer data sets. Elastic SCAD SVM was the only method providing robust classifiers in sparse and nonsparse situations.
Conclusions
The proposed Elastic SCAD SVM algorithm provides the advantages of the SCAD penalty and at the same time avoids sparsity limitations for nonsparse data. We were first to demonstrate that the integration of the interval search algorithm and penalized SVM classification techniques provides fast solutions on the optimization of tuning parameters.
The penalized SVM classification algorithms as well as fixed grid and interval search for finding appropriate tuning parameters were implemented in our freely available R package 'penalizedSVM'.
We conclude that the Elastic SCAD SVM is a flexible and robust tool for classification and feature selection tasks for highdimensional data such as microarray data sets.
Background
Classification and prediction methods play important roles in data analysis for a wide range of applications. Frequently, classification is performed on highdimensional data, where the number of features is much larger compared to the number of samples ('large p small n' problem) [1]. In those cases, classification by Support Vector Machines (SVM), originally developed by Vapnik [2], is one of the most powerful techniques. The SVM classifier aims to separate the samples from different classes by a hyperplane with largest margin.
Often we do not only require a prediction rule but also need to identify relevant components of the classifier. Thus, it would be useful to combine feature selection methods with SVM classification. Feature selection methods aim at finding the features most relevant for prediction. In this context, the objective of feature selection is threefold: (i) improving the prediction performance of the predictors, (ii) providing faster and more costeffective predictors, and (iii) gaining a deeper insight into the underlying processes that generated the data.
Three main groups of feature selection methods exist: filter, wrapper and embedded methods [1,36]. Filter methods simply rank individual features by independently assigning a score to each feature. These methods ignore redundancy and inevitably fail in situations where only a combination of features is predictive. Also, if there is a preset limit on the number of features to be chosen (e.g. top 10 features), this limit is arbitrary and may not include all informative features. Because of these drawbacks, the filter methods are not included in this work.
Connecting filtering with a prediction procedure, wrapper methods wrap feature selection around a particular learning algorithm. Thereby, prediction performance of a given learning method assesses only the usefulness of subsets of variables. After a subset with lowest prediction error is estimated, the final model with reduced number of features is built [5]. However, wrapper methods have the drawback of high computational load, making them less applicable when the dimensionality increases. Wrapper methods also share the arbitrariness of filter methods in feature selection.
The third group of feature selection procedures are embedded methods, which perform feature selection within learning classifiers to achieve better computational efficiency and better performance than wrapper methods. The embedded methods are less computationally expensive and less prone to overfitting than the wrappers [7].
Guyon [1] proposed the recursive feature elimination (RFE) method, which belongs to the wrapper methods. RFE iteratively keeps a subset of features which are ranked by their contribution to the classifier. This approach is computationally expensive and selecting features based only on their ranks may not derive acceptable prediction rules.
An alternative to SVM with RFE is to use penalized SVM with appropriate penalty functions. Penalized SVM belongs to embedded methods and provides an automatic feature selection. The investigation of the widely used family of penalization functions such as LASSO, SCAD, Elastic Net [810] and a novel proposed penalty Elastic SCAD in combination with SVM classification, is the objective of the paper. The ridge penalty [4] corresponds to the ordinary SVM, which does not provide any feature selection, is used as reference with respect to prediction accuracy.
Although feature selection methods can be applied to any highdimensional data, we illustrate the use of these methods on microarray gene expression data due to their relevance in cancer research. Data from microarray experiments are usually stored as large matrices of expression levels of genes in rows and different experimental conditions in columns. Microarray technology allows to screen thousand of genes simultaneously. Detailed reviews on the technology and statistical methods often used in microarray analyses are presented in [1113].
Since SVM is extremely sensitive to the choice of tuning parameters, the search for optimal parameters becomes an essential part of the classification algorithm [14]. The problem of choosing appropriate tuning parameters is discussed and an interval search technique from Froehlich and Zell [15] is proposed to use for SVM classification.
In this paper, we investigate the behaviour of feature selection SVM classifier techniques including commonly used penalization methods together with a novel penalization method, the Elastic SCAD. We compare them to SVM classification with and without recursive feature elimination (RFE [1]) for situations of 'large p small n' problems.
The RFE SVM is chosen as as a stateoftheart representative of feature selection methods in applications [16,17].
A simulation study is designed to investigate the behaviour of different penalization approaches. Publicly available microarray data sets are chosen for illustration purposes as applications on real highdimensional data.
Methods
Support Vector Machines
Suppose a training data set with input data vector x_{i }∈ ℝ^{p }and corresponding class labels y_{i }∈ {1, 1}, i = 1,..., n is given. The SVM finds a maximal margin hyperplane such that it maximises the distance between classes. A linear hyperplane can always perfectly separate n samples in n + 1 dimensions. Since we can assume that highdimensional data with p ≫ n is generally linear separable [6], increasing complexity by using nonlinear kernels is usually not needed. Thus, we use a linear SVM model throughout the paper.
The linear SVM separates classes by a linear boundary
where w = (w_{1}, w_{2},..., w_{p}) is a unique vector of coefficients of the hyperplane with w_{2 }= 1 and b denotes the intercept of the hyperplane. We use '·' to denote the inner product operator. The class assignment for a test data vector x_{test }∈ R^{p }is given by y_{test }= sign [f (x_{test})].
Soft margin SVM
Soft margin SVM allows some data points to be on the wrong side of the margin. To account for erroneous decisions, slack variables ξ_{i }≥ 0, i = 1,..., n are defined as the distance between a misclassified data point and the corresponding margin. For data points on the correct side of the margin ξ_{i }= 0, for data points inside the margin 0 < ξ_{i }≤ 1 and for misclassified data points ξ_{i }> 1. The sum of nonzero ξ_{i }is penalized with a cost parameter C and then added to the optimisation function penalty in the minimisation problem:
The optimisation problem (2) is called the soft margin SVM. The cost parameter C is a data dependent tuning parameter that controls the balance between minimizing the coefficients of the hyperplane and correct classification of the training data set. C is often chosen by cross validation. Problem (2) can be solved by using convex optimisation techniques, namely by the method of Lagrange multipliers [4]. Convex optimisation techniques provide a unique solution for hyperplane parameters w and b
where α_{i }≥ 0, i = 1,..., n are Lagrange multipliers. The data points with positive α_{i}, are called support vectors (SVs). All data points lying on the correct side of their margin have α_{i }= 0. Thus, they do not have any impact on the hyperplane, and we can rewrite Eq. (3) as
where the set of indices of the support vectors S is determined by S := {i : α_{i }> 0}.
The coefficient can be calculated from for any i with α_{i }> 0. In praxis, an average of all solutions for is used for numerical stability.
SVM as a penalization method
Hastie et al. [4] showed that the SVM optimisation problem is equivalent to a penalization problem which has the "loss and penalty" form
where the loss term is described by a sum of the hinge loss functions l (y_{i}, f (x_{i})) = [1  y_{i }f (x_{i})]_{+ }= max(1  y_{i }f (x_{i}), 0) for each sample vector x_{i}, i = 1,..., n. The penalty term is denoted as pen_{λ }(w) and can have different forms:
Ridge penalty
The penalty term for ordinary SVM uses the L_{2 }norm:
The L_{2 }penalty shrinks the coefficients to control their variance. However, the ridge penalty provides no shrinkage of the coefficients to zero and hence no feature selection is performed.
LASSO
The use of a L_{1 }penalization function is originally proposed by Tibshirani [8] for generalized linear models. The technique for parameter estimation with constraints is called LASSO ('least absolute shrinkage and selection operator'). Later, Bradley [18] adapted the L_{1}regularisation to SVM. Then, the penalty term has the form
As a result of singularity of the L_{1 }penalty function, L_{1 }SVM automatically selects features by shrinking coefficients of the hyperplane to zero.
However, the L_{1 }norm penalty has two limitations. First, the number of selected features is bounded by the number of samples. Second, it tends to select only one feature from a group of correlated features and drops the others.
Fung and Mangasarian [19] have published a fast L_{1 }SVM modification, the Newton Linear Programming Support Vector Machine (NLPSVM), which we use in our analyses.
Smoothly clipped absolute deviation penalty (SCAD)
The SCAD penalty is a nonconvex penalty function first proposed by Fan and Li [20]. Later, Zhang et al. [10] combined the SVM technique with the SCAD penalty for feature selection. The SCAD penalty function for a single coefficient w_{j }is defined as
where w_{j }, j = 1,..., p are the coefficients defining the hyperplane and a > 2 and λ > 0 are tuning parameters. Fan and Li [21] showed that SCAD prediction is not sensitive to selection of the tuning parameter a. Their suggested value a = 3.7 is therefore used in our analyses.
The penalty term for SCAD SVM has the form
The SCAD penalty corresponds to a quadratic spline function with knots λ at and aλ. For small coefficients w_{j}, j = 1,..., p, SCAD yields the same behaviour as L_{1}. For large coefficients, however, SCAD applies a constant penalty, in contrast to L_{1}. This reduces the estimation bias. Furthermore, the SCAD penalty holds better theoretical properties than the L_{1 }penalty [21].
Elastic Net
To overcome the limitations of LASSO, Zou and Hastie [9] proposed a linear combination of L_{1 }and L_{2 }penalties which they called Elastic Net:
The Elastic Net penalty provides automatic feature selection similar to L_{1}, but is no longer bounded by the sample size. Moreover, at the same time this penalty manages to select highly correlated features (grouping effect). Increasing λ_{1 }reduces the number of features of the classifier whereas for large λ_{2 }one observes better control of the grouping effect. Wang [22] adapted the Elastic Net penalty to SVM classification problems. Therefore, the Elastic Net SVM optimisation problem can be written as
where λ_{1}, λ_{2 }≥ 0 are the corresponding tuning parameters.
Elastic SCAD
Fan and Li [21] demonstrated the advantages of the SCAD penalty over the L_{1 }penalty. However, using the SCAD penalty might be too strict in selecting features for nonsparse data. A modification of the SCAD penalty analogously to Elastic Net could keep the advantages of the SCAD penalty, and, at the same time, avoid too restrictive sparsity limitations for nonsparse data.
We therefore propose a combination of the SCAD and the L_{2 }penalties. The new penalty term has the form
λ_{1}, λ_{2 }≥ 0 are the tuning parameters. We expect that the Elastic SCAD will improve the SCAD method for less sparse data. According to the nature of the SCAD and L_{2 }penalties, the Elastic SCAD should show good prediction accuracy for both, sparse and nonsparse data.
It can be shown that the combined penalty provides sparsity, continuity, and asymptotic normality when the tuning parameter for the ridge penalty converges to zero, i.e. λ_{2 }→ 0. The asymptotic normality and sparsity of Elastic SCAD leads to the oracle property in the sense of Fan and Li [21].
The Elastic SCAD SVM optimisation problem has the form
where λ_{1}, λ_{2 }≥ 0 are the tuning parameters.
Elastic SCAD SVM: Algorithm
By solving Eq. (9) the same problems as for SCAD SVM occur: the hinge loss function is not differentiable at zero and the SCAD penalty is not convex in w. The Elastic SCAD SVM objective function can be locally approximated by a quadratic function and the minimisation problem can be solved iteratively similar to the SCAD approach [10,21].
For simplicity, we rename the SCAD penalty from to . Accordingly, the firstorder derivative of the penalty is denoted by . Denote the penalized objective function in Eq. (9) by
For each i (with respect to the fact that ) the loss term can be split according to
Given an initial value (b_{0}; w_{0}) close to the minimum of A(b, w), we consider the following local quadratic approximations:
When w_{j0 }is close to zero, set ; otherwise use the approximation for the SCAD penalty
where due to symmetrical nature of the SCAD penalty w_{j} is used instead of w_{j}.
It can be shown that both approximations and their original functions have the same gradient at the point (b_{0}, w_{0}). Therefore, the solution of the local quadratic function corresponds approximately to the solution of the original problem.
The local quadratic approximation of A(b, w) has the form
By minimisation of A(b, w) with respect to w and b, terms without optimisation parameters w and b can be dropped (due to derivatives of constants):
To write the equations in matrix form we define:
Moreover, we define the matrix X = [1, x_{1},..., x_{p}], where 1 is the vector of 1s with length n and x_{j }is the jth input vector. Set
Minimizing A (b, w) is then equivalent to minimizing the quadratic function
The solution to Eq. (10) satisfies the linear equation system
The Elastic SCAD SVM can be implemented by the following iterative algorithm.
Step 1 Set k = 1 and specify the initial value (b^{(1)}, w^{(1)}) by standard L_{2 }SVM according to Zhang et al. [10].
Step 2 Store the solution of the kth iteration: (b_{0}, w_{0}) = (b^{(k)}, w^{(k)}).
Step 3 Minimize Ã (b, w) by solving Eq. (11), and denote the solution as (b^{(k+1)}, w^{(k+1)}).
Step 4 Let k = k + 1. Go to step 2 until convergence.
If elements are close to zero, for instance, smaller than 10^{4}, then the jth variable is considered to be redundant and in the next step will be removed from the model. The algorithm stops after convergence of (b^{(k)}, w^{(k)}).
Choosing tuning parameters
All SVM problems with or without feature selection use one or two tuning parameters which balance the tradeoff between data fit and model complexity. Since these parameters are data dependent, finding optimal tuning parameters is part of the classification task.
Fixed grid search
Tuning parameters are usually determined by a grid search. The grid search method calculates a target value, e.g. the misclassification rate, at each point over a fixed grid of parameter values. This method may offer some protection against local minima but it is not very efficient. The density of the grid plays a critical role in finding global optima. For very sparse grids, it is very likely to find local optimal points. By increasing the density of the grid, the computation cost increases rapidly with no guaranty of finding global optima. The major disadvantage of the fixed grid approach lies in the systematic check of the misclassification rates in each point of the grid. There is no possibility to skip redundant points or to add new ones.
When more parameters are included in the model, the computation complexity is increased. Thus, the fixed grid search is only suitable for tuning of very few parameters.
Interval search
Froehlich and Zell [15] suggested an efficient algorithm of finding a global optimum on the tuning parameter space using a method called EPSGO (Efficient Parameter Selection via Global Optimisation).
The main idea of the EPSGO algorithm is to treat the search for an optimal tuning parameter as a global optimisation problem. For that purpose, the Gaussian Process model is learned from the points in the parameter space which have been already visited. Thereby, training and testing of the GP is very efficient in comparison to the calculation of the original SVM models. New points in the parameter space are sampled by using the expected improvement criterion as described in the EGO algorithm [23], which avoids stacking in local minima. The stopping criteria of the EPSGO algorithm is either convergence of the algorithm or no change of the optimum during the last ten iterations.
Stratified cross validation
Using kfold cross validation, the data set is randomly split into k disjoint parts of roughly equal size, usually k = 5 or k = 10. In addition, the data is often split in a way that each fold contains approximately the same distribution of class labels as the whole data set, denoted by stratified cross validation. For each subset, one fits the model using the other k  1 parts and calculates the prediction error of the selected kth part of the data.
The case k = n is called leave one out cross validation (LOO CV). The choice of k determines a tradeoff between bias and variance of the prediction error. Kohavi [24] showed that tenfold stratified cross validation showed better performance in terms of bias and variance compared to 10 < k < n. Hastie et al. [4] recommended to perform five or tenfold cross validation as a good compromise between variance and bias. We used both, five and tenfold stratified cross validation for simulation study and real applications, respectively.
In the next two sections the application of penalized SVM classification methods are compared. We used simulated and publicly available data to investigate the behaviour of different feature selection SVMs. For all comparisons the R packages "penalizedSVM" [25] and "e1071" [26] were used which are freely available from the CRAN http://cran.rproject.org/ webcite, R version 2.10.1. The R package "e1071" is a wrapper for the wellknown LIBSVM software [27]. We used five and tenfold stratified cross validation in combination with interval search for tuning parameters as described above.
Results and Discussion
Simulation study
Simulation design
A comprehensive simulation study evaluating the performance of four feature selection SVM classifiers, L_{1 }SVM, SCAD SVM, Elastic Net SVM and Elastic SCAD SVM, was performed. We used the ordinary L_{2 }SVM algorithm with a liner kernel as a reference for prediction accuracy.
Two independent data sets are simulated: a training set for building the classifier and a test set for estimating of the prediction errors of classifiers. First, the training data is generated, and the optimal tuning parameters are found using fivefold stratified cross validation according to the interval search approach [15]. Then, the classification hyperplane is computed using the estimated tuning parameters. Finally, application of the classification rule to the test data provides the prediction characteristics such as misclassification error, sensitivity and specificity.
Training and test input data are represented by a data matrix X = {x_{i}}, i = 1,..., n, where x_{i }∈ ℝ^{p }describes feature patterns for the ith sample. The input data X follows a multivariate normal distribution with mean μ and covariance matrix Σ. The class labels Y = {Y_{i}}, i = 1,..., n are generated according to a logistic regression model
where β = {β _{1},..., β _{p}} is a vector of coefficients of a classifier and u_{i }are realisations of a variable following a U 0[1] distribution.
In our simulations the percentage of relevant features varies between 1% and 20%. Coefficients β _{j}, j = 1,..., p are always defined as
with equal numbers of positive and negative coefficients. The intersect β _{0 }is set to zero.
We also consider to have 'clumps' of correlated features. The clumpy dependency is supposed to describe the most common type of dependency in microarray studies [28]. We define "clumps" of correlated features as blocks of one relevant and four redundant features with a covariance matrix Σ*^{(k)}, where k is the number of the current block. The diagonal elements of Σ*^{(k) }for each block are equal to one and the offdiagonal elements are equal to ρ = 0.8. In total, we design five blocks of correlated features and therefore the covariance matrix has the form
where
Due to clumping blocks, the vector of β has a more complex form
with
where r denotes the number of relevant features. Using correlated blocks we investigate the ability of selecting correlated features, the so called grouping effect.
Optimal tuning parameters are found by an interval search in tuning parameter space using five fold cross validation. We select a large tuning parameter interval to be certain not to stick in local optima. The tuning parameter space for L_{1 }and SCAD SVM is onedimensional with λ_{1 }∈ [λ_{1,min}, λ_{1,max}]. Elastic SCAD has two tuning parameters λ_{1}, λ_{2 }∈ [λ_{l, min}, λ_{l, max}], l = 1, 2. Elastic Net applies LARS paths. for fixed λ_{2 }a λ_{1 }path is calculated and the optimal λ_{1 }is identified (for details refer to [17]). Thus, the optimal pair of parameter for Elastic Net was found in the twodimensional space ℝ × [λ_{l, min}, λ_{l, max}] We set the search interval for both parameters to [λ_{l, min}, λ_{l, max}] = [2^{10}, 2^{10}], l = 1, 2.
The performance of classifiers is characterised by the Youden index. The Youden index describes as equally weighted sum of true positive results ("sensitivity") and false positive results ("1  specificity"):
The maximal Youden index is one, when the true positive rate is one and the false positive rate is zero. For a random classifier the expected Youden index is zero. The sensitivity and specificity have equal weights in this index. Most often the costs and consequences of true positives and false positives will differ greatly. Therefore, Gu and Pepe [29] recommend reporting the two measures separately. For our simulated data, we consider the Youden index to be an appropriate indicator for feature selection methods performance, since we weight errors equally.
It is worth to mention, that for discrete classier the Youden index and the area under the curve (AUC) provide the same message due to their linear relationship. According to Greiner et al. [30], if there is only one point in the ROC plot, the ROC curve is estimated by connecting the three points. the point corresponds to the classifier, the (0, 0) and (1, 1) edges of the plot. Then geometrically, the estimated AUC corresponds to the average of estimated sensitivity and specificity. Thus, the Youden index and the AUC have a linear relationship. AUC = (sensitivity + specificity)/2 = (Youden index +1)/2. Optimizing the AUC will lead to the same results as optimizing the Youden index when dealing with discrete classifiers. Nevertheless, for real data application, the AUC values are presented in a separate column due to higher level of familiarity in bioinformatics.
Finally, the misclassification rate, size of the classifiers and frequencies of the selected features within 100 simulation runs are computed.
Simulation results
The performance of the feature selection methods applied to simulated data using p = 1000 features and n = 500 samples for training and testing is presented in the next section. The percentage of relevant features varies between 1% and 20% in four steps, i.e. r = 10, 50, 100, 200. We assume to have correlated blocks of features as described in the design section. The optimal tuning parameters were chosen as described above. Multiple comparisons in performance measures between the proposed prediction methods and the best method (the MCB test) for each simulation step will be done according to Hsu [31] based on 100 simulation runs. We used a noninferiority margin of a procedure to distinguish methods with similar performance.
Misclassification rate
Table 1 summarises the average misclassification rates depending on the number of relevant features. The numbers in parentheses are the standard errors of the estimates. For very sparse models (10 out of 1000 features are relevant) SCAD showed the lowest misclassification error (18%), followed by Elastic Net and Elastic SCAD (19.4% and 20.8% respectively), where both lie in indifference zone for best methods if the noninferiority margin was set to Δ = 0.05. For less sparse to nonsparse models (r = 50 and r = 100) Elastic Net showed the best performance. For r = 200 relevant features L_{1 }and Elastic Net showed nearly the same results (32.9% and 33.1% respectively). The same was observed for SCAD (34.7%) and Elastic SCAD (34.2%). For r ≥ 50 the misclassification rate was indistinguishable for all feature selection methods with exception of the L_{1 }SVM. The L_{2 }SVM classifiers showed larger misclassification errors for sparse models (r = 10 and r = 50) than all other feature selection methods. For less sparse models differences in misclassification error levelled out.
Table 1. Mean misclassification rate of feature selection methods applied to simulated test data
Youden index
The average Youden index for very sparse models (r = 10) was considerably high for all feature selection methods: 0.96 for SCAD, 0.95 for Elastic Net, 0.92 for Elastic SCAD, and 0.81 for L_{1 }SVM (Table 2). By increasing number of informative features, the Elastic Net SVM showed the best Youden index (0.71%  0.27%) among all feature selection methods, closely followed by the Elastic SCAD SVM (0.67%  0.27%), both being indistinguishable.
Table 2. Average Youden index for classifiers applied to simulated test data
All methods except the L1 SVM provided significantly comparable Youden indexes at the level α = 0.05 and a relevant difference Δ = 0.10 for r = 10. By increasing model complexity, the Elastic Net SVM showed the best Youden Index among all feature selection methods, closely followed by the Elastic SCAD SVM. Starting from r > 100 the is no significant difference between Elastic Net and Elastic SCAD SVMs. With increasing number of relevant features, the Youden index decreases from 0.9 to 0.27 for 'elastic' methods to 0.14 for the L_{1 }SVM and to 0.16 for the SCAD SVM. respectively.
Sparsity of the classifier
The SCAD SVM provided the most sparse classifier (in terms of selecting the smallest number of features) for r = 10 and r = 50 out of 1000 features (cf. Table 3). It selected 12 and 61 features, respectively. For less sparse models the Elastic Net and the Elastic SCAD SVMs had similar performance, selecting the smallest number of features.
Table 3. Median number of features selected
Selection Frequencies
A frequencies plot for the simulation study is represented in 'Additional file 1  Frequencies plot'. With increasing number of relevant features (r), a decrease of the proportion of true positives (in red) and an increase of the proportion of false positives (in blue) for all feature selection models was observed, respectively. At the same time we observed an increase of the false positives, which are correlated with the true positives (in green) in classifiers.
Additional file 1. Frequencies plot. Frequencies of selected features in the classifiers after 100 runs. In xaxis: features, yaxis: frequency of appearing of each features in classifiers after 100 runs. Features: true positives or nonzero (in red), zero features correlated with true positives (in green) and true negatives or zero (in blue). Algorithms from left to right: SCAD SVM, 1norm (L_{1}) SVM, Elastic Net SVM and Elastic SCAD SVM. Number of features: from top to bottom from very sparse till nonsparse models, r: 10, 50, 100, 200 out of 1000 features are relevant.
Format: PDF Size: 566KB Download file
This file can be viewed with: Adobe Acrobat Reader
The percentage of true positives in the classifiers is shown in Table S1 (Additional file 2  Tables S1, S2, S3). For r = 10 relevant features the Elastic Net SVM found almost all true positives (99.8%), followed by the Elastic SCAD SVM with 97.6%. For r = 50 the Elastic SCAD SVM achieved the sparsest solution followed by the L_{1 }SVM. In less sparse models, the L_{1 }SVM showed highest true positive rates of 84.5% and 86%.
Additional file 2. Tables S1, S2, S3. Table S1: Mean frequency percentages for nonzero features in the classifier. Mean frequency percentages for nonzero features in the classifier (true positives) after 100 runs. Standard deviations in parentheses. Table S2: Mean frequency percentages for zero features, high correlated with nonzero features in the classifier. Mean frequency percentages for zero features, high correlated with nonzero features in the classifier after 100 runs. Standard deviations in parentheses. Table S3: Mean frequency percentages for independent nonzero features in the classifier (false positives). Mean frequency percentages for independent nonzero features in the classifier (false positives) after 100 runs. Standard deviations in parentheses.
Format: PDF Size: 72KB Download file
This file can be viewed with: Adobe Acrobat Reader
Grouping effect
We further evaluated the ability of feature selection methods to select correlated features of true positives. Although for all scenarios L_{1 }SVM has found the largest percentage of correlated features, which increases with increasing number of relevant features (23.6  62.5%), the level of correlated features is comparable to the level of nonrelevant features (Table S2).
Comparing Tables S1, S2 and S3 one can observe that the SCAD and the L_{1 }SVMs failed to find features highly correlated with true positives more often than with independent false positives. The Elastic Net and the Elastic SCAD SVMs managed to discover correlated features (in green) more often than the independent false positives (in blue), at least for sparse models (r = 10 and r = 50).
The false positive rate
For very sparse models, the false positive rate (FPR) was the smallest for the SCAD SVM, followed by the Elastic Net and the Elastic SCAD SVMs (Table S3). For other less sparse models the Elastic Net SVM selected fewer false positives than the remaining methods. The second best method is the Elastic SCAD SVM.
Conclusions
• As expected from theory the SCAD SVM and the L_{1 }SVM produced classifiers with low prediction error for very sparse situations.
• For less sparse and nonsparse models, the Elastic Net and the Elastic SCAD SVM showed better results than the L_{1 }and the L_{2 }SVMs with respect to accuracy, Youden index and sparsity of classifiers.
• The SCAD SVM and the L_{1 }SVM were not able to find correlated features. The Elastic Net and the Elastic SCAD SVMs found correlated features more frequently than one would expected under random selection. Although the grouping effect strength weakens with increasing number of relevant features, the Elastic Net and Elastic SCAD SVMs still managed the grouping effects.
• In general, the Elastic Net and the Elastic SCAD SVMs showed similar performance. Additionally, the Elastic SCAD SVM provided more sparse classifiers than the Elastic Net SVM.
Applications
NKI breast cancer data set
Two studies on breast cancer from the Netherlands Cancer Institute (NKI) were published by the van't Veer group [32], [33]. In the first paper, a set of 78 lymph node negative patients with preselected 4919 clones were used to find a predictor for distant metastases. The classifier was trained and validated on patients who developed distant metastases within five years after surgery and patients being metastasisfree for at least the first five years. The resulting predictor was a 70gene signature also known as MammaPrint(R). We will use the MammaPrint(R) signature as reference in the analysis of the NKI breast cancer data set. The signature was generated based on genewise correlations between the gene expression and metastasis occurrence. The data set was taken from http://www.rii.com/publications/2002/vantveer.html webcite.
In a subsequent validation study, data from 295 patients (which partially included patients from the first study) were used to validate the signature [33]. Among the patients, 151 were lymph node negative and 144 had lymph node positive disease. The preprocessed data containing 4919 clones is available from http://www.rii.com/publications/2002/nejm.html webcite.
After excluding patients being identical to the training set and 10 patients with no metastasis information, 253 patients remained. Among the 253 patients there are 114 lymph node negative and 139 lymph node positive patients.
In this paper, we combined the 78 lymph node negative sample set from the first publication with 114 lymph node negative patients from the validation study. In total, a data set with 192 lymph node negative samples was used. The estimation of classifier performance was computed by a tenfold stratified crossvalidation.
Results on NKI breast cancer data set
Table 4 shows the misclassification error, sensitivity, specificity, Youden index and AUC value of four feature selection methods, RFE SVM and standard L_{2 }SVM based on tenfold stratified cross validation.
Table 4. Summary of classifiers for the NKI data set with distant metastasis as endpoint
RFE SVM was used according to Guyon's approach [1], where at each iteration half of features with lowest ranks are eliminated. To increase the classifier's stability, RFE SVM with fivefold stratified cross validation was repeated five times. According to the average cross validation error the optimal number of features was 2^{8 }= 256. Optimal tuning parameters for penalized SVM methods were found by the interval search on the tuning parameter space as described in the method section using tenfold stratified cross validation.
The SCAD SVM reduced the number of features from 4919 to 476, L_{1 }SVM selected 1573 features, Elastic Net 109 features, and the Elastic SCAD had 459 features in the classifier. For the NKI data set the best predictor with respect to misclassification error was L_{1 }SVM. Elastic Net and Elastic SCAD SVMs provided similar results, followed by SCAD SVM, which was slightly worse.
The relationship between the true positive rate (TPR, sensitivity) and the false positive rate (FPR, 1specificity) for each classifier is depicted as a point in the ROC plot (Figure 1). Isolines with constant Youden index are plotted as dashed lines. Taking the Youden index as an additional criterion, one could prioritise L_{1 }SVM. RFE SVM and both 'elastic' methods lay clustered in the ROC plot with clear distance to the L_{1 }classifier. The L_{2 }was placed inbetween L_{1 }and this cluster, being not far from the cluster.
Figure 1. ROC plot for the NKI breast data set. The characteristics for the different feature selection methods were derived using tenfold stratified cross validation. TPR and FPR values are presented as points (x axis: 1 specificity = FPR, y axis. sensitivity = TPR). RFE_256 is RFE SVM with 256 top ranked features, ENet is Elastic Net SVM, ESCAD is Elastic SCAD SVM. '70_sign' stands for the 70gene signature classifier. Gray dashed lines depict isolines of the Youden index.
Interestingly, the MammaPrint(R) signature ("70_sign") neither showed good test accuracy nor a reliable sensitivity or specificity. L_{2 }SVM and the feature selection methods outperformed the published signature.
Conclusions
For the two data sets from van't Veer group feature selection methods produced signatures with similar prediction accuracy, but being different in size. L_{1 }SVM with a nonsparse classifier provided the best sensitivity and specificity, followed by more sparse predictors from Elastic Net SVM and Elastic SCAD SVM.
MAQCII breast cancer data set
This data set is part of the MicroArray Quality Control (MAQC)II project, which has been designed to investigate numerous data analysis methods and to reach consensus on the "best practices" for development and validation of microarraybased classifiers for clinical and preclinical applications. One biological endpoint is estrogen receptor (ER) status. Out of 230 patients in total, 89 patients have negative ER status and 141 patients positive ER status. A clinical endpoint is pathological complete response (pCR) to preoperative chemotherapy. Among the 230 patients in the data set, 182 patients showed no pCR and 48 had a pCR.
The preprocessed data contains 22283 features and is available from GEO database, accession number GSE20194.
Results on MAQCII breast cancer data set
The feature selection methods SCAD SVM, L_{1 }SVM, Elastic Net SVM and Elastic SCAD SVM with internal tenfold stratified cross validation were applied to build classifiers. Additionally, the L_{2 }SVM and the RFE SVM were used as reference models. To achieve performance measurements tenfold stratified cross validation was used.
pCR prediction
Based on the minimal average misclassification error, the optimal number of features of RFE SVM classifier was obtained to be 2^{11 }= 2048 (Table 5). The penalized SVM methods provided moderately sparse models, Elastic SCAD SVM with 148 features, Elastic Net SVM with 398 features and dense models, L_{1}, SCAD and RFE SVMs with more than 1000 features.
Table 5. Summary of classifiers for the MAQCII data set with pCR status as endpoint
The misclassification error rate was similar for all methods with the Elastic SCAD classifier showing the lowest error rate of 15%. With nearly equally high specificity (9194%), we observed large variations in sensitivity of different feature selection methods as shown in the corresponding ROC plot (Figure 2). The Elastic SCAD outperformed all methods with sensitivity of 52%. Interestingly, the Elastic Net showed the smallest sensitivity of 15% resulting in a small Youden index of 0.06.
Figure 2. ROC plot for MAQCII breast data set with pCR as endpoint. The characteristics for the different feature selection methods were derived using tenfold statrifierd cross validation. TPR and FPR values are presented as points (x axis: 1 specificity = FPR, y axis. sensitivity = TPR). RFE_256 is RFE SVM with 1024 top ranked features, ENet is Elastic Net SVM, ESCAD is Elastic SCAD SVM. Gray dashed lines depict isolines of the Youden index.
Overall, Elastic SCAD showed better classification characteristics than other methods. Moreover, the higher specificity of the Elastic SCAD classifier is of clinical importance. The patients that did not respond to the therapy were recognized with higher probability.
ER status
We also used the MAQCII data set to predict the ER status. Here, the L_{1 }SVM failed to derive a sparse solution, whereas SCAD, Elastic Net and Elastic SCAD SVM classifiers were similar (Table 6). Moreover, Elastic SCAD showed the smallest error rate and highest sensitivity over all methods.
Table 6. Summary of classifiers for the MAQCII data set with ER status as endpoint
All classification methods provided small misclassification errors, high sensitivity and high specificity. The ROC plot in Figure 3 demonstrates this performance of predictors. As presented in Table 6 the Elastic Net, SCAD and Elastic SCAD SVMs selected small numbers of features, 3, 32 and 59 out of 22283, respectively. The extreme sparseness of the Elastic Net SVM was paid by lower sensitivity and specificity compared to other methods. The misclassification test error was similar for all methods (714%). The Elastic SCAD SVM classifier showed the smallest error rate of 7%.
Figure 3. ROC plot for MAQCII breast data set with ER as endpoint. The characteristics for the different feature selection methods were derived using tenfold stratified cross validation. TPR and FPR values are presented as points (x axis: 1 specificity = FPR, y axis. sensitivity = TPR). RFE_256 is RFE SVM with 1024 top ranked features, ENet is Elastic Net SVM, ESCAD is Elastic SCAD SVM. Gray dashed lines depict isolines of the Youden index.
For this classification task, the sparse classifier Elastic SCAD and SCAD showed the best characteristics.
Screening on two additional breast cancer data sets
These data sets were recently analysed and published by Johannes et. al. [34]. The first data set, the Mainz cohort, contains of 154 lymph nodenegative, relapse free patients and 46 lymph nodenegative patients that suffered a relapse (GEO acession number GSE11121). The relapse is defined as appearance of distant metastasis within five years after the treatment. The second data set, the Rotterdam cohort, represents 286 lymph nodenegative breast cancer samples including 107 relapse events (GSE2034). Both data sets were generated using the Affymetrix HGU133A platform, normalized with the same methods and relapse as the primary classification endpoint. We trained the feature selection classifiers on the whole cohort, Mainz data or Rotterdam data, and used the other cohort as an independent validation data set, respectively as presented in Tables 7 and 8.
Table 7. Summary of classifiers for Mainz cohort, validated on Rotterdam cohort with relapse as endpoint
Table 8. Summary of classifiers for Rotterdam cohort, validated on Mainz cohort with relapse as endpoint
We can see that all feature selection methods had lower misclassification test error than the L_{2 }SVM containing all features for breast cancer data sets. The classifiers perform different for each data set. The Elastic Net SVM had small error rate for the Rotterdam cohort, but failed to classify the Mainz samples adequately. The L_{2 }SVM classifier including all features had the second best Youden index for the Mainz set, however for Rotterdam data showed the worst Youden index. Using both, the test error and AUC value as a combined measure of sensitivity and the specificity, one would conclude that the L_{1}, SCAD and Elastic SCAD SVMs provide reasonable and robust solutions with respect to the combined analysis of the two breast cancer data sets.
Altogether, Elastic SCAD seems to provide an overall acceptable compromise for sparse and nonsparse data.
Conclusions
In highdimensional prediction tasks, feature selection plays an important role. In this paper, we proposed a novel feature selection method for SVM classification using a combination of two penalties, SCAD and L_{2}. The commonly used penalty functions L_{1}, SCAD and Elastic Net were investigated in parallel with the new method on simulated and public data. To address the problem of finding optimal tuning parameters for SVM classification the efficient parameter search algorithm from Froehlich and Zell [15] was implemented.
In almost all cases, the four feature selection classifies outperformed ordinary Support Vector Classification using the L_{2 }penalty. From the simulation study we concluded that for sufficiently large sample sizes, feature selection methods with combined penalties are more robust to changes of the model complexity than using single penalties alone.
The SCAD SVM followed by the L_{1 }SVM, as expected, showed very good performance in terms of prediction accuracy for very sparse models, but failed for less sparse models. Combined penalty functions in combination with the SVM algorithm, Elastic Net and Elastic SCAD, performed well for sparse and less sparse models.
Comparisons with commonly used penalty functions in the simulation study illustrated that the Elastic SCAD and the Elastic Net SVMs showed similar performance with respect to prediction accuracy. Both 'elastic' methods were able to consider correlation structures in the input data (grouping effect). However, the Elastic SCAD SVM in general provides more sparse classifiers than the Elastic Net SVM.
Finally, applied to publicly available breast cancer data sets, the Elastic SCAD SVM performed very flexible and robust in sparse and nonsparse situations. Results from the simulation study and real data application render Elastic SCAD SVM with automatic feature selection a promising classification method for highdimensional applications.
Authors' contributions
NB and AB contributed to the design of the simulation study and to theoretical investigations of problems. NB performed simulations, analyses and wrote the manuscript. GT and AB participated in the preparation of the manuscript. AB and PL supervised the work. All authors read and approved the manuscript.
Acknowledgements
We would like to thank the editor and the anonymous reviewers for their constructive comments, which significantly improved the quality of this manuscript.
This work was supported by Grant 01GS0883 from the German Federal Ministry of Education and Research (BMBF) within the National Genome Research Network NGFNplus.
References

Guyon I, Weston J, Barnhill S, Vapnik V: Gene Selection for Cancer Classification using Support Vector Machines.

Vapnik V: The Nature of Statistical Learning Theory. New York: Springer; 1995.

Hastie T, Tibshirani R, Friedman J: The elements of statistical learning: data mining inference and prediction. New York: Springer; 2001.

Inza I, Sierra B, Blanco R, Larranaga P: Gene selection by sequential search wrapper approaches in microarray cancer class prediction.

Markowetz F, Spang R: Molecular diagnosis: classification, model selection and performance evaluation.
Methods Inf Med 2005, 44(3):438443. PubMed Abstract  Publisher Full Text

Guyon I, Elisseff A: An Introduction to Variable and Feature Selection.
Journal of Machine Learning Research 2003, 3:11571182. Publisher Full Text

Tibshirani R: Regression shrinkage and selection via the lasso.
Journal of the Royal Statistical Society, Series B 1996, 58:267288.

Zou H, Hastie T: Regularization and variable selection via the elastic net.
Journal of the Royal Statistical Society, Series B 2005, 67(2):301320. Publisher Full Text

Zhang HH, Ahn J, Lin X, Park C: Gene selection using support vector machines with nonconvex penalty.
Bioinformatics 2006, 22(1):8895. PubMed Abstract  Publisher Full Text

Allison DB, Cui X, Page GP, Sabripour M: Microarray data analysis: from disarray to consolidation and consensus.
Nature Reviews Genetics 2006, 7:5565. PubMed Abstract  Publisher Full Text

Hoheisel JD: Microarray technology: beyond transcript profiling and genotype analysis.
Nature Reviews Genetics 2006, 7:200210. PubMed Abstract  Publisher Full Text

Quackenbush J: Computational analysis of microarray data.
Nature Review Genetics 2001, 2:418427. Publisher Full Text

Li X, Xu R: High Dimensional Data Analysis in Oncology. New York: Springer; 2008.

Froehlich H, Zell A: Efficient parameter selection for support vector machines in classification and regression via modelbased global optimization.

Liu Q, Sung AH, Chen Z, Liu J, Huang X, Deng Y: Feature selection and classification of MAQCII breast cancer and multiple myeloma microarray gene expression data.
PLoS ONE 2009, 4:e8250. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Wang L, Zhu J, Zou H: Hybrid huberized support vector machines for microarray classification and gene selection.
Bioinformatics 2008, 24(3):412419. PubMed Abstract  Publisher Full Text

Bradley PS, Mangasarian OL: Feature selection via concave minimization and support vector machines.
Machine Learning Proceedings of the Fifteenth International Conference 1998, 8290.

Fung G, Mangasarian OL: A feature selection newton method for support vector machine classification.
Computational Optimization and Applications Journal 2004, 28(2):185202.

Comments onWavelets in Statistics. A Review by Antoniadis
1997.

Fan J, Li R: Variable selection via nonconcave penalized likelihood and its oracle properties.
Journal of American Statistical Association 2001, 96:13481360. Publisher Full Text

Wang L, Zhu J, Zou H: The double regularized support vector machine.

Jones D, Schonlau M, Welch W: Efficient global optimization of expensive blackbox functions.
Journal of Global Optimization 1998, 13:455492. Publisher Full Text

Kohavi R: A study of crossvalidation and bootstrap for accuracy estimation and model selection.
Proceedings of the Fourteenth International Joint Conference on Artificial Intelligence 1995, 2:11371143.

Becker N, Werft W, Toedt G, Lichter P, Benner A: penalizedSVM: a Rpackage for feature selection SVM classification.
Bioinformatics 2009, 25(13):17111712. PubMed Abstract  Publisher Full Text

Dimitriadou E, Hornik K, Leisch F, Meyer D, Weingessel A: [http://CRAN.Rproject.org/package=e1071] webcite
e1071: Misc Functions of the Department of Statistics (e1071), TU Wien. 2009.
[R package version 1.522]

Chang C, Lin C: LIBSVM: a Library for Support Vector Machines. [http://www.csie.ntu.edu.tw/~cjlin/libsvm] webcite
2001.

Storey JD, Tibshirani R: Estimating false discovery rates under dependence, with applications to DNA microarrays. Tech rep., Stanford University: Stanford Technical Report; 2001.

Gu W, Pepe M: Measures to Summarize and Compare the Predictive Capacity of Markers. [http://www.bepress.com/uwbiostat/paper342] webcite
Working paper 342., UW Biostatistics Working Paper Series 2009.

Greiner M, Pfeiffer D, Smith RD: Principles and practical application of the receiveroperating characteristic analysis for diagnostic tests. [http:/ / www.sciencedirect.com/ science/ article/ B6TBK408BJCN3/ 2/ 3a4753dc80ec448666ef990ee4c33078] webcite
Preventive Veterinary Medicine 2000, 45(12):2341. PubMed Abstract  Publisher Full Text

Hsu JC: Multiple Comparisons Theory and Methods. Chapman & Hall; 1996.

van't Veer LJ, Dai H, van de Vijver MJ, He YD, Hart AA, Mao HL, Mand Peterse, van der Kooy K, Marton MJ, Witteveen AT, Schreiber GJ, Kerkhoven RM, Roberts PS, Cand Linsley, Bernards R, H FS: Gene expression profiling predicts clinical outcome of breast cancer.
Nature 2002, 415:530536. PubMed Abstract  Publisher Full Text

van de Vijver MJ, He YD, van't Veer LJ, Dai H, Hart AA, Voskuil DW, Schreiber GJ, Peterse HL, Roberts C, Marton MJ, Parrish M, Atsma D, Witteveen A, Glas A, Delahaye L, van der Velde T, Bartelink H, Rodenhuis S, Rutgers ET, Friend SH, Bernards R: A geneexpression signature as a predictor of survival in breast cancer.
N Engl J Med 2002, 347:19992009. PubMed Abstract  Publisher Full Text

Johannes M, Brase JC, Froehlich H, Gade S, Gehrmann M, Faelth M, Sueltmann H, Beissbarth T: Integration Of Pathway Knowledge Into A Reweighted Recursive Feature Elimination Approach For Risk Stratification Of Cancer Patients.
Bioinformatics 2010, 26(17):21362144. PubMed Abstract  Publisher Full Text