ch9.Rmd

---
title: "9.7 Exercises"
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)
```


## Conceptual

(1) This problem involves hyperplanes in two dimensions.

![](ch9_files/exc1.png) 
(2) We have seen that in $p = 2$ dimensions, a linear decision boundary takes the form $β_0+β_1X_1+β_2X_2 = 0$. We now investigate a non-linear decision boundary.

![](ch9_files/exc-2.png)
![](ch9_files/exc2_c.png)  

All observations except (-1,1) are assigned to the blue class. 

(d) Argue that while the decision boundary in (c) is not linear in terms of $X_1$ and $X_2$, it is linear in terms of $X_1, X_1^2, X_2$ and $X_2^2$.

We can see this by expanding the equation of the boundary. Doing this we arrive get the following expression, which is clearly linear in terms of $X_1, X_1^2, X_2$ and $X_2^2$. 

![](ch9_files/exc2_d.png) 

(3)

![](ch9_files/exec3_a_b.png)

(c) Describe the classification rule for the maximal margin classifier.

A: Classify to "Red" if $-(1/2) + X_1 - X_2 < 0$ and classify to "Blue" otherwise.

![](ch9_files/exec3_d_e.png)
Support vectors are circled.

(f) Argue that a slight movement of the seventh observation would not affect the maximal margin hyperplane.

The seventh observation (4,1) is not a support vector and also is relatively far away of the margin (as can be seen in the figure above), so a small movement of it would not touch the margin, and thus would not move the maximal margin hyperplane.

![](ch9_files/exec3_g.png)

![](ch9_files/exec3_h.png)

## Applied

(4) Generate a simulated two-class data set with 100 observations and two features in which there is a visible but non-linear separation between the two classes. Show that in this setting, a support vector machine with a polynomial kernel (with degree greater than 1) or a radial kernel will outperform a support vector classifier on the training data. Which technique performs best on the test data? Make plots and report training and test error rates in order to back up your assertions.

```{r}
set.seed(1991)
sim_data <- 
  tibble(
    x1 = c(rnorm(30, mean = 1),
           rnorm(30, mean = 3),
           rnorm(40, mean = 2)),
    x2 = c(rnorm(30, mean = 2),
           rnorm(30, mean = -1),
           rnorm(40, mean = -4)),
    class = c(rep("red", 30), rep("blue", 30), rep("red", 40))
  ) %>%
  mutate(class = as.factor(class))

sim_data %>% 
  ggplot(aes(x1, x2, color = class)) +
  geom_point() +
  scale_color_identity()
```

Splitting between test and train:
```{r}
set.seed(1989)
sim_data_train <- 
  sim_data %>% 
  sample_frac(0.5)

sim_data_test <- 
  sim_data %>% 
  anti_join(sim_data_train)
```

Training 3 models: linear, polynomial, and radial.
```{r}
svm_lineal <- svm(class ~ .,
                  data = sim_data_train,
                  kernel = "linear",
                  cost = 10)

svm_polynomial <- svm(class ~ .,
                  data = sim_data_train,
                  kernel = "polynomial",
                  degree = 3,
                  cost = 10)

svm_radial <- svm(class ~ .,
                  data = sim_data_train,
                  kernel = "radial",
                  gamma = 1,
                  cost = 10)
```

```{r}
sim_data_test <- 
  sim_data_test %>% 
  modelr::add_predictions(svm_lineal, var = "pred_linear") %>% 
  modelr::add_predictions(svm_polynomial, var = "pred_poly") %>% 
  modelr::add_predictions(svm_radial, var = "pred_radial")

sim_data_test
```
Performance of linear kernel:
```{r}
sim_data_test %>% 
  select(class, pred_linear) %>% 
  table()
```
Performance of polynomial kernel:
```{r}
sim_data_test %>% 
  select(class, pred_poly) %>% 
  table()
```
Performance of radial kernel:
```{r}
sim_data_test %>% 
  select(class, pred_radial) %>% 
  table()
