Data Mining Algorithms In R/Classification/penalizedSVM

Introduction
Classifiers are some of the most common data analysis tools. There are a lot of implemented techniques, but we may point SVM (Support Vector Machine) as one of the most powerful, especially in high-dimension data. The most well known SVM algorithm was created by Vladimir Vapnik.

The standard SVM implementation SVM takes a input dataset and, for each given input, predicts which of two possible classes the input set belongs to. That's most common use the algorithm to predict if the input belongs to certain dichotomy, or not. Because of this characteristic, SVM is a called a non-probabilistic binary linear classifier.

On Machine Learning-based algorithms such as SVM, the input data has to be separated on two sets: a training set and a test set. The difference between the training and the test set is that, on the training the examples' classes are known beforehand. The test set contains the examples that should have their classes predicted. Given a set of training examples, an SVM algorithm builds a model that predicts what are the categories of the test set's examples.

Representing the examples in a d-dimension space, the model built by the algorithm is a hyper-plane that separates the examples belonging to the dichotomy from the ones that aren't. New examples are then mapped into that space so we can predict the categories they belong, based on the side of the gap it falls on.

The technique described here is a variation of the standard SVM using penalty functions. The technique is implemented on the R-package called penalized SVM, that has smoothly clipped absolute deviation (SCAD), 'L1-norm', 'Elastic Net' ('L1-norm' and'L2-norm') and 'Elastic SCAD' (SCAD and 'L2-norm') as available penalties.

Algorithm

 * $$\begin{align}&\quad \; \text{Given a training dataset }(x_1,y_1),...,(x_n,y_n),\text{ where }x_i\text{ is a }d\text{-tuple of the }d\text{ input}\\&\text{parameters }i\text{, and }y_i\in -1, 1\text{ where }y_i=1\text{ means }i\text{ belongs to the dichotomy, and }y_i=-1,\\&\text{the opposite. The SVM divides the space by a linear boundary:}\\\\&\qquad\qquad\qquad\qquad\qquad\qquad\qquad\ \ f(x)=\begin{matrix}\sum_{j=1}^d w_jh_j+b\end{matrix}\\ \\&\quad\ \;\text{where }w=(w_1,...,w_d)\text{ are the coefficients of the hyper-plan and }b\text{ denotes the in-}\\&\text{tercept of the hyperplane. The output for an example from the test set }x_{test}\text{ would be}\\ &y_{test}=sign[f(x_{test)}].\text{ So, the example belongs to the dichotomy if }f(x_{test})>0.\end{align}$$

Finding the best model to classify new examples is, as a matter of fact, similar to the problem of finding the optimal hyperplane. The quality of a model built by the algorithm is measured by the margin of its hyperplane. The margin is the distance between the hyperplanes wx+b = 1 and wx+b = -1,being so that the closest points to wx+b = 0, in each side of that hyperplane are crossed by the other two hyperplanes. The best hyperplanes are those with the biggest margins.

The problem of finding the optimal hyperplane with maximal margin is solved by convex optimization. Maximizing the margin can be achieved by solving:

$$\min_{b,w} \sum \left[1-y_if\left(x_i\right)\right]_++\operatorname{pen}_\lambda\left(w\right)$$

The term

$$pen_\lambda(w)=\lambda\begin{Vmatrix} w \end{Vmatrix}^2_2$$

for SVM implementation has the form of L2 norm ('ridge penalty'). This penalty causes the reduction of the coefficients, but not always to values greater than zero.

The quality of a model can also be measured by the error computed for the predictions on the training set. A low prediction is required, but reducing it to zero may cause a problem know as over fitting. This means that the model built classifies very well the examples on training set, but is not appropriate to classify the ones on the test set. The predictions made in this case, would usually be wrong for the new examples.

Not only the low prediction defines a model's quality. It is possible to improve it further by identifying covariates that plays important roles on discrimination and assessing their contribution to the classifier. This can be achieved by applying a feature selection method.

Feature selection methods are divided into two classes: filter and wrapper methods. Filter methods drop irrelevant features before the model is built by the algorithm. Wrapper methods increases the prediction power by providing the selection within the optimization procedure.

The R package 'penalizedSVM' provides some feature selection methods. One of them is the wrapper feature selection SCAD (Smoothly Clipped Absolute Deviation). SCAD is a non-convex penalty function first proposed by Fan and discussed in .On SVM is combined with SCAD for feature selection. The penalization term for SCAD SVM has the form:


 * $$pen_\lambda(w) = \textstyle \sum_{j=1}^d p_\lambda(w_j)$$

