生存模型的構(gòu)建方法包括：
1. Lasso回歸；
2. Cox多因素回歸扔役；
3. 隨機(jī)森林帆喇；
4. 支持向量機(jī)

可以把log_rank_test或cox篩選出的基因單獨(dú)做模型構(gòu)建，也可以取交集之后再做模型構(gòu)建亿胸。

1. Lasso回歸（機(jī)器學(xué)習(xí)算法）

目的：從若干個(gè)基因中挑選真正對(duì)生存有影響的基因坯钦。
Lasso回歸可以對(duì)這些基因進(jìn)行統(tǒng)計(jì)和打分，從而挑出關(guān)鍵基因侈玄。

1.1 準(zhǔn)備輸入數(shù)據(jù)（表達(dá)矩陣數(shù)據(jù)和臨床信息數(shù)據(jù)）

load("TCGA-KIRC_sur_model.Rdata")
ls()
exprSet[1:4,1:4]
meta[1:4,1:4]

表達(dá)矩陣數(shù)據(jù)

臨床信息數(shù)據(jù)

1.2 構(gòu)建lasso回歸模型

輸入數(shù)據(jù)是表達(dá)矩陣（僅含tumor樣本）和每個(gè)病人對(duì)應(yīng)的生死（順序必須一致）婉刀。

x=t(exprSet) #轉(zhuǎn)換成基因在列
y=meta$event #結(jié)局
library(glmnet)

• 1.2.1 挑選合適的λ值

Lambda 是構(gòu)建模型的重要參數(shù)。他的大小關(guān)系著模型選擇的基因個(gè)數(shù)

#調(diào)優(yōu)參數(shù)
set.seed(1006) 
#??不設(shè)置的話序仙，每次運(yùn)行時(shí)后面的結(jié)果（選出的基因）會(huì)不斷變化突颊。這也就是說(shuō)，如果文獻(xiàn)中使用了Lasso回歸潘悼，是沒(méi)有辦法復(fù)現(xiàn)結(jié)果的律秃。
cv_fit <- cv.glmnet(x=x, y=y) #cv.glmnet()就是一個(gè)調(diào)優(yōu)參數(shù)的過(guò)程
plot(cv_fit)

圖的橫軸是選擇的λ值，圖最上方的數(shù)字代表選擇對(duì)應(yīng)的λ值時(shí)治唤，模型所使用的基因數(shù)棒动。
兩條虛線分別指示了兩個(gè)特殊的λ值,一個(gè)是lambda.min，一個(gè)是lambda.1se宾添，這兩個(gè)值之間的lambda都認(rèn)為是可用的船惨。lambda.1se構(gòu)建的模型最簡(jiǎn)單，即使用的基因數(shù)量少缕陕，而lambda.min則準(zhǔn)確率更高一點(diǎn)粱锐，使用的基因數(shù)量更多一點(diǎn)。

#系數(shù)圖
fit <- glmnet(x=x, y=y) #glmnet是構(gòu)建模型的
plot(fit,xvar = "lambda")

這張圖中的每一條線代表一個(gè)基因榄檬。橫坐標(biāo)和上面那張圖一樣是log Lambda卜范，縱坐標(biāo)是系數(shù)。從左往右看鹿榜，系數(shù)一直在0附近的，不管λ是幾都不會(huì)被選到锦爵。
這兩張圖就是Lasso回歸最經(jīng)典的兩張圖舱殿。

• 1.2.2 用這兩個(gè)λ值重新建模

兩個(gè)都構(gòu)建一下，再比較哪個(gè)模型更好险掀。實(shí)際操作也可直接選擇lambda.min到lambda.1se中的任意一個(gè)λ值構(gòu)建模型沪袭，λ值的選擇不是固定的。

model_lasso_min <- glmnet(x=x, y=y,lambda=cv_fit$lambda.min) # 得到的模型是一個(gè)列表
model_lasso_1se <- glmnet(x=x, y=y,lambda=cv_fit$lambda.1se) # 得到的模型是一個(gè)列表
View(model_lasso_1se) #從中選一個(gè)查看一下格式