```
We see that the classifier with radial kernel has the best performance on test data (just 2 of 50 observations missclasified, versus 14 and 11 missclassifications when we use other kernels).

Plot of radial kernel:
```{r}
plot(svm_radial, data = sim_data_train)
```

```{r}
plot(svm_polynomial, data = sim_data_train)
```

```{r}
plot(svm_lineal, data = sim_data_train)

```

(5) We have seen that we can fit an SVM with a non-linear kernel in order to perform classification using a non-linear decision boundary. We will now see that we can also obtain a non-linear decision boundary by performing logistic regression using non-linear transformations of the features

(a) Generate a data set with $n = 500$ and $p = 2$, such that the observations belong to two classes with a quadratic decision boundary between them.
```{r}
sim_data2 <- 
  tibble(
    x1 = runif(500) - 0.5,
    x2 = runif(500) - 0.5,
    y = 1*(x1^2-x2^2 > 0)
  )

sim_data2
```
(b) Plot the observations
```{r}
ggplot(sim_data2,
       aes(x1, x2, color = factor(y))) +
  geom_point()
```

(c) Fit a logistic regression model to the data, using $X_1$ and $X_2$ as predictors. 

```{r}
lreg_sim2 <- glm(y ~ x1 + x2, data = sim_data2, family = "binomial")

summary(lreg_sim2)
```

(d) Apply this model to the training data in order to obtain a predicted class label for each training observation. Plot the observations, colored according to the predicted class labels. The decision boundary should be linear.

```{r}
sim_data2 %>% 
  modelr::add_predictions(lreg_sim2, var = "pred_lreg", 
                          type = "response") %>% 
  mutate(pred_lreg_class = ifelse(pred_lreg > 0.5, 1, 0)) %>% 
  ggplot(aes(x1, x2, color = factor(pred_lreg_class))) +
  geom_point() +
  labs(color = "predicted class")
```

(e) Now fit a logistic regression model to the data using non-linear functions of X1 and X2 as predictors (e.g. $X^2$, $X_1 \times X_2$, and so forth).
```{r}
lreg2_sim2 <- 
  glm(y ~ x1 * x2 + I(x1^2) * I(x2^2) + 
        x1:I(x1^2) + x2:I(x2^2), data = sim_data2, 
      family = "binomial")

summary(lreg2_sim2)
```

(f) Apply this model to the training data in order to obtain a predicted class label for each training observation. Plot the observations, colored according to the predicted class labels.

```{r}
sim_data2 %>% 
  modelr::add_predictions(lreg2_sim2, var = "pred_lreg", 
                          type = "response") %>% 
  mutate(pred_lreg_class = ifelse(pred_lreg > 0.5, 1, 0)) %>% 
  ggplot(aes(x1, x2, color = factor(pred_lreg_class))) +
  geom_point() +
  labs(color = "predicted class")
```

(g) Fit a support vector classifier to the data with X1 and X2 as predictors. Obtain a class prediction for each training observation. Plot the observations, colored according to the predicted class labels

```{r}
svm_sim2 <- 
  svm(factor(y) ~ ., data = sim_data2,
    kernel = "linear",
    cost = 1)

sim_data2 %>% 
  modelr::add_predictions(svm_sim2, var = "pred_svm") %>% 
  ggplot(aes(x1, x2, color = pred_svm)) +
  geom_point() +
  labs(color = "predicted class")
```

(h) Fit a SVM using a non-linear kernel to the data. Obtain a class prediction for each training observation. Plot the observations, colored according to the predicted class labels.
```{r}
svm_nl_sim2 <- 
  svm(factor(y) ~ ., data = sim_data2,
    kernel = "polynomial",
    degree = 2,
    cost = 1)

sim_data2 %>% 
  modelr::add_predictions(svm_nl_sim2, var = "pred_svm") %>% 
  ggplot(aes(x1, x2, color = pred_svm)) +
  geom_point() +
