Assignment 8
Matthew Smith
2022-11-17
Applied
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 as follows:
library(tidyverse)set.seed(409)
x1 <- runif(500) - 0.5
x2 <- runif(500) - 0.5
y <- as_factor(1 * (x1^2 - x2^2 > 0))
df <- tibble(x1 = x1, x2 = x2, y = y)head(df)## # A tibble: 6 × 3
## x1 x2 y
## <dbl> <dbl> <fct>
## 1 0.223 0.260 0
## 2 -0.373 -0.245 1
## 3 0.369 -0.153 1
## 4 0.203 0.244 0
## 5 0.213 -0.320 0
## 6 0.181 0.0381 1b) Plot the observations, colored according to their class labels. Your plot should should display \(X_{1}\) in the \(x\)-axis and \(X_{2}\) on the \(y\)-axis.
ggplot(data = df, mapping = aes(x = x1, y = x2, color = y)) +
geom_point() +
theme_minimal()
c) Fit a logistic regression model to the data, using \(X_{1}\) and \(X_{2}\) as predictors.
glm_fit <- glm(y ~ x1 + x2, family = binomial)
summary(glm_fit)##
## Call:
## glm(formula = y ~ x1 + x2, family = binomial)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.330 -1.176 1.031 1.151 1.314
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.007496 0.089825 0.083 0.9335
## x1 -0.527708 0.314042 -1.680 0.0929 .
## x2 -0.235026 0.313903 -0.749 0.4540
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 693.14 on 499 degrees of freedom
## Residual deviance: 689.59 on 497 degrees of freedom
## AIC: 695.59
##
## Number of Fisher Scoring iterations: 3According to the summary of the logistic regression model, neither \(X_{1}\) or \(X_{2}\) are statistically significant for predicting the response variable, \(y\).
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.
glm_prob <- predict(glm_fit, df, type = "response")
df <- mutate(df, pred = as_factor(ifelse(glm_prob > 0.50, 1, 0)))
ggplot(data = df, mapping = aes(x = x1, y = x2, color = pred)) +
geom_point() +
theme_minimal()
The decision boundary appears to be linear.
e) Now fit a logistic regression model to the data using non-linear functions of \(X_{1}\) and \(X_{2}\) as predictors (e.g. \(X_{1}^2\), \(X_{1}\) x \(X_{2}\). \(log(X_{2})\), and so forth).
I am using the square of \(X_{1}\) and \(X_{2}\) as well as the interaction term.
glm_fit2 <- glm(y ~ x1 + x2 + I(x1 ^ 2) + I(x2 ^ 2) + I(x1 * x2), data = df, family = binomial)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 the predicted class labels are obviously non-linear.
glm_prob2 <- predict(glm_fit2, df, type = "response")
df2 <- mutate(df, nonlin_pred = as_factor(ifelse(glm_prob2 > 0.5, 1, 0)))
ggplot(data = df2, mapping = aes(x = x1, y = x2, color = nonlin_pred)) +
geom_point() +
theme_minimal()
The decision boundary is obviously non-linear and is similar to the true decision boundary seen in part b).
g) Fit a support vector classifier to the data with \(X_{1}\) and \(X_{2}\) 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(y ~ x1 + x2, df, kernel = "linear", cost = 0.1)
df3 <- mutate(df, svm_pred = as_factor(predict(svm_fit, df)))
ggplot(data = df3, mapping = aes(x = x1, y = x2, color = svm_pred)) +
geom_point() +
theme_minimal()
The SVM with a linear kernel appears to split the points right down the middle, classifying roughly half of the points as 1 and the other half as 0.
h) Fit an 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_fit2 <- svm(y ~ x1 + x2, df, gamma = 1)
df4 <- mutate(df, svm_pred2 = as_factor(predict(svm_fit2, df)))
ggplot(data = df4, mapping = aes(x = x1, y = x2, color = svm_pred2)) +
geom_point() +
theme_minimal()
The SVM with a non-linear kernel does a much better job than the SVM with a linear kernel. It resembles the true decision boundary seen in part b).
i) Comment on your results.
Both linear models, the base logistic regression and the SVM with a linear kernel, did a poor job of predicting the correct class. This shouldn’t come as a surprise as the decision boundary for the two classes is quite non-linear. It is interesting to see the different ways the two linear models tried to interpret the non-linearity of the data. The two non-linear models on the other hand, the logistic regression with squared and interaction terms and the SVM with a non-linear kernel, did a great job as seen by how closely the plots resemble the actual decision boundary in part b). The SVM with a non-linear kernel appears to be the better option, since manually tweaking the various non-linear and interaction terms can be time consuming, especially as the number of predictors increases.
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
library(ISLR2)
auto <- tibble(Auto) %>%
mutate(mpglevel = as_factor(ifelse(mpg > median(mpg), 1, 0)))
head(auto)## # A tibble: 6 × 10
## mpg cylinders displacement horse…¹ weight accel…² year origin name mpgle…³
## <dbl> <int> <dbl> <int> <int> <dbl> <int> <int> <fct> <fct>
## 1 18 8 307 130 3504 12 70 1 chev… 0
## 2 15 8 350 165 3693 11.5 70 1 buic… 0
## 3 18 8 318 150 3436 11 70 1 plym… 0
## 4 16 8 304 150 3433 12 70 1 amc … 0
## 5 17 8 302 140 3449 10.5 70 1 ford… 0
## 6 15 8 429 198 4341 10 70 1 ford… 0
## # … with abbreviated variable names ¹horsepower, ²acceleration, ³mpglevelb) 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(72)
tune_out <- tune(svm, mpglevel ~. -mpg, 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
## 0.01
##
## - best performance: 0.08935897
##
## - Detailed performance results:
## cost error dispersion
## 1 1e-02 0.08935897 0.04413022
## 2 1e-01 0.09192308 0.04407224
## 3 1e+00 0.09192308 0.04238230
## 4 5e+00 0.10730769 0.04671255
## 5 1e+01 0.09455128 0.04389999
## 6 1e+02 0.10730769 0.05266580svm_linear <- tune_out$best.modelThe cross-validation error is minimized for
cost = 0.01.
c) Now repeat b), but this time using SVMs with radial and
polynomial basis kernels, with different values of gamma
and degree and cost. Comment on your
results.
set.seed(72)
tune_out2 <- tune(svm, mpglevel ~. -mpg, data = auto, kernel = "polynomial",
ranges = list(cost = c(0.1, 1, 5, 10), degree = c(2, 3, 4)))
summary(tune_out2)##
## Parameter tuning of 'svm':
##
## - sampling method: 10-fold cross validation
##
## - best parameters:
## cost degree
## 10 2
##
## - best performance: 0.5054487
##
## - Detailed performance results:
## cost degree error dispersion
## 1 0.1 2 0.5360897 0.05818879
## 2 1.0 2 0.5360897 0.05818879
## 3 5.0 2 0.5360897 0.05818879
## 4 10.0 2 0.5054487 0.09351422
## 5 0.1 3 0.5385897 0.05188816
## 6 1.0 3 0.5385897 0.05188816
## 7 5.0 3 0.5385897 0.05188816
## 8 10.0 3 0.5385897 0.05188816
## 9 0.1 4 0.5435897 0.04099113
## 10 1.0 4 0.5435897 0.04099113
## 11 5.0 4 0.5435897 0.04099113
## 12 10.0 4 0.5435897 0.04099113svm_poly <- tune_out2$best.modelThe lowest cross-validation error for the polynomial kernel is
obtained when cost = 10 and degree = 2.
set.seed(72)
tune_out3 <- tune(svm, mpglevel ~. -mpg, 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_out3)##
## Parameter tuning of 'svm':
##
## - sampling method: 10-fold cross validation
##
## - best parameters:
## cost gamma
## 5 0.1
##
## - best performance: 0.07410256
##
## - Detailed performance results:
## cost gamma error dispersion
## 1 0.1 1e-02 0.11493590 0.06556038
## 2 1.0 1e-02 0.08935897 0.04413022
## 3 5.0 1e-02 0.08942308 0.04581353
## 4 10.0 1e-02 0.08942308 0.05587127
## 5 0.1 1e-01 0.08942308 0.04738124
## 6 1.0 1e-01 0.08679487 0.05165746
## 7 5.0 1e-01 0.07410256 0.05055413
## 8 10.0 1e-01 0.07923077 0.04604708
## 9 0.1 1e+00 0.50608974 0.14648749
## 10 1.0 1e+00 0.08173077 0.05389129
## 11 5.0 1e+00 0.08673077 0.04708817
## 12 10.0 1e+00 0.08673077 0.04708817
## 13 0.1 5e+00 0.54108974 0.04608499
## 14 1.0 5e+00 0.46179487 0.06245986
## 15 5.0 5e+00 0.45673077 0.04672139
## 16 10.0 5e+00 0.45673077 0.04672139
## 17 0.1 1e+01 0.54108974 0.04608499
## 18 1.0 1e+01 0.48737179 0.06299764
## 19 5.0 1e+01 0.48480769 0.06236539
## 20 10.0 1e+01 0.48480769 0.06236539
## 21 0.1 1e+02 0.54608974 0.03690147
## 22 1.0 1e+02 0.54608974 0.03690147
## 23 5.0 1e+02 0.54608974 0.03690147
## 24 10.0 1e+02 0.54608974 0.03690147svm_radial <- tune_out3$best.modelThe lowest cross-validation error for the radial basis kernel is
obtained when cost = 5 and gamme = 0.1.
d) Make some plots to back up your assertions in b) and c).
svm_linear <- svm(mpglevel ~., data = auto, kernel = "linear", cost = 0.01)
svm_poly <- svm(mpglevel ~., data = auto, kernel = "polynomial", cost = 10, degree = 2)
svm_radial <- svm(mpglevel ~., data = auto, kernel = "radial", cost = 5, gamma = 0.1)plot_pairs <- function(fit){
for(name in names(auto)[!(names(auto)) %in% c("mpg", "name", "mpglevel")]){
plot(fit, auto, as.formula(paste("mpg~", name, sep = " ")))
}
}plot_pairs(svm_linear)
plot_pairs(svm_poly)
plot_pairs(svm_radial)
8. This problem involves the OJ data set which is part
of the ISLR2 package.
a) Create a training set containing a random sample of 800 observations, and a test set containing the remaining observations.
set.seed(2323)
in_train <- sample(nrow(OJ), 800)
train_oj <- OJ[in_train, ]
test_oj <- OJ[-in_train, ]b) Fit a support vector classifier to the training data set
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_oj_lin <- svm(Purchase ~., kernel = "linear", data = train_oj, cost = 0.01)
summary(svm_oj_lin)##
## Call:
## svm(formula = Purchase ~ ., data = train_oj, kernel = "linear", cost = 0.01)
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: linear
## cost: 0.01
##
## Number of Support Vectors: 423
##
## ( 211 212 )
##
##
## Number of Classes: 2
##
## Levels:
## CH MMThe support vector classifier created 423 support vectors out of the
800 data points in the training set. There is a nearly even split with
211 belonging to the level CH and the other 212 belonging
to the level MM.
c) What are the training and test error rates?
train_oj2 <- train_oj %>%
as_tibble() %>%
mutate(train_pred = predict(svm_oj_lin, newdata = .)) %>%
summarise('Train Error Rate' = mean(Purchase != train_pred))
train_oj2## # A tibble: 1 × 1
## `Train Error Rate`
## <dbl>
## 1 0.16test_oj2 <- test_oj %>%
as_tibble() %>%
mutate(test_pred = predict(svm_oj_lin, newdata = .)) %>%
summarise('Test Error Rate' = mean(Purchase != test_pred))
test_oj2## # A tibble: 1 × 1
## `Test Error Rate`
## <dbl>
## 1 0.2The training error rate is 16% and the test error rate is 20%.
d) Use the tune() function to select an optimal
cost. Consider values in the range 0.01 to 10.
set.seed(2323)
tune_svm <- tune(svm, Purchase ~., data = train_oj, kernel = "linear",
ranges = list(cost = 10^seq(-2, 1, by = 0.25)))
summary(tune_svm)##
## Parameter tuning of 'svm':
##
## - sampling method: 10-fold cross validation
##
## - best parameters:
## cost
## 0.03162278
##
## - best performance: 0.15875
##
## - Detailed performance results:
## cost error dispersion
## 1 0.01000000 0.16250 0.03385016
## 2 0.01778279 0.16250 0.03726780
## 3 0.03162278 0.15875 0.03775377
## 4 0.05623413 0.16125 0.03701070
## 5 0.10000000 0.16250 0.03679900
## 6 0.17782794 0.16000 0.03855011
## 7 0.31622777 0.15875 0.04126894
## 8 0.56234133 0.16125 0.03928617
## 9 1.00000000 0.16125 0.03972562
## 10 1.77827941 0.16375 0.04016027
## 11 3.16227766 0.16250 0.04082483
## 12 5.62341325 0.16500 0.03855011
## 13 10.00000000 0.16000 0.03987829best_tune <- tune_svm$best.parameters$costThe optimal cost is 0.03162278 according to the tuning
done above.
e) Compute the training and test error rates using the new
value for cost.
lin_svm <- svm(Purchase ~., kernel = "linear", data = train_oj, cost = best_tune)
train_oj3 <- train_oj %>%
as_tibble() %>%
mutate(train_pred3 = predict(lin_svm, newdata = .)) %>%
summarise('Train Error Rate' = mean(Purchase != train_pred3))
train_oj3## # A tibble: 1 × 1
## `Train Error Rate`
## <dbl>
## 1 0.156test_oj3 <- test_oj %>%
as_tibble() %>%
mutate(test_pred3 = predict(lin_svm, newdata = .)) %>%
summarise('Test Error Rate' = mean(Purchase != test_pred3))
test_oj3## # A tibble: 1 × 1
## `Test Error Rate`
## <dbl>
## 1 0.2The training error decrease to 15.625%, while the test error rate stays the same at 20%.
f) Repeat parts b) through e) using a support vector machine
with radial kernel. Use the default value for
gamma.
set.seed(2323)
svm_oj_rad <- svm(Purchase ~., kernel = "radial", data = train_oj)
summary(svm_oj_rad)##
## Call:
## svm(formula = Purchase ~ ., data = train_oj, kernel = "radial")
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: radial
## cost: 1
##
## Number of Support Vectors: 361
##
## ( 180 181 )
##
##
## Number of Classes: 2
##
## Levels:
## CH MMThe support vector classifier created 361 support vectors out of the
800 data points in the training set. There is a nearly even split with
180 belonging to the level CH and the other 181 belonging
to the level MM.
train_oj4 <- train_oj %>%
as_tibble() %>%
mutate(train_pred4 = predict(svm_oj_rad, newdata = .)) %>%
summarise('Train Error Rate' = mean(Purchase != train_pred4))
train_oj4## # A tibble: 1 × 1
## `Train Error Rate`
## <dbl>
## 1 0.145test_oj4 <- test_oj %>%
as_tibble() %>%
mutate(test_pred4 = predict(svm_oj_rad, newdata = .)) %>%
summarise('Test Error Rate' = mean(Purchase != test_pred4))
test_oj4## # A tibble: 1 × 1
## `Test Error Rate`
## <dbl>
## 1 0.215The training error decrease to 14.5%, while the test error rate increases to 21.48148%. This is not an improvement over the SVM with a linear kernel.
g) Repeat parts b) through e) using a support vector machine
with a polynomial kernel. Set degree = 2.
set.seed(2323)
svm_oj_poly <- svm(Purchase ~., kernel = "polynomial", data = train_oj)
summary(svm_oj_poly)##
## Call:
## svm(formula = Purchase ~ ., data = train_oj, kernel = "polynomial")
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: polynomial
## cost: 1
## degree: 3
## coef.0: 0
##
## Number of Support Vectors: 417
##
## ( 211 206 )
##
##
## Number of Classes: 2
##
## Levels:
## CH MMThe support vector classifier created 417 support vectors out of the
800 data points in the training set. The split is 211 belonging to the
level CH and the other 206 belonging to the level
MM.
train_oj5 <- train_oj %>%
as_tibble() %>%
mutate(train_pred5 = predict(svm_oj_poly, newdata = .)) %>%
summarise('Train Error Rate' = mean(Purchase != train_pred5))
train_oj5## # A tibble: 1 × 1
## `Train Error Rate`
## <dbl>
## 1 0.155test_oj5 <- test_oj %>%
as_tibble() %>%
mutate(test_pred5 = predict(svm_oj_poly, newdata = .)) %>%
summarise('Test Error Rate' = mean(Purchase != test_pred5))
test_oj5## # A tibble: 1 × 1
## `Test Error Rate`
## <dbl>
## 1 0.215The training error rate is 15.5%, a slight increase over the SVM with a radial basis kernel. The test error rate is 21.48148, which is the same as the radial basis kernel SVM.
h) Overall, which approach seems to give the best results on this data?
error_results <- tibble("SVM Kernel" = c("Linear", "Radial", "Polynomial"),
"Training Error" = c(train_oj3[[1]], train_oj4[[1]], train_oj5[[1]]),
"Test Error" = c(test_oj3[[1]], test_oj4[[1]], test_oj5[[1]]))
error_results## # A tibble: 3 × 3
## `SVM Kernel` `Training Error` `Test Error`
## <chr> <dbl> <dbl>
## 1 Linear 0.156 0.2
## 2 Radial 0.145 0.215
## 3 Polynomial 0.155 0.215Overall, the approach that seems to give the best results in this data is the SVM with a linear basis kernel. This is based on the fact that the linear basis kernel produces the lowest missclassification error rate on the test data set.