Homework 8 - Chapter 9
Valeria De La Rocha
Problem 5
We can use logistic regression with transformed predictors to create a non-linear decision boundary and compare it with support vector methods.
(a) Generate the data
set.seed(1)
x1 <- runif(500) - 0.5
x2 <- runif(500) - 0.5
y <- factor(1 * (x1^2 - x2^2 > 0))
sim_data <- data.frame(x1, x2, y)
head(sim_data)## x1 x2 y
## 1 -0.23449134 0.05417706 1
## 2 -0.12787610 0.18827524 0
## 3 0.07285336 0.15805755 0
## 4 0.40820779 0.16334273 1
## 5 -0.29831807 -0.02776580 1
## 6 0.39838968 0.46952817 0(b) Plot the observations by actual class
ggplot(sim_data, aes(x = x1, y = x2, color = y)) +
geom_point() +
labs(title = "Actual Class Labels", x = "X1", y = "X2", color = "Class") +
theme_minimal()
The two classes are separated by a curved, quadratic boundary. The pattern looks similar to crossing diagonal regions rather than a straight line.
(c) Fit logistic regression using X1 and X2
logistic_linear <- glm(y ~ x1 + x2,
data = sim_data,
family = binomial)
summary(logistic_linear)##
## Call:
## glm(formula = y ~ x1 + x2, family = binomial, data = sim_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) Predict classes using the linear logistic model
linear_prob <- predict(logistic_linear, type = "response")
sim_data$pred_logistic_linear <- factor(ifelse(linear_prob > 0.5, 1, 0))
mean(sim_data$pred_logistic_linear != sim_data$y)## [1] 0.43ggplot(sim_data, aes(x = x1, y = x2, color = pred_logistic_linear)) +
geom_point() +
labs(title = "Predicted Classes: Linear Logistic Regression",
x = "X1", y = "X2", color = "Predicted") +
theme_minimal()
The linear logistic regression does not capture the true pattern well because it can only create a straight decision boundary from X1 and X2.
(e) Fit logistic regression using non-linear predictors
logistic_nonlinear <- glm(
y ~ x1 + x2 + I(x1^2) + I(x2^2) + I(x1 * x2),
data = sim_data,
family = binomial
)
summary(logistic_nonlinear)##
## Call:
## glm(formula = y ~ x1 + x2 + I(x1^2) + I(x2^2) + I(x1 * x2), family = binomial,
## data = sim_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) Predict classes using the non-linear logistic model
nonlinear_prob <- predict(logistic_nonlinear, type = "response")
sim_data$pred_logistic_nonlinear <- factor(ifelse(nonlinear_prob > 0.5, 1, 0))
mean(sim_data$pred_logistic_nonlinear != sim_data$y)## [1] 0ggplot(sim_data, aes(x = x1, y = x2, color = pred_logistic_nonlinear)) +
geom_point() +
labs(title = "Predicted Classes: Non-Linear Logistic Regression",
x = "X1", y = "X2", color = "Predicted") +
theme_minimal()
The transformed predictors allow logistic regression to produce a curved boundary. This model follows the original class pattern much better than the linear model.
(g) Fit a support vector classifier
svc_fit <- svm(y ~ x1 + x2,
data = sim_data,
kernel = "linear",
cost = 1,
scale = TRUE)
sim_data$pred_svc <- predict(svc_fit, sim_data)
mean(sim_data$pred_svc != sim_data$y)## [1] 0.478ggplot(sim_data, aes(x = x1, y = x2, color = pred_svc)) +
geom_point() +
labs(title = "Predicted Classes: Linear Support Vector Classifier",
x = "X1", y = "X2", color = "Predicted") +
theme_minimal()
The support vector classifier also uses a linear boundary, so it has the same basic problem as the first logistic regression model.
(h) Fit an SVM with a non-linear kernel
svm_radial_fit <- svm(y ~ x1 + x2,
data = sim_data,
kernel = "radial",
cost = 10,
gamma = 1,
scale = TRUE)
sim_data$pred_svm_radial <- predict(svm_radial_fit, sim_data)
mean(sim_data$pred_svm_radial != sim_data$y)## [1] 0.014ggplot(sim_data, aes(x = x1, y = x2, color = pred_svm_radial)) +
geom_point() +
labs(title = "Predicted Classes: Radial SVM",
x = "X1", y = "X2", color = "Predicted") +
theme_minimal()
(i) Comments on the results
The linear logistic regression had a training error rate of 0.430, and the linear support vector classifier had a training error rate of 0.478. These models did not perform well because they can only create a straight decision boundary, while the actual boundary is quadratic.
The logistic regression model with squared terms had a training error rate of 0.000. It correctly captured the quadratic relationship because it included (X_1^2) and (X_2^2). The very large coefficient estimates occurred because the transformed predictors almost perfectly separated the two classes.
The radial SVM also performed well, with a training error rate of 0.014. It created a flexible non-linear boundary without requiring me to manually add squared predictors.
Overall, the non-linear logistic regression performed best on the training data, followed closely by the radial SVM. The linear logistic regression and linear support vector classifier were not appropriate for this data because the relationship between the predictors and the class labels was not linear.
Problem 7
The goal is to predict whether a car has high or low gas mileage
using the Auto data set.
(a) Create the binary mileage variable
data(Auto)
Auto2 <- Auto
Auto2$mpg01 <- factor(ifelse(Auto2$mpg > median(Auto2$mpg), 1, 0))
# Remove mpg and the car name from the predictor set
Auto_svm <- Auto2 %>%
select(-mpg, -name)
table(Auto_svm$mpg01)##
## 0 1
## 196 196Cars above the median MPG are labeled 1, and cars at or below the median are labeled 0. The original MPG variable is removed so that it does not directly give away the answer.
(b) Linear support vector classifier with different cost values
set.seed(1)
linear_tune_auto <- tune(
svm,
mpg01 ~ .,
data = Auto_svm,
kernel = "linear",
ranges = list(cost = c(0.01, 0.1, 1, 10, 100)),
tunecontrol = tune.control(cross = 10)
)
linear_tune_auto$performances## cost error dispersion
## 1 1e-02 0.08923077 0.04698309
## 2 1e-01 0.09185897 0.04393409
## 3 1e+00 0.08435897 0.03662670
## 4 1e+01 0.08948718 0.03898410
## 5 1e+02 0.08692308 0.03887151linear_tune_auto$best.parameters## cost
## 3 1linear_tune_auto$best.performance## [1] 0.08435897I compared costs from 0.01 through 100 using 10-fold cross-validation. I selected the cost with the smallest cross-validation error. A very small cost allows more classification errors, while a large cost penalizes errors more strongly. The table above shows whether increasing the cost improves the result or begins to overfit the data.
(c) Radial and polynomial SVMs
set.seed(1)
radial_tune_auto <- tune(
svm,
mpg01 ~ .,
data = Auto_svm,
kernel = "radial",
ranges = list(
cost = c(0.1, 1, 10),
gamma = c(0.01, 0.1, 1)
),
tunecontrol = tune.control(cross = 10)
)
radial_tune_auto$performances## cost gamma error dispersion
## 1 0.1 0.01 0.11224359 0.03836937
## 2 1.0 0.01 0.08673077 0.04551036
## 3 10.0 0.01 0.08673077 0.04040897
## 4 0.1 0.10 0.08923077 0.04698309
## 5 1.0 0.10 0.08923077 0.04376306
## 6 10.0 0.10 0.08166667 0.04149504
## 7 0.1 1.00 0.08673077 0.04535158
## 8 1.0 1.00 0.06634615 0.03244101
## 9 10.0 1.00 0.08923077 0.02732003radial_tune_auto$best.parameters## cost gamma
## 8 1 1radial_tune_auto$best.performance## [1] 0.06634615set.seed(1)
polynomial_tune_auto <- tune(
svm,
mpg01 ~ .,
data = Auto_svm,
kernel = "polynomial",
ranges = list(
cost = c(0.1, 1, 10),
gamma = c(0.01, 0.1, 1),
degree = c(2, 3)
),
tunecontrol = tune.control(cross = 10)
)
polynomial_tune_auto$performances## cost gamma degree error dispersion
## 1 0.1 0.01 2 0.55115385 0.04366593
## 2 1.0 0.01 2 0.51807692 0.07574405
## 3 10.0 0.01 2 0.31673077 0.09564273
## 4 0.1 0.10 2 0.31673077 0.09564273
## 5 1.0 0.10 2 0.28871795 0.10630818
## 6 10.0 0.10 2 0.17621795 0.04945319
## 7 0.1 1.00 2 0.17621795 0.04945319
## 8 1.0 1.00 2 0.18128205 0.06917126
## 9 10.0 1.00 2 0.17608974 0.05616208
## 10 0.1 0.01 3 0.55115385 0.04366593
## 11 1.0 0.01 3 0.55115385 0.04366593
## 12 10.0 0.01 3 0.29108974 0.09527924
## 13 0.1 0.10 3 0.25794872 0.09147506
## 14 1.0 0.10 3 0.09185897 0.03993389
## 15 10.0 0.10 3 0.08679487 0.04533461
## 16 0.1 1.00 3 0.08179487 0.03391477
## 17 1.0 1.00 3 0.09442308 0.04647330
## 18 10.0 1.00 3 0.12487179 0.04202595polynomial_tune_auto$best.parameters## cost gamma degree
## 16 0.1 1 3polynomial_tune_auto$best.performance## [1] 0.08179487For the radial model, I compared several combinations of cost and gamma. Gamma controls how local or flexible the boundary becomes. For the polynomial model, I also compared degrees 2 and 3. I used the lowest cross-validation error to choose the best settings.
auto_model_comparison <- data.frame(
Model = c("Linear", "Radial", "Polynomial"),
Best_CV_Error = c(
linear_tune_auto$best.performance,
radial_tune_auto$best.performance,
polynomial_tune_auto$best.performance
)
)
auto_model_comparison## Model Best_CV_Error
## 1 Linear 0.08435897
## 2 Radial 0.06634615
## 3 Polynomial 0.08179487The model with the lowest value in the comparison table is the best model based on cross-validation. In general, the radial model can improve the result when the relationship between the car variables and mileage group is not completely linear.
(d) Plots supporting the results
best_linear_auto <- linear_tune_auto$best.model
best_radial_auto <- radial_tune_auto$best.model
best_polynomial_auto <- polynomial_tune_auto$best.model
plot(best_linear_auto, Auto_svm, horsepower ~ weight,
main = "Linear SVC: Horsepower and Weight")
plot(best_radial_auto, Auto_svm, horsepower ~ weight,
main = "Radial SVM: Horsepower and Weight")
plot(best_polynomial_auto, Auto_svm, horsepower ~ weight,
main = "Polynomial SVM: Horsepower and Weight")
plot(best_radial_auto, Auto_svm, displacement ~ year,
main = "Radial SVM: Displacement and Year")
The plots show only two predictors at a time because the complete models use more than two predictors. Weight and horsepower show a strong separation between high- and low-mileage cars. The radial model can create a more flexible boundary than the linear model.
Problem 8
This problem compares linear, radial, and polynomial support vector
models using the OJ data set.
(a) Create the training and test sets
data(OJ)
set.seed(1)
train_index <- sample(seq_len(nrow(OJ)), 800)
OJ_train <- OJ[train_index, ]
OJ_test <- OJ[-train_index, ]
nrow(OJ_train)## [1] 800nrow(OJ_test)## [1] 270table(OJ_train$Purchase)##
## CH MM
## 485 315table(OJ_test$Purchase)##
## CH MM
## 168 102The training set contains 800 observations, and the test set contains the remaining 270 observations.
Helper function for error rates
error_rate <- function(actual, predicted) {
mean(actual != predicted)
}(b) Linear support vector classifier with cost = 0.01
linear_001 <- svm(
Purchase ~ .,
data = OJ_train,
kernel = "linear",
cost = 0.01,
scale = TRUE
)
summary(linear_001)##
## Call:
## svm(formula = Purchase ~ ., data = OJ_train, 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 MMThe summary reports the SVM type, kernel, cost, number of classes, and number of support vectors. The response has two classes, CH and MM. The support vectors are the training observations that help define the decision boundary.
(c) Training and test error rates
linear_001_train_pred <- predict(linear_001, OJ_train)
linear_001_test_pred <- predict(linear_001, OJ_test)
linear_001_train_error <- error_rate(OJ_train$Purchase, linear_001_train_pred)
linear_001_test_error <- error_rate(OJ_test$Purchase, linear_001_test_pred)
linear_001_train_error## [1] 0.175linear_001_test_error## [1] 0.1777778These values are the proportions of observations classified incorrectly in the training and test sets.
(d) Tune the cost for the linear model
set.seed(1)
linear_tune_oj <- tune(
svm,
Purchase ~ .,
data = OJ_train,
kernel = "linear",
ranges = list(cost = c(0.01, 0.1, 1, 5, 10)),
tunecontrol = tune.control(cross = 10)
)
linear_tune_oj$performances## cost error dispersion
## 1 0.01 0.17625 0.02853482
## 2 0.10 0.17250 0.03162278
## 3 1.00 0.17500 0.02946278
## 4 5.00 0.17250 0.03162278
## 5 10.00 0.17375 0.03197764linear_tune_oj$best.parameters## cost
## 2 0.1linear_tune_oj$best.performance## [1] 0.1725I used 10-fold cross-validation and selected the cost with the lowest cross-validation error. This gives a fair comparison because each cost is evaluated on held-out folds from the training data.
(e) Error rates using the selected linear cost
best_linear_oj <- linear_tune_oj$best.model
best_linear_train_pred <- predict(best_linear_oj, OJ_train)
best_linear_test_pred <- predict(best_linear_oj, OJ_test)
best_linear_train_error <- error_rate(OJ_train$Purchase, best_linear_train_pred)
best_linear_test_error <- error_rate(OJ_test$Purchase, best_linear_test_pred)
best_linear_train_error## [1] 0.165best_linear_test_error## [1] 0.162963(f) Radial SVM
First, I fit the radial model with cost = 0.01 and the default gamma.
radial_001 <- svm(
Purchase ~ .,
data = OJ_train,
kernel = "radial",
cost = 0.01,
scale = TRUE
)
summary(radial_001)##
## Call:
## svm(formula = Purchase ~ ., data = OJ_train, kernel = "radial", cost = 0.01,
## scale = TRUE)
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: radial
## cost: 0.01
##
## Number of Support Vectors: 634
##
## ( 319 315 )
##
##
## Number of Classes: 2
##
## Levels:
## CH MMradial_001_train_pred <- predict(radial_001, OJ_train)
radial_001_test_pred <- predict(radial_001, OJ_test)
radial_001_train_error <- error_rate(OJ_train$Purchase, radial_001_train_pred)
radial_001_test_error <- error_rate(OJ_test$Purchase, radial_001_test_pred)
radial_001_train_error## [1] 0.39375radial_001_test_error## [1] 0.3777778Next, I tuned the cost while keeping the default gamma.
set.seed(1)
radial_tune_oj <- tune(
svm,
Purchase ~ .,
data = OJ_train,
kernel = "radial",
ranges = list(cost = c(0.01, 0.1, 1, 5, 10)),
tunecontrol = tune.control(cross = 10)
)
radial_tune_oj$performances## cost error dispersion
## 1 0.01 0.39375 0.04007372
## 2 0.10 0.18625 0.02853482
## 3 1.00 0.17125 0.02128673
## 4 5.00 0.18000 0.02220485
## 5 10.00 0.18625 0.02853482radial_tune_oj$best.parameters## cost
## 3 1radial_tune_oj$best.performance## [1] 0.17125best_radial_oj <- radial_tune_oj$best.model
best_radial_train_pred <- predict(best_radial_oj, OJ_train)
best_radial_test_pred <- predict(best_radial_oj, OJ_test)
best_radial_train_error <- error_rate(OJ_train$Purchase, best_radial_train_pred)
best_radial_test_error <- error_rate(OJ_test$Purchase, best_radial_test_pred)
best_radial_train_error## [1] 0.15125best_radial_test_error## [1] 0.1851852(g) Polynomial SVM with degree = 2
First, I fit the polynomial model using cost = 0.01 and degree = 2.
polynomial_001 <- svm(
Purchase ~ .,
data = OJ_train,
kernel = "polynomial",
degree = 2,
cost = 0.01,
scale = TRUE
)
summary(polynomial_001)##
## Call:
## svm(formula = Purchase ~ ., data = OJ_train, kernel = "polynomial",
## degree = 2, cost = 0.01, scale = TRUE)
##
##
## 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 MMpolynomial_001_train_pred <- predict(polynomial_001, OJ_train)
polynomial_001_test_pred <- predict(polynomial_001, OJ_test)
polynomial_001_train_error <- error_rate(OJ_train$Purchase, polynomial_001_train_pred)
polynomial_001_test_error <- error_rate(OJ_test$Purchase, polynomial_001_test_pred)
polynomial_001_train_error## [1] 0.3725polynomial_001_test_error## [1] 0.3666667Next, I tuned the cost while keeping degree = 2.
set.seed(1)
polynomial_tune_oj <- tune(
svm,
Purchase ~ .,
data = OJ_train,
kernel = "polynomial",
degree = 2,
ranges = list(cost = c(0.01, 0.1, 1, 5, 10)),
tunecontrol = tune.control(cross = 10)
)
polynomial_tune_oj$performances## cost error dispersion
## 1 0.01 0.39125 0.04210189
## 2 0.10 0.32125 0.05001736
## 3 1.00 0.20250 0.04116363
## 4 5.00 0.18250 0.03496029
## 5 10.00 0.18125 0.02779513polynomial_tune_oj$best.parameters## cost
## 5 10polynomial_tune_oj$best.performance## [1] 0.18125best_polynomial_oj <- polynomial_tune_oj$best.model
best_polynomial_train_pred <- predict(best_polynomial_oj, OJ_train)
best_polynomial_test_pred <- predict(best_polynomial_oj, OJ_test)
best_polynomial_train_error <- error_rate(OJ_train$Purchase, best_polynomial_train_pred)
best_polynomial_test_error <- error_rate(OJ_test$Purchase, best_polynomial_test_pred)
best_polynomial_train_error## [1] 0.15best_polynomial_test_error## [1] 0.1888889(h) Overall comparison and model-selection strategy
oj_results <- data.frame(
Model = c(
"Linear, cost 0.01",
"Tuned linear",
"Radial, cost 0.01",
"Tuned radial",
"Polynomial degree 2, cost 0.01",
"Tuned polynomial degree 2"
),
Training_Error = c(
linear_001_train_error,
best_linear_train_error,
radial_001_train_error,
best_radial_train_error,
polynomial_001_train_error,
best_polynomial_train_error
),
Test_Error = c(
linear_001_test_error,
best_linear_test_error,
radial_001_test_error,
best_radial_test_error,
polynomial_001_test_error,
best_polynomial_test_error
)
)
oj_results %>% arrange(Test_Error)## Model Training_Error Test_Error
## 1 Tuned linear 0.16500 0.1629630
## 2 Linear, cost 0.01 0.17500 0.1777778
## 3 Tuned radial 0.15125 0.1851852
## 4 Tuned polynomial degree 2 0.15000 0.1888889
## 5 Polynomial degree 2, cost 0.01 0.37250 0.3666667
## 6 Radial, cost 0.01 0.39375 0.3777778I fit six models: the original and tuned versions of the linear, radial, and polynomial SVMs. For each kernel, I used 10-fold cross-validation on the training set to select the cost with the lowest cross-validation error. The test set was not used while tuning the cost values.The tuned linear model had a test error rate of 0.163, which was the lowest test error. The original linear model had a test error rate of 0.178. The tuned radial model had a test error rate of 0.185, and the tuned polynomial model had a test error rate of 0.189.
Although the radial and polynomial models had slightly lower training errors, their test errors were higher. This suggests that the more flexible models did not improve performance on new observations. Overall, I would select the tuned linear support vector classifier with cost = 0.1. It had the lowest test error and was also simpler than the radial and polynomial models.