labs(color = "predicted class")
```

(i) Comment on your results.

A: We see that the linear support vector classifier and the regression with the original feature space yield to similar results: a decision border that is just a line (an hyperplane in $p=2$). Similarly, both the SVM with polynomial kernel and the regression with non-linear functions of the predictors succeed in learning the true data generating process (because both of them use an enlarged feature space that allows a non-linear boundary in $p=2$).

(6) At the end of Section 9.6.1, it is claimed that in the case of data that is just barely linearly separable, a support vector classifier with a small value of cost that misclassifies a couple of training observations may perform better on test data than one with a huge value of cost that does not misclassify any training observations. You will now investigate this claim.

(a) Generate two-class data with $p = 2$ in such a way that the classes are just barely linearly separable.

```{r}
set.seed(1989)
sim_data3 <- 
  tibble(
    x1 = runif(500, -0.5, 0.5),
    x2 = runif(500, -0.5, 0.5),
    class = factor(ifelse(x1-2*x2 > 0.1, 1, -1))
  )

sim_data3 %>% 
  ggplot(aes(x1,x2, color = class)) +
  geom_point()
```

(b) Compute the cross-validation error rates for support vector classifiers with a range of cost values. How many training errors are misclassified for each value of cost considered, and how does this relate to the cross-validation errors obtained?

```{r}
cost_range <- c(0.01, 0.1, 1, 5, 10, 100, 
                       1000, 1e+04, 1e+05, 1e+09)

tune_out <- tune(svm,
                 class ~ .,
                 data = sim_data3,
                 kernel = "linear",
                 ranges = list(
                   cost = cost_range
                 ))

summary(tune_out)
```
```{r}
qplot(factor(cost), error, data = tune_out$performances, geom = "col")
```

Looking for misclassified observations in training data:
```{r}
train_svm <- function(cost) {
  svm(class ~ .,
      data = sim_data3,
      kernel = "linear",
      cost = cost)
}

count_train_errors <- function(model) {
  sim_data3 %>% 
  modelr::add_predictions(model) %>% 
  mutate(error = class != pred) %>% 
  pull(error) %>% 
  sum()
}


svm_by_cost <- 
  tibble(
  cost_range = cost_range,
  model = map(cost_range, train_svm),
  n_train_errors = map_dbl(model, count_train_errors)
)

svm_by_cost %>% 
  ggplot(aes(factor(cost_range), n_train_errors)) +
  geom_col() + 
  geom_label(aes(label = n_train_errors)) +
  labs(x = "cost", y = "N. of training errors")
```

The number of training observations misclassified always goes downwards as we increase the cost parameter, but the cross validation starts going up after certain cost value.

(c) Generate an appropriate test data set, and compute the test errors corresponding to each of the values of cost considered. Which value of cost leads to the fewest test errors, and how does this compare to the values of cost that yield the fewest training errors and the fewest cross-validation errors?

```{r}
#New dataset with the same data-generating process but another seed
set.seed(2020)
sim_data4 <- 
  tibble(
    x1 = runif(500, -0.5, 0.5),
    x2 = runif(500, -0.5, 0.5),
    class = factor(ifelse(x1-2*x2 > 0.1, 1, -1))
  )

count_test_errors <- function(model) {
  sim_data4 %>% 
  modelr::add_predictions(model) %>% 
  mutate(error = class != pred) %>% 
  pull(error) %>% 
  sum()
}

svm_by_cost %>% 
  mutate(n_test_errors = map_dbl(model, count_test_errors)) %>% 
  ggplot(aes(factor(cost_range), n_test_errors)) +
  geom_col() + 
  geom_label(aes(label = n_test_errors)) +
  labs(x = "cost", y = "N. of test errors")
```

`cost = 1000` leads to the fewest test errors. This is close to what we get with CV errors, where values 10, 100 and 1000 all minimized the CV error. 

(7) In this problem, you will use support vector approaches in order to predict whether a given car gets high or low gas mileage based on the `Auto` data set.

(a) Create a binary variable that takes on a 1 for cars with gas mileage above the median, and a 0 for cars with gas mileage below the median.

```{r}
auto <- ISLR::Auto %>% 
  as_tibble() %>% 
  mutate(high_mileage = factor(ifelse(mpg > median(mpg), 1, 0))) %>%
  select(-mpg)
