MachineLearning 32. 机器学习之支持向量机递归特征消除的特征筛选 (mSVM-RFE)

简介

本文提出了一种新的特征选择方法，该方法使用类似于支持向量机递归特征消除 (SVM-RFE)的反向消除过程。与 SVM-RFE 方法不同的是，在每一步中，该方法通过对原始训练数据的子样本上训练的多个线性支持向量机的权重向量进行统计分析来计算特征排序得分。我们在四个用于癌症分类的基因表达数据集上测试了所提出的方法。结果表明，所提出的特征选择方法比原 SVM-RFE 方法选择了更好的基因子集，提高了分类精度。基于基因本体的相似性评估表明，选择的子集在功能上是多样的，进一步验证了我们的基因选择方法。本研究还表明，对于基于基因表达的癌症分类，可以推荐来自多个训练集和测试集的平均测试误差作为性能质量的参考。

一种从一组初始特征开始反向工作的迭代算法。在每一轮中，

1)拟合一个简单的线性支持向量机，

2)根据SVM解决方案中的权重对特征进行排序，

3)消除权重最低的特征。

软件包安装

这个 mRFE 依赖软件包 SVM-RFE，进一步优化而得到了多重 SVM-RFE

if(!require(SVM-RFE))
  devtools::install_github("llrs/SVM-RFE")
install.packages('e1071')

加载 mRFE，这个软件包需要在github上下载源代码 https://codeload.github.com/johncolby/SVM-RFE/zip/refs/heads/master ，然后 source(),即可使用了。

library(mRFE)
# > Loading required package: e1071 basic example code Set up R environment
library(e1071)
source("./SVM-RFE-master/msvmRFE.R")

输入数据应该格式化为数据框，每一行显示一个观测值，每一列显示一个特征。第一列应该包含真正的类标签。因素特征应该编码为多个数字“虚拟”特征。

BreastCancer <- read.csv("wisc_bc_data.csv", stringsAsFactors = FALSE)
str(BreastCancer)
## 'data.frame':	568 obs. of  32 variables:
##  $ id                     : int  842517 84300903 84348301 84358402 843786 844359 84458202 844981 84501001 845636 ...
##  $ diagnosis              : chr  "M" "M" "M" "M" ...
##  $ radius_mean            : num  20.6 19.7 11.4 20.3 12.4 ...
##  $ texture_mean           : num  17.8 21.2 20.4 14.3 15.7 ...
##  $ perimeter_mean         : num  132.9 130 77.6 135.1 82.6 ...
##  $ area_mean              : num  1326 1203 386 1297 477 ...
##  $ smoothness_mean        : num  0.0847 0.1096 0.1425 0.1003 0.1278 ...
##  $ compactne_mean         : num  0.0786 0.1599 0.2839 0.1328 0.17 ...
##  $ concavity_mean         : num  0.0869 0.1974 0.2414 0.198 0.1578 ...
##  $ concave_points_mean    : num  0.0702 0.1279 0.1052 0.1043 0.0809 ...
##  $ symmetry_mean          : num  0.181 0.207 0.26 0.181 0.209 ...
##  $ fractal_dimension_mean : num  0.0567 0.06 0.0974 0.0588 0.0761 ...
##  $ radius_se              : num  0.543 0.746 0.496 0.757 0.335 ...
##  $ texture_se             : num  0.734 0.787 1.156 0.781 0.89 ...
##  $ perimeter_se           : num  3.4 4.58 3.44 5.44 2.22 ...
##  $ area_se                : num  74.1 94 27.2 94.4 27.2 ...
##  $ smoothness_se          : num  0.00522 0.00615 0.00911 0.01149 0.00751 ...
##  $ compactne_se           : num  0.0131 0.0401 0.0746 0.0246 0.0335 ...
##  $ concavity_se           : num  0.0186 0.0383 0.0566 0.0569 0.0367 ...
##  $ concave_points_se      : num  0.0134 0.0206 0.0187 0.0188 0.0114 ...
##  $ symmetry_se            : num  0.0139 0.0225 0.0596 0.0176 0.0216 ...
##  $ fractal_dimension_se   : num  0.00353 0.00457 0.00921 0.00511 0.00508 ...
##  $ radius_worst           : num  25 23.6 14.9 22.5 15.5 ...
##  $ texture_worst          : num  23.4 25.5 26.5 16.7 23.8 ...
##  $ perimeter_worst        : num  158.8 152.5 98.9 152.2 103.4 ...
##  $ area_worst             : num  1956 1709 568 1575 742 ...
##  $ smoothness_worst       : num  0.124 0.144 0.21 0.137 0.179 ...
##  $ compactne_worst        : num  0.187 0.424 0.866 0.205 0.525 ...
##  $ concavity_worst        : num  0.242 0.45 0.687 0.4 0.535 ...
##  $ concave_points_worst   : num  0.186 0.243 0.258 0.163 0.174 ...
##  $ symmetry_worst         : num  0.275 0.361 0.664 0.236 0.399 ...
##  $ fractal_dimension_worst: num  0.089 0.0876 0.173 0.0768 0.1244 ...
dim(BreastCancer)
## [1] 568  32
table(BreastCancer$diagnosis)
##   B   M 
## 357 211
rownames(BreastCancer) = BreastCancer[, 1]
input = na.omit(BreastCancer[, -1])
sum(is.na(input))  # 检测数据是否有缺失
## [1] 0
input$diagnosis = as.factor(input$diagnosis)

