Chapter_9_Support_Vector_Machines.Rmd

---
title: "Chapter_9_Support_Vector_Machines"
author: "Russ Conte"
date: "8/7/2021"
output: html_document
---
```{r}
library(e1071)
```

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.

We now 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. We begin by generating the observations, which belong to two classes, and checking whether the classes are linearly separable.

```{r}
set.seed(1)
x = matrix(rnorm(20*2), ncol = 2)
y = c(rep(-1, 10), rep(1, 10))
x[y == 1, ] = x[y == 1,] + 1
plot(x, col = (3 - y))
```

They are not. 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 variable. We now create a data frame with the response coded as a factor.

```{r}
dat = data.frame( x = x, y = as.factor(y))
library(e1071)
svmfit = svm(y ~., data = dat, 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, one might prefer to use scale=TRUE.
We can now plot the support vector classifier obtained:

```{r}
plot(svmfit, dat)
```

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:

```{r}
svmfit$index
```

We can obtain some basic information about the support vector classifier fit using the summary() command:

```{r}
summary(svmfit)

```

This 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?

```{r}
svmfit = svm(y~., data = dat, kernel = "linear", cost = 0.1, scale = FALSE)
plot(svmfit, dat)
svmfit$index

```

Now 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. Unfor- tunately, 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.

```{r}
set.seed((1))
tune.out = tune(svm, y~., data = dat, 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:

summary(tune.out)

```

We see that cost=0.1 results in the lowest cross-validation error rate. The tune() function stores the best model obtained, which can be accessed as follows:

```{r}
bestmod = tune.out$best.model
summary(bestmod)
```

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. We begin by generating a test data set.

```{r}
xtest = matrix(rnorm(20*2), ncol = 2)
ytest = sample(c(-1, 1), 20, rep = TRUE)
xtest[ytest == 1,] = xtest[ytest == 1,] + 1
testdat = data.frame(x = xtest, y = as.factor(ytest))
ypred = predict(bestmod, testdat)
table(predict = ypred, truth = testdat$y)
```

Now consider a situation in which the two classes are linearly separable.
Then we can find a separating hyperplane using the svm() function. We first further separate the two classes in our simulated data so that they are linearly separable:

```{r}
x[y == 1,] = x[y==1,] + 0.5
plot(x, col = (y + 5)/2, pch = 19)
```

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.

```{r}
dat = data.frame(x = x, y = as.factor(y))
svmfit = svm(y~., data = dat, kernel = "linear", cost = 1e5)
summary(svmfit)
```

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. We now try a smaller value of cost:

```{r}
svmfit = svm(y~., data = dat, kernel = "linear", cost = 1)
summary(svmfit)
plot(svmfit, dat)
```

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.

## 9.6.2 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 (this is d in (9.22)), and in the latter case we use gamma to specify a value of γ for the radial basis kernel (9.24).

Generate some data with a non-linear class boundry

```{r}
set.seed(1)
x = matrix(rnorm(200*2), ncol = 2)
x[1:100,] = x[1:100,] + 2
x[101:150,] = x[101:150,] - 2
y = c(rep(1, 150), rep(2, 50))
dat = data.frame(x = x, y = as.factor(y))

#Plotting the data makes it clear that the class boundary is indeed non- linear:

plot(x, col = y)
```

The data is randomly split into training and testing groups. We then fit the training data using the svm() function with a radial kernel and γ = 1:

```{r}
train = sample(200, 100)
svmfit = svm(y~., data = dat[train,], kernel = "radial", gamma = 1, cost = 1)
plot(svmfit, dat[train,])
```

The plot shows that the resulting SVM has a decidedly non-linear boundary. The summary() function can be used to obtain some information about the SVM fit:

```{r}
summary(svmfit)
```

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. However, this comes at the price of a more irregular decision boundary that seems to be at risk of overfitting the data.

```{r}
svmfit = svm(y~., data = dat[train,], kernel = "radial", gamma = 1, cost = 1e5)

plot(svmfit, dat[train,])
```


We can perform cross-validation using tune() to select the best choice of γ and cost for an SVM with a radial kernel:

```{r}
set.seed(1)
tune.out = tune(svm, y~., data = dat[train,], kernel = "radial", ranges = list(cost = c(0.1, 1, 10, 100, 1000), gamma = c(0.5, 1, 2, 3, 4)))
summary(tune.out)
```

Therefore, the best choice of parameters involves cost=1 and gamma=0.5. We can view the test set predictions for this model by applying the predict() function to the data. Notice that to do this we subset the dataframe dat using -train as an index set.

```{r}
table(true = dat[-train, "y"], pred = predict(tune.out$best.model, newdata = dat[-train,]))
(67+21) / (67+21+2+10) # 12% misclassification error rate

```

## 9.6.3 ROC Curves

The ROCR package can be used to produce ROC curves such as those in Figures 9.10 and 9.11. 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.

```{r}
library(ROCR)

rocplot=function(pred, truth, ...){
  predob = prediction(pred, truth)
  perf = performance(predob , "tpr", "fpr")
  plot(perf ,...)}

svmfit.opt = svm(y~., data = dat[train,], kernel = "radial", gamma = 2, cost = 1, decision.values = T)
fitted = attributes(predict(svmfit.opt, dat[train,], decision.values = TRUE))$decision.values

# Now we can produce the ROC plot.

par(mfrow = c(1,2))
rocplot(fitted, dat[train, "y"], main = "Training Data")

```

When we compute the ROC curves on the test data, the model with γ = 2 appears to provide the most accurate results.

```{r}
fitted = attributes(predict(svmfit.opt, dat[-train], decision.values = T))$decision.values
rocplot(fitted, dat[-train, "y"], main = "Test Data")
fitted = attributed(predict(svmfit.flex, dat[-train,], decision.values = T))$decision.values
rocplot(fitted, dat[-train, "y"], add = "T", col = "red")

```


## 9.6.4 SVM with Multiple Classes

If the response is a factor containing more than two levels, then the svm() function will perform multi-class classification using the one-versus-one ap- proach. We explore that setting here by generating a third class of obser- vations.

```{r}
library(e1071)
set.seed (1)
x=rbind(x, matrix(rnorm(50*2), ncol=2))
y=c(y, rep(0,50))
x[y==0,2]=x[y==0,2]+2
dat=data.frame(x=x, y=as.factor(y))
par(mfrow=c(1,1))
plot(x,col=(y+1))

# we now fit an SVM to the data

svmfit=svm(y∼., data=dat, kernel="radial", cost=10, gamma=1)
plot(svmfit , dat)

```

## 9.6.5 Application to Gene Expression Data

We now examine the Khan data set, which consists of a number of tissue samples corresponding to four distinct types of small round blue cell tu- mors. For each tissue sample, gene expression measurements are available. The data set consists of training data, xtrain and ytrain, and testing data, xtest and ytest.
We examine the dimension of the data:

```{r}
library(ISLR)
names(Khan)
```

This data set consists of expression measurements for 2,308 genes.
The training and test sets consist of 63 and 20 observations respectively.

```{r}
table(Khan$ytrain)
table(Khan$ytest)
```

We will use a support vector approach to predict cancer subtype using gene expression measurements. In this data set, there are a very large number of features relative to the number of observations. This suggests that we should use a linear kernel, because the additional flexibility that will result from using a polynomial or radial kernel is unnecessary.

```{r}
dat = data.frame(x = Khan$xtrain, y = as.factor(Khan$ytrain))
out = svm(y~., data = dat, kernel = "linear", cost = 10)
summary(out) # not rendering the same as the book

```

```{r}
dat.te = data.frame(x = Khan$xtest, y = as.factor(Khan$ytest))
pred.te = predict(out, newdata = dat.te)
table(pred.te, dat.te$y)
```