選中的基因與系數(shù)存放于模型的子集beta中樟氢，beta是一個(gè)稀疏矩陣冈绊，用到的基因有一個(gè)s0值（系數(shù)值）侠鳄，沒(méi)用的基因只記錄了“.”（系數(shù)是0），所以可以用下面代碼挑出用到的基因死宣。

head(model_lasso_min$beta,20)
choose_gene_min=rownames(model_lasso_min$beta)[as.numeric(model_lasso_min$beta)!=0] #as.numeric不等于0的基因就是有系數(shù)的也就是被選中了的基因伟恶。
choose_gene_1se=rownames(model_lasso_1se$beta)[as.numeric(model_lasso_1se$beta)!=0]
length(choose_gene_min)
# [1] 35 ##選到了35個(gè)基因
length(choose_gene_1se)
# [1] 11 ##選到了11個(gè)基因
save(choose_gene_min,file = "lasso_choose_gene_min.Rdata")

• 1.2.3 模型預(yù)測(cè)
  使用predict函數(shù)進(jìn)行模型預(yù)測(cè)（提供模型和數(shù)據(jù)就可以得到預(yù)測(cè)結(jié)果）
  newx參數(shù)是預(yù)測(cè)對(duì)象。輸出結(jié)果lasso.prob是一個(gè)矩陣毅该，第一列是min的預(yù)測(cè)結(jié)果博秫，第二列是1se的預(yù)測(cè)結(jié)果，預(yù)測(cè)結(jié)果是概率眶掌，或者說(shuō)百分比挡育，不是絕對(duì)的0和1。
  將每個(gè)樣本的生死和預(yù)測(cè)結(jié)果放在一起朴爬，直接cbind即可即寒。

lasso.prob <- predict(cv_fit, newx=x, s=c(cv_fit$lambda.min,cv_fit$lambda.1se) ) 
re=cbind(y ,lasso.prob)
head(re)
#                              y         1         2
# TCGA-A3-3307-01A-01T-0860-13 0 0.1151638 0.2300976
# TCGA-A3-3308-01A-02R-1324-13 0 0.4020126 0.3591050
# TCGA-A3-3311-01A-02R-1324-13 1 0.2906333 0.2954458
# TCGA-A3-3313-01A-02R-1324-13 1 0.3831590 0.3456060
# TCGA-A3-3316-01A-01T-0860-13 0 0.3284941 0.3113726
# TCGA-A3-3317-01A-01T-0860-13 0 0.3832764 0.3342198
## y是每個(gè)病人，1和2是給的lamda值是min和1se的時(shí)候召噩，模型給每個(gè)病人計(jì)算出來(lái)的預(yù)測(cè)值評(píng)分是多少
re=as.data.frame(re)
colnames(re)=c('event','prob_min','prob_1se')
#將真實(shí)結(jié)局event和預(yù)測(cè)值合并在一起
re$event=as.factor(re$event)

• 1.2.4 模型評(píng)估（繪制ROC曲線）

計(jì)算AUC取值范圍在0.5-1之間母赵，越接近于1越好◎汲＃可以根據(jù)預(yù)測(cè)結(jié)果繪制ROC曲線市咽。

library(ROCR)
library(caret)
#min
pred_min <- prediction(re[,2], re[,1])
auc_min = performance(pred_min,"auc")@y.values[[1]]
perf_min <- performance(pred_min,"tpr","fpr")
plot(perf_min,colorize=FALSE, col="blue") 
lines(c(0,1),c(0,1),col = "gray", lty = 4 )
text(0.8,0.2, labels = paste0("AUC = ",round(auc_min,3)))
#1se
pred_1se <- prediction(re[,3], re[,1])
auc_1se = performance(pred_1se,"auc")@y.values[[1]]
perf_1se <- performance(pred_1se,"tpr","fpr")
plot(perf_1se,colorize=FALSE, col="red") 
lines(c(0,1),c(0,1),col = "gray", lty = 4 )
text(0.8,0.2, labels = paste0("AUC = ",round(auc_1se,3)))

把兩個(gè)模型的圖畫一起