这里我们已经表明，我们希望k=10作为mSVM-RFE的“多重”部分的k-fold交叉验证。要使用标准SVM-RFE，可以使用k=1。还要注意这个二分之一。以上参数。这允许你每轮将特征切成两半，而不是一个接一个。这对于具有许多特征的数据集非常有用。这里是对半。Above =100，所以每轮的特征将被切成两半，直到剩下的少于100个。输出是一个特征指数的向量，现在从最有用到最没用排序。请注意，由于CV抽签的随机性，排名分数接近的特征，以及可能包含的无用特征，这些排名可能会发生变化。但是，这个演示的输出应该是相同的，因为我们都将随机种子重置为相同的值。

要进行特征排序，使用svmRFE()函数:

set.seed(12345)
svmRFE(input, k = 10, halve.above = 100)
## Scaling data...Done!
##  [1] 21 30 22  7 20 18  6 17 14  2 12 24 15  5 29 27 26 16 11  1  3  8 23 19  4
## [26] 13 25 28 10  9

基本上，方法是将整个特征选择和泛化误差估计过程包装在外部交叉验证的顶层循环中。对于10倍CV，我们首先定义哪些观察值在哪些折叠中。

在R语言中，许多并行函数喜欢在列表上操作，所以接下来我们将把折叠重新格式化为列表，其中每个列表元素都是包含该折叠的测试集索引的向量。

