Support Vector Machines

Question 5

We have seen that we can fit an SVM with a non-liner kernel order to perform classification using 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 behind to two classes with a quadratic decision boundary between them.

library(caret)
library(ggplot2)
library(e1071)

set.seed(1)

x1 <- runif(500) - 0.5
x2 <- runif(500) - 0.5

y <- 1 * (x1^2 - x2^2 > 0)

data <- data.frame(
  X1 = x1,
  X2 = x2,
  y = factor(y)
)

head(data)

summary(data)

       X1                  X2            y      
 Min.   :-0.498163   Min.   :-0.498685   0:261  
 1st Qu.:-0.241871   1st Qu.:-0.242053   1:239  
 Median :-0.023730   Median :-0.000780          
 Mean   :-0.004345   Mean   : 0.003728          
 3rd Qu.: 0.234146   3rd Qu.: 0.251293          
 Max.   : 0.496077   Max.   : 0.499931          

b.)

Plot the observations, colored according to their class labels.

(x1 on x axis, x2 on y axis)

plot_boundary <- function(model, data, title){
  
  x1_grid <- seq(min(data$x1), max(data$x1), length.out = 200)
  x2_grid <- seq(min(data$x2), max(data$x2), length.out = 200)
  
  grid <- expand.grid(
    x1 = x1_grid,
    x2 = x2_grid
  )
  
  grid$pred <- predict(model, newdata = grid)
  
  ggplot(data, aes(x=x1, y=x2, color=y)) +
    geom_point(size=2) +
    geom_contour(
      data = grid,
      aes(z=as.numeric(pred)),
      breaks=c(1.5),
      color="black",
      linewidth=1
    ) +
    labs(
      title=title,
      x="X1",
      y="X2"
    ) +
    theme_minimal()
}

ggplot(data, aes(x=x1, y=x2, color=y)) +
  geom_point(size=2) +
  labs(
    title="Original Data",
    x="X1",
    y="X2"
  ) +
  theme_minimal()

c.)

Fit a logistic regression model to the data, using X1 and X2 as predictors.

logistic_linear <- glm(
  y ~ x1 + x2,
  data=data,
  family=binomial
)

summary(logistic_linear)

Call:
glm(formula = y ~ x1 + x2, family = binomial, data = data)

Coefficients:
             Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.087260   0.089579  -0.974    0.330
x1           0.196199   0.316864   0.619    0.536
x2          -0.002854   0.305712  -0.009    0.993

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 692.18  on 499  degrees of freedom
Residual deviance: 691.79  on 497  degrees of freedom
AIC: 697.79

Number of Fisher Scoring iterations: 3

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.

linear_prob <- predict(
  logistic_linear,
  type="response"
)

linear_pred <- ifelse(
  linear_prob > .5,
  1,
  0
)

data$linear_pred <- factor(linear_pred)


ggplot(data, aes(x=x1, y=x2, color=linear_pred)) +
  geom_point(size=2) +
  labs(
    title="Linear Logistic Regression Predictions",
    x="X1",
    y="X2"
  ) +
  theme_minimal()

e.)

Now fit a logistic regression model to the data using non-linear functions of X1 and X2 as predictors.

logistic_quad <- glm(
  y ~ x1 + x2 + I(x1^2) + I(x2^2) + I(x1*x2),
  data=data,
  family=binomial
)

summary(logistic_quad)

Call:
glm(formula = y ~ x1 + x2 + I(x1^2) + I(x2^2) + I(x1 * x2), family = binomial, 
    data = data)

Coefficients:
             Estimate Std. Error z value Pr(>|z|)
(Intercept)    -10.16     713.54  -0.014    0.989
x1              42.10   15492.58   0.003    0.998
x2             -66.81   14788.95  -0.005    0.996
I(x1^2)      16757.98  519013.02   0.032    0.974
I(x2^2)     -16671.65  508668.89  -0.033    0.974
I(x1 * x2)    -206.38   41802.81  -0.005    0.996

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 6.9218e+02  on 499  degrees of freedom
Residual deviance: 3.5810e-06  on 494  degrees of freedom
AIC: 12