plot(perf_min,colorize=FALSE, col="blue") 
plot(perf_1se,colorize=FALSE, col="red",add = T) 
lines(c(0,1),c(0,1),col = "gray", lty = 4 )
text(0.8,0.3, labels = paste0("AUC = ",round(auc_min,3)),col = "blue")
text(0.8,0.2, labels = paste0("AUC = ",round(auc_1se,3)),col = "red")

一般畫成正方形，可以自己調(diào)整

ggplot繪圖（和上面的圖含義一樣抵蚊，ggplot2畫的圖更好看施绎，繪圖也更靈活）

tpr_min = performance(pred_min,"tpr")@y.values[[1]]
tpr_1se = performance(pred_1se,"tpr")@y.values[[1]]
dat = data.frame(tpr_min = perf_min@y.values[[1]],
                 fpr_min = perf_min@x.values[[1]],
                 tpr_1se = perf_1se@y.values[[1]],
                 fpr_1se = perf_1se@x.values[[1]])
library(ggplot2)
ggplot() + 
  geom_line(data = dat,aes(x = fpr_min, y = tpr_min),color = "blue") + 
  geom_line(aes(x=c(0,1),y=c(0,1)),color = "grey")+
  theme_bw()+
  annotate("text",x = .75, y = .25,
           label = paste("AUC of min = ",round(auc_min,2)),color = "blue")+
  scale_x_continuous(name  = "fpr")+
  scale_y_continuous(name = "tpr")


ggplot() + 
  geom_line(data = dat,aes(x = fpr_min, y = tpr_min),color = "blue") + 
  geom_line(data = dat,aes(x = fpr_1se, y = tpr_1se),color = "red")+
  geom_line(aes(x=c(0,1),y=c(0,1)),color = "grey")+
  theme_bw()+
  annotate("text",x = .75, y = .25,
           label = paste("AUC of min = ",round(auc_min,2)),color = "blue")+
  annotate("text",x = .75, y = .15,label = paste("AUC of 1se = ",round(auc_1se,2)),color = "red")+
  scale_x_continuous(name  = "fpr")+
  scale_y_continuous(name = "tpr")

• 1.2.5 切割數(shù)據(jù)構(gòu)建模型并預(yù)測(cè)

對(duì)數(shù)據(jù)進(jìn)行切割，一組做訓(xùn)練集贞绳，使用表達(dá)矩陣和臨床信息用來(lái)構(gòu)建模型谷醉。一組做測(cè)試集，輸入表達(dá)矩陣來(lái)對(duì)結(jié)局進(jìn)行預(yù)測(cè)冈闭，測(cè)序構(gòu)建的模型預(yù)測(cè)的結(jié)果俱尼。

用R包c(diǎn)aret切割數(shù)據(jù)，生成的結(jié)果是一組代表列數(shù)的數(shù)字萎攒，用這些數(shù)字來(lái)給表達(dá)矩陣和meta取子集即可遇八。

library(caret) #機(jī)器學(xué)習(xí)R包，可以比較科學(xué)的拆分?jǐn)?shù)據(jù)耍休。
set.seed(12345679)
sam<- createDataPartition(meta$event, p = .5,list = FALSE) #根據(jù)event來(lái)切分刃永，一半一半和7比3都可以
head(sam)

train <- exprSet[,sam] #切分表達(dá)矩陣，sam是一組羊精，非sam是另一組斯够。
test <- exprSet[,-sam]
train_meta <- meta[sam,] #切分臨床信息
test_meta <- meta[-sam,]

#看一下，不要讓臨床信息差的太多
prop.table(table(train_meta$stage))
prop.table(table(test_meta$stage)) 
prop.table(table(test_meta$race)) 
prop.table(table(train_meta$race)) 

切割后的train數(shù)據(jù)集建模

#計(jì)算lambda
x = t(train)
y = train_meta$event
cv_fit <- cv.glmnet(x=x, y=y)
plot(cv_fit)
#構(gòu)建模型
model_lasso_min <- glmnet(x=x, y=y,lambda=cv_fit$lambda.min)
model_lasso_1se <- glmnet(x=x, y=y,lambda=cv_fit$lambda.1se)
#挑出基因
head(model_lasso_min$beta)
choose_gene_min=rownames(model_lasso_min$beta)[as.numeric(model_lasso_min$beta)!=0]
choose_gene_1se=rownames(model_lasso_1se$beta)[as.numeric(model_lasso_1se$beta)!=0]
length(choose_gene_min)
# [1] 18
length(choose_gene_1se)
# [1] 3

