Closed ixxmu closed 3 years ago
常规型
Step1.安装加载R包
Step2.数据准备——因人而异
Step.3Lasso回归分析
Bootstrap法
方法由来及作用
代码实现
弹性网络回归
方法由来及作用
代码实现
参考资料
#install.packages("glmnet")
library(glmnet)
library(survival)
library(dplyr)
cliDat临床数据(matrix格式) exprset预测变量数据(matrix格式)
LIHC_cliDat <- data.table::fread("TCGA-LIHC.survival.tsv.gz",header = T)
load("LIHC_exprDat.Rdata")
gene <- c("TREM1","EXO1","COLEC10","C7","FCN3","DEFA6","KLK7")
tumorSample <- colnames(LIHC_exprDat)[as.numeric(substr(colnames(LIHC_exprDat),14,15))<10]
tumorSample <- substr(tumorSample,1,16)
#确保表达数据与临床信息样本名的一致性
LIHC_cliDat <- LIHC_cliDat[LIHC_cliDat$sample%in%tumorSample,]
LIHC_cliDat <- filter(LIHC_cliDat,OS.time>30)
colnames(LIHC_exprDat) <- substr(colnames(LIHC_exprDat),1,16)
exprset <- LIHC_exprDat[gene,LIHC_cliDat$sample]
exprset <- t(as.matrix(exprset))
exprset[1:4,1:4]
colnames(LIHC_cliDat) <- c("sample","fustat","_PATIENT","futime")
cliDat <- data.frame(time=LIHC_cliDat$futime,status=LIHC_cliDat$fustat)
#cliDat$status <- ifelse(cliDat$status=="Alive",0,1)
rownames(cliDat) <- LIHC_cliDat$sample
cliDat <- as.matrix(cliDat)
save(exprset,cliDat,file = "lassoData.Rdata")
说明
肝癌患者是文献中的“感兴趣结局事件”,所以必须排除正常患者 为了保证生存时间是疾病导致,而不是其他非疾病因素(如车祸等)导致,所以选择 OS.time>30
,具体数值可因癌症类型或者主观判断抉择。(比如某些严重恶性肿瘤,本身生存时间比较短,如果筛选时间过长,很可能导致样本数目严重不足)最好将数据保存为Matrix格式,这是lasso回归要求的输入数据格式。当然,也可以在分析前 as.matrix
rm(list = ls())
load("lassoData.Rdata")
x <- as.matrix(exprset)
y <- as.matrix(Surv(cliDat[,1],cliDat[,2]))
fit <- glmnet(x,y, family = "cox",alpha = 1,nlambda = 1000)
说明
x,表示预测变量(协变量)矩阵y, 因变量(生存对象矩阵),注意创建生存对象时time不可以<=0!!!!!!family ,"gaussian"(高斯)意味y为线性数据(连续),"binomial"("multinomial"),意味y为二分类变量(多分类变量)"cox"意味y为生存资料,y用Surv构建(包括"time"和"status"两列)alpha ,α=1表示lasso回归,α=0表示岭回归
nlambda,表示λ值的数量,默认为100,一般而言,但当%Dev值变化不大时,就会停止
print(fit)
# Df %Dev Lambda
# 1 0 0.0000000 0.14030
# 2 1 0.0002849 0.13900
# 3 1 0.0005652 0.13770
# 4 2 0.0008558 0.13650
# 5 2 0.0013460 0.13520
# 6 2 0.0018270 0.13400
# ...
# 260 7 0.0286300 0.01288
# 261 7 0.0286400 0.01276
# 262 7 0.0286400 0.01265
# 263 7 0.0286500 0.01253
# 264 7 0.0286600 0.01241
# 265 7 0.0286700 0.01230
说明
Df自由度(预示预测变量的数量),%Dev解析残差百分比(越接近1越好),Lambda表示λ值
plot(fit,xvar = "lambda",label = T)
plot(fit,xvar = "dev",label = T)
说明
xvar表示以什么作为x轴,label表示是否展示标签
coef(fit,s=0.0125)#预测系数
# 7 x 1 sparse Matrix of class "dgCMatrix"
# 1
# TREM1 0.141542140
# EXO1 0.160796869
# COLEC10 -0.017576378
# C7 -0.016905219
# FCN3 -0.032273735
# DEFA6 0.018806098
# KLK7 0.007648985
cvfit <- cv.glmnet(x,y,
family = "cox",
nfolds = 10, #多少折检验,常用10折交叉检验
nlambda = 100,#λ数目,会影响图片竖杠数目
alpha=1)
plot(cvfit)
说明
注意:每次拟合的结果可能略有不同
plot(fit,xvar = "lambda",label = T)
abline(v = log(cvfit$lambda.min),col="#007dc3",lty=3,lwd=2)
text(log(cvfit$lambda.min),0.15,labels = paste0("λ.min=",round(cvfit$lambda.min,4)),
pos = 4,col ="#007dc3",cex = .7 )
abline(v = log(cvfit$lambda.1se),col="#b30838",lty=3,lwd=2)
text(log(cvfit$lambda.1se),0.01,labels = paste0("λ.1se=",round(cvfit$lambda.1se,4)),
pos = 2,col ="#b30838",cex = .7 )lambda.min(最小λ):表示均方误差最小时的λ值lambda.1se(最小λ+标准误):表示在不显著改变均方误差时尽可能减少待测特征的数量
fitMin <- glmnet(x,y,family = "cox",lambda = cvfit$lambda.min,alpha = 1)
fitLse <- glmnet(x,y,family = "cox",lambda = cvfit$lambda.1se,alpha = 1)
cbind(fitLse$beta,fitMin$beta)#里面的beta等同于coef()的结果。
# 7 x 2 sparse Matrix of class "dgCMatrix"
# s0 s0
# TREM1 0.01138985 0.116068474
# EXO1 0.01896033 0.140059358
# COLEC10 . -0.006842854
# C7 . -0.007127897
# FCN3 . -0.013155785
# DEFA6 . .
# KLK7 . .
lassoGene<- rownames(fitMin$beta)[as.numeric(fitMin$beta)!=0]
predict(fit,newx = x,s = cvfit$lambda.min,type = "link")
说明
验证一下:
sum(exprset[1,lassoGene]*fitMin$beta[as.numeric(fitMin$beta)!=0])#验证风险评分
作用:提高预测变量选择的一致性,使模型更具普适性
常规型
Step1.安装加载R包
Step2.数据准备——因人而异
Step.3Lasso回归分析
Bootstrap法
方法由来及作用
代码实现
弹性网络回归
方法由来及作用
代码实现
参考资料
#install.packages("glmnet")
library(glmnet)
library(survival)
library(dplyr)
cliDat临床数据(matrix格式) exprset预测变量数据(matrix格式)
LIHC_cliDat <- data.table::fread("TCGA-LIHC.survival.tsv.gz",header = T)
load("LIHC_exprDat.Rdata")
gene <- c("TREM1","EXO1","COLEC10","C7","FCN3","DEFA6","KLK7")
tumorSample <- colnames(LIHC_exprDat)[as.numeric(substr(colnames(LIHC_exprDat),14,15))<10]
tumorSample <- substr(tumorSample,1,16)
#确保表达数据与临床信息样本名的一致性
LIHC_cliDat <- LIHC_cliDat[LIHC_cliDat$sample%in%tumorSample,]
LIHC_cliDat <- filter(LIHC_cliDat,OS.time>30)
colnames(LIHC_exprDat) <- substr(colnames(LIHC_exprDat),1,16)
exprset <- LIHC_exprDat[gene,LIHC_cliDat$sample]
exprset <- t(as.matrix(exprset))
exprset[1:4,1:4]
colnames(LIHC_cliDat) <- c("sample","fustat","_PATIENT","futime")
cliDat <- data.frame(time=LIHC_cliDat$futime,status=LIHC_cliDat$fustat)
#cliDat$status <- ifelse(cliDat$status=="Alive",0,1)
rownames(cliDat) <- LIHC_cliDat$sample
cliDat <- as.matrix(cliDat)
save(exprset,cliDat,file = "lassoData.Rdata")
说明
肝癌患者是文献中的“感兴趣结局事件”,所以必须排除正常患者 为了保证生存时间是疾病导致,而不是其他非疾病因素(如车祸等)导致,所以选择 OS.time>30
,具体数值可因癌症类型或者主观判断抉择。(比如某些严重恶性肿瘤,本身生存时间比较短,如果筛选时间过长,很可能导致样本数目严重不足)最好将数据保存为Matrix格式,这是lasso回归要求的输入数据格式。当然,也可以在分析前 as.matrix
rm(list = ls())
load("lassoData.Rdata")
x <- as.matrix(exprset)
y <- as.matrix(Surv(cliDat[,1],cliDat[,2]))
fit <- glmnet(x,y, family = "cox",alpha = 1,nlambda = 1000)
说明
x,表示预测变量(协变量)矩阵y, 因变量(生存对象矩阵),注意创建生存对象时time不可以<=0!!!!!!family ,"gaussian"(高斯)意味y为线性数据(连续),"binomial"("multinomial"),意味y为二分类变量(多分类变量)"cox"意味y为生存资料,y用Surv构建(包括"time"和"status"两列)alpha ,α=1表示lasso回归,α=0表示岭回归
nlambda,表示λ值的数量,默认为100,一般而言,但当%Dev值变化不大时,就会停止
print(fit)
# Df %Dev Lambda
# 1 0 0.0000000 0.14030
# 2 1 0.0002849 0.13900
# 3 1 0.0005652 0.13770
# 4 2 0.0008558 0.13650
# 5 2 0.0013460 0.13520
# 6 2 0.0018270 0.13400
# ...
# 260 7 0.0286300 0.01288
# 261 7 0.0286400 0.01276
# 262 7 0.0286400 0.01265
# 263 7 0.0286500 0.01253
# 264 7 0.0286600 0.01241
# 265 7 0.0286700 0.01230
说明
Df自由度(预示预测变量的数量),%Dev解析残差百分比(越接近1越好),Lambda表示λ值
plot(fit,xvar = "lambda",label = T)
plot(fit,xvar = "dev",label = T)
说明
xvar表示以什么作为x轴,label表示是否展示标签
coef(fit,s=0.0125)#预测系数
# 7 x 1 sparse Matrix of class "dgCMatrix"
# 1
# TREM1 0.141542140
# EXO1 0.160796869
# COLEC10 -0.017576378
# C7 -0.016905219
# FCN3 -0.032273735
# DEFA6 0.018806098
# KLK7 0.007648985
cvfit <- cv.glmnet(x,y,
family = "cox",
nfolds = 10, #多少折检验,常用10折交叉检验
nlambda = 100,#λ数目,会影响图片竖杠数目
alpha=1)
plot(cvfit)
说明
注意:每次拟合的结果可能略有不同
plot(fit,xvar = "lambda",label = T)
abline(v = log(cvfit$lambda.min),col="#007dc3",lty=3,lwd=2)
text(log(cvfit$lambda.min),0.15,labels = paste0("λ.min=",round(cvfit$lambda.min,4)),
pos = 4,col ="#007dc3",cex = .7 )
abline(v = log(cvfit$lambda.1se),col="#b30838",lty=3,lwd=2)
text(log(cvfit$lambda.1se),0.01,labels = paste0("λ.1se=",round(cvfit$lambda.1se,4)),
pos = 2,col ="#b30838",cex = .7 )lambda.min(最小λ):表示均方误差最小时的λ值lambda.1se(最小λ+标准误):表示在不显著改变均方误差时尽可能减少待测特征的数量
fitMin <- glmnet(x,y,family = "cox",lambda = cvfit$lambda.min,alpha = 1)
fitLse <- glmnet(x,y,family = "cox",lambda = cvfit$lambda.1se,alpha = 1)
cbind(fitLse$beta,fitMin$beta)#里面的beta等同于coef()的结果。
# 7 x 2 sparse Matrix of class "dgCMatrix"
# s0 s0
# TREM1 0.01138985 0.116068474
# EXO1 0.01896033 0.140059358
# COLEC10 . -0.006842854
# C7 . -0.007127897
# FCN3 . -0.013155785
# DEFA6 . .
# KLK7 . .
lassoGene<- rownames(fitMin$beta)[as.numeric(fitMin$beta)!=0]
predict(fit,newx = x,s = cvfit$lambda.min,type = "link")
说明
验证一下:
sum(exprset[1,lassoGene]*fitMin$beta[as.numeric(fitMin$beta)!=0])#验证风险评分
作用:提高预测变量选择的一致性,使模型更具普适性
已经在微信付费
第一种和第二种几乎一样,当然第二种我之前已经实现此种方法,第三种弹性网络作者倒入了文章中的source代码,并没有做出过多解释,我会总结之后把详细代买分享出来
微信真是业界毒瘤,md,原来微信的付费文章不支持复制黏贴,fuck,下面就是补全的代码部分,大概调整了下格式,看完应该可以理解,或许之后会补上图片,或许就这了,真的参考价值不大,30大洋不值得
Bootstrap法
方法由来及作用
作用:提高预测变量选择的一致性,使模型更具普适性
image-20210204143017559 代码实现 初阶版:for循环
library(glmnet) library(survival) library(dplyr) library(caret)
load("lassoData.Rdata")
outdata <- data.frame() for (i in 1:500){
trainIndex <- createDataPartition(cliDat[,2], p = 0.75, list =FALSE)
xTrain <- exprset[trainIndex,]
cliTrain <- cliDat[trainIndex,]
yTrain <- as.matrix(Surv(cliTrain[,1],cliTrain[,2]))
cvfit <- cv.glmnet(xTrain,yTrain,
family = "cox",
nfolds = 10, #多少折检验,常用10折交叉检验
nlambda = 500,#λ数目,会影响图片竖杠数目
alpha=1)
fitMin <- glmnet(xTrain,yTrain,family = "cox",lambda = cvfit$lambda.min,alpha = 1)
lassoGene<- rownames(fitMin$beta)[as.numeric(fitMin$beta)!=0]#提取模型基因
index <- cbind(i,fitMin$lambda,fitMin$dev.ratio,lassoGene)
outdata <- rbind(outdata,index)
}
table(outdata$lassoGene)
高阶版:foreach并行操作
load("lassoData.Rdata")
# install.packages("foreach")
# install.packages("doParallel")
library(foreach);library(doParallel)
################################################################## #~~~~~~~~~~~~~~~~~~~~~~~~~~设置并行后端~~~~~~~~~~~~~~~~~~~~~~~~~~
parallel::detectCores()
n.cores <- parallel::detectCores()/2
my.cluster <- parallel::makeCluster(
n.cores, #线程数
type = "PSOCK" #设置并程后端:FORK(Linux、Mac)和PSOCK(All)
)
print(my.cluster) #检查集群
doParallel::registerDoParallel(cl = my.cluster) #注册
foreach::getDoParRegistered() #检查是否注册
foreach::getDoParWorkers() #检查线程 #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ####################################################################
outdata <- foreach (
i=1:500,
.combine = "c" ) %dopar% {
trainIndex <- caret::createDataPartition(cliDat[,2], p = 0.75, list =FALSE)
xTrain <- exprset[trainIndex,]
cliTrain <- cliDat[trainIndex,]
yTrain <- as.matrix(survival::Surv(cliTrain[,1],cliTrain[,2]))
cvfit <- glmnet::cv.glmnet(xTrain,yTrain,
family = "cox",
nfolds = 10, #多少折检验,常用10折交叉检验
nlambda = 500,#λ数目,会影响图片竖杠数目
alpha=1)
fitMin <- glmnet::glmnet(xTrain,yTrain,family = "cox",lambda = cvfit$lambda.min,alpha = 1)
lassoGene<- rownames(fitMin$beta)[as.numeric(fitMin$beta)!=0] #提取模型基因
return(lassoGene)
}
parallel::stopCluster(cl = my.cluster)
table(outdata)
初阶版和高阶版的区别在于,高阶版采用并行运行的方式,大大缩短运行时间!!!
为什么两个版本的结局不一样,其实每次运行同一代码,结果也不一定一样。因为createDataPartition,它是随机分组,每次分组情况存在一定的差异,如果没有差异就没有重复运行的意义了!!!
弹性网络回归
方法由来及作用 image-20210213193337710 作用:通过类似于穷举法的原理,确定α值,使变量的筛选更加灵活,在尽量不损失原始数据的基础上,尽可能去除无关变量(很好地,同时克服lasso回归和岭回归的优缺点)。同时,也保证了模型的拟合优度。
如果模型中含有较多的无关变量时,因lasso回归可以将无关变量排除,故lasso回归比岭回归模型更优,其在不同数据集中的方差更小。 相反,如果模型中大多数变量为相关变量时,因岭回归不会误删一些变量,故岭回归比lasso回归模型更优,其在不同数据集中的方差更小。
代码实现
# install_github("YuLab-SMU/yulab.utils")
# yulab.utils::install_zip_gh("nathanvan/parallelsugar")
# library(parallelsugar) rm(list = ls()) source("Elastic Net(from Cell)_Mac.R") source("Elastic Net(from Cell)_Window.R")
library(survival) library(survminer) library(RColorBrewer) library(glmnet)
load("lassoData.Rdata")
cliDat <- as.data.frame(cliDat)
exprset <- as.data.frame(exprset)
head(exprset);head(cliDat)
system.time(fun.cox.elasticnet(DATA_ORG =exprset, #预测因素数据框
time = cliDat$time, #时间
status = cliDat$status, #生存状态
cores=6, #线程,detectCores()查看()我的是12
nfold=10, #XX折检验
min.elnet=0.5, #alpha最小值
max.elnet=0.6, #alpha最大值,这里可以看出是弹性回归
percentage=0.25, #bootstrap取样比例
REPEATS=50 #检验alpha值次数 ))
结果解读:CVM与平均拟合误差有关
注意区分Mac用户和window用户(我是window用户),window用户需要提前下载好parallelsugar包。
(对于线程数的设置,parallel包的mclapply函数对于window系统不支持多线程运行)
参考资料 glmnet:https://web.stanford.edu/~hastie/glmnet/glmnet_alpha.html#lin foreach:https://www.blasbenito.com/post/02_parallelizing_loops_with_r/ parlapply:https://stackoverflow.com/questions/46343775/mc-cores-1-is-not-support-on-windows parallelsugar:https://github.com/nathanvan/parallelsugar
https://mp.weixin.qq.com/s/XwNVAWw21zy04lXFqjDIcw