比较三种机器学习模型（随机森林，支持向量机，逻辑回归）的分类效果

这是2014年发表在nature communication上的一张图，大意如下：
a：作者想用181个假基因预测癌症亚型，将306个样本分成训练集（223个），验证集（83个）
b: 用训练集的数据（n=223）进行分析，采用五折交叉验证避免过拟合，选用随机森林，支持向量机，逻辑回归三种机器学习方法进行分类。
c: 三种机器学习方法的AUC曲线，发现SVM的效果最好。（AUC=0.962）
d: 用整个训练集（n=223）对验证集（n=83）进行分析，用SVM进行分类，AUC为0.922。
（原文题目：The Pan-Cancer analysis of pseudogene expression reveals biologically and clinically relevant tumour subtypes）

我把这个思路搬运到我的课题里，复现图如下：

相似度：75%。a和b是用AI画的，所以今天分享下c和d的代码

一 差异分析

假设我们有了这样的一个表达矩阵，叫Expr_all。以及表示表型数据的temp (character)。差异分析你可以用DESeq，但本次我就用T检验。

temp

##转置
Expr_all<-t(Expr_all)
##将每列变成numeric
Expr_all[,1]<-as.numeric(Expr_all[,1])
for(i in 2:535){Expr_all[,i]<-as.numeric(Expr_all[,i])}
##T检验
pvalue<-sapply(1:266,function(i){t.test(Expr_all[which(temp=="cancer"),i],Expr_all[which(temp=="normal"),i])[[3]]})
##将T检验得到的P值进行矫正，降低假阳性
pvalue1<-p.adjust(pvalue1)
##挑出P<0.05的基因
t_expr<-Expr_all[,which(pvalue1<0.05)]
t_expr<-as.data.frame(t_expr)
t_expr<-t(t_expr)
t_expr<-cbind(temp,t_expr)

得到的t_expr矩阵如下

t_expr

二 构建三种机器学习模型

1. 划分训练集和验证集

library(randomForest)
library(caret)
library(pROC)
library(caret)
inTrain<-createDataPartition(y=t_expr[,1],p=0.25,list=F)
test<-t_expr[inTrain,]
train<-t_expr[-inTrain,]

2. 构建五折交叉验证RF模型

folds<-createFolds(y=train[,1],k=5)
fc<-as.numeric()
mod_pre<-as.numeric()
for(i in 1:5){
  fold_test<-train[folds[[i]],]
  fold_train<-train[-folds[[i]],]
  model<-randomForest(temp~.,data=fold_train,proximity=T,importance=T)
  model_pre<-predict(model,newdata = fold_test,type="prob")
  fc<-append(fc,as.numeric(fold_test$temp))
  mod_pre<-append(mod_pre,model_pre[,1])
}
df<-cbind(fc,as.numeric(mod_pre))

3. 构建五折交叉验证SVM模型

fc<-as.numeric()
mod_pre<-as.numeric()
for(i in 1:5){
  fold_test<-t_expr[folds[[i]],]
  fold_train<-t_expr[-folds[[i]],]
  model<-svm(temp~.,data=fold_train,probability=T)
  model_pre<-predict(model,newdata = fold_test,decision.values = TRUE, probability = TRUE)
  fc<-append(fc,as.numeric(fold_test$temp))
  mod_pre<-append(mod_pre,as.numeric(attr(model_pre, "probabilities")[,2]))
}
df<-cbind(df,cbind(fc,mod_pre))

4. 构建五折交叉验证LR模型

max=0
num=0
fc<-as.numeric()
mod_pre<-as.numeric()
for(i in 1:5){
  fold_test<-train[folds[[i]],]
  fold_train<-train[-folds[[i]],]
  model<-glm(temp~.,family=binomial(link=logit),data=fold_train)
  model_pre<-predict(model,type='response',newdata=fold_test)
  fc<-append(fc,fold_test$temp)
  mod_pre<-append(mod_pre,as.numeric(model_pre))
}
df<-cbind(df,cbind(fc,mod_pre))

5. 画三者的AUC曲线

pdf("/Users/baiyunfan/desktop/ROC.pdf",height=6,width=6)
mycol <- c("slateblue","seagreen3","dodgerblue","firebrick1","lightgoldenrod","magenta","orange2")
x<-plot.roc(df[,1],df[,2],
            smooth=F,
            lwd=2,
            ylim=c(0,1),
            xlim=c(1,0),
            legacy.axes=T,
            main="",
            col=mycol[2])
x<-plot.roc(df[,3],df[,4],
            smooth=F,
            add=T,
            lwd=2,
            ylim=c(0,1),
            xlim=c(1,0),
            legacy.axes=T,
            main="",
            col=mycol[3])
x<-plot.roc(df[,5],df[,6],
            smooth=F,
            add=T,
            lwd=2,
            ylim=c(0,1),
            xlim=c(1,0),
            legacy.axes=T,
            main="",
            col=mycol[4])

legend.name <- c(paste("RF","AUC",0.9836,sep=" "),paste("SVM","AUC",0.9903,sep=" "),paste("LR","AUC",0.9986,sep=" "))
legend("bottomright", 
       legend=legend.name,
       col = mycol[2:4],
       lwd = 2,
       bty="n")
dev.off()

发现LR的AUC最高，于是单拿出来用整个train对test画

6. 用整个训练集对验证集建立LR模型以及AUC曲线

model<-glm(temp~.,data=train,family=binomial(link=logit))
model_pre<-predict(model,type='response',newdata=test)
pdf("/Users/baiyunfan/desktop/RF.pdf",height=6,width=6)
x<-plot.roc(test[,1],model_pre,
            smooth=F,
            lwd=2,
            ylim=c(0,1),
            xlim=c(1,0),
            legacy.axes=T,
            main="",
            col=mycol[4])
legend.name <- paste("LR","AUC",0.9557,sep=" ")
legend("bottomright", 
       legend=legend.name,
       col = mycol[4],
       lwd = 2,
       bty="n")
dev.off()

AUC掉下来了，0.9557，大家来探讨下为什么吧～