模型預(yù)測(cè)
用訓(xùn)練集構(gòu)建模型，預(yù)測(cè)測(cè)試集的生死读规，注意newx參數(shù)變了抓督。

lasso.prob <- predict(cv_fit, newx=t(test), s=c(cv_fit$lambda.min,cv_fit$lambda.1se) )
re=cbind(event = test_meta$event ,lasso.prob)
re=as.data.frame(re)
colnames(re)=c('event','prob_min','prob_1se')
re$event=as.factor(re$event)
head(re)

再畫ROC曲線

library(ROCR)
library(caret)
# 訓(xùn)練集模型預(yù)測(cè)測(cè)試集
#min
pred_min <- prediction(re[,2], re[,1])
auc_min = performance(pred_min,"auc")@y.values[[1]]
perf_min <- performance(pred_min,"tpr","fpr")

#1se
pred_1se <- prediction(re[,3], re[,1])
auc_1se = performance(pred_1se,"auc")@y.values[[1]]
perf_1se <- performance(pred_1se,"tpr","fpr")

tpr_min = performance(pred_min,"tpr")@y.values[[1]]
tpr_1se = performance(pred_1se,"tpr")@y.values[[1]]
dat = data.frame(tpr_min = perf_min@y.values[[1]],
                 fpr_min = perf_min@x.values[[1]],
                 tpr_1se = perf_1se@y.values[[1]],
                 fpr_1se = perf_1se@x.values[[1]])

ggplot() + 
  geom_line(data = dat,aes(x = fpr_min, y = tpr_min),color = "blue") + 
  geom_line(data = dat,aes(x = fpr_1se, y = tpr_1se),color = "red")+
  geom_line(aes(x=c(0,1),y=c(0,1)),color = "grey")+
  theme_bw()+
  annotate("text",x = .75, y = .25,
           label = paste("AUC of min = ",round(auc_min,2)),color = "blue")+
  annotate("text",x = .75, y = .15,label = paste("AUC of 1se = ",round(auc_1se,2)),color = "red")+
  scale_x_continuous(name  = "fpr")+
  scale_y_continuous(name = "tpr")

2. Cox多因素回歸

如果Lasso回歸挑出的基因數(shù)目還是太多，就可以通過(guò)Cox多因素回歸再進(jìn)行篩選束亏。
使用Lasso回歸挑出的基因作為Cox多因素回歸的輸入數(shù)據(jù)铃在，使用逐步回歸法去挑選可選范圍內(nèi)最好的模型。（通常來(lái)說(shuō)做數(shù)據(jù)挖掘的文章構(gòu)建模型的話最后一步都是多因素Cox）

2.1 準(zhǔn)備輸入數(shù)據(jù)

if(!require(My.stepwise))install.packages("My.stepwise")
load("TCGA-KIRC_sur_model.Rdata") #KIRC整理好的生存數(shù)據(jù)
load("lasso_choose_gene_min.Rdata") #Lasso回歸選出的35個(gè)基因（這里沒(méi)有區(qū)分lnc和mRNA）

2.2 構(gòu)建coxph模型

將用于建模的基因（例如lasso回歸選中的基因）從表達(dá)矩陣中取出來(lái)枪汪，可作為列添加在meta表格的后面,組成的數(shù)據(jù)框賦值給dat涌穆。

library(stringr)
e=t(exprSet[choose_gene_min,]) #從矩陣中取出Lasso回歸選出的35個(gè)基因
colnames(e)= str_replace_all(colnames(e),"-","_") #因?yàn)閙iRNA的名字中有減號(hào)-，而cox的公式中有加號(hào)雀久，會(huì)造成干擾宿稀，因此需要將基因名中的-替換成_。
dat=cbind(meta,e)

