svm_snippet.Rmd

---
title: "SVM - snippet"
author: "Filippo Biscarini"
date: "2023-01-30"
output: html_document
editor_options: 
  markdown: 
    wrap: sentence
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

# Load Packages
library("ISLR") 
library("knitr")
library("e1071")
library("plotly")
library("kernlab")      
library("tidyverse") 
library("RColorBrewer") 
```

## SVM: support vector machines. A super-quick snippet

SVM: developed (in the 1990s) to construct boundaries between groups of observations based on their measurements $\rightarrow$ **classification**.
SVMs are considered to be one of the best **out of the box** classifiers.

SVMs are a generalization of the **maximum margin classifier** (that required the classes to be linearly separable).
SVMs can handle **any number of classes** and sample sizes, and can construct **boundaries of virtually any shape** (linear, polynomial, highly complex/wiggly).

**Examples with synthetic data**: a **merely qualitative** explanation of separating hyperplanes, margins and SVM follows.

```{r synthetic_dataset}
set.seed(223)
set.seed(229)
# Construct sample data set - completely separated
x <- matrix(rnorm(20*2), ncol = 2)
y <- c(rep(-1,10), rep(1,10))
x[y==1,] <- x[y==1,] + 3/2
dat <- data.frame(x=x, y=as.factor(y))

# Plot data
ggplot(data = dat, aes(x = x.2, y = x.1, color = y, shape = y)) + 
  geom_point(size = 2) +
  scale_color_manual(values=c("#000000", "#FF0000")) +
  theme(legend.position = "none")
```

### Separating (hyper)plane and support vectors

We try to find a (hyper)plane that separates the two classes in the feature space:

$$
\beta_0 + \beta_1x_1 + \beta_2x_2 + \ldots + \beta_px_p = 0 
$$

All points for which the equation above holds (is equal to zero) are points that **belong to the hyperplane**.
The separating hyperplane will yield: $\sum_{j=1}^p \beta_j x_j > 0$ for classes labelled `+1`, and $\sum_{j=1}^p \beta_j x_j < 0$ for classes labelled as `-1`.
This can be combined in $y_i(\sum_{j=1}^p \beta_j x_{ij}) > 0$ (where $y_i$ is either $\pm 1$).

We can find **several possible separating (hyper)planes** (see figure above): which shall we pick?

**Maximal margin classifier**: among all separating hyperplanes, find the one that makes the biggest gap or margin between the two classes.
This means **maximising the distance** between the **separating hyperplane** and the **closest points**.

We can use the `svm()` function in the `e1071` *R package* to find this boundary.

```{r pressure, echo=FALSE}
svmfit <- svm(y ~ ., data = dat, kernel = "linear", scale = FALSE)
# Plot Results
plot(svmfit, dat)
```

The "closest points to the boundary" are represented by an `X` in the plot above, and are the **support vectors**.
Only the support vectors do affect the classification line.
The other points (marked with a `o`s) don't affect the calculation of the classification boundary (the separating hyperplane).
The support vectors "support" the hyperplane: if these points change, the separating hyperplane changes as well.

```{r}
# fit model and produce plot
# kernfit <- ksvm(x, y, type = "C-svc", kernel = 'vanilladot')
# plot(kernfit, data = x)
```

```{r}
make.grid=function(x,n=75){
  grange=apply(x,2,range)
  x1=seq(from=grange[1,1],to=grange[2,1],length=n)
  x2=seq(from=grange[1,2],to=grange[2,2],length=n)
  expand.grid(x.1=x1,x.2=x2)
  }
