ch9_lab.Rmd

---
title: "9.6 Lab: Support Vector Machines"
output: 
  github_document:
    md_extensions: -fancy_lists+startnum
  html_notebook: 
    md_extensions: -fancy_lists+startnum
---

```{r setup, message=FALSE, warning=FALSE}
library(tidyverse)
library(e1071)
library(ROCR)
library(ISLR)
```

## 9.6.1 Support Vector Classifier

`e1071::svm` allows to fit a support vector classifier when using `kernel="linear"`. Instead of setting a "budget" for violations of the margins, it uses a `cost` parameter.

### Two-dimensions, not linearly separable

```{r}
set.seed(1989)
x <- matrix(rnorm(20*2), ncol = 2)
y <- c(rep(-1, 10), rep(1, 10))

# mean shifting rows based on the assigned class
x[y==1,] <- x[y==1,] + 1

plot(x, col = (3-y))
```

The classes are not linearly separable.

Fitting the support vector classifier. We must encode the response variable as factor in order to do classification (and not regression).
```{r}
data <- 
  tibble(
  x1 = x[,1],
  x2 = x[,2],
  y = as.factor(y)
)

svmfit <-
  svm(
    y ~ .,
    data = data,
    kernel = "linear",
    cost = 10,
    scale = FALSE
  )
```

`scale` argument is for scaling variables to have mean 0 and std. deviation 1.

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

Support vectors are plotted as crosses and the remaining observations are plotted as circles.

We can get the indexes of the support vectors:
```{r}
svmfit$index
```
We can obtain some basic information about the support vector classifier fit using the `summary()` command:
```{r}
summary(svmfit)
```

This tells that there were 14 support vectors, 7 in each class.

Let's check what happens when we use a smaller value in the cost parameter.
```{r}
svmfit2 <- svm(y ~ ., data = data, kernel = "linear", cost = 0.1,
               scale = FALSE)

plot(svmfit2, data)
```

```{r}
svmfit2$index
```
Now almost all the train observations are support vectors. This is because lowering the cost parameter is equivalent to increasing the violation budget C. We end up with model with higher bias, but lower variance.

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.

We can use `tune()` for performing cross validation using a set of parameters. By default it suses 10-fold CV
```{r}
tune_out <- tune(svm,
                 y ~ .,
                 data = data,
                 kernel = "linear",
                 ranges = list(
                   cost = c(0.001, 0.01, 0.1, 1, 5, 10, 100)
                 ))

summary(tune_out)
```

`cost = 100` results in the lowest CV error rate.

Accessing the best model obtained:
```{r}
bestmod <- tune_out$best.model
summary(bestmod)
```

We can use `predict()` to predict classes for test observations:
```{r}
set.seed(2020)

xtest <- matrix(rnorm(20*2), ncol = 2)
ytest <- sample(c(-1, 1), 20, rep = TRUE)
xtest[ytest==1,] <- xtest[ytest==1,] + 1

testdat <- 
  tibble(
    x1 = xtest[,1],
    x2 = xtest[,2],
    y = ytest
  )

ypred <- predict(bestmod, testdat)

table(predict = ypred, truth = testdat$y)
```

Using the optimal value of `cost`, 17 of the 20 test observations are correctly classified.

What if we used `cost = 0.01`?
```{r}
svmfit3 <- svm(y ~ ., data = data, kernel = "linear",
               cost = .01, scale = FALSE)

ypred2 <- predict(svmfit3, testdat)

table(predict = ypred2, truth = testdat$y)
```

Now only 15 (out of 20) test observations are correctly classified.

### Two dimmensions, linearly separable

```{r}
x_separable <- x
x_separable[y == 1,] <- x_separable[y == 1,] + 2.5

plot(x_separable, col = (y+5)/2, pch = 19)
```

Now we fit the classifier using a very large value of `cost` so that no observations are misclassified (even if we end up with few support vectors and high variance)
```{r}
data_separable <- 
  tibble(
    x1 = x_separable[,1],
    x2 = x_separable[,2],
    y = as.factor(y)
  )

svmfit4 <- svm(y ~ ., data = data_separable,
               kernel = "linear", cost = 1e5, scale = FALSE)

summary(svmfit4)
```

```{r}
plot(svmfit4, data_separable)
```

All the training observations are correctly classified, and we end up with just 3 support vectors (Xs in the plot). It could be that this model will perform poorly on test data.

Trying now with a smaller value of cost:
```{r}
svmfit5 <- svm(y ~ ., data = data_separable,
               kernel = "linear", cost = 0.5, scale = FALSE)

summary(svmfit5)
```

```{r}
plot(svmfit5, data_separable)
```

Now we misclassify one observations, but the margin is wider and we have one additional support vector. This model would probably perform better on test data.

## 9.6.2 Support Vector Machine

To use a non-linear kernel we just have to change the argument `kernel` to `"polynomial"` or `"radial"`. In the former case we also have to specify a `degree` and in the latter, a `gamma`.

```{r data_svm}
set.seed(1989)

x_svm <- matrix(rnorm(200*2), ncol = 2)
x_svm[1:100,] <- x_svm[1:100,] + 2
x_svm[101:150,] <- x_svm[101:150,] - 2
y_svm <- c(rep(1, 150), rep(2, 50))

data_svm <- 
  tibble(
    x1 = x_svm[,1],
    x2 = x_svm[,2],
    y = as.factor(y_svm)
  )

plot(x_svm, col = y_svm)
```

We see that the class boundary is indeed non-linear.

Splitting into test and train sets.
```{r}
set.seed(1989)
data_svm_train <- data_svm %>% 
  sample_frac(0.5)

data_svm_test <- data_svm %>% 
  anti_join(data_svm_train)
```

