library(ISLR)
set.seed(42)

Problem 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. For instance, you can do this as follows:

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

(b) Plot the observations, colored according to their class labels. Your plot should display X1 on the x-axis, and X2 on the yaxis.

plot(x1[y == 0], x2[y == 0], col = "#dc6900", xlab = "X1", ylab = "X2", pch = "+")
points(x1[y == 1], x2[y == 1], col = "steelblue", pch = 4)

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

lm.fit  <-  glm(y ~ x1 + x2, family = binomial)
summary(lm.fit)
## 
## Call:
## glm(formula = y ~ x1 + x2, family = binomial)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -1.301  -1.216   1.075   1.136   1.207  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  0.09978    0.08976   1.112    0.266
## x1          -0.17659    0.30658  -0.576    0.565
## x2          -0.20067    0.30978  -0.648    0.517
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 691.79  on 499  degrees of freedom
## Residual deviance: 691.08  on 497  degrees of freedom
## AIC: 697.08
## 
## 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.

data <- data.frame(x1 = x1, x2 = x2, y = y)
lm.prob <- predict(lm.fit, data, type = "response")
lm.pred <- ifelse(lm.prob > 0.52, 1, 0)
data.pos <- data[lm.pred == 1, ]
data.neg <- data[lm.pred == 0, ]
plot(data.pos$x1, data.pos$x2, col = "steelblue", xlab = "X1", ylab = "X2", pch = "+")
points(data.neg$x1, data.neg$x2, col = "#dc6900", pch = 4)

(e) Now fit a logistic regression model to the data using non-linear functions of X1 and X2 as predictors (e.g. X21 , X1×X2, log(X2), and so forth). 

lm.fit <- glm(y ~ poly(x1, 2) + poly(x2, 2) + I(x1 * x2), data = data, family = binomial)
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