dat$gender=as.numeric(factor(dat$gender)) #把levels轉(zhuǎn)換成數(shù)字作為輸入
dat$stage=as.numeric(factor(dat$stage))
colnames(dat) #dat矩陣行是樣本赖捌，列是感興趣的臨床信息和挑出的35個(gè)基因
#   [1] "ID"             "event"          "death"          "last_followup"  "race"          
#  [6] "age"            "gender"         "stage"          "time"           "age_group"     
# [11] "hsa_mir_101_1"  "hsa_mir_1179"   "hsa_mir_1185_1" "hsa_mir_1248"   "hsa_mir_1269"  
# [16] "hsa_mir_1277"   "hsa_mir_1305"   "hsa_mir_130a"   "hsa_mir_130b"   "hsa_mir_144"   
# [21] "hsa_mir_149"    "hsa_mir_181b_2" "hsa_mir_181c"   "hsa_mir_18a"    "hsa_mir_223"   
# [26] "hsa_mir_2276"   "hsa_mir_27b"    "hsa_mir_28"     "hsa_mir_3149"   "hsa_mir_34c"   
# [31] "hsa_mir_3613"   "hsa_mir_3614"   "hsa_mir_365_2"  "hsa_mir_3651"   "hsa_mir_3667"  
# [36] "hsa_mir_3676"   "hsa_mir_376b"   "hsa_mir_3917"   "hsa_mir_548q"   "hsa_mir_612"   
# [41] "hsa_mir_627"    "hsa_mir_676"    "hsa_mir_9_2"    "hsa_mir_939"    "hsa_mir_95"    

逐步回歸法構(gòu)建最優(yōu)模型（從若干個(gè)因素中挑選出更關(guān)鍵的）

library(survival)
library(survminer)
library(My.stepwise)
vl <- colnames(dat)[c(6,7,8,11:ncol(dat))] #選擇感興趣的輸入因素祝沸，這里選了列名中從age開(kāi)始到最后一列的35個(gè)因素作為輸入
My.stepwise.coxph(Time = "time",
                  Status = "event",
                  variable.list = vl, #從這些因素中挑選更關(guān)鍵的因素，輸出結(jié)果中最后一個(gè)模型是最好的
                  data = dat)

使用輸出結(jié)果里的最后一個(gè)模型（逐步回歸法經(jīng)過(guò)若干次迭代后選出的最好的模型）

# 復(fù)制輸出結(jié)果中的最后一個(gè)公式
model = coxph(formula = Surv(time, event) ~ stage + hsa_mir_223 + age + 
    hsa_mir_34c + hsa_mir_3917 + hsa_mir_3651 + hsa_mir_144 + 
    hsa_mir_2276 + hsa_mir_3667 + hsa_mir_3149 + hsa_mir_627 + 
    hsa_mir_101_1, data = dat)

逐步回歸法最重要的是選出了～后面的這些因素越庇，作為它認(rèn)為的最優(yōu)因素組合來(lái)構(gòu)建模型罩锐，～后的因素也可以通過(guò)其它方式挑選，公式還是一樣的寫法卤唉。

2.3 模型可視化-森林圖

ggforest(model,data = dat)

每一行是一個(gè)因素涩惑，n=516是用來(lái)建模的樣本數(shù)，第三列和右邊的圖都是表示SR值桑驱，以1為界竭恬，小于1的是保護(hù)因素，大于1的是危險(xiǎn)因素熬的。最后一列是p值痊硕。圖中有一些p值不顯著的也被選中了，但這不影響這是逐步回歸法選出的一個(gè)最優(yōu)模型押框。左下角的Events: 158是說(shuō)有158個(gè)人結(jié)局是1岔绸。Global p-value: 3.0657e-35盒揉，這個(gè)是全局的p值道媚，顯示這個(gè)模型還不錯(cuò)最域。AIC是衡量模型精確度和簡(jiǎn)潔度的標(biāo)準(zhǔn)镀脂，這個(gè)值越小越好薄翅。逐步回歸法選擇模型的標(biāo)準(zhǔn)就是AIC鼎天。Concordance index和ROC曲線的AUC值是相同的取值范圍相同的作用但荤，只是意義不同，它和AUC一樣是評(píng)價(jià)模型好壞的一個(gè)指標(biāo)腹躁。模型不好的話Concordance index就更接近于0.5垒手，好的話就更接近于1鳖悠。

