一文学会SVM
一文学会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)。
