Lab of Chapter 9

9.6.1 Support Vector Classifier

For data not linearly seperable

Setup data set:

library(e1071)
set.seed(1)
x <- matrix(rnorm(20 * 2), ncol = 2)
y <- c(rep(-1, 10), rep(1, 10))
x[y == 1, ] <- x[y == 1, ] + 1

这里从 x 数据集中选择了一半观测点, 给这些点的x坐标、y坐标分别加1，相当于将这些点向右上方平移了 \(\sqrt 2\)。

Check whether the classes are linearly separable:

plot(x, col = (3 - y))

Fit the support vector classifier and plot:

dat <- data.frame(x = x, y = as.factor(y))
svmfit <- svm(y ~ ., data = dat, kernel = 'linear', cost = 10, scale = FALSE)
plot(svmfit, dat)

由于这幅图将x2作为横坐标，x1作为纵坐标，直接看与上面原始数据点图对应不上，但只要把原始数据点图坐标变换到上图的方向（需要在图片浏览软件里先旋转，再翻转），就会发现20个数据点都能对应上。

上图中，青色和紫色的边界就是超平面所在位置（由于构建svm时指定了 kernel = 'linear'，所以边界是直线，由于绘图算法的原因，图中是一条有毛刺的直线），叉号代表支持向量，圆圈代表非支持向量，所以叉号覆盖的、与超平面平行的带状区域，就是margin所在区域（margin的含义可参考P342上的图9.3）。 在 (x1 = -1.2, x2 = 1) 位置上的红色叉号处于青色区域内，表示这个点被错误的分类了。

Summary of the model:

svmfit$index

[1]  1  2  5  7 14 16 17

summary(svmfit)

Call:
svm(formula = y ~ ., data = dat, kernel = "linear", cost = 10, scale = FALSE)


Parameters:
   SVM-Type:  C-classification 
 SVM-Kernel:  linear 
       cost:  10 
      gamma:  0.5 

Number of Support Vectors:  7

 ( 4 3 )


Number of Classes:  2 

Levels: 
 -1 1

(4 3) 表明上面的7个支持向量，4个属于 -1 类（黑色，都位于青色区域中），3个属于 1 类（红色，其中2个位于紫色区域中，1个位于青色区域中，表示被错误的分类），与上图中青色和紫色区域中的叉号分布吻合。

Build support vector classifier with smaller cost:

svmfit <- svm(y ~ ., data = dat, kernel = 'linear', cost = 0.1, scale = FALSE)
plot(svmfit, dat)

svmfit$index

 [1]  1  2  3  4  5  7  9 10 12 13 14 15 16 17 18 20

与 cost = 10 的输出图相比，margin 明显变宽了，覆盖了更多的数据点成为 support vector，由于更多的数据点参与了超平面的确定，可以认为超平面对数据变化的敏感度降低了，也就是variance降低了(bias升高了)。

Choose cost value with cross-validation:

set.seed(1)
tune.out <- tune(svm, y ~., data = dat, kernel = 'linear', ranges = list(cost = c(0.001, 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.1

- best performance: 0.1 

- Detailed performance results:
   cost error dispersion
1 1e-03  0.70  0.4216370
2 1e-02  0.70  0.4216370
3 1e-01  0.10  0.2108185
4 1e+00  0.15  0.2415229
5 5e+00  0.15  0.2415229
6 1e+01  0.15  0.2415229
7 1e+02  0.15  0.2415229

Summary the best model:

best.model <- tune.out$best.model
summary(best.model)

Call:
best.tune(method = svm, train.x = y ~ ., data = dat, 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:  0.1 
      gamma:  0.5 

Number of Support Vectors:  16

 ( 8 8 )


Number of Classes:  2 

Levels: 
 -1 1

Calculate the test error. First build the test dataset:

xtest <- matrix(rnorm(20 * 2), ncol = 2)
ytest <- sample(c(1, -1), 20, replace = TRUE)
xtest[ytest == 1, ] = xtest[ytest == 1, ] + 1
testdat <- data.frame(x = xtest, y = as.factor(ytest))

Predict the class on the test dataset:

ypred <- predict(best.model, testdat)
table(predict = ypred, truth = testdat$y)

       truth
predict -1 1
     -1  6 2
     1   3 9

The correct rate is: (7 + 9) / (7 + 4 + 0 + 9) = 80%.

Predict based on the model with cost = 0.01:

svmfit <- svm(y ~ ., data = dat, kernel = 'linear', cost = 0.01, scale = FALSE)
ypred <- predict(svmfit, testdat)
table(predict = ypred, truth = testdat$y)

       truth
predict -1 1
     -1  8 8
     1   1 3

The correct rate is (7 + 4) / (7 + 4 + 0 + 9) = 55%, which is lower than the model with cost 0.1.

For data linearly seperable

Move the half data point more faraway:

x[y == 1, ] = x[y == 1, ] + 0.5
plot(x, col = (y + 5) / 2, pch = 19)

这里生成的数据向右上平移0.5后仍然不可分，所以上面的代码将平移量改成了1。

Build svm:

dat <- data.frame(x = x, y = as.factor(y))
svmfit <- svm(y ~ ., data = dat, kernel = 'linear', scale = FALSE, cost = 1e5)
summary(svmfit)

Call:
svm(formula = y ~ ., data = dat, kernel = "linear", cost = 1e+05, scale = FALSE)


Parameters:
   SVM-Type:  C-classification 
 SVM-Kernel:  linear 
       cost:  1e+05 
      gamma:  0.5 

Number of Support Vectors:  3

 ( 1 2 )


Number of Classes:  2 

Levels: 
 -1 1

plot(svmfit, data = dat)

可以看到高 cost 导致 margin非常窄。

Build the svm with lower cost:

svmfit <- svm(y ~ ., data = dat, kernel = 'linear', scale = FALSE, cost = 1)
summary(svmfit)

Call:
svm(formula = y ~ ., data = dat, kernel = "linear", cost = 1, scale = FALSE)


Parameters:
   SVM-Type:  C-classification 
 SVM-Kernel:  linear 
       cost:  1 
      gamma:  0.5 

Number of Support Vectors:  7

 ( 3 4 )


Number of Classes:  2 

Levels: 
 -1 1

plot(svmfit, data = dat)

降低 cost 后，margin 宽了许多，其中一个数据点 (x1 = -0.5, x2 = 1.5) 被划分到错误的区域中。所以 bias 比 cost = 1e5 时高，相应地 variance 降低了。

9.6.2 Support Vector Machine

Generate some data with a non-linear class boundary:

set.seed(1)
x <- matrix(rnorm(200 * 2), ncol = 2)
x[1:100, ] <- x[1:100, ] + 2
x[101:150, ] <- x[101:150, ] - 2
y <- c(rep(1, 150), rep(2, 50))
dat <- data.frame(x = x, y = as.factor(y))
plot(x, col = y)

前100个点向右上平移，第101 ~ 150个点向左下平移，这两部分标记为1，中间夹的未移动的50个点标记为2。

Split the data into training and test set, and build a SVM (with \(\gamma = 1\)) on the training set:

train <- sample(200, 100)
svmfit <- svm(y ~ ., data = dat[train, ], kernel = 'radial', gamma = 1, cost = 1)
plot(svmfit, data = dat[train, ])

summary(svmfit)

Call:
svm(formula = y ~ ., data = dat[train, ], kernel = "radial", gamma = 1, 
    cost = 1)


Parameters:
   SVM-Type:  C-classification 
 SVM-Kernel:  radial 
       cost:  1 
      gamma:  1 

Number of Support Vectors:  37

 ( 17 20 )


Number of Classes:  2 

Levels: 
 1 2

Fit with a large cost:

svmfit <- svm(y ~ ., data = dat[train, ], kernel = 'radial', gamma = 1, cost = 1e5)
plot(svmfit, data = dat[train, ])

Select best \(\gamma\) and cost with cross-validation:

set.seed(1)
tune.out <- tune(svm, y ~ ., data = dat[train, ], kernel = 'radial', ranges = list(cost = c(0.1, 1, 10, 100, 1000), gamma = c(0.5, 1, 2, 3, 4)))
summary(tune.out)

Parameter tuning of ‘svm’:

- sampling method: 10-fold cross validation 

- best parameters:
 cost gamma
    1     2

- best performance: 0.12 

- Detailed performance results:
    cost gamma error dispersion
1  1e-01   0.5  0.27 0.11595018
2  1e+00   0.5  0.13 0.08232726
3  1e+01   0.5  0.15 0.07071068
4  1e+02   0.5  0.17 0.08232726
5  1e+03   0.5  0.21 0.09944289
6  1e-01   1.0  0.25 0.13540064
7  1e+00   1.0  0.13 0.08232726
8  1e+01   1.0  0.16 0.06992059
9  1e+02   1.0  0.20 0.09428090
10 1e+03   1.0  0.20 0.08164966
11 1e-01   2.0  0.25 0.12692955
12 1e+00   2.0  0.12 0.09189366
13 1e+01   2.0  0.17 0.09486833
14 1e+02   2.0  0.19 0.09944289
15 1e+03   2.0  0.20 0.09428090
16 1e-01   3.0  0.27 0.11595018
17 1e+00   3.0  0.13 0.09486833
18 1e+01   3.0  0.18 0.10327956
19 1e+02   3.0  0.21 0.08755950
20 1e+03   3.0  0.22 0.10327956
21 1e-01   4.0  0.27 0.11595018
22 1e+00   4.0  0.15 0.10801234
23 1e+01   4.0  0.18 0.11352924
24 1e+02   4.0  0.21 0.08755950
25 1e+03   4.0  0.24 0.10749677

Apply the best model on the test set:

table(prediction = predict(tune.out$best.model, newdata = dat[-train, ]), truth = dat[-train, ]$y)

          truth
prediction  1  2
         1 74  7
         2  3 16

So the correct rate is:

(68 + 21) / (68 + 1 + 10 + 21)

[1] 0.89

9.6.3 ROC Curves

Define a ROC plot function:

library(ROCR)

Loading required package: gplots

Attaching package: ‘gplots’

The following object is masked from ‘package:stats’:

    lowess

rocplot <- function(pred, truth, ...) {
  predob <- prediction(pred, truth)
  perf <- performance(predob, 'tpr', 'fpr')
  plot(perf, ...)
}

Decision values 是将计算出的系数代入到观测值向量后得到的值，超平面两侧的观测点各自的 decision values 分别大于和小于零，绝对值越大，离分隔面越远，当这个值为0时，表示观测点在分隔面上。

In order to obtain the fitted values for a given SVM model fit, we use decision.values=TRUE when fitting svm() . Then the predict() function will output the fitted values:

svmfit.opt <- svm(y ~ ., data = dat[train, ], kernel = 'radial', cost = 1, gamma = 2, decision.values = TRUE)
fit <- attributes(predict(svmfit.opt, dat[train, ], decision.values = TRUE))$decision.values
par(mfrow = c(1,2))
rocplot(fit, dat[train, ]$y, main = 'Training Data')
# By increasing $\gamma$ we can produce a more flexible fit and generate further improvements in accuracy:
svmfit.flex <- svm(y ~ ., data = dat[train, ], kernel = 'radial', cost = 1, gamma = 50, decision.values = TRUE)
fit <- attributes(predict(svmfit.flex, dat[train, ], decision.values = TRUE))$decision.values
rocplot(fit, dat[train, ]$y, add = TRUE, col = 'red')

注意这里不能把增加 gamma 后的代码放在新的 code chunk 里，不同 chunk 间的 plot 对象是不共享的。

Predict and plot on test dataset:

fit <- attributes(predict(svmfit.opt, dat[-train, ], decision.values = TRUE))$decision.values
rocplot(fit, dat[-train, ]$y, main = 'Test Data')
fit <- attributes(predict(svmfit.flex, dat[-train, ], decision.values = TRUE))$decision.values
rocplot(fit, dat[-train, ]$y, add = TRUE, col = 'red')

So the model svmfit.opt with gamma = 2 is better than the model svmfit.flex with gamma = 50.

9.6.4 SVM with Multiple Classes

set.seed(1)
x <- rbind(x, matrix(rnorm(50 * 2),  ncol = 2))
y <- c(y, rep(0, 50))
x[y == 0, 2] <- x[y == 0, 2] + 2
dat <- data.frame(x = x, y = as.factor(y))
par(mfrow = c(1,1))
plot(x, col = (y + 1))

可以看到黑色(y = 0)数据点相对于原始位置（绿色点所在区域）向上平移了2个单位（只给y坐标加了2）。

plot() 函数的 col 参数自带8种颜色，超过8后重新开始循环：

lst <- c(1:17)
barplot(lst, axes=FALSE, col=lst)

Fit an SVM to the data:

svmfit <- svm(y ~ ., data = dat, kernel = 'radial', gamma = 1, cost = 10)
plot(svmfit, dat)

再次提醒此图的坐标轴与上面的数据分布图不同，需要旋转90度再水平翻转才能与数据分布图对照。

响应变量的水平数为2时，使用标准 SVM 分类，如果大于2，使用 one-vs-one 方法（见9.4.1节）确定观测点的响应变量类型。

9.6.5 Application to Gene Expression Data

library(ISLR)
names(Khan)

[1] "xtrain" "xtest"  "ytrain" "ytest" 

dim(Khan$xtrain)

[1]   63 2308

dim(Khan$xtest)

[1]   20 2308

table(Khan$ytrain)

 1  2  3  4 
 8 23 12 20 

So there are 63 observations in training set and 20 observations in test set, each with 2308 features. The response variable has 4 levels, marked as 1 ~ 4.

Fit a SVM with a linear kernel on training set:

dat <- data.frame(x = Khan$xtrain, y = as.factor(Khan$ytrain))
out <- svm(y ~ ., data = dat, kernel = 'linear', cost = 10)
summary(out)

Call:
svm(formula = y ~ ., data = dat, kernel = "linear", cost = 10)


Parameters:
   SVM-Type:  C-classification 
 SVM-Kernel:  linear 
       cost:  10 
      gamma:  0.0004332756 

Number of Support Vectors:  58

 ( 20 20 11 7 )


Number of Classes:  4 

Levels: 
 1 2 3 4

table(out$fitted, dat$y)

   
     1  2  3  4
  1  8  0  0  0
  2  0 23  0  0
  3  0  0 12  0
  4  0  0  0 20

非对角线上全部是0，也就是没有 training bias，由于特征数大大高于观测数，这个结果是合理的。

Predict on the test set:

pred <- predict(out, Khan$xtest)
table(pred, Khan$ytest)

    
pred 1 2 3 4
   1 3 0 0 0
   2 0 6 2 0
   3 0 0 4 0
   4 0 0 0 5

非对角线上有2个观测：有两个属于3类的观测被预测为了2类。