一文学会SVM

SVM 背景知识

支持向量机，因其英文名为support vector machine，故一般简称SVM，就是一种二类分类模型，其基本模型定义为特征空间上的间隔最大的线性分类器，其学习策略便是间隔最大化，最终可转化为一个凸二次规划问题的求解。

看起来这个定义是不是专有名词太多呀！其实还有要完全理解SVM原理及算法，还需要理解 线性回归，最小二乘法，逻辑回归，线性分类器，线性可分，核函数，损失函数，

但是不要怕，不具体理解SVM原理及算法，我们仍然是可以使用它，左右不过是一个分类器罢了，就是根据一堆自变量来预测因变量，所以就是变量预测，值得一提的是，SVM通常应用于二元分类变量预测，但是经过一些改进也可以勉强对多元分类变量预测，同时基于SVM的SVR也可以预测连续变量。

通俗的理解，我们想根据年收入来预测某家庭是贫穷还是富有，可以简单的按照年收入50万来进行分类，这个时候就只有一个自变量，就是收入的金额这个数值，因变量也很简单，就是二元分类情况。只不过通常我们要使用SVM的场景，因变量肯定不止一个，阈值也没有那么简单找到。

SVM示例

二元分类变量预测

毫无疑问，生物学领域最经典的二元分类变量就是病人的生死问题啦！

load('~/Documents/Nutstore/github/TCGA-KIRC-miRNA-example/Rdata/TCGA-LUAD-survival_input.Rdata')
## 上面的测试数据大家可以发邮件给我索要，我的邮箱是  jmzeng1314@163.com 

# 首先你会有一个表达矩阵如下，每个病人的每个基因都有表达量。
exprSet[1:4,1:2]
##              TCGA-05-4244-01A-01T-1108-13 TCGA-05-4249-01A-01T-1108-13
## hsa-let-7a-1                         3985                         8916
## hsa-let-7a-2                         7947                        17800
## hsa-let-7a-3                         4128                         9079
## hsa-let-7b                           9756                        32960
# 然后你会有这些病人的临床信息
head(phe)
##                     ID event race age gender stage days age_group
## 52.70.0   TCGA-05-4244     0 <NA>  70   male    iv    0     older
## 52.70.0.2 TCGA-05-4249     0 <NA>  67   male    ib 1158     older
## 52.70.0.3 TCGA-05-4250     1 <NA>  79 female  iiia  121     older
## 58.73.0   TCGA-05-4382     0 <NA>  68   male    ib  607     older
## 58.73.0.1 TCGA-05-4389     0 <NA>  70   male    ia 1369     older
## 58.73.0.2 TCGA-05-4395     1 <NA>  76   male  iiib    0     older
##                time
## 52.70.0    0.000000
## 52.70.0.2 38.600000
## 52.70.0.3  4.033333
## 58.73.0   20.233333
## 58.73.0.1 45.633333
## 58.73.0.2  0.000000
# 当然，我这里举例就只关心 生死这个情况。
table(phe$event)
## 
##   0   1 
## 388 125
# 125个病人去世了，有388活着的。

## 需要进行简单的转换保证病人表达矩阵和生死信息合并。
exprSet=mean(colSums(exprSet))*exprSet/colSums(exprSet)
exprSet=log2(exprSet+1)
data=cbind(t(exprSet),phe$event)
colnames(data)[ncol(data)]='event'
data=as.data.frame(data)
data$event=as.factor(data$event)
# 最后得到的data这个变量就可以进入SVM啦。

我们可以先把上面的TCGA-LUAD的miRNA数据集分成50%的训练集(Train)，50%的测试集(Test)：

## 简单粗暴的区分测试集和训练集
require("mlbench")
## Loading required package: mlbench
library(e1071)
## Warning: package 'e1071' was built under R version 3.4.4
smp.size = floor(0.5*nrow(data)) 
set.seed(123456789)                     
train.ind = sample(seq_len(nrow(data)), smp.size)
train = data[train.ind, ] # 50%
test = data[-train.ind, ] # 50%
table(test$event)
## 
##   0   1 
## 194  63
table(train$event)
## 
##   0   1 
## 194  62

然后就可以使用svm()训练分类模型

model = svm(formula = event ~ .,  # 这里的待预测的变量event是二分类变量，生与死。
            data = train,kernel = "linear")
## 值得注意的是这里默认会选择 kernel = "radial" ，核函数的概念需要理解。
#
summary(model) 
## 
## Call:
## svm(formula = event ~ ., data = train, kernel = "linear")
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  linear 
##        cost:  1 
##       gamma:  0.001485884 
## 
## Number of Support Vectors:  171
## 
##  ( 59 112 )
## 
## 
## Number of Classes:  2 
## 
## Levels: 
##  0 1

共有9 种核函数，常用的为其中的前四个：linear，Polynomial，RBF，Sigmoid 其中 RBF 适用于因变量比较少，而 linear适用于因变量非常多，也就是本例子里面的基因非常多，所以我们现在linear。

学习资料： - https://data-flair.training/blogs/svm-kernel-functions/ - https://www.quora.com/What-are-kernels-in-machine-learning-and-SVM-and-why-do- -we-need-them - https://www.zhihu.com/question/21883548 - https://www.quora.com/How-do-I-select-SVM-kernels

训练好的模型，就可以使用predict()函数去测试集里面看看效果。

train.pred = predict(model, train)
test.pred = predict(model, test)