Number of Fisher Scoring iterations: 25

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. The decision boundary should be obviously non-linear. If it is not then repeat (a)-(e) until you come up with an example in which the predicted class labels are obviously non-linear.

quad_prob <- predict(
  logistic_quad,
  type="response"
)

quad_pred <- ifelse(
  quad_prob > .5,
  1,
  0
)

data$quad_pred <- factor(quad_pred)


ggplot(data, aes(x=x1, y=x2, color=quad_pred)) +
  geom_point(size=2) +
  labs(
    title="Quadratic Logistic Regression Predictions",
    x="X1",
    y="X2"
  ) +
  theme_minimal()

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.

svm_linear <- svm(
  y ~ x1 + x2,
  data=data,
  kernel="linear",
  cost=1
)

summary(svm_linear)

Call:
svm(formula = y ~ x1 + x2, data = data, kernel = "linear", 
    cost = 1)


Parameters:
   SVM-Type:  C-classification 
 SVM-Kernel:  linear 
       cost:  1 

Number of Support Vectors:  480

 ( 239 241 )


Number of Classes:  2 

Levels: 
 0 1

svm_linear_pred <- predict(
  svm_linear,
  data
)

data$svm_linear_pred <- svm_linear_pred


ggplot(data, aes(x=x1, y=x2, color=svm_linear_pred)) +
  geom_point(size=2) +
  labs(
    title="Linear Support Vector Classifier",
    x="X1",
    y="X2"
  ) +
  theme_minimal()

h.)

Fit a SVM vector classifier to the data with X1 and X2 as predictors. Obtain a class prediction for each training observations, colored according to the predicted class labels.

svm_nonlinear <- svm(
  y ~ x1 + x2,
  data=data,
  kernel="radial",
  cost=1,
  gamma=1
)

summary(svm_nonlinear)

Call:
svm(formula = y ~ x1 + x2, data = data, kernel = "radial", 
    cost = 1, gamma = 1)


Parameters:
   SVM-Type:  C-classification 
 SVM-Kernel:  radial 
       cost:  1 

Number of Support Vectors:  147

 ( 73 74 )


Number of Classes:  2 

Levels: 
 0 1

svm_nonlinear_pred <- predict(
  svm_nonlinear,
  data
)

data$svm_nonlinear_pred <- svm_nonlinear_pred


ggplot(data, aes(x=x1, y=x2, color=svm_nonlinear_pred)) +
  geom_point(size=2) +
  labs(
    title="Nonlinear SVM (Radial Kernel)",
    x="X1",
    y="X2"
  ) +
  theme_minimal()

i.) Results

Comment on your results.

The Support Vector Model and linear regression model both perform poorly, indicating a non-linear relationship between X1 and X2. Meanwhile the non-linear SVM produced the most flexible decision boundary since the kernel allows for a better separation of observations.

Question 7

In this problem, you will use support vector approached 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.

library(ISLR2)
library(e1071)
library(caret)
library(ggplot2)

data(Auto)

median.mpg <- median(Auto$mpg)
Auto$mpg01 <- ifelse(Auto$mpg > median.mpg, 1, 0)

Auto$mpg01 <- factor(Auto$mpg01,
                     levels = c(0,1),
                     labels = c("low", "high"))
table(Auto$mpg01)

 low high 
 196  196 

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. Note you will need to fit the classifier without the gas mileage variable to produce sensible results.

set.seed(1)
train_index <- createDataPartition(Auto$mpg01,
                                   p = .7,
                                   list = FALSE)
train <- Auto[train_index,]
test <- Auto[-train_index,]

train.svm <- train[, !names(train) %in% c("mpg")]

