生存模型：
? Cox單因素分析 (前面一篇文章講生存分析的時(shí)候講了)
? Lasso回歸 （本篇文章）
? Cox多因素分析
? 隨機(jī)森林
? 支持向量機(jī)
1、lasso回歸：從一堆的基因中找到幾個(gè)關(guān)鍵的用來(lái)建模或者預(yù)測(cè)
輸入數(shù)據(jù)：表達(dá)矩陣（exprSet）和病人對(duì)應(yīng)的生死(meta$event)
判斷一下是否相等：這一步是至關(guān)重要的，如果樣本和病人ID順序不對(duì)，后面都是錯(cuò)的
> identical(stringr::str_sub(colnames(exprSet),1,12),meta$ID)
[1] TRUE
這個(gè)包我沒(méi)有下載成功收捣。陨收。势告。岸啡。不知道原因是什么，然后我一直搜索赫编。巡蘸。。
我覺(jué)得如果不是網(wǎng)絡(luò)問(wèn)題擂送，肯定就是鏡像問(wèn)題悦荒，不知道清華鏡像最近咋了
然后我就換了厲害無(wú)比的阿里云鏡像：
劃重點(diǎn)！｀诙帧！！??
options("repos" = c(CRAN="https://mirrors.aliyun.com/CRAN/"))
結(jié)果一頓操作猛如虎就成功了枉氮。本來(lái)都打算放棄了允粤。。问芬。
下面正式開(kāi)始：
這里是舉例子悦析，所以只計(jì)算了10個(gè)λ值，解釋一下輸出結(jié)果三列的意思此衅。
- Df 是自由度
- 列%Dev代表了由模型解釋的殘差的比例强戴，對(duì)于線性模型來(lái)說(shuō)就是模型擬合的 R^2(R-squred)亭螟。它在0和1之間，越接近1說(shuō)明模型的表現(xiàn)越好骑歹，如果是0预烙，說(shuō)明模型的預(yù)測(cè)結(jié)果還不如直接把因變量的均值作為預(yù)測(cè)值來(lái)的有效。
- Lambda 是構(gòu)建模型的重要參數(shù)道媚。
解釋的殘差百分比越高越好扁掸，但是構(gòu)建模型使用的基因的數(shù)量也不能太多，需要取一個(gè)折中值衰琐。
#### 2.1挑選合適的λ值
計(jì)算1000個(gè)也糊，畫(huà)圖，篩選表現(xiàn)最好的λ值
cv_fit <- cv.glmnet(x=x, y=y, nlambda = 1000,alpha = 1)
plot(cv_fit)
兩條虛線分別指示了兩個(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)钞馁。
#### 2.2 用這兩個(gè)λ值重新建模
model_lasso_min <- glmnet(x=x, y=y, alpha = 1, lambda=cv_fit$lambda.min)
model_lasso_1se <- glmnet(x=x, y=y, alpha = 1, lambda=cv_fit$lambda.1se)
這兩個(gè)值體現(xiàn)在參數(shù)lambda上。有了模型匿刮，可以將篩選的基因挑出來(lái)了僧凰。所有基因存放于模型的子集beta中，用到的基因有一個(gè)s0值熟丸，沒(méi)用的基因只記錄了“.”训措，所以可以用下面代碼挑出用到的基因。
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)
length(choose_gene_1se)
### 3.模型預(yù)測(cè)和評(píng)估
#### 3.1自己預(yù)測(cè)自己：自我驗(yàn)證
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, #predict是預(yù)測(cè)铐炫，是一個(gè)泛型函數(shù)
newx=x , #這個(gè)是預(yù)測(cè)對(duì)象垒手，cv_fit是預(yù)測(cè)數(shù)據(jù)
s=c(cv_fit$lambda.min,cv_fit$lambda.1se) )
re=cbind(y ,lasso.prob) #cbind是按列合并
head(re)
#### 3.2 箱線圖
對(duì)預(yù)測(cè)結(jié)果進(jìn)行可視化。以實(shí)際的生死作為分組倒信，畫(huà)箱線圖整體上查看預(yù)測(cè)結(jié)果淫奔。
re=as.data.frame(re)
colnames(re)=c('event','prob_min','prob_1se')
#上面一行的代碼意思是把列名一次改為'event','prob_min','prob_1se'
re$event=as.factor(re$event) #作為因子來(lái)變成分類變量
library(ggpubr)
p1 = ggboxplot(re, x = "event", y = "prob_min",
color = "event", palette = "jco",
add = "jitter")+ stat_compare_means()
p2 = ggboxplot(re, x = "event", y = "prob_1se",
color = "event", palette = "jco",
add = "jitter")+ stat_compare_means()
library(patchwork)
p1+p2
可以看到，真實(shí)結(jié)果是死(0)的樣本堤结，預(yù)測(cè)的值就小一點(diǎn)(靠近0)唆迁，真實(shí)結(jié)果是活著(1)的樣本鸭丛，預(yù)測(cè)的值就大一點(diǎn)(靠近1)，整體上趨勢(shì)是對(duì)的唐责，但不是完全準(zhǔn)確鳞溉，模型是可用的.對(duì)比兩個(gè)λ值構(gòu)建的模型，差別不大鼠哥，1se的預(yù)測(cè)值準(zhǔn)確一點(diǎn)熟菲。
下面是ROC曲線：
#### 3.3 ROC曲線
計(jì)算AUC取值范圍在0.5-1之間，越接近于1越好朴恳〕保可以根據(jù)預(yù)測(cè)結(jié)果繪制ROC曲線。
library(ROCR)
library(caret)
# 自己預(yù)測(cè)自己
#min 下面整段復(fù)制運(yùn)行代碼就會(huì)出圖
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)))
#強(qiáng)迫癥選項(xiàng)于颖，想把兩個(gè)模型畫(huà)一起呆贿。
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")
#round()這個(gè)函數(shù)是取幾位小數(shù)，round(656563,3)就是取3位小數(shù)
AUC
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")
最后補(bǔ)充一個(gè)小小知識(shí)點(diǎn)
#設(shè)置隨機(jī)數(shù)種子森渐，這樣的話每次抽出來(lái)的數(shù)就是一樣的了
set.seed(1234)
sample(1:50,7)