(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. 

lm.prob <- predict(lm.fit, data, type = "response")
lm.pred <- ifelse(lm.prob > 0.55, 1, 0)
data.pos <- data[lm.pred == 1, ]
data.neg <- data[lm.pred == 0, ]
plot(data.pos$x1, data.pos$x2, col = "steelblue", xlab = "X1", ylab = "X2", pch = "+")
points(data.neg$x1, data.neg$x2, col = "#dc6900", pch = 4)

  1. As expected the non-linear decision boundary resembles the true decision boundary produced in step (d).

(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.

library(e1071)
svm.fit <- svm(as.factor(y) ~ x1 + x2, data, kernel = "linear", cost = 0.1)
svm.pred <- predict(svm.fit, data)
data.pos <- data[svm.pred == 1, ]
data.neg <- data[svm.pred == 0, ]
plot(data.pos$x1, data.pos$x2, col = "steelblue", xlab = "X1", ylab = "X2", pch = "+")
points(data.neg$x1, data.neg$x2, col = "#dc6900", pch = 4)

  1. Clearly, the linear kernel fails to find the non-linear decision boundary and classifies all points to a single 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.

svm.fit <- svm(as.factor(y) ~ x1 + x2, data, gamma = 1)
svm.pred <- predict(svm.fit, data)
data.pos <- data[svm.pred == 1, ]
data.neg <- data[svm.pred == 0, ]
plot(data.pos$x1, data.pos$x2, col = "steelblue", xlab = "X1", ylab = "X2", pch = "+")
points(data.neg$x1, data.neg$x2, col = "#dc6900", pch = 4)

  1. Once again, the non-linear decision boundary on our predicted labels resembles the true decision boundary.

(i) Comment on your results.

Problem 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.

attach(Auto)

(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. 

gas.med <- median(Auto$mpg)
new.var <- ifelse(Auto$mpg > gas.med, 1, 0)
Auto$mpglevel <- as.factor(new.var)

(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. 

tune.out <- tune(svm, mpglevel ~ ., data = Auto, kernel = "linear", ranges = list(cost = c(0.01, 
    0.1, 1, 5, 10, 100)))
summary(tune.out)
## 
## Parameter tuning of 'svm':
## 
## - sampling method: 10-fold cross validation 
## 
## - best parameters:
##  cost
##     1
## 
## - best performance: 0.007692308 
## 
## - Detailed performance results:
##    cost       error dispersion
## 1 1e-02 0.076538462 0.03813752
## 2 1e-01 0.053589744 0.02824912
## 3 1e+00 0.007692308 0.01238579
## 4 5e+00 0.015384615 0.01324097
## 5 1e+01 0.015384615 0.01324097
## 6 1e+02 0.035833333 0.02165712

Our cross-validation error is minimized for \(\tt{cost}=1\).

(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. 

tune.out <- tune(svm, mpglevel ~ ., data = Auto, kernel = "polynomial", ranges = list(cost = c(0.1, 
    1, 5, 10), degree = c(2, 3, 4)))
summary(tune.out)
## 
## Parameter tuning of 'svm':
## 
## - sampling method: 10-fold cross validation 
## 
## - best parameters:
##  cost degree
##    10      2
## 
## - best performance: 0.5251282 
## 
## - Detailed performance results:
##    cost degree     error dispersion
## 1   0.1      2 0.5535256 0.03684550
## 2   1.0      2 0.5535256 0.03684550
## 3   5.0      2 0.5535256 0.03684550
## 4  10.0      2 0.5251282 0.08581809
## 5   0.1      3 0.5535256 0.03684550
## 6   1.0      3 0.5535256 0.03684550
## 7   5.0      3 0.5535256 0.03684550
## 8  10.0      3 0.5535256 0.03684550
## 9   0.1      4 0.5535256 0.03684550
## 10  1.0      4 0.5535256 0.03684550
## 11  5.0      4 0.5535256 0.03684550
## 12 10.0      4 0.5535256 0.03684550

In this example, the lowest cross-validation error is obtained for \(\tt{cost} = 10\) and \(\tt{degree} = 2\).

tune.out <- tune(svm, mpglevel ~ ., data = Auto, kernel = "radial", ranges = list(cost = c(0.1, 
    1, 5, 10), gamma = c(0.01, 0.1, 1, 5, 10, 100)))
summary(tune.out)
## 
## Parameter tuning of 'svm':
## 
## - sampling method: 10-fold cross validation 
## 
## - best parameters:
##  cost gamma
##    10   0.1
## 
## - best performance: 0.02044872 
## 
## - Detailed performance results:
##    cost gamma      error dispersion
## 1   0.1 1e-02 0.08666667 0.05274083
## 2   1.0 1e-02 0.07397436 0.05590309
## 3   5.0 1e-02 0.05615385 0.04645773
## 4  10.0 1e-02 0.03057692 0.02638876
## 5   0.1 1e-01 0.08166667 0.05769721
## 6   1.0 1e-01 0.05615385 0.03967253
## 7   5.0 1e-01 0.02301282 0.02244393
## 8  10.0 1e-01 0.02044872 0.02646892
## 9   0.1 1e+00 0.54839744 0.03280489
## 10  1.0 1e+00 0.06384615 0.04872938
## 11  5.0 1e+00 0.06378205 0.04228961
## 12 10.0 1e+00 0.06378205 0.04228961
## 13  0.1 5e+00 0.54839744 0.03280489
## 14  1.0 5e+00 0.48455128 0.05079582
## 15  5.0 5e+00 0.48455128 0.05079582
## 16 10.0 5e+00 0.48455128 0.05079582
## 17  0.1 1e+01 0.54839744 0.03280489
## 18  1.0 1e+01 0.51262821 0.04549029
## 19  5.0 1e+01 0.50750000 0.04673390
## 20 10.0 1e+01 0.50750000 0.04673390
## 21  0.1 1e+02 0.54839744 0.03280489
## 22  1.0 1e+02 0.54839744 0.03280489
## 23  5.0 1e+02 0.54839744 0.03280489
## 24 10.0 1e+02 0.54839744 0.03280489

For the radial basis kernel we received \(\tt{cost} = 10\) and \(\tt{gamma} = 0.01\) as the lowest cross-validation error.

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

Hint: In the lab, we used the plot() function for svm objects only in cases with p = 2. When p > 2, you can use the plot() function to create plots displaying pairs of variables at a time. Essentially, instead of typing > plot(svmfit , dat)

where svmfit contains your fitted model and dat is a data frame containing your data, you can type > plot(svmfit , dat , x1∼x4)

in order to plot just the first and fourth variables. However, you must replace x1 and x4 with the correct variable names. To find out more, type ?plot.svm.

svm.linear <- svm(mpglevel ~ ., data = Auto, kernel = "linear", cost = 1)
svm.poly <- svm(mpglevel ~ ., data = Auto, kernel = "polynomial", cost = 10, 
    degree = 2)
svm.radial <- svm(mpglevel ~ ., data = Auto, kernel = "radial", cost = 10, gamma = 0.01)
plotpairs <- function(fit) {
    for (name in names(Auto)[!(names(Auto) %in% c("mpg", "mpglevel", "name"))]) {
        plot(fit, Auto, as.formula(paste("mpg~", name, sep = "")))
    }
}
plotpairs(svm.linear)

plotpairs(svm.poly)

plotpairs(svm.radial)

Problem 8:

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

attach(OJ)

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

train <- sample(dim(OJ)[1], 800)
OJ.train <- OJ[train, ]
OJ.test <- 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.

svm.linear <- svm(Purchase ~ ., kernel = "linear", data = OJ.train, cost = 0.01)
summary(svm.linear)
## 
## Call:
## svm(formula = Purchase ~ ., data = OJ.train, kernel = "linear", cost = 0.01)
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  linear 
##        cost:  0.01 
## 
## Number of Support Vectors:  437
## 
##  ( 219 218 )
## 
## 
## Number of Classes:  2 
## 
## Levels: 
##  CH MM

This support vector classifier created 437 support vectors out of the 800 training points. From these, 219 belong to level \(\tt{CH}\) and remaining 218 belong to level \(\tt{MM}\).

(c) What are the training and test error rates?

train.pred <- predict(svm.linear, OJ.train)
table(OJ.train$Purchase, train.pred)
##     train.pred
##       CH  MM
##   CH 419  55
##   MM  71 255
(73 + 59)/(430 + 59 + 73 + 238)
## [1] 0.165
test.pred = predict(svm.linear, OJ.test)
table(OJ.test$Purchase, test.pred)
##     test.pred
##       CH  MM
##   CH 156  23
##   MM  26  65
(26 + 20)/(144 + 20 + 26 + 80)
## [1] 0.1703704

The error rate for our training set is \(16.5\)% and test error rate is \(17\)%.

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

tune.out <- tune(svm, Purchase ~ ., data = OJ.train, kernel = "linear", ranges = list(cost = 10^seq(-2, 
    1, by = 0.25)))
summary(tune.out)
## 
## Parameter tuning of 'svm':
## 
## - sampling method: 10-fold cross validation 
## 
## - best parameters:
##  cost
##  0.01
## 
## - best performance: 0.16375 
## 
## - Detailed performance results:
##           cost   error dispersion
## 1   0.01000000 0.16375 0.04466309
## 2   0.01778279 0.16375 0.04267529
## 3   0.03162278 0.17250 0.04116363
## 4   0.05623413 0.16625 0.03729108
## 5   0.10000000 0.16625 0.03775377
## 6   0.17782794 0.16750 0.04417453
## 7   0.31622777 0.17125 0.04372023
## 8   0.56234133 0.17125 0.04604120
## 9   1.00000000 0.17250 0.04556741
## 10  1.77827941 0.17500 0.04409586
## 11  3.16227766 0.17750 0.04241004
## 12  5.62341325 0.17375 0.04466309
## 13 10.00000000 0.17375 0.04427267

Tuning shows that optimal \(\tt{cost} = 0.0178\).

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

svm.linear <- svm(Purchase ~ ., kernel = "linear", data = OJ.train, cost = tune.out$best.parameters$cost)
train.pred <- predict(svm.linear, OJ.train)
table(OJ.train$Purchase, train.pred)
##     train.pred
##       CH  MM
##   CH 419  55
##   MM  71 255
(61 + 72)/(428 + 61 + 72 + 239)
## [1] 0.16625
test.pred <- predict(svm.linear, OJ.test)
table(OJ.test$Purchase, test.pred)
##     test.pred
##       CH  MM
##   CH 156  23
##   MM  26  65
(20 + 24)/(144 + 20 + 24 + 82)
## [1] 0.162963

The error rate for our training set is \(16.6\)% and test error rate is \(16.3\)%.

(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 = OJ.train, kernel = "radial")
summary(svm.radial)
## 
## Call:
## svm(formula = Purchase ~ ., data = OJ.train, kernel = "radial")
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  radial 
##        cost:  1 
## 
## Number of Support Vectors:  375
## 
##  ( 186 189 )
## 
## 
## Number of Classes:  2 
## 
## Levels: 
##  CH MM
train.pred <- predict(svm.radial, OJ.train)
table(OJ.train$Purchase, train.pred)
##     train.pred
##       CH  MM
##   CH 427  47
##   MM  73 253
(76 + 46)/(443 + 76 + 46 + 235)
## [1] 0.1525
test.pred <- predict(svm.radial, OJ.test)
table(OJ.test$Purchase, test.pred)
##     test.pred
##       CH  MM
##   CH 158  21
##   MM  25  66
(25 + 17)/(147 + 25 + 17 + 81)
## [1] 0.1555556

Number of Support Vectors: \(377\) \(185\) belong to level \(\tt{CH}\) and the remaining \(192\) belong to level \(\tt{MM}\) The error rate for our training set is \(15.2\)% and test error rate is \(15.6\)%. This model preformed better than the linear kernel.

tune.out <- tune(svm, Purchase ~ ., data = OJ.train, kernel = "radial", ranges = list(cost = 10^seq(-2,
    1, by = 0.25)))
summary(tune.out)
## 
## Parameter tuning of 'svm':
## 
## - sampling method: 10-fold cross validation 
## 
## - best parameters:
##       cost
##  0.1778279
## 
## - best performance: 0.16125 
## 
## - Detailed performance results:
##           cost   error dispersion
## 1   0.01000000 0.40750 0.06072479
## 2   0.01778279 0.40750 0.06072479
## 3   0.03162278 0.30250 0.07425556
## 4   0.05623413 0.19375 0.04218428
## 5   0.10000000 0.17500 0.02500000
## 6   0.17782794 0.16125 0.02239947
## 7   0.31622777 0.16625 0.02638523
## 8   0.56234133 0.16625 0.02045490
## 9   1.00000000 0.17000 0.02443813
## 10  1.77827941 0.16750 0.03073181
## 11  3.16227766 0.16750 0.03238227
## 12  5.62341325 0.18250 0.03641962
## 13 10.00000000 0.18250 0.03961621
svm.radial <- svm(Purchase ~ ., data = OJ.train, kernel = "radial", cost = tune.out$best.parameters$cost)
train.pred <- predict(svm.radial, OJ.train)
table(OJ.train$Purchase, train.pred)
##     train.pred
##       CH  MM
##   CH 427  47
##   MM  72 254
(79 + 50)/(439 + 50 + 79 + 232)
## [1] 0.16125
test.pred <- predict(svm.radial, OJ.test)
table(OJ.test$Purchase, test.pred)
##     test.pred
##       CH  MM
##   CH 158  21
##   MM  29  62
(23 + 16)/(148 + 16 + 23 + 83)
## [1] 0.1444444

The CV tuning slightly increased training error to \(16\)% and slightly decreases test error to \(14.4\)%. Again, this is better than linear kernel.

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

svm.poly <- svm(Purchase ~ ., data = OJ.train, kernel = "poly", degree = 2)
summary(svm.poly)
## 
## Call:
## svm(formula = Purchase ~ ., data = OJ.train, kernel = "poly", degree = 2)
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  polynomial 
##        cost:  1 
##      degree:  2 
##      coef.0:  0 
## 
## Number of Support Vectors:  466
## 
##  ( 229 237 )
## 
## 
## Number of Classes:  2 
## 
## Levels: 
##  CH MM
train.pred <- predict(svm.poly, OJ.train)
table(OJ.train$Purchase, train.pred)
##     train.pred
##       CH  MM
##   CH 441  33
##   MM 113 213
(35 + 121)/(454 + 35 + 121 + 109)
## [1] 0.216968
test.pred <- predict(svm.poly, OJ.test)
table(OJ.test$Purchase, test.pred)
##     test.pred
##       CH  MM
##   CH 166  13
##   MM  32  59
(10 + 34)/(154 + 10 + 34 + 72)
## [1] 0.162963

The polynomial kernel summary show that \(444\) support vectors were produced. Of these, \(218\) belong to level \(\tt{CH}\) and remaining \(226\) belong to level \(\tt{MM}\). This kernel produces a train error of \(21.6\)% and a test error of \(16.3\)% which are slightly higher than the errors produces by radial kernel but lower than the errors produced by linear kernel.

tune.out <- tune(svm, Purchase ~ ., data = OJ.train, kernel = "poly", degree = 2, 
    ranges = list(cost = 10^seq(-2, 1, by = 0.25)))
summary(tune.out)
## 
## Parameter tuning of 'svm':
## 
## - sampling method: 10-fold cross validation 
## 
## - best parameters:
##      cost
##  5.623413
## 
## - best performance: 0.18375 
## 
## - Detailed performance results:
##           cost   error dispersion
## 1   0.01000000 0.40125 0.06832571
## 2   0.01778279 0.38375 0.05653477
## 3   0.03162278 0.37375 0.06136469
## 4   0.05623413 0.35625 0.04832256
## 5   0.10000000 0.34125 0.05138701
## 6   0.17782794 0.26250 0.03280837
## 7   0.31622777 0.22000 0.03184162
## 8   0.56234133 0.21625 0.03488573
## 9   1.00000000 0.20750 0.03184162
## 10  1.77827941 0.20375 0.03175973
## 11  3.16227766 0.18750 0.03004626
## 12  5.62341325 0.18375 0.03064696
## 13 10.00000000 0.18625 0.02389938
svm.poly <- svm(Purchase ~ ., data = OJ.train, kernel = "poly", degree = 2, cost = tune.out$best.parameters$cost)
train.pred <- predict(svm.poly, OJ.train)
table(OJ.train$Purchase, train.pred)
##     train.pred
##       CH  MM
##   CH 432  42
##   MM  81 245
(46 + 82 )/(443 + 46 + 82 + 229)
## [1] 0.16
test.pred <- predict(svm.poly, OJ.test)
table(OJ.test$Purchase, test.pred)
##     test.pred
##       CH  MM
##   CH 161  18
##   MM  29  62
(13 + 28)/(151 + 13 + 28 + 78)
## [1] 0.1518519

Tuning reduced the training error to \(16\)% and test error to \(15.1\)% which is worse than radial kernel but slightly better than linear kernel.

(h) Overall, which approach seems to give the best results on this data? The Radial Basis Kernel seems to produce the minimum misclassification error on both train and test data.