1 Exercise 5 — Non-linear boundaries: logistic regression vs. SVM

1.1 (a) Generate data with a quadratic decision boundary

set.seed(5)
n  <- 500
x1 <- runif(n) - 0.5
x2 <- runif(n) - 0.5
y  <- as.numeric(x1^2 - x2^2 > 0)     # true boundary: |x1| = |x2|
dat <- data.frame(x1, x2, y = as.factor(y))
table(y)
## y
##   0   1 
## 254 246

The true boundary \(x_1^2 - x_2^2 = 0\) is a pair of crossing diagonals, so the classes form an “X” pattern — strongly non-linear.

1.2 (b) Plot the observations, coloured by true class

plot.class <- function(pred, main) {
  plot(x1, x2, col = ifelse(pred == 1, "red", "blue"),
       pch = 19, cex = 0.6, xlab = expression(X[1]), ylab = expression(X[2]),
       main = main)
}
plot.class(y, "True class labels")

The two classes occupy opposite pairs of the four triangular wedges created by the lines \(x_2 = \pm x_1\).

1.3 (c) Logistic regression with linear predictors

glm.lin <- glm(y ~ x1 + x2, data = dat, family = binomial)
summary(glm.lin)$coefficients
##                Estimate Std. Error    z value  Pr(>|z|)
## (Intercept) -0.03150042 0.08949213 -0.3519910 0.7248450
## x1          -0.06176136 0.30506037 -0.2024562 0.8395601
## x2          -0.11508575 0.31085517 -0.3702231 0.7112163

Neither predictor is significant — a linear model cannot detect this symmetric pattern.

1.4 (d) Predicted labels from the linear logistic model

p.lin    <- predict(glm.lin, dat, type = "response")
pred.lin <- as.numeric(p.lin > 0.5)
plot.class(pred.lin, "Logistic regression: linear predictors")

c(training_error = mean(pred.lin != y), n_predicted_class1 = sum(pred.lin == 1))
##     training_error n_predicted_class1 
##              0.548             92.000

The decision boundary is linear: the model can only cut the square with a straight line, so it slices off one corner (92 of the 500 points are assigned to class 1) rather than recovering the “X” pattern. The training error is about 0.548 — worse than random guessing on these balanced classes, confirming that a linear boundary is fundamentally the wrong shape here.

1.5 (e) Logistic regression with non-linear predictors

glm.nl <- glm(y ~ x1 + x2 + I(x1^2) + I(x2^2) + I(x1 * x2),
              data = dat, family = binomial)
summary(glm.nl)$coefficients
##                  Estimate Std. Error    z value Pr(>|z|)
## (Intercept)  3.326437e+14    5463222   60887820        0
## x1           2.028082e+14   10236293   19812661        0
## x2           1.535589e+13   10431559    1472061        0
## I(x1^2)      2.348870e+16   40604997  578468154        0
## I(x2^2)     -2.714208e+16   40267676 -674041375        0
## I(x1 * x2)  -5.006653e+14   34113923  -14676275        0

1.6 (f) Predicted labels from the non-linear logistic model

pred.nl <- as.numeric(predict(glm.nl, dat, type = "response") > 0.5)
plot.class(pred.nl, "Logistic regression: quadratic predictors")

c(training_error = mean(pred.nl != y))
## training_error 
##          0.066

The decision boundary is now obviously non-linear and recovers the “X” shape almost exactly, with a very small training error.

1.7 (g) Support vector classifier (linear kernel)

svc <- svm(y ~ x1 + x2, data = dat, kernel = "linear", cost = 10, scale = FALSE)
pred.svc <- predict(svc, dat)
plot.class(pred.svc, "Support vector classifier (linear kernel)")

c(training_error = mean(pred.svc != y),
  n_predicted_class1 = sum(pred.svc == 1))
##     training_error n_predicted_class1 
##              0.492              0.000

The support vector classifier also produces a linear boundary and fails even more completely than the linear logistic model: it assigns all 500 observations to a single class (0 predicted as class 1), giving a training error equal to the proportion of the minority class.

1.8 (h) SVM with a radial kernel

svm.rad <- svm(y ~ x1 + x2, data = dat, kernel = "radial", gamma = 1, cost = 10)
pred.rad <- predict(svm.rad, dat)
plot.class(pred.rad, "SVM with radial kernel")

c(training_error = mean(pred.rad != y))
## training_error 
##          0.022

1.9 (i) Comments