# Set up cross validation
nfold = 10
nrows = nrow(input)
folds = rep(1:nfold, len = nrows)[sample(nrows)]
folds
##   [1]  2  2  7  4  9  1  3  3  8  9 10 10  6  3  3  1  6  6  7  1  3  5  2 10  2
##  [26]  2  9  4  2  4  6  2  7  8  2  7  5  5  8  7  3  6  7  8  7  4  4  1  4  8
##  [51]  6  3 10 10  3  6  6 10  3  6  3  4  8  6  6  9  3  7  4  3  9  5  1  7  5
##  [76] 10  4  3  5  7  6  9  7 10  7  3  5  1  7  8  3  8  6  5  5 10  4  9  8  9
## [101]  9  3  7  5  4  3  6  6  3  1  7  8  1  4  2  7 10  2  1  5  1  3  8  5 10
## [126]  5  3  2  2  5  6 10  7  3  7  3  1  9  2  8  5  3  7  6  9  5  3  7  6  1
## [151] 10 10  7  3  5  5  7  5  6 10  6  4  1  3  3 10  4  4 10  7  2  1  2  1  1
## [176]  1  7  4  3 10  2  7  9  6  4  9  5  6  9  6  7  8  8  4  2  6  5  1  6  9
## [201]  2  6  3  2  4  2  1  7  2 10  2  9  1  6  6  5 10  8  4  1  7  8  3 10  3
## [226]  8  2 10  1  6 10  8  9  1  7 10 10  7  6  7  6  9  2  9  7 10  9  4  5  7
## [251]  3  4  5  9  6  9  4  8  4  1  7  8  8  4  7  4  9  4  8  8 10  5  2  8  8
## [276]  7  5  5  1  8  5  9  2 10  4  7 10  6  1  8  2  8  3 10  5  9 10  7  8  9
## [301]  2  1  5  5  6 10  9 10  3  6  7  2 10  2  5  1  3  1  1  6  5  3  1  7  9
## [326]  5  2 10  3  2  9  9  1  8  1  7  3  4  4  6  5  1 10  8  9  1  9  1  1  9
## [351]  6 10  1  1  4  4  8  2  5  1  6  3  7  7  9  8  7  9  6  7  4  3  8  4 10
## [376]  9 10  4  8  4  4  6  7  8  3  9  1  1  8  9  8  4  2  8  1  9 10  9  8  2
## [401]  6  9 10 10  8  7  8  1  9  9  2  7  3  2  2  2  5  6  6  2  5  3  1  6  9
## [426]  4  4  1  2  5  1  4  3  5  4  5  7  4  3  4  1  2  4  4  8 10  1 10  2  9
## [451]  5  3  5  4  3  9 10  6  1  4  4  2  9  2  9  7  2  2  3  6  6  5  8  2  5
## [476]  9  3  8  5  9  1  8  3  4  5  5  6  2  7  4  7  5  6  6  6  6  7  8  9  2
## [501]  7  1  5 10  9  3  5  1  5  8  1  8  1 10  8  5  2  9  8 10 10  1  3  3  3
## [526]  3 10  5 10 10 10  5 10  2  6  4  7  2  8  2  8  9  4  8  4  4  6  9  9  2
## [551]  5  8  8  6  4 10  2  3  7  7  8  4  2  1  3  5  6 10
folds = lapply(1:nfold, function(x) which(folds == x))

对所有训练集进行特征排序

使用lapply或它的一个通用的并行表。现在可以对所有10个训练集执行特征排序。

# Perform feature ranking on all training sets
results = lapply(folds, svmRFE.wrap, input, k = 10, halve.above = 100)
## Scaling data...Done!
length(results)
## [1] 10

获得all folds 的 top 特征.

如果我们要将这些发现应用到某个地方的最终测试集，我们仍然需要所有这些训练数据的最佳特征。

# Obtain top features across ALL folds
top.features = WriteFeatures(results, input, save = F)
head(top.features)
##               FeatureName FeatureID AvgRank
## 1              area_worst        24     2.9
## 2           texture_worst        22     3.0
## 3 fractal_dimension_worst        30     6.1
## 4    fractal_dimension_se        20     6.9
## 5                 area_se        14     7.5
## 6       concave_points_se        18     8.0

使用不同数量的top特征估计泛化误差

现在我们有了10个训练集的特征排名，最后一步是估计如果我们在这些特征上训练最终分类器并将其应用于新的测试集，我们可以预期的泛化误差。在这里，径向基函数核支持向量机对每个训练集进行独立的调优。这包括使用网格搜索对SVM超参数(Cost和Gamma)的每个组合进行内部10倍CV误差估计。然后使用最优参数在整个训练集上训练支持向量机。最后，通过预测相应的测试集来确定泛化误差。这是对外部10倍CV中的每一倍进行的，并且所有10个这些泛化误差估计都是平均的，以提供更多的稳定性。这个过程是重复的，同时改变作为输入的顶级功能的数量，通常会有一个“最佳点”，即不会有太多或太少的功能。概述一下，这个过程看起来像:

为了在前5个功能中实现它，我们这样做:

featsweep = lapply(1:5, FeatSweep.wrap, results, input)
saveRDS(featsweep,file = 'featsweep.rds') #保存 rds

概化误差与特征数量的关系图

为了直观地显示这些结果，我们可以绘制平均泛化误差与用作输入的顶级特征数量的关系。为供参考，给出机会错误率是有用的。通常，如果我们总是简单地选择数据集中流行率更高的类标签，这等于我们得到的“无信息”率。

