Support Vector Machines with R
Anirban Ghatak
1 Support vector Machine
The advent of computers brought on rapid advances in the field of statistical classification, one of which is the Support Vector Machine, or SVM. The goal of an SVM is to take groups of observations and construct boundaries to predict which group future observations belong to based on their measurements. The different groups that must be separated will be called “classes”. SVMs can handle any number of classes, as well as observations of any dimension. SVMs can take almost any shape (including linear, radial, and polynomial, among others), and are generally flexible enough to be used in almost any classification endeavor that the user chooses to undertake.
This tutorial on Support Vector Machines in R is an abbreviated version of p. 359-366 of “Introduction to Statistical Learning with Applications in R” by Gareth James, Daniela Witten, Trevor Hastie and Robert Tibshirani. It was re-implemented in Fall 2016 in tidyverse format by Amelia McNamara and R. Jordan Crouser at Smith College.
1.1 Packages
In this lab, we’ll use the e1071 library in R to demonstrate the support vector classifier
and the SVM. Another option is the LiblineaR library, which is particularly useful for
very large linear problems.
1.2 Support Vector Classifier
The e1071 library contains implementations for a number of statistical
learning methods. In particular, the svm() function can be used to fit a
support vector classifier when the argument kernel="linear" is used. The cost argument allows us to specify the cost of
a violation to the margin. When the cost argument is small, then the margins
will be wide and many support vectors will be on the margin or will
violate the margin. When the cost argument is large, then the margins will
be narrow and there will be few support vectors on the margin or violating
the margin.
We can use the svm() function to fit the support vector classifier for a
given value of the cost parameter. Here we demonstrate the use of this
function on a two-dimensional example so that we can plot the resulting
decision boundary. Let’s start by generating a set of observations, which belong
to two classes:
set.seed(1)
x = matrix(rnorm(20*2), ncol=2)
class = c(rep(-1,10), rep(1,10))
x[class == 1,] = x[class == 1,] + 1Let’s plot the data to see whether the classes are linearly separable:
Nope; not linear. Next, we fit the support vector classifier. Note that in order
for the svm() function to perform classification (as opposed to SVM-based
regression), we must encode the response as a factor:
training_data = data.frame(x = x, class = as.factor(class))
library(e1071)
svmfit = svm(class~.,
data = training_data,
kernel = "linear",
cost = 10,
scale = FALSE)The argument scale = FALSE tells the svm() function not to scale each feature
to have mean zero or standard deviation one; depending on the application,
we might prefer to use scale = TRUE1. We can now plot the support vector classifier by calling the plot() function on the output of the call to svm(), as well as the data used in the call to svm():
The region of feature space that will be assigned to the -1 class is shown in
yellow, and the region that will be assigned to the +1 class is shown in
red. The decision boundary between the two classes is linear (because we
used the argument kernel = "linear"), though due to the way in which the
plotting function is implemented in this library the decision boundary looks
somewhat jagged in the plot. We see that in this case only one observation
is misclassified. (Note also that the second feature is plotted on the x-axis
and the first feature is plotted on the y-axis, in contrast to the behavior of
the usual plot() function in R.)
The support vectors are plotted as crosses and the remaining observations are plotted as circles; we see here that there are seven support vectors. We can determine their identities as follows:
## [1] 1 2 5 7 14 16 17A long sidenote
svm can use C-classification or nu-classification. By default it uses C-classification. C refers to the cost in this case.
The problem with the parameter C is:
- that it can take any positive value
- that it has no direct interpretation.
It is therefore hard to choose correctly and one has to resort to cross validation or direct experimentation to find a suitable value.
In response Schölkopf et al. (2000) reformulated SVM to take a new regularization parameter \(\nu\). This parameter is:
- bounded between 0 and 1
- has a direct interpretation
Interpretation of \(\nu\) (nu):
The parameter nu is an upper bound on the fraction of margin errors and a lower bound of the fraction of support vectors relative to the total number of training examples. For example, if you set it to 0.05 you are guaranteed to find at most 5% of your training examples being misclassified (at the cost of a small margin, though) and at least 5% of your training examples being support vectors.
Relationship between C and nu:
The relation between C and nu is governed by the following formula:
\[\nu = A+B/C\]
A and B are constants which are unfortunately not that easy to calculate.
The takeaway message is that C and \(\nu\) SVM are equivalent regarding their classification power. The regularization in terms of \(\nu\) is easier to interpret compared to C, but the \(\nu\) SVM is usually harder to optimize and runtime doesn’t scale as well as the C variant with number of input samples.
More details (including formulas for A and B) can be found here: Chang CC, Lin CJ - “Training nu-support vector classifiers: theory and algorithms”
We can obtain some basic information about the support vector classifier
fit using the summary() command.
##
## Call:
## svm(formula = class ~ ., data = training_data, kernel = "linear",
## cost = 10, scale = FALSE)
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: linear
## cost: 10
##
## Number of Support Vectors: 7
##
## ( 4 3 )
##
##
## Number of Classes: 2
##
## Levels:
## -1 1This tells us, for instance, that a linear kernel was used with cost = 10, and
that there were seven support vectors, four in one class and three in the
other. What if we instead used a smaller value of the cost parameter?
svmfit = svm(class~.,
data = training_data,
kernel = "linear",
cost = 0.1,
scale = FALSE)
plot(svmfit, training_data)
## [1] 1 2 3 4 5 7 9 10 12 13 14 15 16 17 18 20Now that a smaller value of the cost parameter is being used, we obtain a
larger number of support vectors, because the margin is now wider. Unfortunately,
the svm() function does not explicitly output the coefficients of
the linear decision boundary obtained when the support vector classifier is
fit, nor does it output the width of the margin.
The e1071 library includes a built-in function, tune(), to perform cross-validation. By default, tune() performs ten-fold cross-validation on a set
of models of interest. In order to use this function, we pass in relevant
information about the set of models that are under consideration. The
following command indicates that we want to compare SVMs with a linear
kernel, using a range of values of the cost parameter:
set.seed(1)
tune_out = tune(svm,
class~.,
data = training_data,
kernel = "linear",
ranges = list(cost = c(0.001, 0.01, 0.1, 1,5,10,100)))We can easily access the cross-validation errors for each of these models
using the summary() command:
##
## Parameter tuning of 'svm':
##
## - sampling method: 10-fold cross validation
##
## - best parameters:
## cost
## 0.1
##
## - best performance: 0.05
##
## - Detailed performance results:
## cost error dispersion
## 1 1e-03 0.55 0.4377975
## 2 1e-02 0.55 0.4377975
## 3 1e-01 0.05 0.1581139
## 4 1e+00 0.15 0.2415229
## 5 5e+00 0.15 0.2415229
## 6 1e+01 0.15 0.2415229
## 7 1e+02 0.15 0.2415229The tune() function stores the best model obtained, which can be accessed as
follows:
##
## Call:
## best.tune(method = svm, train.x = class ~ ., data = training_data,
## ranges = list(cost = c(0.001, 0.01, 0.1, 1, 5, 10, 100)), kernel = "linear")
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: linear
## cost: 0.1
##
## Number of Support Vectors: 16
##
## ( 8 8 )
##
##
## Number of Classes: 2
##
## Levels:
## -1 1As usual, the predict() function can be used to predict the class label on a set of
test observations, at any given value of the cost parameter. Let’s
generate a test data set:
xtest = matrix(rnorm(20*2), ncol = 2)
ytest = sample(c(-1,1), 20, rep = TRUE)
xtest[ytest == 1,] = xtest[ytest == 1,] + 1
test_data = data.frame(x = xtest, class = as.factor(ytest))Now we predict the class labels of these test observations. Here we use the best model obtained through cross-validation in order to make predictions:
## true
## predicted -1 1
## -1 9 1
## 1 2 8Thus, with this value of cost, 19 of the test observations are correctly
classified. What if we had instead used cost = 0.01?
svmfit = svm(class~., data = training_data, kernel = "linear", cost = .01, scale = FALSE)
class_pred = predict(svmfit, test_data)
table(predicted = class_pred, true = test_data$class)## true
## predicted -1 1
## -1 11 6
## 1 0 3In this case one additional observation is misclassified.
Now consider a situation in which the two classes are linearly separable.
Then we can find a separating hyperplane using the svm() function. First we’ll give our simulated data a little nudge so that they are linearly separable:
x[class == 1,] = x[class == 1,] + 0.5
ggplot(data.frame(x), aes(X1, X2, colour = factor(class))) +
geom_point()
Now the observations are just barely linearly separable. We fit the support
vector classifier and plot the resulting hyperplane, using a very large value
of cost so that no observations are misclassified.
training_data2 = data.frame(x = x, class = as.factor(class))
svmfit = svm(class~., data = training_data2, kernel = "linear", cost = 1e5)
summary(svmfit)##
## Call:
## svm(formula = class ~ ., data = training_data2, kernel = "linear",
## cost = 1e+05)
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: linear
## cost: 1e+05
##
## Number of Support Vectors: 3
##
## ( 1 2 )
##
##
## Number of Classes: 2
##
## Levels:
## -1 1
No training errors were made and only three support vectors were used.
However, we can see from the figure that the margin is very narrow (because
the observations that are not support vectors, indicated as circles, are very close to the decision boundary). It seems likely that this model will perform
poorly on test data. Let’s try a smaller value of cost:
##
## Call:
## svm(formula = class ~ ., data = training_data2, kernel = "linear",
## cost = 1)
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: linear
## cost: 1
##
## Number of Support Vectors: 7
##
## ( 4 3 )
##
##
## Number of Classes: 2
##
## Levels:
## -1 1
Using cost = 1, we misclassify a training observation, but we also obtain
a much wider margin and make use of seven support vectors. It seems
likely that this model will perform better on test data than the model with
cost = 1e5.
1.3 Figuring out the importance of variables in SVM
Unlike many other standard methods, unfortunately, SVM doesn’t have many straightforward methods to calculate variable importance in R. Hence, rather than depending on a method written by someone else, let us jump into the task of writing a code of our own for the task at hand! Remember, the only condition for the following code is, in your data set, the predicted variable should be the last one in the data frame.
set.seed(1)
varimp_df <- NULL # df with results
ptm1 <- proc.time() # Start the clock!
for(i in 1:(ncol(training_data)-1)) { # the last var is the dep var, hence the -1
#smp_size <- floor(0.70 * nrow(training_data)) # 70/30 split
#train_ind <- sample(seq_len(nrow(training_data)), size = smp_size)
training <- training_data[, -c(i)] # receives all the df less 1 var
testing <- test_data[, -c(i)]
tune.out.linear <- tune(svm, class ~ .,
data = training,
kernel = "linear",
ranges = list(cost =10^seq(1, 3, by = 0.5))) # you can choose any range you see fit
svm.linear <- svm(class ~ .,
kernel = "linear",
data = training,
cost = tune.out.linear[["best.parameters"]][["cost"]])
train.pred.linear <- predict(svm.linear, testing)
testing_y <- as.factor(testing$class)
conf.matrix.svm.linear <- caret::confusionMatrix(train.pred.linear, testing_y)
sens=conf.matrix.svm.linear[["byClass"]][["Sensitivity"]]
spec=conf.matrix.svm.linear[["byClass"]][["Specificity"]]
varimp_df <- rbind(varimp_df,data.frame(
var_no=i,
variable=colnames(training_data[i]),
cost_param=tune.out.linear[["best.parameters"]][["cost"]],
accuracy=conf.matrix.svm.linear[["overall"]][["Accuracy"]],
kappa=conf.matrix.svm.linear[["overall"]][["Kappa"]],
sensitivity=sens,
specificity=spec,
DOR=(sens*spec)/((1-sens)*(1-spec))))
runtime1 <- as.data.frame(t(data.matrix(proc.time() - ptm1)))$elapsed # time for running this loop
runtime1 # divide by 60 and you get minutes, /3600 you get hours
}
varimp_df## var_no variable cost_param accuracy kappa sensitivity specificity DOR
## 1 1 x.1 31.62278 0.6 0.2079208 0.5454545 0.6666667 2.4
## 2 2 x.2 10.00000 0.8 0.5876289 0.9090909 0.6666667 20.0Make sure that the predicted variable is last in the data frame.
I mainly needed Diagnostic odds ratio values from the models and by removing one predictor at a time, I can evaluate the importance of each predictor, i.e. by removing a predictor, the smallest DOR of the model(less predictor number i) means that the predictor has the most importance. You need to know on what indicator you will attribute importance.
You can also add another for loop inside to change between kernels, i.e. linear, polynomial, radial, but you might have to account for the other parameters,e.g. gamma. Change class with your target variable and training_data with your data frame.
1.4 Support Vector Machine
In order to fit an SVM using a non-linear kernel, we once again use the svm()
function. However, now we use a different value of the parameter kernel.
To fit an SVM with a polynomial kernel we use kernel="polynomial", and
to fit an SVM with a radial kernel we use kernel="radial". In the former
case we also use the degree argument to specify a degree for the polynomial
kernel, and in the latter case we use gamma to specify a
value of \(\gamma\) for the radial basis kernel.
The gamma parameter defines how far the influence of a single training example reaches, with low values meaning ‘far’ and high values meaning ‘close’. The gamma parameters can be seen as the inverse of the radius of influence of samples selected by the model as support vectors. If gamma is too large, the radius of the area of influence of the support vectors only includes the support vector itself and no amount of regularization with C will be able to prevent overfitting. When gamma is very small, the model is too constrained and cannot capture the complexity or “shape” of the data. The region of influence of any selected support vector would include the whole training set.
Let’s generate some data with a non-linear class boundary:
set.seed(1)
x = matrix(rnorm(200*2), ncol = 2)
x[1:100,] = x[1:100,]+2
x[101:150,] = x[101:150,]-2
class = c(rep(1,150),rep(2,50))
nonlinear_data = data.frame(x = x, class = as.factor(class))
ggplot(nonlinear_data, aes(x.1, x.2, colour = factor(class))) +
geom_point()
See how one class is kind of stuck in the middle of another class? This suggests that we might want to use a radial kernel in our SVM. Now let’s randomly split this data into training and testing groups, and then fit
the training data using the svm() function with a radial kernel and \(\gamma = 1\):
##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionnonlinear_train = nonlinear_data %>%
sample_frac(0.5)
nonlinear_test = nonlinear_data %>%
setdiff(nonlinear_train)
svmfit = svm(class~., data = nonlinear_train, kernel = "radial", gamma = 1, cost = 1)
plot(svmfit, nonlinear_train)
Not too shabby! The plot shows that the resulting SVM has a decidedly non-linear boundary. We can see from the figure that there are a fair number of training errors in this SVM fit. If we increase the value of cost, we can reduce the number of training errors:
svmfit = svm(class~., data = nonlinear_train, kernel = "radial", gamma = 1, cost = 1e5)
plot(svmfit, nonlinear_train)
However, this comes at the price of a more irregular decision boundary that seems to be at risk of overfitting the data. We can perform cross-validation using tune() to select the best choice of
\(\gamma\) and cost for an SVM with a radial kernel:
set.seed(1)
tune_out = tune(svm, class~., data = nonlinear_train, kernel = "radial",
ranges = list(cost = c(0.1,1,10,100,1000), gamma = c(0.5,1,2,3,4)))
bestmod = tune_out$best.model
summary(bestmod)##
## Call:
## best.tune(method = svm, train.x = class ~ ., data = nonlinear_train,
## ranges = list(cost = c(0.1, 1, 10, 100, 1000), gamma = c(0.5,
## 1, 2, 3, 4)), kernel = "radial")
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: radial
## cost: 1
##
## Number of Support Vectors: 30
##
## ( 15 15 )
##
##
## Number of Classes: 2
##
## Levels:
## 1 2Therefore, the best choice of parameters involves cost = 1 and gamma = 2. We
can plot the resulting fit using the plot() function, and view the test set predictions for this model by applying the predict()
function to the test data.
## pred
## true 1 2
## 1 67 10
## 2 2 2190% of test observations are correctly classified by this SVM. Not bad!
1.5 ROC Curves
The ROCR package can be used to produce ROC curves such as those we saw in earlier lectures. We first write a short function to plot an ROC curve given a vector containing a numerical score for each observation, pred, and a vector containing the class label for each observation, truth:
library(ROCR)
rocplot = function(pred, truth, ...){
predob = prediction(pred, truth)
perf = performance(predob, "tpr", "fpr")
plot(perf,...)}SVMs and support vector classifiers output class labels for each observation. However, it is also possible to obtain fitted values for each observation, which are the numerical scores used to obtain the class labels. For instance, in the case of a support vector classifier, the fitted value for an observation \(X = (X_1,X_2, . . .,X_p)^T\) takes the form \(\hat\beta_0 + \hat\beta_1X_1 + \hat\beta_2X_2 + . . . + \hat\beta_pX_p\).
For an SVM with a non-linear kernel, the equation that yields the fitted value is given in (9.23) on p. 352 of the ISLR book. In essence, the sign of the fitted value determines on which side of the decision boundary the observation lies. Therefore, the relationship between the fitted value and the class prediction for a given observation is simple: if the fitted value exceeds zero then the observation is assigned to one class, and if it is less than zero than it is assigned to the other.
In order to obtain the fitted values for a given SVM model fit, we
use decision.values=TRUE when fitting svm(). Then the predict() function
will output the fitted values. Let’s fit models using the \(\gamma\) selected by cross-validation, and a higher value, which will produce a more flexible fit:
svmfit_opt = svm(class~., data = nonlinear_train, kernel = "radial",
gamma = 2, cost = 1, decision.values = TRUE)
svmfit_flex = svm(class~., data = nonlinear_train, kernel = "radial",
gamma = 50, cost = 1, decision.values = TRUE)Now we can produce the ROC plot to see how the models perform on both the training and the test data:
par(mfrow = c(1,2))
# Plot optimal parameter model's performance on training data
fitted_opt_train = attributes(predict(svmfit_opt, nonlinear_train,
decision.values = TRUE))$decision.values
rocplot(fitted_opt_train, nonlinear_train$class, main = "Training Data")
# Add more flexible model's performance to the plot
fitted_flex_train = attributes(predict(svmfit_flex, nonlinear_train,
decision.values = TRUE))$decision.values
rocplot(fitted_flex_train, nonlinear_train$class, add = TRUE, col = "red")
# Plot optimal parameter model's performance on test data
fitted_opt_test = attributes(predict(svmfit_opt, nonlinear_test,
decision.values = TRUE))$decision.values
rocplot(fitted_opt_test, nonlinear_test$class, main = "Test Data")
# Add more flexible model's performance to the plot
fitted_flex_test = attributes(predict(svmfit_flex, nonlinear_test,
decision.values = TRUE))$decision.values
rocplot(fitted_flex_test, nonlinear_test$class, add = TRUE, col = "red")
1.6 SVM with Multiple Classes
If the response is a factor containing more than two levels, then the \({\tt svm()}\) function will perform multi-class classification using the one-versus-one approach. We explore that setting here by generating a third class of observations:
x = rbind(x, matrix(rnorm(50*2), ncol = 2))
class = c(class, rep(0,50))
x[class == 0,2] = x[class == 0,2]+2
data_3_classes = data.frame(x = x, class = as.factor(class))
ggplot(data_3_classes, aes(x.1, x.2, colour = factor(class))) +
geom_point()
Fitting an SVM to multiclass data uses identical syntax to fitting a simple two-class model:
svmfit = svm(class~., data = data_3_classes, kernel = "radial", cost = 10, gamma = 1)
plot(svmfit, data_3_classes)
The e1071 library can also be used to perform support vector regression,
if the response vector that is passed in to svm() is numerical rather than a
factor.
SVM tries to maximize the distance between the separating plane and the support vectors. If one feature (i.e. one dimension in this space) has very large values, it will dominate the other features when calculating the distance. If you rescale all features (e.g. to [0, 1]), they all have the same influence on the distance metric. So if the feature values are extremely disproportionate to each other, it is a good practice to rescale the features.↩︎