where the SCAD penalty function for each coefficient wj is defined as:


 * $$p_\lambda(w_j)=\begin{cases} \lambda\left\vert w_j\right\vert & \text{if }\left\vert w_j\right\vert \le \lambda \\ -\frac{\left\vert w_j\right\vert^2-2a\alpha\left\vert w_j\right\vert+\lambda^2}{2(a-1)} & \text{if }\left\vert w_j\right\vert \le a\lambda \\ \frac{a+1/\lambda^2}{2} & \text{if }\left\vert w_j\right\vert > a\lambda \end{cases}$$

with tuning parameters a > 2 (in the package, a = 3.7) and λ > 0. pλ (w) corresponds to a quadratic spline function with knots at λ and aλ.

For small coefficients, the SCAD has the same behavior as the L1. For large coefficients, however, the SCAD applies a constant penalty, in contrast to the L1 penalty, which increases linearly as the coefficient increases. This absolute maximum of the SCAD penalty, which is independent from the input data, decreases the possible bias for estimating large coefficients.

Implementation
The package described in this section is penalizedSVM. This package provides feature selection SVM using penalty functions. The smoothly clipped absolute deviation (SCAD), 'L1-norm', 'Elastic Net' ('L1-norm' and 'L2-norm') and 'Elastic SCAD'(SCAD and 'L2-norm') penalties are available. The tuning parameters can be found using either a fixed grid or a interval search.

This package has several dependencies. The packages that need to be installed and theirs descriptions are listed below:

vector machines, shortest path computation, bagged clustering, naive Bayes classifier, ... edition). separate shrinkage for variances and correlations. The details of the method are explained in Sch\"afer and Strimmer (2005) and Opgen-Rhein and Strimmer (2007). analysis, mixed linear models, heteroscedastic regression, Tweedie family generalized linear models, the inverse-Gaussian distribution and Gauss quadrature. with jumps to the limiting linear model (LLM). diagnostic plots. Contact the maintainer for a package version that implements sensitivity analysis functionality.
 * e1071 Functions for latent class analysis, short time Fourier transform, fuzzy clustering, support
 * MASSFunctions and datasets to support Venables and Ripley, 'Modern Applied Statistics with S' (4th
 * corpcorThis package implements a James-Stein-type shrinkage estimator for the covariance matrix, with
 * stadmod Various statistical modeling functions including growth curve comparisons, limiting dilution
 * tgp Bayesian nonstationary, semiparametric nonlinear regression and design by treed Gaussian processes
 * mlepg Maximum likelihood Gaussian process modeling for univariate and multi-dimensional outputs with
 * lhs This package provides a number of methods for creating and augmenting Latin Hypercube Samples

The major function on penalizedSVM is svm.fs. Its use can be described as follow:

svm.fs( x, y, fs.method = c("scad", "1norm", "scad+L2", "DrHSVM"), grid.search=c("interval","discrete"), lambda1.set=NULL, lambda2.set=NULL, bounds=NULL, parms.coding= c("log2","none"), maxevals=500, inner.val.method = c("cv", "gacv"), cross.inner= 5, show= c("none", "final"), calc.class.weights=FALSE, class.weights=NULL, seed=123, maxIter=700, verbose=TRUE, ...)
 * 1) Default S3 method:
 * 1) tuning parameter settings
 * 2) chose the search method for tuning lambda1,2: 'interval' or 'discrete'
 * 1) fixed grid for lambda1, lambda2
 * 1) define range for lambda1,2 for interval search
 * 1) parms.coding="none" or "log2"
 * 1) internal parameter for DIRECT
 * 1) valuidation settings
 * 2) fot nested validation, 'cross.outer'-fold cv
 * 3) cross.outer= 0,
 * 4) method for the inner validation: cross validation, gacv
 * 1) 'cross.inner'-fold cv
 * 1) show plots in Direct?
 * 1) other settings
 * 2) internal parameter for svm
 * 1) seed
 * 1) max iterations for the feature selection svm method
 * 1) verbose?

where the arguments are:

sizes. Not all factor levels have to be supplied (default weight: 1). All components have to be named.
 * x: input matrix with genes in columns and samples in rows!
 * y: numerical vector of class labels, -1, 1
 * fs.method: feature selection method. Available ’scad’, ’1norm’ for 1-norm, "DrHSVM" for            Elastic Net and "scad+L2" for Elastic SCAD
 * grid.search: chose the search method for tuning lambda1,2: ’interval’ or ’discrete’, default: ’interval’
 * lambda1.set: for fixed grid search: fixed grid for lambda1, default: NULL
 * lambda2.set: for fixed grid search: fixed grid for lambda2, default: NULL
 * bounds: for interval grid search: fixed grid for lambda2, default: NULL
 * parms.coding: for interval grid search: parms.coding: none or log2, default: log2
 * maxevals: the maximum number of DIRECT function evaluations, default: 500.
 * cross.outer: fold of outer cross validation, default is 0, no cv.
 * calc.class.weights: calculate class.weights for SVM, default: FALSE
 * class.weights: a named vector of weights for the different classes, used for asymmetric class
 * inner.val.method: method for the inner validation: cross validation, gacv, default cv
 * cross.inner: ’cross.inner’-fold cv, default: 5
 * show: for interval search: show plots of DIRECT algorithm: none, final iteration, all iterations. Default: none
 * seed: seed
 * maxIter: maximal iteration, default: 700
 * verbose: verbose?, default: TRUE
 * ... additional: argument(s)