森林圖繪圖更多參考：http://www.reibang.com/p/58c90b2f3910

2.4 模型預(yù)測(cè)

fp <- predict(model,newdata = dat)
library(Hmisc) #使用Hmisc包中的函數(shù)計(jì)算C Index
options(scipen=200)
with(dat,rcorr.cens(fp,Surv(time, event))) #用1-返回的C Index值是我們所說(shuō)的C Index，也就是上面森林圖中標(biāo)出的0.81
#         C Index             Dxy            S.D.               n         missing      uncensored 
#      0.19234306     -0.61531388      0.03361968    516.00000000      0.00000000    158.00000000 
#  Relevant Pairs      Concordant       Uncertain 
#  89174.00000000  17152.00000000 176550.00000000 

C-index用于計(jì)算生存分析中的COX模型預(yù)測(cè)值與真實(shí)之間的區(qū)分度（discrimination）反砌，也稱為Harrell's concordanceindex宴树。C-index在0.5-1之間蠢莺。0.5為完全不一致,說(shuō)明該模型沒(méi)有預(yù)測(cè)作用,1為完全一致,說(shuō)明該模型預(yù)測(cè)結(jié)果與實(shí)際完全一致。

2.5 切割數(shù)據(jù)構(gòu)建模型并預(yù)測(cè)

• 2.5.1 切割數(shù)據(jù)
  用R包c(diǎn)aret切割數(shù)據(jù)辕翰，生成的結(jié)果是一組代表列數(shù)的數(shù)字夺衍，用這些數(shù)字來(lái)給表達(dá)矩陣和meta取子集即可。

library(caret)
set.seed(12345679)
sam<- createDataPartition(meta$event, p = .5,list = FALSE)
train <- exprSet[,sam]
test <- exprSet[,-sam]
train_meta <- meta[sam,]
test_meta <- meta[-sam,]

• 2.5.2 切割后的train數(shù)據(jù)集建模
  和上面的建模方法一樣喜命。

e=t(train[choose_gene_min,])
colnames(e)= str_replace_all(colnames(e),"-","_")
dat=cbind(train_meta,e)

dat$gender=as.numeric(factor(dat$gender))
dat$stage=as.numeric(factor(dat$stage))
colnames(dat)

#install.packages("My.stepwise")
# library(My.stepwise)
# vl <- colnames(dat)[c(6,7,8,11:ncol(dat))]
# My.stepwise.coxph(Time = "time",
#                   Status = "event",
#                   variable.list = vl,
#                   data = dat)

model = coxph(formula = Surv(time, event) ~ stage + hsa_mir_223 + age + 
    hsa_mir_34c + hsa_mir_181b_2 + hsa_mir_3614, data = dat)

• 2.5.3 模型可視化

ggforest(model, data =dat)

• 2.5.4 用切割后的數(shù)據(jù)test數(shù)據(jù)集驗(yàn)證模型

e=t(test[choose_gene_min,])
colnames(e)= str_replace_all(colnames(e),"-","_")
test_dat=cbind(test_meta,e)
test_dat$gender=as.numeric(factor(test_dat$gender))
test_dat$stage=as.numeric(factor(test_dat$stage))

fp <- predict(model,newdata = test_dat)
library(Hmisc)
with(test_dat,rcorr.cens(fp,Surv(time, event)))
#        C Index            Dxy           S.D.              n        missing 
#     0.22565299    -0.54869403     0.04919064   258.00000000     0.00000000 
#     uncensored Relevant Pairs     Concordant      Uncertain 
#    75.00000000 21440.00000000  4838.00000000 44858.00000000 

3 隨機(jī)森林

3.1 準(zhǔn)備輸入數(shù)據(jù)

輸入數(shù)據(jù)是腫瘤樣本表達(dá)矩陣exprSet和臨床信息meta

load("TCGA-KIRC_sur_model.Rdata")
library(randomForest)
library(ROCR)
library(genefilter)
library(Hmisc)