Error in .rs.exprMutatesPackageLibrary(part) : 
  argument "part" is missing, with no default

test.svm <- test[, !names(test) %in% c("mpg")]

Error in .rs.exprMutatesPackageLibrary(part) : 
  argument "part" is missing, with no default

cost_values <- c(0.001,0.01,0.1,1,10,100)

linear_results <- data.frame(
  Cost = cost_values,
  CV_Error = NA
)


for(i in 1:length(cost_values)){
  
  svm_model <- svm(
    mpg01 ~ .,
    data=train.svm,
    kernel="linear",
    cost=cost_values[i],
    cross=10
  )
  
  linear_results$CV_Error[i] <- svm_model$tot.accuracy
}


linear_results

linear_results$CV_Error <- 100 - linear_results$CV_Error

linear_results

c.)

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

radial_results <- data.frame()


for(cost in c(0.01,0.1,1,10,100)){
  
  for(gamma in c(0.001,0.01,0.1,1)){
    
    svm_model <- svm(
      mpg01 ~ .,
      data=train.svm,
      kernel="radial",
      cost=cost,
      gamma=gamma,
      cross=10
    )
    
    radial_results <- rbind(
      radial_results,
      data.frame(
        Cost=cost,
        Gamma=gamma,
        Accuracy=svm_model$tot.accuracy,
        Error=100-svm_model$tot.accuracy
      )
    )
  }
}


radial_results

radial_results[which.min(radial_results$Error),]

Error in .rs.exprMutatesPackageLibrary(expr) : 
  argument "part" is missing, with no default

poly_results <- data.frame()


for(cost in c(0.01,0.1,1,10,100)){
  
  for(degree in c(2,3,4)){
    
    svm_model <- svm(
      mpg01 ~ .,
      data=train.svm,
      kernel="polynomial",
      cost=cost,
      degree=degree,
      cross=10
    )
    
    
    poly_results <- rbind(
      poly_results,
      data.frame(
        Cost=cost,
        Degree=degree,
        Accuracy=svm_model$tot.accuracy,
        Error=100-svm_model$tot.accuracy
      )
    )
  }
}


poly_results

poly_results[which.min(poly_results$Error),]

Error in .rs.exprMutatesPackageLibrary(expr) : 
  argument "part" is missing, with no default

d.)

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

best_linear <- linear_results[
which.min(linear_results$CV_Error),]

best_radial <- radial_results[
which.min(radial_results$Error),]

best_poly <- poly_results[
which.min(poly_results$Error),]


best_linear

best_radial

best_poly

ggplot(linear_results,
       aes(x=Cost,
           y=CV_Error))+
  geom_line()+
  geom_point()+
  scale_x_log10()+
  labs(
    title="Linear SVM Cross Validation Error",
    x="Cost",
    y="CV Error (%)"
  )

ggplot(radial_results,
       aes(x=Gamma,
           y=Error,
           color=factor(Cost)))+
  geom_line()+
  geom_point()+
  scale_x_log10()+
  labs(
    title="Radial SVM Error by Gamma and Cost",
    x="Gamma",
    y="CV Error (%)",
    color="Cost"
  )

ggplot(poly_results,
       aes(x=Degree,
           y=Error,
           color=factor(Cost)))+
  geom_point(size=3)+
  geom_line()+
  labs(
    title="Polynomial SVM Error by Degree and Cost",
    x="Polynomial Degree",
    y="CV Error (%)",
    color="Cost"
  )

Question 8

This problem involves the OJ data set which is a 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.

library(ISLR2)
library(e1071)
library(caret)

data(OJ)

set.seed(1)

train <- sample(1:nrow(OJ), 800)

train.oj <- OJ[train,]
test.oj <- OJ[-train,]

b.)

Fit a support vector classifier to the remaining data using cost = .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.

svm.linear.oj <- svm(Purchase ~ ., 
                     data = train.oj,
                     kernel = "linear",
                     cost = .01,
                     scale = TRUE)