data.frame(
  model = c("Logistic (linear)", "Logistic (quadratic)",
            "SVC (linear kernel)", "SVM (radial kernel)"),
  training_error = round(c(mean(pred.lin != y), mean(pred.nl != y),
                           mean(pred.svc != y), mean(pred.rad != y)), 4)
)
##                  model training_error
## 1    Logistic (linear)          0.548
## 2 Logistic (quadratic)          0.066
## 3  SVC (linear kernel)          0.492
## 4  SVM (radial kernel)          0.022

The results show a clear pattern. Linear methods fail: both linear logistic regression and the linear support vector classifier produce a straight-line boundary and cannot separate classes whose true boundary is the pair of lines \(x_2 = \pm x_1\). Both have a training error near \(50\%\) — no better than guessing on these balanced classes. The SVC degenerates completely (every point assigned to one class), while the linear logistic model merely lops off a corner of the square; neither is usable.

Non-linear methods succeed: logistic regression with quadratic terms and the SVM with a radial kernel both recover the “X”-shaped boundary with very small training error. The key insight is that we do not need an SVM to get a non-linear boundary — logistic regression achieves the same thing if we manually supply the right non-linear transformations. The practical advantage of the SVM is that the kernel does this automatically: the radial kernel finds the non-linear boundary without our having to guess which transformations (\(x_1^2\), \(x_1x_2\), …) to include.

2 Exercise 7 — SVMs on the Auto data set

2.1 (a) Binary response: above/below median gas mileage

Auto <- na.omit(Auto)
Auto$mpg01 <- as.factor(ifelse(Auto$mpg > median(Auto$mpg), 1, 0))
Auto2 <- subset(Auto, select = -c(mpg, name))   # drop mpg and name
table(Auto$mpg01)
## 
##   0   1 
## 196 196

2.2 (b) Support vector classifier with various values of cost

set.seed(1)
tune.lin <- tune(svm, mpg01 ~ ., data = Auto2, kernel = "linear",
                 ranges = list(cost = c(0.01, 0.1, 1, 5, 10, 100)))
summary(tune.lin)
## 
## Parameter tuning of 'svm':
## 
## - sampling method: 10-fold cross validation 
## 
## - best parameters:
##  cost
##     1
## 
## - best performance: 0.08435897 
## 
## - Detailed performance results:
##    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 5e+00 0.08948718 0.03898410
## 5 1e+01 0.08948718 0.03898410
## 6 1e+02 0.08692308 0.03887151
lin.best <- tune.lin$best.parameters$cost
c(best_cost = lin.best, cv_error = round(min(tune.lin$performances$error), 4))
## best_cost  cv_error 
##    1.0000    0.0844

The CV error reaches its minimum at cost = 1, with a CV error of about 0.0844 (roughly 8.4%). What stands out is how insensitive the linear classifier is to cost here: across the whole range \(0.01\)–\(100\) the CV error only moves between about 0.084 and 0.092, a spread of well under two percentage points, with a shallow minimum at cost = 1. This tells us the two mileage groups are close to linearly separable, so the amount of slack allowed by cost barely matters.

2.3 (c) Radial and polynomial kernels

2.3.1 Radial kernel: tuning cost and gamma

set.seed(1)
tune.rad <- tune(svm, mpg01 ~ ., data = Auto2, kernel = "radial",
                 ranges = list(cost  = c(0.1, 1, 5, 10, 100),
                               gamma = c(0.01, 0.1, 1, 5)))
tune.rad$best.parameters
##    cost gamma
## 12    1     1
round(min(tune.rad$performances$error), 4)
## [1] 0.0663

2.3.2 Polynomial kernel: tuning cost and degree

set.seed(1)
tune.poly <- tune(svm, mpg01 ~ ., data = Auto2, kernel = "polynomial",
                  ranges = list(cost   = c(0.1, 1, 5, 10, 100),
                                degree = c(2, 3, 4)))
tune.poly$best.parameters
##    cost degree
## 10  100      3
round(min(tune.poly$performances$error), 4)
## [1] 0.0842
data.frame(
  kernel   = c("linear", "radial", "polynomial"),
  cv_error = round(c(min(tune.lin$performances$error),
                     min(tune.rad$performances$error),
                     min(tune.poly$performances$error)), 4)
)
##       kernel cv_error
## 1     linear   0.0844
## 2     radial   0.0663
## 3 polynomial   0.0842