```

(b) Fit a support vector classifier to the data with various values of cost, in order to predict whether a car gets high or low gas mileage. Report the cross-validation errors associated with different values of this parameter. Comment on your results.

```{r}
tune_out_auto <- 
  tune(
  svm,
  high_mileage ~ .,
  data = auto,
  kernel = "linear",
  ranges = list(cost = c(1e-09, 1e-06, 1e-04, cost_range))
)

tune_out_auto$performances
```

`cost = 1e-02` achieves the minimum CV error.

(c) Now repeat (b), this time using SVMs with radial and polynomial basis kernels, with different values of gamma and degree and cost. Comment on your results.

```{r}
tune_out_auto_poly <- 
  tune(
  svm,
  high_mileage ~ .,
  data = auto,
  kernel = "polynomial",
  ranges = list(cost = c(1e-09, 1e-06, 1e-04, cost_range),
                degree = 2:10)
)

tune_out_auto_poly$performances %>% 
  arrange(error)
```


```{r}
tune_out_auto_radial <- 
  tune(
  svm,
  high_mileage ~ .,
  data = auto,
  kernel = "radial",
  ranges = list(cost = c(1e-09, 1e-06, 1e-04, cost_range),
                gamma = c(0.001, 0.01, 0.1, 1, 10, 100, 1000))
)

tune_out_auto_radial$performances %>% 
  arrange(error)
```
Here the lowest CV error is achieved with `cost = 5` and `gamma = 1e-01`. This error value is lower than the minimum errors obtained with linear and polynomial kernels.

(d) Make some plots to back up your assertions in (b) and (c).

Linear kernel. CV error for different `cost` values:
```{r}
tune_out_auto$performances %>% 
  ggplot(aes(factor(cost), error, group = 1)) +
  geom_line() +
  geom_point() +
  geom_vline(xintercept = which.min(tune_out_auto$performances$error),
             color = "red")
```

Polynomial kernel. CV error for different `degree` and `cost` values:
```{r}
tune_out_auto_poly$performances %>% 
  ggplot(aes(factor(cost), factor(degree), fill = error)) +
  geom_tile() +
  scale_fill_viridis_b(direction = -1) +
  labs(x = "cost",
       y = "degree") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
```

Radial kernel. CV error for different `gamma` and `cost` values:
```{r}
tune_out_auto_radial$performances %>% 
  ggplot(aes(factor(cost), factor(gamma), fill = error)) +
  geom_tile() +
  scale_fill_viridis_b(direction = -1) +
  labs(x = "cost",
       y = "gamma") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
```

Comparing the minimum CV error for each kernel
```{r}
get_minimum_error <- function(performances_df) {
  performances_df %>% 
    pull(error) %>% 
    min()
}

tibble(
  kernel = c("linear", "polynomial", "radial"),
  performances = list(tune_out_auto$performances,
                      tune_out_auto_poly$performances,
                      tune_out_auto_radial$performances),
  minimum_error = map_dbl(performances, get_minimum_error)
) %>% 
  ggplot(aes(kernel, minimum_error, fill = minimum_error)) +
  scale_fill_viridis_b(direction = -1) +
  geom_col()
```

The radial kernel allows us to obtain the minimum cross-validation error.

(8) This problem involves the OJ data set which is part of the ISLR package.

(a) Create a training set containing a random sample of 800 observations, and a test set containing the remaining observations.

```{r}
set.seed(1989)
oj_train <- ISLR::OJ %>% 
  as_tibble() %>% 
  sample_n(800)

oj_test <- ISLR::OJ %>% 
  as_tibble() %>% 
  anti_join(oj_train)
```

(b) Fit a support vector classifier to the training data using `cost=0.01`, with `Purchase` as the response and the other variables as predictors. Use the `summary()` function to produce summary statistics, and describe the results obtained.

```{r}
svm_oj_linear <- 
  svm(Purchase ~ .,
    data = oj_train,
    kernel = "linear",
    cost = 0.01,
    scale = TRUE)