xgrid=make.grid(x)
ygrid=predict(svmfit,xgrid)
plot(xgrid,col=c("red","blue")[as.numeric(ygrid)],pch=20,cex=.2)
points(x,col=y+3,pch=19)
points(x[svmfit$index,],pch=5,cex=2)
```

```{r}
beta=drop(t(svmfit$coefs)%*%x[svmfit$index,])
beta0=svmfit$rho
plot(xgrid,col=c("red","blue")[as.numeric(ygrid)],pch=20,cex=.2)
points(x,col=y+3,pch=19)
points(x[svmfit$index,],pch=5,cex=2)
abline(beta0/beta[2],-beta[1]/beta[2])
abline((beta0-0.5)/beta[2],-beta[1]/beta[2],lty=2)
abline((beta0+0.35)/beta[2],-beta[1]/beta[2],lty=2)
```

#### Support Vector Classifiers

In most cases however, the two classes are not perfectly separable by a hyperplane (e.g. plot below), therefore we need to get creative, e.g.:

-   [we **soften** what we mean by "separates"]{style="color:red"}
-   we **transform** (enrich, enlarge) the feature space so that separation is possible.

Moving away from a perfectly-separating hyperplane will also reduce problems with generalization (remember the bias/variance trade-off and the training/test data).

```{r}
set.seed(129)
# Construct sample data set - not completely separated
x <- matrix(rnorm(20*2), ncol = 2)
y <- c(rep(-1,10), rep(1,10))
x[y==1,] <- x[y==1,] + 1
dat <- data.frame(x=x, y=as.factor(y))

# Plot data set
ggplot(data = dat, aes(x = x.2, y = x.1, color = y, shape = y)) + 
  geom_point(size = 2) +
  scale_color_manual(values=c("#000000", "#FF0000")) +
  theme(legend.position = "none")
```

The **support vector classifier** will allow some observations (data points) to **violate the margin**: some points will be allowed to be inside the margin or even on the other side of the hyperplane (misclassifications) $\rightarrow$ **soft margin** (the support vector classifier is also called *soft margin classifier*).

To do this, **slack (support) variables** $\epsilon_i$ are introduced for each example $i$ in $[1,n]$:

-   $\epsilon_i=0$ if example $i$ is on the correct side of the margin;
-   $\epsilon_i$ is between 0 and 1, if example $i$ is on the wrong side of the margin, BUT on the correct side of the hyperplane
-   $\epsilon_i > 1$ if example $i$ is on the wrong side of the margin, AND on the wrong side of the hyperplane

[**QUESTION**: How many observations/examples (data points) will be allowed to violate the margin/hyperplane?]{style="color:red"}

This is controlled by the **cost** $C$:

$$
\sum_{i=1}^n \epsilon_i \leq C
$$

We use the same function `svm()` to fit the **support vector classifier**; you can see that we now use the **cost** argument of the function to specify the extent to which misclassifications (or simple violations of the margin) are allowed:

```{r}
# Fit Support Vector Classifier model to data set
svmfit <- svm(y~., data = dat, kernel = "linear", cost = 100)
# Plot Results
plot(svmfit, dat)
```

```{r}
# Fit Support Vector Machine model to data set
kernfit <- ksvm(x,y, type = "C-svc", kernel = 'vanilladot', C = 100)
# Plot results
plot(kernfit, data = x)
```

[**QUESTION**: How do we decide how costly these misclassifications actually are?]{style="color:red"}

The cost of the support vector classifier is a **hyperparameter** (tuning parameter) of the model: we already know that to decide on the value for hypeparameters of machine learning models we use **cross-validation**.

We can use the `tune()` function of the `e1071` *R package*:

```{r}
# find optimal cost of misclassification
tune.out <- tune(svm, y~., data = dat, kernel = "linear",
                 ranges = list(cost = c(0.001, 0.01, 0.1, 1, 5, 10, 100)))