Comments. The radial kernel is the clear winner, cutting the CV error to about 0.0663 versus roughly \(0.084\) for the other two — a reduction of about two percentage points. Its best setting is an intermediate gamma = 1 with cost = 1, and the tuning grid shows why intermediate is right: the best error at each gamma is 0.0842, 0.0791, 0.0663, 0.0917 for gamma = 0.01, 0.1, 1, 5. Too small a gamma gives an almost-linear boundary (no better than the linear kernel), while too large a gamma (\(5\)) makes the kernel very local and overfits, giving the worst error of all.

The polynomial kernel needs a large cost (100) and degree 3 to become competitive, and even then it only ties the linear classifier rather than beating it. So the ordering is radial \(\ll\) polynomial \(\approx\) linear: the classes are nearly linearly separable (hence the strong linear baseline), but a moderately flexible radial boundary still extracts a genuine further improvement.

2.4 (d) Plots supporting the conclusions

par(mfrow = c(1, 2))
# linear: CV error vs cost
plot(tune.lin$performances$cost, tune.lin$performances$error, type = "b", log = "x",
     xlab = "cost (log scale)", ylab = "CV error", main = "Linear kernel")
# radial: CV error vs cost, one line per gamma
perf <- tune.rad$performances
gams <- sort(unique(perf$gamma))
plot(range(perf$cost), range(perf$error), type = "n", log = "x",
     xlab = "cost (log scale)", ylab = "CV error", main = "Radial kernel")
for (i in seq_along(gams)) {
  sub <- perf[perf$gamma == gams[i], ]
  lines(sub$cost, sub$error, type = "b", col = i, pch = 19)
}
legend("topright", legend = paste("gamma =", gams), col = seq_along(gams), lty = 1, cex = 0.8)

par(mfrow = c(1, 1))

The left panel shows how flat the linear classifier’s CV error is across four orders of magnitude of cost — the vertical axis spans less than two percentage points, so no value of cost is meaningfully better than another. The right panel shows the radial kernel’s behaviour in gamma: the gamma = 1 curve sits lowest, the very small gamma = 0.01 curve is higher (the kernel is so smooth that it behaves almost like the linear classifier), and the large gamma = 5 curve is worst at every cost — a very local kernel that overfits the training data. This is the classic U-shaped flexibility trade-off, with the optimum at an intermediate gamma.

svm.rad.best <- svm(mpg01 ~ ., data = Auto2, kernel = "radial",
                    cost = tune.rad$best.parameters$cost,
                    gamma = tune.rad$best.parameters$gamma)
par(mfrow = c(1, 2))
plot(Auto2$weight, Auto2$displacement, col = ifelse(Auto2$mpg01 == 1, "red", "blue"),
     pch = 19, cex = 0.6, xlab = "weight", ylab = "displacement",
     main = "True classes")
plot(Auto2$weight, Auto2$displacement,
     col = ifelse(predict(svm.rad.best, Auto2) == 1, "red", "blue"),
     pch = 19, cex = 0.6, xlab = "weight", ylab = "displacement",
     main = "Radial SVM predictions")

par(mfrow = c(1, 1))

Plotting displacement against weight (two of the most informative predictors) shows that high-mileage cars (red) are light with small engines, and the tuned radial SVM reproduces that separation almost perfectly.

3 Exercise 8 — SVMs on the OJ data set

3.1 (a) Training set of 800 observations

set.seed(1)
train <- sample(nrow(OJ), 800)
OJ.train <- OJ[train, ]
OJ.test  <- OJ[-train, ]
err.rates <- function(fit) {
  c(train = mean(predict(fit, OJ.train) != OJ.train$Purchase),
    test  = mean(predict(fit, OJ.test)  != OJ.test$Purchase))
}

3.2 (b) Support vector classifier with cost = 0.01

svc <- svm(Purchase ~ ., data = OJ.train, kernel = "linear", cost = 0.01)
summary(svc)
## 
## 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:  435
## 
##  ( 219 216 )
## 
## 
## Number of Classes:  2 
## 
## Levels: 
##  CH MM

The fit uses 435 support vectors (219 in one class and 216 in the other) — a large fraction of the 800 training observations, which is expected with such a small cost (a wide, soft margin).

3.3 (c) Training and test error rates

round(err.rates(svc), 4)
##  train   test 
## 0.1750 0.1778

3.4 (d) Cross-validation to select cost

set.seed(1)
tune.lin <- tune(svm, Purchase ~ ., data = OJ.train, kernel = "linear",
                 ranges = list(cost = c(0.01, 0.05, 0.1, 0.5, 1, 5, 10)))