summary(svm.linear.oj)

Call:
svm(formula = Purchase ~ ., data = train.oj, kernel = "linear", 
    cost = 0.01, scale = TRUE)


Parameters:
   SVM-Type:  C-classification 
 SVM-Kernel:  linear 
       cost:  0.01 

Number of Support Vectors:  435

 ( 219 216 )


Number of Classes:  2 

Levels: 
 CH MM

c.)

What are the training and test error rates?

train.pred <- predict(svm.linear.oj, train.oj)

test.pred <- predict(svm.linear.oj, test.oj)

confusionMatrix(train.pred, train.oj$Purchase)

Confusion Matrix and Statistics

          Reference
Prediction  CH  MM
        CH 420  75
        MM  65 240
                                          
               Accuracy : 0.825           
                 95% CI : (0.7969, 0.8507)
    No Information Rate : 0.6062          
    P-Value [Acc > NIR] : <2e-16          
                                          
                  Kappa : 0.6314          
                                          
 Mcnemar's Test P-Value : 0.4469          
                                          
            Sensitivity : 0.8660          
            Specificity : 0.7619          
         Pos Pred Value : 0.8485          
         Neg Pred Value : 0.7869          
             Prevalence : 0.6062          
         Detection Rate : 0.5250          
   Detection Prevalence : 0.6188          
      Balanced Accuracy : 0.8139          
                                          
       'Positive' Class : CH              
                                          

confusionMatrix(test.pred, test.oj$Purchase)

Confusion Matrix and Statistics

          Reference
Prediction  CH  MM
        CH 153  33
        MM  15  69
                                          
               Accuracy : 0.8222          
                 95% CI : (0.7713, 0.8659)
    No Information Rate : 0.6222          
    P-Value [Acc > NIR] : 6.769e-13       
                                          
                  Kappa : 0.6083          
                                          
 Mcnemar's Test P-Value : 0.01414         
                                          
            Sensitivity : 0.9107          
            Specificity : 0.6765          
         Pos Pred Value : 0.8226          
         Neg Pred Value : 0.8214          
             Prevalence : 0.6222          
         Detection Rate : 0.5667          
   Detection Prevalence : 0.6889          
      Balanced Accuracy : 0.7936          
                                          
       'Positive' Class : CH              
                                          

train.err <- mean(train.pred != train.oj$Purchase)
train.err

[1] 0.175

d.)

Use the tune() function to select an optimal cost. Consider values in the range 0.01 to 10.

set.seed(1)

tune.linear <- tune(svm,
                    Purchase ~ .,
                    data = train.oj,
                    kernel = "linear",
                    ranges = list(cost = c(0.01,
                                           0.1,
                                           1,
                                           5,
                                           10)))

summary(tune.linear)

Parameter tuning of ‘svm’:

- sampling method: 10-fold cross validation 

- best parameters:

- best performance: 0.1725 

- Detailed performance results:

best.linear <- tune.linear$best.model

summary(best.linear)

Call:
best.tune(METHOD = svm, train.x = Purchase ~ ., data = train.oj, 
    ranges = list(cost = c(0.01, 0.1, 1, 5, 10)), kernel = "linear")


Parameters:
   SVM-Type:  C-classification 
 SVM-Kernel:  linear 
       cost:  0.1 

Number of Support Vectors:  342

 ( 171 171 )


Number of Classes:  2 

Levels: 
 CH MM

e.)

Compute the training and test error rates using this new value for cost.

train.pred.best <- predict(best.linear, train.oj)
test.pred.best <- predict(best.linear, test.oj)

train.err.best <- mean(train.pred.best != train.oj$Purchase)
test.err.best <- mean(test.pred.best != test.oj$Purchase)

cat("Train Error:", train.err.best, "\n")

Train Error: 0.165 

cat("Test Error:", test.err.best, "\n")