# extract the best model
(bestmod <- tune.out$best.model)
```

```{r}
# Create a confusion matrix of predictions vs observations
ypred <- predict(bestmod, dat)
(misclass <- table(predict = ypred, truth = dat$y))
```

Important: only observations that either lie on the margin or violate the margin will affect the hyperplane!
$\rightarrow$ **support vectors**

If $C$ is large there are many support vectors (many observations violating the margin), and viceversa.

### Support Vector Machines

We said that to overcome the limitations of the *maximum margin classifier*, we need to be creative in two ways.
We already saw how softening the margin leads to the *support vector classifier*.
We now look at the second extension of the method, which deals with transforming the feature space to accommodate non-linear separation boundaries:

-   we **soften** what we mean by "separates"
-   [we **transform** (enrich, enlarge) the feature space so that separation is possible]{style="color:red"}

This transformation of the feature space can be obtained by using **quadratic**, **cubic**, or **higher-order polynomial** functions of the original features (or other more complicated functions of the features can be used) $\rightarrow$ this is known as the **kernel trick** (mathematically very efficient).

We are not going too much into the mathematical details, we just say that there are several types of kernels: `linear` ($\rightarrow$ support vector classifier), `polynomial`, `radial`, `sigmoid` etc.

We create a dataset with clearly a **non-linear boundary** between the classes:

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

# Plot data
ggplot(data = dat, aes(x = x.2, y = x.1, color = y, shape = y)) + 
  geom_point(size = 2) +
  scale_color_manual(values=c("#000000", "#FF0000")) +
  theme(legend.position = "none")
```

------------------------------------------------------------------------

**Intuition: feature space enlargement**

We now manually enlarge the feature space to provide some intuition behind SVM.
We add a third variable to the dataset `z`, that has positive values for observations from one class and negative values for observations from the other class.

```{r}

temp <- dat

temp$z <- NA
n1 = table(temp$y)[[1]]
n2 = table(temp$y)[[2]]

temp[temp$y == 1,"z"] <-  -1*abs(rnorm(n = n1, mean = 0, sd = 1)/10)
temp[temp$y == 2,"z"] <-  +1*abs(rnorm(n = n2, mean = 0, sd = 1)/10)


p <- plot_ly(data = temp, 
             x = ~x.1, y = ~x.2, z = ~z,
             type = "scatter3d",
             color = ~y) %>%
  add_markers() %>%
  layout(scene = list(xaxis = list(title = 'x1'),
                      yaxis = list(title = 'x2'),
                      zaxis = list(title = 'z')),
         annotations = list(
           x = 0.005,
           y = 0.01,
           text = 'y',
           xref = 'paper',
           yref = 'paper',
           showarrow = FALSE
         ))

p
```

------------------------------------------------------------------------

Now, with the same `svm()` function we fit a **support vector machine** to the data, specifying the `kernel` argument:

```{r}
# set pseudorandom number generator
set.seed(123)
# sample training data and fit model
train <- base::sample(200,100, replace = FALSE)
svmfit <- svm(y~., data = dat[train,], kernel = "radial", gamma = 1, cost = 1)
# plot classifier
plot(svmfit, dat)
```

```{r}
# Fit radial-based SVM in kernlab
kernfit <- ksvm(x[train,],y[train], type = "C-svc", kernel = 'rbfdot', C = 1, scaled = c())
# Plot training data
plot(kernfit, data = x[train,])
```

### SVM: multiclass classification