3.2 構(gòu)建隨機(jī)森林模型

輸入數(shù)據(jù)是表達(dá)矩陣(僅含tumor樣本)和對(duì)應(yīng)的生死沟沙。（和Lasso回歸一樣）

x=t(exprSet)
y=meta$event
#構(gòu)建模型，一個(gè)叫randomForest的函數(shù),運(yùn)行時(shí)間很長(zhǎng),存Rdata跳過(guò)
tmp_rf="TCGA_KIRC_miRNA_rf_output.Rdata"
if(!file.exists(tmp_rf)){
  rf_output=randomForest(x=x, y=y,importance = TRUE, ntree = 10001, proximity=TRUE )
  save(rf_output,file = tmp_rf)
}
load(file = tmp_rf)
#top30的基因
varImpPlot(rf_output, type=2, n.var=30, scale=FALSE, 
           main="Variable Importance (Gini) for top 30 predictors",cex =.7)

隨機(jī)森林算法篩選基因不像Lasso回歸那么絕對(duì)壁榕，不是直接分是用到了還是沒(méi)用到矛紫，隨機(jī)森林只能幫我們把基因的重要性進(jìn)行排名 ，要選前多少個(gè)基因是我們自己定的牌里。（這里選了前30颊咬，如上圖）

rf_importances=importance(rf_output, scale=FALSE)
head(rf_importances)
#                   %IncMSE IncNodePurity
# hsa-let-7a-1 1.852761e-04     0.1787383
# hsa-let-7a-2 2.167420e-04     0.1916623
# hsa-let-7a-3 2.218169e-04     0.1858544
# hsa-let-7b   7.399404e-05     0.1628646
# hsa-let-7c   7.658155e-05     0.1635053
# hsa-let-7d   1.974099e-04     0.2382185
#解釋量top30的基因,和圖上是一樣的，從下往上看牡辽。
choose_gene=rownames(tail(rf_importances[order(rf_importances[,2]),],30)) #選前30個(gè)基因喳篇，這個(gè)數(shù)值可變

3.3 模型預(yù)測(cè)和評(píng)估

rf.prob <- predict(rf_output, x) #rf.prob是預(yù)測(cè)值，越接近于0就表示模型預(yù)測(cè)這個(gè)人活著态辛，越接近于1就表示模型預(yù)測(cè)這個(gè)人死了麸澜。
re=cbind(y ,rf.prob) #預(yù)測(cè)值和真實(shí)值放在一起查看比較
head(re)
#                              y   rf.prob
# TCGA-A3-3307-01A-01T-0860-13 0 0.1364447
# TCGA-A3-3308-01A-02R-1324-13 0 0.1793771
# TCGA-A3-3311-01A-02R-1324-13 1 0.6709712
# TCGA-A3-3313-01A-02R-1324-13 1 0.7742376
# TCGA-A3-3316-01A-01T-0860-13 0 0.2035863
# TCGA-A3-3317-01A-01T-0860-13 0 0.1619938

ROC曲線（隨機(jī)森林如果不做另外一個(gè)數(shù)據(jù)的驗(yàn)證，它的AUC值是1奏黑，所以繪制ROC曲線沒(méi)有什么意義）

library(ROCR)
#library(caret)

pred <- prediction(re[,2], re[,1])
auc = performance(pred,"auc")@y.values[[1]]
auc
# [1] 1

3.4 切割數(shù)據(jù)構(gòu)建模型并預(yù)測(cè)

• 3.4.1 切割數(shù)據(jù)

用R包c(diǎn)aret切割數(shù)據(jù)炊邦，生成的結(jié)果是一組代表列數(shù)的數(shù)字编矾，用這些數(shù)字來(lái)給表達(dá)矩陣和meta取子集即可。

library(caret)
set.seed(12345679)
sam<- createDataPartition(meta$event, p = .5,list = FALSE)

train <- exprSet[,sam]
test <- exprSet[,-sam]
train_meta <- meta[sam,]
test_meta <- meta[-sam,]

• 3.4.2 切割后的train數(shù)據(jù)集建模