Test Error: 0.162963 

f.)

Repeat parts (b) through (e) using a support vector machine with a radial kernel. Use the default value for gamma.

svm.radial <- svm(
  Purchase ~ .,
  data = train.oj,
  kernel = "radial",
  cost = 0.01
)

summary(svm.radial)

Call:
svm(formula = Purchase ~ ., data = train.oj, kernel = "radial", 
    cost = 0.01)


Parameters:
   SVM-Type:  C-classification 
 SVM-Kernel:  radial 
       cost:  0.01 

Number of Support Vectors:  634

 ( 319 315 )


Number of Classes:  2 

Levels: 
 CH MM

# Errors

train.pred.radial <- predict(svm.radial, train.oj)

test.pred.radial <- predict(svm.radial, test.oj)

train.err.rad <- mean(train.pred.radial != train.oj$Purchase)

test.err.rad <- mean(test.pred.radial != test.oj$Purchase)

cat("Train Error:", train.err.rad, "\n")

Train Error: 0.39375 

cat("Test Error:", test.err.rad, "\n")

Test Error: 0.3777778 

#tuning

set.seed(1)

tune.radial <- tune(
  svm,
  Purchase ~ .,
  data = train.oj,
  kernel = "radial",
  ranges = list(cost = c(0.01,
                         0.1,
                         1,
                         5,
                         10))
)

summary(tune.radial)

Parameter tuning of ‘svm’:

- sampling method: 10-fold cross validation 

- best parameters:

- best performance: 0.17125 

- Detailed performance results:

best.radial <- tune.radial$best.model

train.error.radial <-
  mean(predict(best.radial, train.oj) != train.oj$Purchase)

test.error.radial <-
  mean(predict(best.radial, test.oj) != test.oj$Purchase)

cat("Radial Training Error:", train.error.radial, "\n")

Radial Training Error: 0.15125 

cat("Radial Test Error:", test.error.radial, "\n")

Radial Test Error: 0.1851852 

g.)

Repeat parts (b) through (e) using a support vector machine with a polynomial kernel. Set degree = 2.

svm.poly <- svm(
  Purchase ~ .,
  data = train.oj,
  kernel = "polynomial",
  degree = 2,
  cost = 0.01
)

summary(svm.poly)

Call:
svm(formula = Purchase ~ ., data = train.oj, kernel = "polynomial", 
    degree = 2, cost = 0.01)


Parameters:
   SVM-Type:  C-classification 
 SVM-Kernel:  polynomial 
       cost:  0.01 
     degree:  2 
     coef.0:  0 

Number of Support Vectors:  636

 ( 321 315 )


Number of Classes:  2 

Levels: 
 CH MM

# Errors

train.pred.poly <- predict(svm.poly, train.oj)

test.pred.poly <- predict(svm.poly, test.oj)

mean(train.pred.poly != train.oj$Purchase)

[1] 0.3725

mean(test.pred.poly != test.oj$Purchase)

[1] 0.3666667

#Tuning:

set.seed(1)

tune.poly <- tune(
  svm,
  Purchase ~ .,
  data = train.oj,
  kernel = "polynomial",
  degree = 2,
  ranges = list(cost = c(0.01,
                         0.1,
                         1,
                         5,
                         10))
)

summary(tune.poly)

Parameter tuning of ‘svm’:

- sampling method: 10-fold cross validation 

- best parameters:

- best performance: 0.18125 

- Detailed performance results:

best.poly <- tune.poly$best.model

train.error.poly <-
  mean(predict(best.poly, train.oj) != train.oj$Purchase)

test.error.poly <-
  mean(predict(best.poly, test.oj) != test.oj$Purchase)

train.error.poly

[1] 0.15

test.error.poly

[1] 0.1888889

h.)

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

The tuned SVM classifier with a linear kernel performed the best when compared to the other SVM classifiers. It achieved both the lowest training and testing errors.