机器学习之Lasso回归

生存模型：
• Cox单因素分析 (前面一篇文章讲生存分析的时候讲了)
• Lasso回归（本篇文章）
• Cox多因素分析
• 随机森林
• 支持向量机

1、lasso回归：从一堆的基因中找到几个关键的用来建模或者预测

输入数据：表达矩阵（exprSet）和病人对应的生死(meta$event)

image.png

病人的ID

判断一下是否相等：这一步是至关重要的，如果样本和病人ID顺序不对，后面都是错的

image.png

> identical(stringr::str_sub(colnames(exprSet),1,12),meta$ID)
[1] TRUE

image.png

这个包我没有下载成功。。。。不知道原因是什么，然后我一直搜索。。。
我觉得如果不是网络问题，肯定就是镜像问题，不知道清华镜像最近咋了
然后我就换了厉害无比的阿里云镜像：

划重点！！！！😄

options("repos" = c(CRAN="https://mirrors.aliyun.com/CRAN/"))

结果一顿操作猛如虎就成功了。本来都打算放弃了。。。

image.png

下面正式开始：

这里是举例子，所以只计算了10个λ值，解释一下输出结果三列的意思。
- Df 是自由度
- 列%Dev代表了由模型解释的残差的比例，对于线性模型来说就是模型拟合的 R^2(R-squred)。它在0和1之间，越接近1说明模型的表现越好，如果是0，说明模型的预测结果还不如直接把因变量的均值作为预测值来的有效。
- Lambda 是构建模型的重要参数。
解释的残差百分比越高越好，但是构建模型使用的基因的数量也不能太多，需要取一个折中值。

image.png

#### 2.1挑选合适的λ值
计算1000个，画图，筛选表现最好的λ值
cv_fit <- cv.glmnet(x=x, y=y, nlambda = 1000,alpha = 1)
plot(cv_fit)

image.png

两条虚线分别指示了两个特殊的λ值,一个是lambda.min,一个是lambda.1se,这两个值之间的lambda都认为是合适的。lambda.1se构建的模型最简单，即使用的基因数量少，而lambda.min则准确率更高一点，使用的基因数量更多一点。

#### 2.2 用这两个λ值重新建模
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)

这两个值体现在参数lambda上。有了模型，可以将筛选的基因挑出来了。所有基因存放于模型的子集beta中，用到的基因有一个s0值，没用的基因只记录了“.”，所以可以用下面代码挑出用到的基因。

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)

image.png

### 3.模型预测和评估
#### 3.1自己预测自己：自我验证
newx参数是预测对象。输出结果lasso.prob是一个矩阵，第一列是min的预测结果，第二列是1se的预测结果，预测结果是概率，或者说百分比，不是绝对的0和1。
将每个样本的生死和预测结果放在一起，直接cbind即可。
lasso.prob <- predict(cv_fit,  #predict是预测，是一个泛型函数
                      newx=x ,  #这个是预测对象，cv_fit是预测数据
                      s=c(cv_fit$lambda.min,cv_fit$lambda.1se) )
re=cbind(y ,lasso.prob)  #cbind是按列合并
head(re)

image.png

#### 3.2 箱线图
对预测结果进行可视化。以实际的生死作为分组，画箱线图整体上查看预测结果。
re=as.data.frame(re)
colnames(re)=c('event','prob_min','prob_1se')
#上面一行的代码意思是把列名一次改为'event','prob_min','prob_1se'
re$event=as.factor(re$event)  #作为因子来变成分类变量
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

image.png

可以看到，真实结果是死(0)的样本，预测的值就小一点(靠近0)，真实结果是活着(1)的样本，预测的值就大一点(靠近1)，整体上趋势是对的，但不是完全准确，模型是可用的.对比两个λ值构建的模型，差别不大，1se的预测值准确一点。

下面是ROC曲线：

#### 3.3 ROC曲线
计算AUC取值范围在0.5-1之间，越接近于1越好。可以根据预测结果绘制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)))

AUC几乎接近0.9了已经很不错了

#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)))

image.png

#强迫症选项，想把两个模型画一起。
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()这个函数是取几位小数，round(656563,3)就是取3位小数

image.png

AUC

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")

#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")

image.png

最后补充一个小小知识点

#设置随机数种子，这样的话每次抽出来的数就是一样的了
set.seed(1234)
sample(1:50,7)