featsweep <- readRDS("featsweep.rds")  # 读取 rd
# Make plot
no.info = min(prop.table(table(input[, 1])))
errors = sapply(featsweep, function(x) ifelse(is.null(x), NA, x$error))
top.features[1:which.min(errors), "FeatureName"]
## [1] "area_worst"              "texture_worst"          
## [3] "fractal_dimension_worst" "fractal_dimension_se"
PlotErrors(errors, no.info = no.info)

Lasso回归筛选特征

同样的数据，我们再次使用Lasso回归看下最后筛选出来的特征变量与SVM-RFE筛选出来的特征变量有何异同。首先构建 LASSO回归模型，记得alpha=1为Lasso回归，alpha=0为logistics回归，然后是二分类因变量，family=“binomial”，然后交叉验证cv.glmnet(),最后 cvfit$lambda.min找到最优lambda，大于lambda.min都为最后筛选的特征变量。

######## Lasso
library(glmnet)
library(VennDiagram)
x = as.matrix(input[, -1])
y <- ifelse(input$diagnosis == "B", 0, 1)
fit = glmnet(x, y, family = "binomial", alpha = 1)
## Call:  glmnet(x = x, y = y, family = "binomial", alpha = 1) 
##     Df  %Dev  Lambda
## 1    0  0.00 0.38330
## 2    2  8.19 0.34920
## 3    2 15.55 0.31820
## 4    2 21.82 0.28990
## 5    2 27.26 0.26420
## 6    2 32.05 0.24070
## 7    3 36.36 0.21930
## 8    3 40.25 0.19980
## 9    3 43.78 0.18210
## 10   3 46.97 0.16590
## 11   2 49.88 0.15120
## 12   2 52.53 0.13770
## 13   3 54.95 0.12550
## 14   3 57.18 0.11440
## 15   3 59.22 0.10420
## 16   4 61.40 0.09494
## 17   4 63.44 0.08650
## 18   4 65.32 0.07882
## 19   4 67.06 0.07182
## 20   4 68.68 0.06544
## 21   4 70.17 0.05962
## 22   4 71.55 0.05433
## 23   4 72.84 0.04950
## 24   4 74.03 0.04510
## 25   5 75.15 0.04110
## 26   5 76.22 0.03745
## 27   5 77.22 0.03412
## 28   6 78.18 0.03109
## 29   7 79.11 0.02833
## 30   7 80.00 0.02581
## 31   7 80.82 0.02352
## 32   7 81.59 0.02143
## 33   8 82.31 0.01952
## 34   8 82.99 0.01779
## 35   8 83.62 0.01621
## 36   8 84.20 0.01477
## 37   8 84.75 0.01346
## 38   9 85.26 0.01226
## 39   9 85.74 0.01117
## 40   9 86.18 0.01018
## 41  10 86.66 0.00928
## 42  10 87.13 0.00845
## 43  10 87.56 0.00770
## 44  10 87.95 0.00702
## 45  10 88.31 0.00639
## 46  10 88.64 0.00583
## 47  10 88.94 0.00531
## 48  11 89.24 0.00484
## 49  12 89.53 0.00441
## 50  13 89.83 0.00402
## 51  13 90.11 0.00366
## 52  13 90.36 0.00333
## 53  15 90.60 0.00304
## 54  15 90.85 0.00277
## 55  17 91.08 0.00252
## 56  17 91.32 0.00230
## 57  17 91.53 0.00209
## 58  16 91.71 0.00191
## 59  16 91.87 0.00174
## 60  17 92.03 0.00158
## 61  17 92.18 0.00144
## 62  17 92.31 0.00132
## 63  17 92.44 0.00120
## 64  17 92.56 0.00109
## 65  15 92.67 0.00099
## 66  16 92.76 0.00091
## 67  16 92.87 0.00083
## 68  17 92.97 0.00075
## 69  18 93.07 0.00069
## 70  19 93.18 0.00062
## 71  19 93.29 0.00057
## 72  19 93.40 0.00052
## 73  21 93.51 0.00047
## 74  21 93.64 0.00043
## 75  22 93.75 0.00039
## 76  22 93.93 0.00036
## 77  23 94.08 0.00033
## 78  23 94.20 0.00030
## 79  23 94.32 0.00027
## 80  23 94.43 0.00025
## 81  23 94.53 0.00022
## 82  23 94.64 0.00020
## 83  23 94.75 0.00019
## 84  24 94.86 0.00017
## 85  24 94.96 0.00015
## 86  23 95.06 0.00014
## 87  23 95.15 0.00013
## 88  24 95.23 0.00012
## 89  25 95.31 0.00011
## 90  25 95.39 0.00010
## 91  25 95.46 0.00009
## 92  26 95.52 0.00008
## 93  26 95.58 0.00007
## 94  26 95.63 0.00007
## 95  26 95.67 0.00006
## 96  26 95.72 0.00006
## 97  27 95.75 0.00005
## 98  27 95.79 0.00005
## 99  27 95.83 0.00004
## 100 27 95.86 0.00004
par(mfrow = c(1, 2))
plot(fit, label = TRUE)
# 绘制LASSO回归10折交叉验证图
cvfit = cv.glmnet(x, y, nfold = 10, family = "binomial", type.measure = "class")
plot(cvfit)