[example from: <https://uc-r.github.io/svm>]

```{r}
# construct data set
x <- rbind(x, matrix(rnorm(50*2), ncol = 2))
y <- c(y, rep(0,50))
x[y==0,2] <- x[y==0,2] + 2.5
dat <- data.frame(x=x, y=as.factor(y))

# plot data set
ggplot(data = dat, aes(x = x.2, y = x.1, color = y, shape = y)) + 
  geom_point(size = 2) +
  scale_color_manual(values=c("#000000","#FF0000","#00BA00")) +
  theme(legend.position = "none")
```

```{r}
# fit model
svmfit <- svm(y~., data = dat, kernel = "radial", cost = 10, gamma = 1)
# plot results
plot(svmfit, dat)
```

```{r}
#construct table
ypred <- predict(svmfit, dat)
(misclass <- table(predict = ypred, truth = dat$y))
```

```{r}
# fit and plot
kernfit <- ksvm(as.matrix(dat[,2:1]),dat$y, type = "C-svc", kernel = 'rbfdot', 
                C = 100, scaled = c())

# Create a fine grid of the feature space
x.1 <- seq(from = min(dat$x.1), to = max(dat$x.1), length = 100)
x.2 <- seq(from = min(dat$x.2), to = max(dat$x.2), length = 100)
x.grid <- expand.grid(x.2, x.1)

# Get class predictions over grid
pred <- predict(kernfit, newdata = x.grid)

# Plot the results
cols <- brewer.pal(3, "Set1")
plot(x.grid, pch = 19, col = adjustcolor(cols[pred], alpha.f = 0.05))

classes <- matrix(pred, nrow = 100, ncol = 100)
contour(x = x.2, y = x.1, z = classes, levels = 1:3, labels = "", add = TRUE)

points(dat[, 2:1], pch = 19, col = cols[predict(kernfit)])
```

## Application to a classification problem

Now we use a clinical dataset on breast cancer ($n=699$).
Examples are classified as `benign` or `malignant` based on results from a tissue biopsy:

1.  get the data

```{r}
library("mlbench")
data(BreastCancer)
head(BreastCancer)
```

2.  Split the data $\rightarrow$ Training and test sets:

```{r}
BreastCancer = na.omit(BreastCancer)
BreastCancer <- select(BreastCancer, -Id)
BreastCancer$Class = factor(BreastCancer$Class)
prop=0.8 ## proportion training
vec_train <- sample(nrow(BreastCancer),prop*nrow(BreastCancer))

train_set = BreastCancer[vec_train,]
test_set = BreastCancer[-vec_train,]
```

3.  Fit the SVM model to the training data:

```{r}
svmfit = svm(factor(Class) ~ ., data=train_set, scale=TRUE, kernel="radial", cost=5, decision.values=TRUE)
```

4.  Measure performance on the training set:

```{r}
table(svmfit$fitted, train_set$Class)
```

5.  Measure performance on the **test set**:

```{r}
svm.pred <- predict(svmfit, test_set)
table(svm.pred, test_set$Class)
```

## SVM for regression problems (SVR)

In most regression methods the objective is to minimize the error: e.g.
OLS (Ordinary Least Squares), or extensions of OLS like Ridge Regression, Lasso-penalised regression etc.

SVR changes the perspective: aim for an **acceptable range for the errors, not the minimum**, and **minimise the model coefficients** subject to the **constraint that the errors fall within the range**:

-   $min\frac{1}{2}\sum_{j=1}^p(w_j)^2$
-   $|y_i-\mathbf{w}_i\mathbf{x}_i| \leq \epsilon$

As it turns out, just like with SVM for classification, also in SVR we need to **soften the margin** and allow some errors to exceed the acceptable margin (subject to their overall sum being within a given cost $C$ -tunable hyperparameter).

-   $min \left( \frac{1}{2}\sum_{j=1}^p(w_j)^2 + C\sum_{i=1}^n|\xi| \right)$
-   $|y_i-\mathbf{w}_i\mathbf{x}_i| \leq \epsilon + |\xi_i|$

We use data on chick weights as a function of time and diet:

1.  Get the data:

```{r}
data("ChickWeight")
head(ChickWeight)
ChickWeight <- ChickWeight |> mutate(Diet = as.factor(Diet))
```

```{r}
p <- ggplot(ChickWeight, aes(x = Time, y = weight, color = Diet)) + geom_jitter(aes(color=Diet), width=0.2, alpha=0.3) 
p <- p + geom_smooth(se=FALSE)
p
```

2.  A little preprocessing:

```{r}
y = ChickWeight$weight
diet <- model.matrix(~-1+Diet, data=ChickWeight)
X = bind_cols("Time"=ChickWeight[,2], diet)
```

3.  Split the data into the training and test sets:

```{r}
prop=0.8 ## proportion training
vec_train <- sample(nrow(X),prop*nrow(X))

train_X = X[vec_train,]
train_y = y[vec_train]

test_X = X[-vec_train,]
test_y = y[-vec_train]
```

4.  Fit the SVR model:

```{r}
svmfit = svm(x = train_X, y = train_y, type = "eps-regression")
print(svmfit)
```

5.  Fine-tuning of the hyperparameters:

```{r}
# find optimal cost of misclassification
tune.out <- tune(svm, train.x = train_X, train.y = train_y, kernel="radial",
                 ranges = list(epsilon = seq(0.1,1.1,0.1), cost = 2^(seq(0.5,10,.5))))
# extract the best model
(bestmod <- tune.out$best.model)
```

```{r}
#Predict using SVM regression
# predYsvm = predict(svmfit, test_X)
preds = predict(bestmod, test_X)

df = data.frame("y"=test_y, "y_hat"=preds)
ggplot(df, aes(y, y_hat)) + geom_point()
```

6.  Measure model performance on the **test data**:

```{r}
rmse = sqrt(mean((preds-test_y)^2))
print(paste("RMSE is:", rmse))
```

Is this a **large error or not**?

```{r}
print(paste("The RMSE is ", round(rmse/mean(ChickWeight$weight),3)*100, "% of the average value for the response variable (chick weight)"))
print(paste("The RMSE is", round(rmse/sd(ChickWeight$weight),2), "standard deviations of the response"))
```

Distribution of prediction errors (test data):

```{r}
hist(abs(preds-test_y), xlab = "errors", breaks = 20, main = "Distribution of errors")
```

## SVM with tidymodels

```{r}
library("tidymodels")
```

We first define the **SVM model**:

```{r 'svm_with_tidymodels'}
svm_mod <-
  svm_rbf(cost = tune(), rbf_sigma = tune()) %>%
  set_mode("classification") %>%
  set_engine("kernlab")
```

Then the **preprocessing recipe**:

```{r}
breast_rec <-
  recipe(Class ~ ., data = BreastCancer)  %>%
  # remove any zero variance predictors
  step_zv(all_predictors()) %>% 
  # remove any linear combinations
  step_lincomb(all_numeric())
```

And the **cross-validation** partitioning for the **fine-tuning of the hyperparameters**:

```{r}
breast_rs <- vfold_cv(BreastCancer, v = 5, repeats = 5)
```

We choose a target metric, AUC:

```{r}
roc_vals <- metric_set(roc_auc)
```

We now specify that we want to save the predictions from fine-tuning (cross-validation):

```{r}
ctrl <- control_grid(verbose = FALSE, save_pred = TRUE)
```

Now we do the fine tuning with default grid values:

```{r}
formula_res <-
  svm_mod %>% 
  tune_grid(
    Class ~ .,
    resamples = breast_rs,
    metrics = roc_vals,
    control = ctrl
  )
formula_res
```

```{r}
# formula_res %>% 
#   select(.metrics) %>% 
#   dplyr::slice(1) %>% 
#   pull(1)
```

From the `tidymodels` object, we can extract the AUC values corresponding to the tested values for the hyperparameters (averages over 25 replicates: 5-fold CV repeated 5 times):

```{r}
estimates <- collect_metrics(formula_res)
estimates
```

The function `show_best()` shows the best `n` models (default, the top 5 modfels):

```{r}
top_models <- show_best(formula_res, metric = "roc_auc")
top_models
```

We can also collect all predictions from cross-validation (fine-tuning):

```{r}
collect_predictions(formula_res)
```

And plot their distribution:

```{r}
augment(formula_res) %>%
  ggplot(aes(.pred_malignant, color = Class)) +
  geom_boxplot()
```
Now, we focus on the **best model** from corss-validation:

```{r}
svm_best <- 
  formula_res %>% 
  collect_metrics() %>% 
  arrange(mean) |>
  dplyr::slice(1)

svm_best <- top_models |> dplyr::slice(1)
svm_best
```

And plot the ROC curve for the best model:

```{r}
svm_auc <- 
  formula_res %>% 
  collect_predictions(parameters = svm_best) %>% 
  roc_curve(Class, .pred_benign) %>% 
  mutate(model = "SVM")

autoplot(svm_auc)
```

Finally, we can look at predictions (confusion matrix) from one repetition and one validation fold:

```{r}
formula_res %>% 
  collect_predictions(parameters = svm_best) %>% 
  filter(id == "Repeat1", id2 == "Fold1") |>
  mutate(pred = factor(ifelse(.pred_benign >= 0.5, "benign", "malignant"))) |>
  conf_mat(Class, pred)
```