View
To visualize the result of the algorithm, you can use the function show, as exampled above.

Case Study
In this section, we illustrate a case study with penalizedSVM:

Scenario
The objective is to tell if a picture belongs to the same category of another one selected previously. For that purpose, we use the difference between the values of those characteristics. The attributes for the data used as input represents the difference from the value of the characteristic in the picture selected and the ones we want to classify.

Dataset
The training set and the test set were generated using the commands shown below: > train<-sim.data(n = 200, ng = 100, nsg = 10, corr=FALSE, seed=seed ) > print(str(train)) List of 3 $ x  : num [1:100, 1:200] -0.5605 2.1988 -0.0736 1.074 0.3563 ... ..- attr(*, "dimnames")=List of 2 .. ..$ : chr [1:100] "pos1" "pos2" "pos3" "pos4" ... .. ..$ : chr [1:200] "1" "2" "3" "4" ... $ y   : Named num [1:200] 1 -1 1 1 1 -1 -1 -1 1 -1 ...  ..- attr(*, "names")= chr [1:200] "1" "2" "3" "4" ... $ seed: num 123 NULL > test<-sim.data(n =20, ng = 100, nsg = 10, corr=FALSE, seed=seed+1 ) > print(str(test)) List of 3 $ x  : num [1:100, 1:20] -1.3851 -1.1036 -0.2677 0.2836 -0.0951 ... ..- attr(*, "dimnames")=List of 2 .. ..$ : chr [1:100] "pos1" "pos2" "pos3" "pos4" ... .. ..$ : chr [1:20] "1" "2" "3" "4" ... $ y   : Named num [1:20] -1 1 -1 1 1 1 1 1 1 -1 ...  ..- attr(*, "names")= chr [1:20] "1" "2" "3" "4" ... $ seed: num 124 NULL

Execution
To build the model, the following command were used:

> bounds=t(data.frame(log2lambda1=c(-10, 10))) >colnames(bounds)<-c("lower", "upper") >system.time( scad<- svm.fs(t(train$x)[,1:100], y=train$y, fs.method="scad", bounds=bounds, + cross.outer= 0, grid.search = "interval", maxIter = 10, + inner.val.method = "cv", cross.inner= 5, maxevals=500, + seed=seed, parms.coding = "log2", show="none", verbose=FALSE ) )
 * 1) computation intensive; for demostration reasons only for the first 100 features
 * 2) and only for 10 Iterations maxIter=10, default maxIter=700

Output
To see the model created:

> print(str(scad$model)) List of 11 $ w      : Named num [1:23] 0.625 0.616 0.353 0.258 0.959 ...  ..- attr(*, "names")= chr [1:23] "pos1" "pos2" "pos3" "pos4" ... $ b      : num -0.115 $ xind   : int [1:23] 1 2 3 4 5 6 7 8 9 10 ... $ index   : int [1:83] 3 4 9 14 17 18 22 35 37 40 ... $ fitted  : num [1:200] 2.6 1.24 0.65 1 1.15 ... $ type   : num 0 $ lambda1 : num 0.126 $ lambda2 : NULL $ iter   : num 10 $ q.val  : num 0.195 $ fit.info:List of 13 ..$ fmin      : num 0.195 ..$ xmin      : Named num -2.99 .. ..- attr(*, "names")= chr "log2lambda1" ..$ iter      : num 26 ..$ neval     : num 46 ..$ maxevals  : num 500 ..$ seed      : num 123 ..$ bounds    : num [1, 1:2] -10 10 .. ..- attr(*, "dimnames")=List of 2 .. .. ..$ : chr "log2lambda1" .. .. ..$ : chr [1:2] "lower" "upper" ..$ Q.func    : chr ".calc.scad" ..$ points.fmin:'data.frame':	1 obs. of 2 variables: .. ..$ log2lambda1: num -2.99 .. ..$ f         : num 0.195 ..$ Xtrain    : num [1:46, 1] -7.52 -2.26 -1.34 9.99 9.03 ... .. ..- attr(*, "dimnames")=List of 2 .. .. ..$ : NULL .. .. ..$ : chr "log2lambda1" ..$ Ytrain    : num [1:46] 3.65e-01 3.20e-01 4.60e-01 1.00e+16 1.00e+16 ...  ..$ gp.seed    : num [1:25] 123 124 125 126 127 128 129 130 131 132 ... ..$ model.list :List of 1 .. ..$ model:List of 10 .. .. ..$ w              : Named num [1:23] 0.625 0.616 0.353 0.258 0.959 ...  .. .. .. ..- attr(*, "names")= chr [1:23] "pos1" "pos2" "pos3" "pos4" ... .. .. ..$ b              : num -0.115 .. .. ..$ xind           : int [1:23] 1 2 3 4 5 6 7 8 9 10 ...  .. .. ..$ index          : int [1:83] 3 4 9 14 17 18 22 35 37 40 ...  .. .. ..$ fitted         : num [1:200] 2.6 1.24 0.65 1 1.15 ...  .. .. ..$ type           : num 0 .. .. ..$ lambda1        : num 0.126 .. .. ..$ iter           : num 10 .. .. ..$ q.val          : num 0.195 .. .. ..$ inner.val.method: chr "cv" NULL To predict a class for the examples on the test set:

