R语言glmnet包分别拟合cox比例风险回归模型的lasso回归和岭回归

在上一期的分享中，我们演示了

R语言glmnet包分别拟合二分类logistic模型的lasso回归和岭回归的简单示例，今天我们继续分享生存分析中glmnet包拟合cox比例风险回归模型的lasso回归和岭回归的简单操作。

常规数据在表示时，只需要一个值，比如患者的血压，性别等数据，不是连续型就是离散型；生存数据则有两个值，第一个是生存时间，可以看做是一个连续型的变量，第二个是生存事件，可以看做是离散型的变量。

生存分析是研究生存时间的分布规律，以及生存时间和相关因素之间关系的一种统计分析方法。在生存分析中，研究的主要对象是寿命超过某一时间的概率。还可以描述其他一些事情发生的概率，例如产品的失效、出狱犯人第一次犯罪、失业人员第一次找到工作等等。

Cox模型‌，也称为Cox回归模型或比例风险模型（Proportional Hazards Model），是由英国统计学家David Cox于1972年提出的一种半参数回归模型。该模型主要用于生存分析，能够同时分析多个风险因素对生存时间的影响，适用于医学随访研究和其他需要分析生存数据的领域。

Cox模型在医学研究中有广泛的应用，主要用于以下几个方面：

‌生存分析‌：研究患者的生存时间及其相关因素。

‌风险因素分析‌：确定与特定事件（如死亡、疾病复发）相关的风险因素。

‌生存率比较‌：比较不同组群之间的生存率。

‌预后评估‌：评估影响患者预后的因素，为治疗决策提供依据。

‌调整分析‌：调整其他变量的影响，更准确地评估特定因素对生存时间的影响。

今天我们仍以熟悉的示例数据集为例，演示一下R语言glmnet包分别拟合cox比例风险回归模型的lasso回归和岭回归的简单示例。

我们先整理一下Rstudio的环境并加载数据及常用包

rm(list=ls()) #移除所有变量数据

install.packages("") #安装包

library() #加载包

#1.数据及包准备

#读取Excel数据

library(readxl) #加载包

data <- read_excel("C:/Users/hyy/Desktop/示例数据.xlsx")

data

library(glmnet)

library(Matrix)

library(survival)

library(survminer)

library(ggplot2)

library(ggpubr)

library(Hmisc)

#Lasso回归前的生存分析

#根据变量在表中定位设置自变量、因变量、时间变量

lambdas <- seq(0,5, length.out = 1000)

X <- as.matrix(data[,3:9])

Y <- as.matrix(data[,2])

Time <- as.matrix(data[,10])

fit_sur <- Surv(Time,Y)

fit_sur

plot(fit_sur)

#根据分类指标分别绘制生存曲线

fit_sur2 <- survfit(Surv(Time,Y)~指标8, data=data)

fit_sur2

plot(fit_sur2)

ggsurvplot(fit_sur2, data = data)

#累积风险曲线

ggsurvplot(fit_sur2,

           data = data,

           fun = "cumhaz", 

           conf.int = TRUE, # 可信区间

           palette = "lancet", # 支持ggsci配色，自定义颜色，brewer palettes中的配色，等

           ggtheme = theme_bw() # 支持ggplot2及其扩展包的主题

)

#累积事件曲线：

ggsurvplot(fit_sur2,

           data = data,

           fun = "event", 

           conf.int = TRUE, # 可信区间

           palette = "grey",

           ggtheme = theme_pubclean() 

)

#或者根据变量名设置自变量、因变量、时间变量,拟合cox模型

data$status <- as.matrix(data$结局)

survfit <- coxph(Surv(时间,status)~指标1+指标2+指标3+指标4+指标5+指标6+指标7,data=data)

summary(survfit)

## 绘制Schoenfeld残差图

p <- cox.zph(survfit)

p$table

p$residuals$schoenfeld %>% ggplot(aes(p$time,p$x)) + 

  geom_line() + 

  labs(title = "Schoenfeld Residuals Over Time")

#或者使用数据集合进行拟合cox模型开始进行Lasso回归

data$status <- as.numeric(data$结局)

x<-as.matrix(X)#指定自变量为矩阵

y<-as.matrix(Surv(data$时间,data$status))#指定因变量为矩阵

set.seed(12345)

coxfit<-glmnet(x,y,family = "cox",alpha = 1)

plot(coxfit,label = T)

plot(coxfit, xvar = "lambda", label = TRUE)

plot(coxfit, xvar = "dev", label = TRUE)

print(coxfit) 

#交叉验证

set.seed(12345)

lscoxmodel <- cv.glmnet(x,y,  family = "cox",alpha = 1)

plot(lscoxmodel)

plot(lscoxmodel$glmnet.fit, "lambda", label = T)

#查看具体筛选的变量系数值和Lambda值(lambda.min处)

coef<-coef(lscoxmodel,s="lambda.min")

coef

#查出lambda.min处的Lambda值

Lambda<-lscoxmodel$lambda.min

Lambda

#lambda.min处保留的变量的个数，也就是系数不为0

bl.index<-which(coef!=0)

bl.index

#lambda.min处保留的各变量的系数

bl.coef<-coef[bl.index]

bl.coef

#列出lambda.min处保留的各变量的变量名

rownames(coef)[bl.index]

#岭回归

set.seed(12345)

ridgecox <- glmnet(x, y, family = "cox",alpha = 0)

plot(ridgecox,label = T)

plot(ridgecox, xvar = "lambda", label = TRUE)

plot(ridgecox, xvar = "dev", label = TRUE)

#交叉验证

set.seed(12345)

rgcox_model <- cv.glmnet(x, y, family = "cox",alpha = 0,lambda = lambdas,nfolds =3)

plot(rgcox_model)

plot(rgcox_model$glmnet.fit, "lambda", label = T)

#查看具体筛选的变量系数值和Lambda值(lambda.min处)

coefridge<-coef(rgcox_model,s="lambda.min")

coefridge

#查出lambda.min处的Lambda值

Lambdaridge<-rgcox_model$lambda.min

Lambdaridge

#lambda.min处保留的变量的个数，也就是系数不为0

ridge.index<-which(coefridge!=0)

ridge.index

#lambda.min处保留的各变量的系数

ridge.coef<-coef[ridge.index]

ridge.coef

#列出lambda.min处保留的各变量的变量名

rownames(coef)[ridge.index]