summary(svm_oj_linear)
```

About half the training observations ended up as support vectors.

(c) What are the training and test error rates?
```{r}
oj_train_pred <- oj_train %>% 
  modelr::add_predictions(svm_oj_linear) %>% 
  mutate(error = pred != Purchase)

oj_test_pred <- oj_test %>% 
  modelr::add_predictions(svm_oj_linear) %>% 
  mutate(error = pred != Purchase)
```

Train error:
```{r}
oj_train_pred$error %>% mean()
```

Test error:
```{r}
oj_test_pred$error %>% mean()
```
(d) Use the `tune()` function to select an optimal `cost`. Consider values in the range 0.01 to 10.
```{r}
tune_linear_oj <- tune(
  svm,
  Purchase ~ .,
  data = oj_train,
  kernel = "linear",
  scale = TRUE,
  ranges = list(cost = c(0.01, 0.05, 0.1, 0.3, 0.5, 1, 3, 5, 7, 10))
)

tune_linear_oj$performances %>% 
  arrange(error)
```
(e) Compute the training and test error rates using this new value for cost.
```{r}
oj_train_pred <- oj_train %>% 
  modelr::add_predictions(tune_linear_oj$best.model,
                          var = "pred_tune_linear") %>% 
  mutate(error_tune_linear = pred_tune_linear != Purchase)

oj_test_pred <- oj_test %>% 
  modelr::add_predictions(tune_linear_oj$best.model,
                          var = "pred_tune_linear") %>% 
  mutate(error_tune_linear = pred_tune_linear != Purchase)
```

Train error:
```{r}
oj_train_pred$error_tune_linear %>% mean()
```
Test error:
```{r}
oj_test_pred$error_tune_linear %>% mean()
```
(f) Repeat parts (b) through (e) using a support vector machine with a radial kernel. Use the default value for `gamma`.

```{r}
tune_radial_oj <- tune(
  svm,
  Purchase ~ .,
  data = oj_train,
  kernel = "radial",
  scale = TRUE,
  ranges = list(cost = c(0.01, 0.05, 0.1, 0.3, 0.5, 1, 3, 5, 7, 10))
)

tune_radial_oj$performances %>% 
  arrange(error)
```
```{r}
oj_train_pred <- oj_train %>% 
  modelr::add_predictions(tune_radial_oj$best.model,
                          var = "pred_tune_radial") %>% 
  mutate(error_tune_radial = pred_tune_radial != Purchase)

oj_test_pred <- oj_test %>% 
  modelr::add_predictions(tune_radial_oj$best.model,
                          var = "pred_tune_radial") %>% 
  mutate(error_tune_radial = pred_tune_radial != Purchase)
```

Train error:
```{r}
oj_train_pred$error_tune_radial %>% mean()
```
Test error:
```{r}
oj_test_pred$error_tune_radial %>% mean()
```
(g) Repeat parts (b) through (e) using a support vector machine with a polynomial kernel. Set `degree=2`.


```{r}
tune_poly_oj <- tune(
  svm,
  Purchase ~ .,
  data = oj_train,
  kernel = "polynomial",
  scale = TRUE,
  degree = 2,
  ranges = list(cost = c(0.01, 0.05, 0.1, 0.3, 0.5, 1, 3, 5, 7, 10))
)

tune_poly_oj$performances %>% 
  arrange(error)
```
```{r}
oj_train_pred <- oj_train %>% 
  modelr::add_predictions(tune_poly_oj$best.model,
                          var = "pred_tune_poly") %>% 
  mutate(error_tune_poly = pred_tune_poly != Purchase)

oj_test_pred <- oj_test %>% 
  modelr::add_predictions(tune_poly_oj$best.model,
                          var = "pred_tune_poly") %>% 
  mutate(error_tune_poly = pred_tune_poly != Purchase)
```

Train error:
```{r}
oj_train_pred$error_tune_poly %>% mean()
```
Test error:
```{r}
oj_test_pred$error_tune_poly %>% mean()
```

(h) Overall, which approach seems to give the best results on this
data?

SVM with radial kernel and `cost = 0.5` leads to the lower CV error rate.