library(caret)
## Warning: package 'caret' was built under R version 4.3.3
## Loading required package: ggplot2
## Loading required package: lattice
## Warning: package 'lattice' was built under R version 4.3.3
library(ggplot2)
library(kernlab)
## Warning: package 'kernlab' was built under R version 4.3.3
##
## Attaching package: 'kernlab'
## The following object is masked from 'package:ggplot2':
##
## alpha
library(tidyverse)
## Warning: package 'tidyverse' was built under R version 4.3.3
## Warning: package 'tibble' was built under R version 4.3.3
## Warning: package 'tidyr' was built under R version 4.3.3
## Warning: package 'readr' was built under R version 4.3.3
## Warning: package 'dplyr' was built under R version 4.3.3
## Warning: package 'lubridate' was built under R version 4.3.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.1 ✔ stringr 1.5.1
## ✔ lubridate 1.9.4 ✔ tibble 3.2.1
## ✔ purrr 1.0.2 ✔ tidyr 1.3.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ kernlab::alpha() masks ggplot2::alpha()
## ✖ purrr::cross() masks kernlab::cross()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ purrr::lift() masks caret::lift()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(e1071)
## Warning: package 'e1071' was built under R version 4.3.3
##
## Attaching package: 'e1071'
##
## The following object is masked from 'package:ggplot2':
##
## element
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.
n <- 500
x1 <- runif(n) - 0.5
x2 <- runif(n) - 0.5
y <- as.factor(ifelse(x1^2 - x2^2 > 0, "Class1", "Class0"))
data <- data.frame(x1 = x1, x2 = x2, y = y)
ggplot(data, aes(x = x1, y = x2, color = y)) +
geom_point(alpha = 0.8, size = 2) +
theme_minimal()
log_model <- glm(y ~ x1 + x2, data = data, family = "binomial")
summary(log_model)
##
## Call:
## glm(formula = y ~ x1 + x2, family = "binomial", data = data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.05321 0.09035 0.589 0.556
## x1 -0.04597 0.30935 -0.149 0.882
## x2 -0.22692 0.31437 -0.722 0.470
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 692.86 on 499 degrees of freedom
## Residual deviance: 692.31 on 497 degrees of freedom
## AIC: 698.31
##
## Number of Fisher Scoring iterations: 3
pred.log <- predict(log_model, type = "response")
pred.glm <- as.factor(as.integer(pred.log > 0.5))
data$pred_linear <- pred.glm
ggplot(data, aes(x = x1, y = x2, color = pred_linear)) +
geom_point(size = 1.5) +
theme_minimal()
nl.glm <- 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
summary(nl.glm)
##
## Call:
## glm(formula = y ~ poly(x1, 2) + poly(x2, 2) + I(x1 * x2), family = binomial,
## data = data)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 489.8 2933.5 0.167 0.867
## poly(x1, 2)1 -39701.5 221643.6 -0.179 0.858
## poly(x1, 2)2 284717.8 1609246.8 0.177 0.860
## poly(x2, 2)1 -31953.1 186391.4 -0.171 0.864
## poly(x2, 2)2 -277457.4 1567379.9 -0.177 0.859
## I(x1 * x2) -1408.2 14496.7 -0.097 0.923
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 6.9286e+02 on 499 degrees of freedom
## Residual deviance: 2.4065e-04 on 494 degrees of freedom
## AIC: 12
##
## Number of Fisher Scoring iterations: 25
glm.probs.nl <- predict(nl.glm, type = "response")
glm.pred.nl <- as.factor((glm.probs.nl > 0.5))
data$pred_nonlinear <- glm.pred.nl
ggplot(data, aes(x = x1, y = x2, color = pred_nonlinear)) +
geom_point(size = 1.5) +
labs(title = "Logistic regression",
color = "Predicted class") +
theme_minimal()
svm.linear.fit <- train(
y ~ x1 + x2,
data = data,
method = "svmLinear",
tuneGrid = data.frame(C = 1)
)
svm.linear.pred <- predict(svm.linear.fit, newdata = data)
data$pred_svm_linear <- svm.linear.pred
ggplot(data, aes(x = x1, y = x2, color = pred_svm_linear)) +
geom_point(size = 1.5) +
labs(title = "Support vector classifier",
color = "Predicted class") +
theme_minimal()
svm.radial.fit <- train(
y ~ x1 + x2,
data = data,
method = "svmRadial",
tuneGrid = data.frame(C = 1, sigma = 1)
)
svm.radial.pred <- predict(svm.radial.fit, newdata = data)
data$pred_svm_radial <- svm.radial.pred
ggplot(data, aes(x = x1, y = x2, color = pred_svm_radial)) +
geom_point(size = 1.5) +
labs(title = "SVM predicted",
color = "Predicted class") +
theme_minimal()
mean(svm.radial.pred == data$y)
## [1] 0.97
Without knowing the true shape of the boundary, SVM non linear would be the best choice since it is a more robust choice than guessing polynomial terms.
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.
library(ISLR2)
## Warning: package 'ISLR2' was built under R version 4.3.3
data(Auto)
head(Auto)
## mpg cylinders displacement horsepower weight acceleration year origin
## 1 18 8 307 130 3504 12.0 70 1
## 2 15 8 350 165 3693 11.5 70 1
## 3 18 8 318 150 3436 11.0 70 1
## 4 16 8 304 150 3433 12.0 70 1
## 5 17 8 302 140 3449 10.5 70 1
## 6 15 8 429 198 4341 10.0 70 1
## name
## 1 chevrolet chevelle malibu
## 2 buick skylark 320
## 3 plymouth satellite
## 4 amc rebel sst
## 5 ford torino
## 6 ford galaxie 500
Auto$mpg01 <- as.factor(ifelse(Auto$mpg > median(Auto$mpg), 1, 0))
Auto_svm <- subset(Auto, select = -c(mpg, name))
set.seed(42)
tune_linear <- tune(svm, mpg01 ~ ., data = Auto_svm, kernel = "linear",
ranges = list(cost = c(0.001, 0.01, 0.1, 1, 5, 10, 100)))
summary(tune_linear)
##
## Parameter tuning of 'svm':
##
## - sampling method: 10-fold cross validation
##
## - best parameters:
## cost
## 5
##
## - best performance: 0.08653846
##
## - Detailed performance results:
## cost error dispersion
## 1 1e-03 0.13288462 0.07351673
## 2 1e-02 0.08916667 0.05258186
## 3 1e-01 0.09673077 0.05699840
## 4 1e+00 0.09423077 0.04632467
## 5 5e+00 0.08653846 0.03776796
## 6 1e+01 0.08653846 0.03776796
## 7 1e+02 0.08653846 0.03776796
best_model <- tune_linear$best.model
summary(best_model)
##
## Call:
## best.tune(METHOD = svm, train.x = mpg01 ~ ., data = Auto_svm, ranges = list(cost = c(0.001,
## 0.01, 0.1, 1, 5, 10, 100)), kernel = "linear")
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: linear
## cost: 5
##
## Number of Support Vectors: 83
##
## ( 41 42 )
##
##
## Number of Classes: 2
##
## Levels:
## 0 1
set.seed(42)
tune_radial <- tune(svm, mpg01 ~ ., data = Auto_svm, kernel = "radial",
ranges = list(cost = c(0.1, 1, 5, 10, 100),
gamma = c(0.001, 0.01, 0.1, 1, 5)))
summary(tune_radial)
##
## Parameter tuning of 'svm':
##
## - sampling method: 10-fold cross validation
##
## - best parameters:
## cost gamma
## 1 1
##
## - best performance: 0.07878205
##
## - Detailed performance results:
## cost gamma error dispersion
## 1 0.1 0.001 0.59679487 0.05312225
## 2 1.0 0.001 0.11750000 0.05956907
## 3 5.0 0.001 0.08916667 0.05258186
## 4 10.0 0.001 0.09173077 0.05525432
## 5 100.0 0.001 0.08653846 0.05996931
## 6 0.1 0.010 0.11237179 0.05439623
## 7 1.0 0.010 0.08660256 0.05519479
## 8 5.0 0.010 0.09416667 0.05728285
## 9 10.0 0.010 0.08660256 0.06001684
## 10 100.0 0.010 0.08647436 0.05341235
## 11 0.1 0.100 0.09173077 0.05254364
## 12 1.0 0.100 0.08916667 0.05258186
## 13 5.0 0.100 0.08897436 0.04597723
## 14 10.0 0.100 0.09653846 0.05358794
## 15 100.0 0.100 0.08891026 0.04725465
## 16 0.1 1.000 0.09673077 0.05699840
## 17 1.0 1.000 0.07878205 0.04472958
## 18 5.0 1.000 0.08634615 0.04391746
## 19 10.0 1.000 0.09403846 0.04383004
## 20 100.0 1.000 0.10692308 0.04869339
## 21 0.1 5.000 0.59679487 0.05312225
## 22 1.0 5.000 0.09653846 0.05220694
## 23 5.0 5.000 0.11185897 0.05276112
## 24 10.0 5.000 0.10679487 0.05240247
## 25 100.0 5.000 0.10679487 0.05240247
bestmodel.radial <- tune_radial$best.model
summary(bestmodel.radial)
##
## Call:
## best.tune(METHOD = svm, train.x = mpg01 ~ ., data = Auto_svm, ranges = list(cost = c(0.1,
## 1, 5, 10, 100), gamma = c(0.001, 0.01, 0.1, 1, 5)), kernel = "radial")
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: radial
## cost: 1
##
## Number of Support Vectors: 184
##
## ( 92 92 )
##
##
## Number of Classes: 2
##
## Levels:
## 0 1
Poly
set.seed(42)
tune_poly <- tune(svm, mpg01 ~ ., data = Auto_svm, kernel = "polynomial",
ranges = list(cost = c(0.1, 1, 5, 10, 100),
degree = c(2, 3, 4)))
summary(tune_poly)
##
## Parameter tuning of 'svm':
##
## - sampling method: 10-fold cross validation
##
## - best parameters:
## cost degree
## 10 3
##
## - best performance: 0.07621795
##
## - Detailed performance results:
## cost degree error dispersion
## 1 0.1 2 0.28108974 0.06318810
## 2 1.0 2 0.27025641 0.09718216
## 3 5.0 2 0.18570513 0.07190928
## 4 10.0 2 0.19083333 0.07899168
## 5 100.0 2 0.18326923 0.06521573
## 6 0.1 3 0.19429487 0.11454046
## 7 1.0 3 0.09679487 0.05172882
## 8 5.0 3 0.08141026 0.05172618
## 9 10.0 3 0.07621795 0.05015358
## 10 100.0 3 0.09416667 0.03322969
## 11 0.1 4 0.27339744 0.08041984
## 12 1.0 4 0.26493590 0.10509261
## 13 5.0 4 0.19365385 0.04717381
## 14 10.0 4 0.17076923 0.04911761
## 15 100.0 4 0.14006410 0.05044478
bestmodel.poly <- tune_poly$best.model
summary(bestmodel.poly)
##
## Call:
## best.tune(METHOD = svm, train.x = mpg01 ~ ., data = Auto_svm, ranges = list(cost = c(0.1,
## 1, 5, 10, 100), degree = c(2, 3, 4)), kernel = "polynomial")
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: polynomial
## cost: 10
## degree: 3
## coef.0: 0
##
## Number of Support Vectors: 99
##
## ( 49 50 )
##
##
## Number of Classes: 2
##
## Levels:
## 0 1
svm_linear_final <- svm(mpg01 ~ ., data = Auto_svm, kernel = "linear",
cost = tune_linear$best.parameters$cost)
svm_radial_final <- svm(mpg01 ~ ., data = Auto_svm, kernel = "radial",
cost = tune_radial$best.parameters$cost,
gamma = tune_radial$best.parameters$gamma)
svm_poly_final <- svm(mpg01 ~ ., data = Auto_svm, kernel = "polynomial",
cost = tune_poly$best.parameters$cost,
degree = tune_poly$best.parameters$degree)
plot(svm_linear_final, Auto_svm, weight ~ horsepower)
plot(svm_radial_final, Auto_svm, weight ~ horsepower)
plot(svm_poly_final, Auto_svm, weight ~ horsepower)
plot(svm_linear_final, Auto_svm, displacement ~ weight)
plot(svm_radial_final, Auto_svm, displacement ~ weight)
This problem involves the OJ data set which is part of the ISLP package.
data(OJ)
head(OJ)
## Purchase WeekofPurchase StoreID PriceCH PriceMM DiscCH DiscMM SpecialCH
## 1 CH 237 1 1.75 1.99 0.00 0.0 0
## 2 CH 239 1 1.75 1.99 0.00 0.3 0
## 3 CH 245 1 1.86 2.09 0.17 0.0 0
## 4 MM 227 1 1.69 1.69 0.00 0.0 0
## 5 CH 228 7 1.69 1.69 0.00 0.0 0
## 6 CH 230 7 1.69 1.99 0.00 0.0 0
## SpecialMM LoyalCH SalePriceMM SalePriceCH PriceDiff Store7 PctDiscMM
## 1 0 0.500000 1.99 1.75 0.24 No 0.000000
## 2 1 0.600000 1.69 1.75 -0.06 No 0.150754
## 3 0 0.680000 2.09 1.69 0.40 No 0.000000
## 4 0 0.400000 1.69 1.69 0.00 No 0.000000
## 5 0 0.956535 1.69 1.69 0.00 Yes 0.000000
## 6 1 0.965228 1.99 1.69 0.30 Yes 0.000000
## PctDiscCH ListPriceDiff STORE
## 1 0.000000 0.24 1
## 2 0.000000 0.24 1
## 3 0.091398 0.23 1
## 4 0.000000 0.00 1
## 5 0.000000 0.00 0
## 6 0.000000 0.30 0
train <- sample(1:nrow(OJ), 800)
OJ.train <- OJ[train, ]
OJ.test <- OJ[-train, ]
set.seed(42)
svm.linear.oj <- svm(Purchase ~ ., data = OJ.train, kernel = "linear",
cost = 0.01)
summary(svm.linear.oj)
##
## 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: 438
##
## ( 218 220 )
##
##
## Number of Classes: 2
##
## Levels:
## CH MM
set.seed(42)
train.pred <- predict(svm.linear.oj, OJ.train)
train.err <- mean(train.pred != OJ.train$Purchase)
train.err
## [1] 0.16875
test.pred <- predict(svm.linear.oj, OJ.test)
test.err <- mean(test.pred != OJ.test$Purchase)
test.err
## [1] 0.1703704
set.seed(42)
tune.linear.oj <- tune(svm, Purchase ~ ., data = OJ.train, kernel = "linear",
ranges = list(cost = c(0.01, 0.1, 1, 5, 10)))
summary(tune.linear.oj)
##
## Parameter tuning of 'svm':
##
## - sampling method: 10-fold cross validation
##
## - best parameters:
## cost
## 5
##
## - best performance: 0.17125
##
## - Detailed performance results:
## cost error dispersion
## 1 0.01 0.17875 0.03438447
## 2 0.10 0.17375 0.04226652
## 3 1.00 0.17250 0.04556741
## 4 5.00 0.17125 0.04210189
## 5 10.00 0.17375 0.03793727
oj.linear <- tune.linear.oj$best.model
summary(oj.linear)
##
## Call:
## best.tune(METHOD = svm, train.x = Purchase ~ ., data = OJ.train,
## ranges = list(cost = c(0.01, 0.1, 1, 5, 10)), kernel = "linear")
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: linear
## cost: 5
##
## Number of Support Vectors: 335
##
## ( 169 166 )
##
##
## Number of Classes: 2
##
## Levels:
## CH MM
set.seed(42)
train.pred1 <- predict(oj.linear, OJ.train)
train.err.opt <- mean(train.pred1 != OJ.train$Purchase)
train.err.opt
## [1] 0.16125
test.pred1 <- predict(oj.linear, OJ.test)
test.err.opt <- mean(test.pred1 != OJ.test$Purchase)
test.err.opt
## [1] 0.162963
svm.radial.oj <- svm(Purchase ~ ., data = OJ.train, kernel = "radial", cost = 0.01)
summary(svm.radial.oj)
##
## Call:
## svm(formula = Purchase ~ ., data = OJ.train, kernel = "radial", cost = 0.01)
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: radial
## cost: 0.01
##
## Number of Support Vectors: 618
##
## ( 308 310 )
##
##
## Number of Classes: 2
##
## Levels:
## CH MM
train.pred2 <- predict(svm.radial.oj, OJ.train)
mean(train.pred2 != OJ.train$Purchase)
## [1] 0.385
test.pred2 <- predict(svm.radial.oj, OJ.test)
mean(test.pred2 != OJ.test$Purchase)
## [1] 0.4037037
set.seed(42)
tune.radial.oj <- tune(svm, Purchase ~ ., data = OJ.train, kernel = "radial",
ranges = list(cost = c(0.01, 0.1, 1, 5, 10)))
summary(tune.radial.oj)
##
## Parameter tuning of 'svm':
##
## - sampling method: 10-fold cross validation
##
## - best parameters:
## cost
## 1
##
## - best performance: 0.1825
##
## - Detailed performance results:
## cost error dispersion
## 1 0.01 0.3850 0.06368324
## 2 0.10 0.1950 0.02581989
## 3 1.00 0.1825 0.04495368
## 4 5.00 0.1850 0.04923018
## 5 10.00 0.1850 0.04706674
oj.radial <- tune.radial.oj$best.model
set.seed(42)
train.pred3 <- predict(oj.radial, OJ.train)
mean(train.pred3 != OJ.train$Purchase) # ~0.14-0.15
## [1] 0.15375
test.pred3 <- predict(oj.radial, OJ.test)
mean(test.pred3 != OJ.test$Purchase)
## [1] 0.1481481
svm.poly <- svm(Purchase ~ ., data = OJ.train, kernel = "polynomial",
degree = 2, cost = 0.01)
summary(svm.poly)
##
## Call:
## svm(formula = Purchase ~ ., data = OJ.train, 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: 622
##
## ( 308 314 )
##
##
## Number of Classes: 2
##
## Levels:
## CH MM
set.seed(42)
train.pred4 <- predict(svm.poly, OJ.train)
mean(train.pred4 != OJ.train$Purchase)
## [1] 0.36375
test.pred4 <- predict(svm.poly, OJ.test)
mean(test.pred4 != OJ.test$Purchase)
## [1] 0.3925926
set.seed(42)
tune.poly <- tune(svm, Purchase ~ ., data = OJ.train, 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:
## cost
## 10
##
## - best performance: 0.17875
##
## - Detailed performance results:
## cost error dispersion
## 1 0.01 0.38500 0.06368324
## 2 0.10 0.32375 0.04980866
## 3 1.00 0.19875 0.03653860
## 4 5.00 0.18500 0.04594683
## 5 10.00 0.17875 0.04966904
tune_poly <- tune.poly$best.model
set.seed(42)
train.pred5 <- predict(tune_poly, OJ.train)
mean(train.pred5 != OJ.train$Purchase)
## [1] 0.15625
test.pred5 <- predict(tune_poly, OJ.test)
mean(test.pred5 != OJ.test$Purchase)
## [1] 0.1777778
linear support vector classifier and the radial-kernel SVM with tuned cost tend to perform about equally well and best overall/