## 训练集的混淆矩阵
table(real=train$event, predict=train.pred) 
##     predict
## real   0   1
##    0 194   0
##    1   0  62
confus.matrix = table(real=train$event, predict=train.pred)
sum(diag(confus.matrix))/sum(confus.matrix)
## [1] 1
# 可以很明显的看到过拟合现象，在训练集预测100%。

## 测试集的混淆矩阵
table(real=test$event, predict=test.pred)
##     predict
## real   0   1
##    0 162  32
##    1  47  16
confus.matrix = table(real=test$event, predict=test.pred)
sum(diag(confus.matrix))/sum(confus.matrix)
## [1] 0.692607
# 可以看到准确率并不高。

关于模型的各种其它检验指标，建议直接阅读我在生信技能树分享的文章：https://vip.biotrainee.com/d/812-

多元分类变量预测

对于多元分类变量预测例子，我们可以使用mlbench中的数据集Glass来简单演示,数据集里面的type变量是6个分类，看下面代码及输出结果。

简单探索一下数据集Glass，里面记录着214个观测值，10个变量。

require("mlbench")
library(e1071)
data(Glass, package="mlbench")
data = Glass

在e1071的包里面，我们可以直接使用svm()这个函数建模。

同样的我们可以先把Glass数据集分成80%的训练集(Train)，20%的测试集(Test)：

smp.size = floor(0.8*nrow(data)) 
set.seed(123456789)                     
train.ind = sample(seq_len(nrow(data)), smp.size)
train = data[train.ind, ] # 80%
test = data[-train.ind, ] # 20%

然后就可以使用svm()训练分类模型

model = svm(formula = Type ~ .,  # 这里的待预测的变量Type必须是Factor
            data = train)
## 本例子变量并不多，所以实验默认的kernel ="radial" 来建模。
summary(model) 
## 
## Call:
## svm(formula = Type ~ ., data = train)
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  radial 
##        cost:  1 
##       gamma:  0.1111111 
## 
## Number of Support Vectors:  148
## 
##  ( 16 53 11 48 11 9 )
## 
## 
## Number of Classes:  6 
## 
## Levels: 
##  1 2 3 5 6 7

训练好的模型，就可以使用predict()函数去测试集里面看看效果。

train.pred = predict(model, train)
test.pred = predict(model, test)

## 训练集的混淆矩阵
table(real=train$Type, predict=train.pred) 
##     predict
## real  1  2  3  5  6  7
##    1 47  8  0  0  0  0
##    2 10 48  0  0  0  0
##    3  9  7  0  0  0  0
##    5  0  0  0 11  0  0
##    6  0  1  0  0  8  0
##    7  1  0  0  0  0 21
confus.matrix = table(real=train$Type, predict=train.pred)
sum(diag(confus.matrix))/sum(confus.matrix)
## [1] 0.7894737
## 测试集的混淆矩阵
table(real=test$Type, predict=test.pred)
##     predict
## real  1  2  3  5  6  7
##    1  9  6  0  0  0  0
##    2  5 13  0  0  0  0
##    3  0  1  0  0  0  0
##    5  0  0  0  2  0  0
##    6  0  0  0  0  0  0
##    7  0  3  0  0  0  4
confus.matrix = table(real=test$Type, predict=test.pred)
sum(diag(confus.matrix))/sum(confus.matrix)
## [1] 0.6511628

SVR预测连续变量

这里直接摘抄：https://rpubs.com/skydome20/R-Note14-SVM-SVR 的笔记，他讲解的非常好！

SVR是本來SVM的延伸型態，能夠處理連續的預測問題。

在e1071套件裡面，並沒有一個函式叫做svr()，而是一樣用svm()。

差別只在於：

• 依變數的型態是factor時，svm()會建立SVM的超平面，來處理分類問題
• 依變數的型態是numeric，svm()會轉為SVR，進行連續值的預測。

這裡簡單手動建立一筆資料：

data = data.frame(x=1:20,
                  y=c(3,4,8,2,6,10,12,13,15,14,17,18,20,17,21,22,25,30,29,31))

# 資料的原始值
plot(data$x, data$y, pch=16, xlab="X", ylab="Y")

image

我們可以先拉一條簡單的線性迴歸：

model <- lm(y ~ x , data) 

# lm預測
lm.pred = predict(model, data)

# 資料的原始值(黑點)
plot(data$x, data$y, pch=16, xlab="X", ylab="Y")
# lm的預測值(紅三角形)
points(lm.pred, pch=2, col="red")
abline(model, col="red")

image

我們可以直接用SVR來建模、預測：

model <- svm(y ~ x , data) # 依變數的型態要是numeric

# 預測
svr.pred = predict(model, data)

# 資料的原始值(黑點)
plot(data$x, data$y, pch=16, xlab="X", ylab="Y")
# SVR的預測值(藍叉)
points(svr.pred, pch=4, col="blue")

image

最後比較一下線性迴歸和SVR的表現，同時計算RMSE(root mean square error)：

# 資料的原始值(黑點)
plot(data$x, data$y, pch=16, xlab="X", ylab="Y")
# lm的預測值(紅三角形)
points(lm.pred, pch=2, col="red")
# SVR的預測值(藍叉)
points(svr.pred, pch=4, col="blue")

image

# (lm, SVR) in RMSE
c(sqrt(mean((data$y - lm.pred)^2)),
  sqrt(mean((data$y - svr.pred)^2))
)
## [1] 1.914203 1.795094

可以發現到，在這個例子，SVR比lm的效果還要好一些些(1.79 < 1.91)。

---