和上面的建模方法一樣铣耘。

#建模
x = t(train)
y = train_meta$event
tmp_rf="TCGA_KIRC_miRNA_t_rf_output.Rdata"
if(!file.exists(tmp_rf)){
  rf_output=randomForest(x=x, y=y,importance = TRUE, ntree = 10001, proximity=TRUE )
  save(rf_output,file = tmp_rf)
}
load(file = tmp_rf)

choose_gene=rownames(tail(rf_importances[order(rf_importances[,2]),],30))
head(choose_gene)
# [1] "hsa-mir-511-1"  "hsa-mir-155"    "hsa-mir-409"    "hsa-mir-1185-1"
# [5] "hsa-mir-1277"   "hsa-mir-149"

3.5 模型預(yù)測(cè)

用訓(xùn)練集構(gòu)建模型洽沟，預(yù)測(cè)測(cè)試集的生死。

x = t(test)
y = test_meta$event
rf.prob <- predict(rf_output, x)
re=cbind(y ,rf.prob)
re=as.data.frame(re)
colnames(re)=c('event','prob')
re$event=as.factor(re$event)

再看AUC值蜗细。

library(ROCR)
# 訓(xùn)練集模型預(yù)測(cè)測(cè)試集
pred <- prediction(re[,2], re[,1])
auc= performance(pred,"auc")@y.values[[1]]
auc
# [1] 0.7121311

模型還可以。

4 支持向量機(jī)（Support Vector Machine, SVM）

前三種預(yù)測(cè)方法都是可以把建模時(shí)所 使用的基因挑出來(lái)怒详，也可以寫出公式炉媒，SVM比前三種方法要簡(jiǎn)單很多。

4.1 準(zhǔn)備輸入數(shù)據(jù)

load("TCGA-KIRC_sur_model.Rdata")
library(ROCR)
library(genefilter)
library(Hmisc)
library(e1071) #SVM的包名

4.2 構(gòu)建支持向量機(jī)模型

• 4.2.1 切割數(shù)據(jù)
  用R包c(diǎn)aret切割數(shù)據(jù)昆烁，生成的結(jié)果是一組代表列數(shù)的數(shù)字吊骤，用這些數(shù)字來(lái)給表達(dá)矩陣和meta取子集即可。

library(caret)
set.seed(12345679)
sam<- createDataPartition(meta$event, p = .5,list = FALSE)

train <- exprSet[,sam]
test <- exprSet[,-sam]
train_meta <- meta[sam,]
test_meta <- meta[-sam,]

• 4.2.2 train數(shù)據(jù)集建模

x=t(train)
y=as.factor(train_meta$event)
model = svm(x,y,kernel = "linear")
summary(model) 
# 
# Call:
# svm.default(x = x, y = y, kernel = "linear")
# 
# 
# Parameters:
#    SVM-Type:  C-classification 
#  SVM-Kernel:  linear 
#        cost:  1 
# 
# Number of Support Vectors:  176
# 
#  ( 108 68 )
# 
# 
# Number of Classes:  2 
# 
# Levels: 
#  0 1

• 4.2.3 模型預(yù)測(cè)
  用訓(xùn)練集構(gòu)建模型静尼，預(yù)測(cè)測(cè)試集的生死白粉。不同于其他模型，這個(gè)預(yù)測(cè)結(jié)果是分類變量鼠渺，直接預(yù)測(cè)生死鸭巴，而不是prob。

x=t(test)
y=as.factor(test_meta$event)
pred = predict(model, x)
table(pred,y)
#     y
# pred   0   1
#    0 142  47
#    1  41  28

上面的結(jié)果列的0和1是預(yù)測(cè)的存活和死亡拦盹，行的0和1是真實(shí)的存活和死亡鹃祖。預(yù)測(cè)值是0的有142+41個(gè)，預(yù)測(cè)值是1的有47+28個(gè)普舆。真實(shí)是0的有142+47個(gè)恬口，真實(shí)是1的有41+28個(gè)。41（false negative）和47（false positive）就屬于誤判沼侣。

代碼來(lái)自2021生信技能樹(shù)數(shù)據(jù)挖掘課