# 查看最佳lambda
cvfit$lambda.min
## [1] 0.001314779
# 获取LASSO选出来的特征
myCoefs <- coef(cvfit, s = "lambda.min")
lasso <- myCoefs@Dimnames[[1]][which(myCoefs != 0)]
lasso <- lasso[-1]
lasso
##  [1] "compactne_mean"         "concavity_mean"         "concave_points_mean"   
##  [4] "fractal_dimension_mean" "radius_se"              "texture_se"            
##  [7] "smoothness_se"          "compactne_se"           "symmetry_se"           
## [10] "fractal_dimension_se"   "radius_worst"           "texture_worst"         
## [13] "area_worst"             "smoothness_worst"       "concavity_worst"       
## [16] "concave_points_worst"   "symmetry_worst"

Lasso 回归和SVM-RFE比较

最后绘制Venn图比较两种方法筛选的特征变量，从图中发现，LASSO回归筛选的变量有12个，而SVM-RFE筛选的只有4个，对于我们临床上做标志物，找marker来说，数量缺少越好，方便后续的验证工作，因此，机器学习的方法很多，在筛选特征变量或者基因的时候可以选择多种方法，最后择优。

# 比较lasso和SVM-REF方法一找出的特征变量，画Venn图
overlap <- intersect(lasso, top.features[1:which.min(errors), "FeatureName"])  #交集
overlap
## [1] "fractal_dimension_se" "texture_worst"        "area_worst"
# 统计交叉基因有多少个
summary(lasso %in% top.features[1:which.min(errors), "FeatureName"])
##    Mode   FALSE    TRUE 
## logical      14       3
# 绘制venn图
grid.newpage()
venn.plot <- venn.diagram(list(LASSO = lasso, SVM_RFE = as.character(top.features[1:which.min(errors),
    "FeatureName"])), NULL, fill = c("#E64B35FF", "#4DBBD5FF"), alpha = c(0.5, 0.5),
    cex = 4, cat.fontface = 3, category.names = c("LASSO", "SVM_RFE"), main = "Overlap")
grid.draw(venn.plot)

Reference

Guyon, I., Weston, J., Barnhill, S. et al. Gene Selection for Cancer Classification using Support Vector Machines. Machine Learning 46, 389–422 (2002). https://doi.org/10.1023/A:1012487302797
Duan KB, Rajapakse JC, Wang H, Azuaje F. Multiple SVM-RFE for gene selection in cancer classification with expression data. IEEE Trans Nanobioscience. 2005 Sep;4(3):228-34. doi: 10.1109/tnb.2005.853657. PMID: 16220686.

桓峰基因，铸造成功的您！

未来桓峰基因公众号将不间断的推出单细胞系列生信分析教程，

敬请期待！！

桓峰基因官网正式上线，请大家多多关注，还有很多不足之处，大家多多指正！http://www.kyohogene.com/

桓峰基因和投必得合作，文章润色优惠85折，需要文章润色的老师可以直接到网站输入领取桓峰基因专属优惠券码：KYOHOGENE，然后上传，付款时选择桓峰基因优惠券即可享受85折优惠哦！https://www.topeditsci.com/