summary(tune.lin)
## 
## Parameter tuning of 'svm':
## 
## - sampling method: 10-fold cross validation 
## 
## - best parameters:
##  cost
##   0.5
## 
## - best performance: 0.16875 
## 
## - Detailed performance results:
##    cost   error dispersion
## 1  0.01 0.17625 0.02853482
## 2  0.05 0.17625 0.02853482
## 3  0.10 0.17250 0.03162278
## 4  0.50 0.16875 0.02651650
## 5  1.00 0.17500 0.02946278
## 6  5.00 0.17250 0.03162278
## 7 10.00 0.17375 0.03197764
lin.cost <- tune.lin$best.parameters$cost
lin.cost
## [1] 0.5

3.5 (e) Error rates with the optimal cost

svc.best <- svm(Purchase ~ ., data = OJ.train, kernel = "linear", cost = lin.cost)
lin.err  <- err.rates(svc.best)
round(lin.err, 4)
##  train   test 
## 0.1650 0.1556

Tuning cost to 0.5 gives a modest but real improvement over cost = 0.01: the test error falls from 0.1778 to 0.1556 (about two percentage points), and the training error drops as well.

3.6 (f) SVM with a radial kernel (default gamma)

svm.rad <- svm(Purchase ~ ., data = OJ.train, kernel = "radial")
c(support_vectors = svm.rad$tot.nSV)
## support_vectors 
##             373
round(err.rates(svm.rad), 4)
##  train   test 
## 0.1512 0.1852
set.seed(1)
tune.rad <- tune(svm, Purchase ~ ., data = OJ.train, kernel = "radial",
                 ranges = list(cost = c(0.01, 0.05, 0.1, 0.5, 1, 5, 10)))
rad.cost <- tune.rad$best.parameters$cost
svm.rad.best <- svm(Purchase ~ ., data = OJ.train, kernel = "radial", cost = rad.cost)
rad.err <- err.rates(svm.rad.best)
c(best_cost = rad.cost)
## best_cost 
##       0.5
round(rad.err, 4)
##  train   test 
## 0.1475 0.1778

3.7 (g) SVM with a polynomial kernel (degree = 2)

svm.poly <- svm(Purchase ~ ., data = OJ.train, kernel = "polynomial", degree = 2)
c(support_vectors = svm.poly$tot.nSV)
## support_vectors 
##             447
round(err.rates(svm.poly), 4)
##  train   test 
## 0.1825 0.2222
set.seed(1)
tune.poly <- tune(svm, Purchase ~ ., data = OJ.train, kernel = "polynomial", degree = 2,
                  ranges = list(cost = c(0.01, 0.05, 0.1, 0.5, 1, 5, 10)))
poly.cost <- tune.poly$best.parameters$cost
svm.poly.best <- svm(Purchase ~ ., data = OJ.train, kernel = "polynomial",
                     degree = 2, cost = poly.cost)
poly.err <- err.rates(svm.poly.best)
c(best_cost = poly.cost)
## best_cost 
##        10
round(poly.err, 4)
##  train   test 
## 0.1500 0.1889

3.8 (h) Which approach is best?

results <- data.frame(
  kernel     = c("Linear", "Radial", "Polynomial (deg 2)"),
  best_cost  = c(lin.cost, rad.cost, poly.cost),
  train_error = round(c(lin.err["train"], rad.err["train"], poly.err["train"]), 4),
  test_error  = round(c(lin.err["test"],  rad.err["test"],  poly.err["test"]),  4)
)
results
##               kernel best_cost train_error test_error
## 1             Linear       0.5      0.1650     0.1556
## 2             Radial       0.5      0.1475     0.1778
## 3 Polynomial (deg 2)      10.0      0.1500     0.1889
best <- results$kernel[which.min(results$test_error)]
best
## [1] "Linear"

Comparing the tuned models on the held-out test set, the Linear kernel gives the lowest test error rate (0.1556), with the radial and polynomial kernels about two and three percentage points behind respectively.

The most instructive feature of the table is that the training and test rankings are exactly reversed. By training error the flexible kernels look best (radial 0.1475 < polynomial 0.15 < linear 0.165), but by test error the ordering flips (linear 0.1556 < radial 0.1778 < polynomial 0.1889). The extra flexibility of the radial and polynomial kernels buys a better fit to the training data but does not generalise — a textbook illustration of overfitting, and precisely why model selection must rest on cross-validation or held-out data rather than training error. Since the simplest model also generalises best, the linear support vector classifier is the right choice for these data.