```{r}
svmfit6 <- svm(y ~ ., data = data_svm_train,
               kernel = "radial", gamma = 1,
               cost = 1)

plot(svmfit6, data_svm_train)
```

```{r}
summary(svmfit6)
```

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}
svmfit7 <- svm(y ~ ., data = data_svm_train,
               kernel = "radial", gamma = 1,
               cost = 1e5)

plot(svmfit7, data_svm_train)
```

We can also increase `gamma` to add flexibility:
```{r}
svmfit8 <- svm(y ~ ., data = data_svm_train,
               kernel = "radial", gamma = 20,
               cost = 1e5)

plot(svmfit8, data_svm_train)
```

Using `tune()` for cross-validation of the parameters:
```{r}
set.seed(1989)

tune_out_svm <- 
  tune(svm, y ~ .,
       data = data_svm_train, kernel = "radial",
       ranges = list(
         cost = c(0.1, 1, 10, 100, 1000),
         gamma = c(0.5, 1, 2, 3, 4)
       ))

summary(tune_out_svm)
```
In this case, the best choice involves `cost = 100` and `gamma = 0.5`. We can see how well it performs on test data:
```{r}
table(
  true = data_svm_test$y,
  pred = predict(tune_out_svm$best.model, newdata = data_svm_test)
)
```
85% of test observations are correctly classified :)

## 9.6.3 ROC Curves

Creating function to plot ROC curve from vectors with predicted score and actual values:
```{r}
rocplot <- function(pred, truth, ...) {
  predob <- prediction(pred, truth)
  perf <- performance(predob, "tpr", "fpr")
  plot(perf, ...)
}
```

To get the fitted values for a given SVM model fit we use `decision.values=TRUE` when fitting `svm()`. Then `predict()` will output the fitted values:
```{r}
svmfit_opt <- svm(
  y ~ .,
  data = data_svm_train,
  kernel = "radial",
  gamma = 0.5,
  cost = 100,
  decision.values = TRUE
)

fitted_svm <- attributes(
  predict(svmfit_opt, data_svm_train, decision.values = TRUE)
)$decision.values
```

Now we can produce the ROC plot:
```{r}
rocplot(fitted_svm, data_svm_train$y, main = "Training Data")
```

Trying to increase flexibility:
```{r}
svmfit_flex <- svm(
  y ~ .,
  data = data_svm_train,
  kernel = "radial",
  gamma = 50,
  cost = 1,
  decision.values = TRUE
)

fitted_svm_flex <- 
  attributes(
    predict(svmfit_flex, data_svm_train, decision.values = TRUE)
  )$decision.values

rocplot(fitted_svm, data_svm_train$y, main = "Training Data")
rocplot(fitted_svm_flex, data_svm_train$y, add = TRUE, col = "red")
```
 
 This looks like it has an overfitting problem. Let's see what happens on test data:
```{r}
fitted_svm_test <- attributes(
  predict(svmfit_opt, data_svm_test, decision.values = TRUE)
)$decision.values

fitted_svm_flex_test <- attributes(
  predict(svmfit_flex, data_svm_test, decision.values = TRUE)
)$decision.values

rocplot(fitted_svm_test, data_svm_test$y, main = "Test Data")
rocplot(fitted_svm_flex_test, data_svm_test$y,
        add = TRUE, col="red")
```
 
there are no big differences between the models.

### 9.6.4 SVM with Multiple Classes

If there are multiple classes, `svm()` will perform classification using the one-vs-one approach.

```{r}
set.seed(1989)
x_mc <- rbind(x, matrix(rnorm(50*2), ncol = 2)) #mc stands for multi-class
y_mc <- c(y, rep(0, 50))
x_mc[y==0,2] <- x_mc[y==0,2] + 2

data_mc <- 
  tibble(
    x1 = x_mc[,1],
    x2 = x_mc[,2],
    y = factor(y_mc)
  )

plot(x_mc, col = (y_mc+2))
```
```{r}
svmfit_mc <- 
  svm(
    y ~ .,
    data = data_mc,
    kernel = "radial",
    cost = 10,
    gamma = 5
  )

plot(svmfit_mc, data_mc)
```
### 9.6.5 Application to Gene Expression Data

```{r}
names(ISLR::Khan)
```
```{r}
dim(Khan$xtrain)
```

```{r}
length(Khan$ytrain)
```
```{r}
dim(Khan$xtest)
```
```{r}
length(Khan$ytest)
```
Values that the response variable can take:
```{r}
table(Khan$ytrain)
```

There are 63 train and 20 test observations (tissue samples) with expression measurements for 2308 genes.

We will use a support vector approach to predict cancer subtype using gene expression measurements.

Since there are a very large number of features relative to the number of observations we should use a linear kernel, because the additional flexibility that will result from using a polynomial or radial kernel is unnecessary.

```{r}
data_genes <- 
  as_tibble(Khan$xtrain) %>% 
  mutate(
    y = as.factor(Khan$ytrain)
  )

data_genes
```

```{r}
svm_genes <- 
  svm(
    y ~ ., 
    data = data_genes,
    kernel = "linear",
    cost = 10
  )

summary(svm_genes)
```

```{r}
table(predicted = svm_genes$fitted, truth = data_genes$y)
```

There are no training errors, which is not surprising since there is high dimensionality. Let's check how the model performs on test data:
```{r}
data_test_genes <- 
  as_tibble(Khan$xtest) %>% 
  mutate(y = as.factor(Khan$ytest))

pred_test_genes <- predict(svm_genes, data_test_genes)

table(predict = pred_test_genes, truth = data_test_genes$y)
```

We see that using `cost=10` yields two test set errors on this data.