Matlab基于 支持向量机 递归特征消除 ( SVM _ RFE )的回归数据特征选择算法，matlab 代码，输出为选择的特征序号（Matlab完整程序和数据） Matlab基于 支持向量机 递归特征消除 ( SVM _ RFE )的回归数据特征选择算法，matlab 代码，输出为选择的特征序号（Matlab完整程序和数据） Matlab基于 支持向量机 递归特征消除 ( SVM _ RFE )的回归数据特征选择算法，matlab 代码，输出为选择的特征序号（Matlab完整程序和数据）经过特征选择后，保留特征的序号为： 126 160 161 163 165 166 237 239 240 370 -----------------------误差计算-------------------------- 评价结果如下所示：平均绝对误差MAE为：0.27933 均方误差MSE为： 0.15813 均方根误差RMSEP为： 0.39765 决定系数R^2为： 0.93392 剩余预测残差RPD为： 4.2631 平均绝对百分比误差MAPE为： 0.00 32 299 支持向量机 递归特征消除 (下文简称 SVM - RFE )是由Guyon等人在对癌症分类时提出来的，最初只能对两类数据进行特征提取。它是一种基于Embedded方法。 支持向量机 支持向量机 广泛用于基于 SVM - RFE 支持向量机 递归特征消除的回归数据特征选择算法，输出为选择的特征序号（Matlab完整程序和数据） Chinese: Options：可用的选项即表示的涵义如下　　-s svm 类型： SVM 设置类型(默认0) 　　0 -- C-SVC 　　1 --v-SVC 　　2 – 一类 SVM 3 -- e -SVR 　　4 -- v-SVR 　　-t 核函数类型：核函数设置类型(默认2) 　　0 – 线性：u'v 　　1 – 多项式：(r*u'v + coef0)^degree 　　2 – RBF函数：exp(-r|u-v|^2) 　　3 –sigmoid：tanh(r*u'v + coef0) 经过特征选择后，保留特征的序号为： 126 160 161 163 165 166 237 239 240 370 评价结果如下所示：平均绝对误差MAE为：0.27933 均方误差MSE为： 0.15813 均方根误差RMSEP为： 0.39765 决定系数R^2为： 0.93392 剩余预测残差RPD为： 4.2631 平均绝对百分比误差MAPE为： 0.00 32 299 上文复现了DEG和WGCNA之后，与原文章有出入，最后我借用sangerbox工具取交集，选出210个gene进行LASSO和 SVM 分析，继续筛选核心基因，其实可以使用更多的 机器学习 办法，如随机森林和Xgboost。LASSO回归结果和原文差别不小，我注意到原文四个基因都在那210个基因中，但是LASSO只选出了6个，其中只有2个基因和原文一致，LASSO具体参数尚不清楚。取前38个基因和LASSO取交集，只有SOCS2和FOSB和原文一致， SVM 估计参数也不一样，恳请指正。 SVM - RFE 与 SVM - RFE CV都是用于对特征进行缩减，用于找到数目最优的特征数。 RFE CV基于 RFE 基础上，添加了交叉验证，使得在每个step中，都可以对现有的特征数目进行评估，以确定比较好的数目。 RFE 中需要指定n_features_to_select，而 RFE CV中需要指定min_features_to_select，二者虽然都要运行到指定挑选特征数目时才停止筛选，但 RFE CV的每步交叉验... 文献标题:Gene Selection for Cancer Classifification using Support Vector Machine s 1. SVM 递归特征消除（ SVM - RFE ） SVM - RFE 是一种以权重大小作为排序标准的 REFE 应用。 1.1 算法流程输入：训练样本集幸存特征子集：特征排序列表：开始循环直到 s = []：限制训练数据集从而得到一个良好的特征指标采用 SVM 进行分类（对偶方式）获得值计算... 基于 支持向量机 递归特征消除 ( SVM _ RFE )的回归数据特征选择算法，matlab 代码，输出为选择的特征序号。评价指标包括:R2、MAE、MSE、RMSE和MAPE等，代码质量极高，方便学习和替换数据。

简 介