>(scad.5cv.test<-predict.penSVM(scad, t(test$x)[,1:100], newdata.labels=test$y) ) $pred.class [1] -1 1 -1 -1 -1 1  1  1  1  -1 -1 -1 1  1  1  -1 1  1  -1 -1 Levels: -1 1

$fitted [,1] 1 -2.5344366 2   2.3440943 3  -1.3972349 4  -0.3613470 5  -2.1187284 6   1.1287477 7   2.5584662 8   1.9155333 9   1.5543941 10 -0.7128084 11 -1.6944994 12 -0.2943272 13  1.8497781 14  2.7800572 15  0.8927699 16 -0.1289518 17  2.4560094 18  0.8756835 19 -2.2114729 20 -1.7342811

$tab newdata.labels pred.class -1 1 -1 6 4        1   1 9

$error [1] 0.25

$sensitivity [1] 0.6923077

$specificity [1] 0.8571429

> test<-sim.data(n = 20, ng = 100, nsg = 10, corr=FALSE, seed=seed+1 ) > print(str(test)) List of 3 $ x  : num [1:100, 1:20] -1.3851 -1.1036 -0.2677 0.2836 -0.0951 ... ..- attr(*, "dimnames")=List of 2 .. ..$ : chr [1:100] "pos1" "pos2" "pos3" "pos4" ... .. ..$ : chr [1:20] "1" "2" "3" "4" ... $ y   : Named num [1:20] -1 1 -1 1 1 1 1 1 1 -1 ...  ..- attr(*, "names")= chr [1:20] "1" "2" "3" "4" ... $ seed: num 124 NULL > (scad.5cv.test<-predict.penSVM(scad, t(test$x)[,1:100], newdata.labels=test$y)) $pred.class [1] -1 1 -1 -1 -1 1  1  1  1  -1 -1 -1 1  1  1  -1 1  1  -1 -1 Levels: -1 1

$fitted [,1] 1 -2.5344366 2   2.3440943 3  -1.3972349 4  -0.3613470 5  -2.1187284 6   1.1287477 7   2.5584662 8   1.9155333 9   1.5543941 10 -0.7128084 11 -1.6944994 12 -0.2943272 13  1.8497781 14  2.7800572 15  0.8927699 16 -0.1289518 17  2.4560094 18  0.8756835 19 -2.2114729 20 -1.7342811

$tab newdata.labels pred.class -1 1 -1 6 4        1   1 9

$error [1] 0.25

$sensitivity [1] 0.6923077

$specificity [1] 0.8571429

Analysis
To analyse the results, the follow commands can be used:

> print(paste("minimal 5-fold cv error:", scad$model$fit.info$fmin, + "by log2(lambda1)=", scad$model$fit.info$xmin)) [1] "minimal 5-fold cv error: 0.195 by log2(lambda1)= -2.99093721912059" > print(" all lambdas with the same minimum? ") [1] " all lambdas with the same minimum? " > print(scad$model$fit.info$ points.fmin) log2lambda1    f 36   -2.990937 0.195 > print(paste(scad$model$fit.info$neval, "visited points")) [1] "46 visited points" >.plot.EPSGO.parms (scad$model$fit.info$Xtrain, scad$model$fit.info$Ytrain, + bound=bounds, Ytrain.exclude=10^16, plot.name=NULL )
 * 1) Plot the results



= References =