【原】绘制Cox回归模型的决策曲线

Memo_Cleon 2024-02-18 发布于上海

展开全文

决策曲线分析（Decision Curve Analysis，DCA）的基本概念参见《二分类logistic回归模型决策曲线的绘制与解读》。此次笔记绘制Cox回归的决策曲线。

生存分析中结局事件的发生会随着时间的变化而变化，这个我们在绘制生存模型的校准曲线和生存模型的ROC曲线时都已提及。生存结局这种时间特异性意味着我们在进行生存模型的评估时，需要先确定一个评估时间点，生存模型的DCA绘制也不例外。绘制生存模型DCA的关键是要计算每一个个体在不同时间点上的事件发生率。这里的事件发生率指的是结局事件的发生概率（Event-probability），或者叫失效事件发生概率（Probability of faililure），有时候也称事件发生风险（Risk of an event）。注意此风险（Risk）非彼风险（Hazard）。

在非竞争风险模型中，结局事件的发生概率=1-生存概率，所以预知结局事件的发生率，只需要先求得每一个个体的在不同时间点上的生存概率就可以了。

Cox回归模型基本形式：

这是Cox风险比例模型的风险函数（hazard function）的形式，其可以转换为生存函数（survival function）形式：

简单地表示为：

X为Cox模型中的协方差矩阵，β为Cox模型的参数估计值向量，S₀(t)是到时间t时基线生存概率。具有特征X的某个个体在t时的生存函数（生存概率，S(t|X)）等于t时刻对应基线生存函数S₀(t)的exp(Xβ)次幂。基线生存函数S₀(t)表示t时间时协变量X全为零时的生存概率，所有个体在某一时刻的S₀(t)都是一样的，因此在某一时刻个体之间的生存率差别仅在于协变量的不同，S₀(t)和β值都可以通过计算来获得，知道了S₀(t)和β值，我们就可以预测任何一个个体在任意时刻的生存概率了。

示例：数据来自survival包，德国乳腺癌研究学组（GBSG）在1984-1989年主导的一项临床试验，含720例淋巴结阳性乳腺癌患者，686例预后变量数据完整。pid：患者编号；age：年龄(年)；meno：绝经状态(0=绝经前，1=绝经后)；size：肿瘤大小(mm)；grade：肿瘤分级；nodes：阳性淋巴结数量；pgr：孕酮受体(fmol/l)；er：雌激素受体(fmol/l)；hormon：激素疗法(0=no，1=yes)；rfstime：无复发时间（首次复发、死亡或最后一次随访的天数）；status：结局状态(0=存活或无复发，1=复发或死亡)。

【1】数据导入，数据集拆分

library(survival)

GBSG<-gbsg

GBSG$grade<-factor(GBSG$grade)

library(caret)

library(ggplot2);library(lattice)

set.seed(20240121)

split<-createDataPartition(GBSG$status,p=0.7,list=FALSE)

split<-as.vector(split)

trainset<-GBSG[split,]

testset<-GBSG[-split,]

trainset$Subset<-0

testset$Subset<-1

BC<-rbind(trainset,testset)

BC$Subset<-factor(BC$Subset,levels=c(0,1),labels=c("trainset","testset"))

【2】建立Cox回归模型

library(StepReg)

stepwiseCox(Surv(rfstime,status==1)~age+meno+size+grade+nodes+pgr+er+hormon,data=trainset,selection="score",select="AIC",method="efron",sle=0.05,sls=0.10,best=1)

全子集回归获得4个建模变量：size、nodes、pgr、hormon，利用这些变量建立第一个模型BC.cox[实际上建模时一般不考虑治疗因素，此处仅演示]。剩余变量建立第二个模型BC.cox2。设置x=TRUE是因为有的函数在计算生存率时需要。

library(survival)

BC.cox<-coxph(Surv(rfstime,status==1)~size+nodes+pgr+hormon,data=trainset,x=TRUE,y=TRUE)

BC.cox2<-coxph(Surv(rfstime,status==1)~age+meno+grade+er,data=trainset,x=TRUE,y=TRUE)

【3】绘制决策曲线

此次笔记示例dcurves::dac和ggDCA::dac来绘制Cox回归模型的决策曲线。

3.1 dcurves:Decision Curve Analysis for Model Evaluation

采用dcurves::dac绘制决策曲线最关键的在于计算事件发生率。

3.1.1 计算事件发生率

3.1.1.1 predict::survival间接计算事件发生率。需要先将生存时间调整为感兴趣的时间点，再通过predict获得感兴趣时间点的生存率。survival包中predict.coxph函数可以获得Cox回归的预测值，这些预测值包括线性预测值（linear predictor，"lp"）、风险得分（risk score=exp(lp)，"risk"）、给定协变量和随访时间的预期事件数（expected number of events，"expected"）, 每个预测变量的参数估计值（terms of the linear predictor，"terms"）, 生存概率（survival probability=exp(-expected)，"survival"）。其中通过参数type="survival"预测的生存概率是通过模型预测获得的受试者在他的生存时间上的生存概率。比如在该示例中，乳腺癌复发或死亡是结局事件，1号患者（pid=1）在第1337天时发生死亡或者复发，而第132号患者（pid=132）在随访时间1838天时没有死亡也没有复发，用我们后面建立的模型BC.cox进行预测，获得的两位患者的生存概率是分别是0.6422和0.3970，其含义是根据患者的特征用我们的模型进行预测，1号患者存活至1337天或到1337天无复发的概率是0.6422（实际该患者确实存活），而132号患者存活着1838天或到1838天无复发的概率是0.3970（实际该患者确实死亡）。简单来说，直接通过predict进行预测，我们只能获得每位患者在一个时间点（其生存时间）的生存概率，因此如果想获得任意时间点的生存概率，需要将生存时间调整为感兴趣的时间点。

#在整个数据集上计算每个个体1年的结局事件的发生率，再从中提取训练集和测试集

BC.s1Y<-BC

BC.s1Y$rfstime=365.25#将生存时间调整为感兴趣的时间点，此处以1年为例

Pr.s1Y<-predict(BC.cox,newdata=BC.s1Y,type="survival") #预测所有患者在1年时的生存概率

BC$Pr.e1Y=1-Pr.s1Y #在数据集BC中添加名称为Pr.e1Y的新列，其值为1年时结局事件发生率，等于1-1年的生存概率。也可通过dplyr包中的mutate函数（创建、修改、删除列）来实现：library(dplyr); BC<-mutate(BC,Pr.e1Y=1-Pr.s1Y)

trainset<-BC[BC$Subset=="trainset",] #从数据集中提取训练集

testset<-BC[BC$Subset=="testset",]#从数据集中提取测试集

#计算两个不同模型在训练集上每个个体3年的结局事件的发生率

trainset.s3Y<-trainset

trainset.s3Y$rfstime=365.25*3

Pr.s3Y<-predict(BC.cox,newdata=trainset.s3Y,type="survival") #获取BC.cox模型的预测的生存率

trainset$Pr.e3Y=1-Pr.s3Y #library(dplyr);trainset<-mutate(trainset,Pr.e3Y=1-Pr.s3Y)

Pr2.s3Y<-predict(BC.cox2,newdata=trainset.s3Y,type="survival") #获取BC.cox2模型的预测的生存率

trainset$Pr2.e3Y=1-Pr2.s3Y

#单独计算测试集上每个个体5年的结局事件的发生率

testet.s5Y<-testset

testet.s5Y$rfstime=365.25*5

Pr.s5Y<-predict(BC.cox,newdata=testet.s5Y,type="survival")

testset$Pr.e5Y=1-Pr.s5Y

3.1.1.2 summary&survfit::survival计算事件发生率。网络教程基本上都是用的这个方法。survfit.coxph可以通过Cox比例风险计算预测的生存函数。summary.survfit可以返回由times指定的时间点上survfit对象中众多结果，其中surv部分是生存率，选择每个个案的第一行[1,]可将数据格式调整为每个个案数据占一行。

#训练集每个个体在1年、3年和5年的结局事件的发生率

trainset$Pr.e1Y<-1-summary(survfit(BC.cox,newdata=trainset),times=365.25)$surv[1,]

trainset$Pr.e3Y<-1-summary(survfit(BC.cox,newdata=trainset),times=365.25*3)$surv[1,]

trainset$Pr.e5Y<-1-summary(survfit(BC.cox,newdata=trainset),times=365.25*5)$surv[1,]

#注：如果不指定newdata参数，survfit帮助文件指出此时会生成一个没有意义的“伪”个体的生存率，这个“伪”个体的协方差取值等于拟合模型的成分的均值。帮助文件中说这个“伪”个体“almost never make sense”，但笔者觉得其实它的实际意义也非常明确，即该值刚好是基线生存函数S₀(t)。因为在survival包中在线性预测值等于其均值中心化后的值与相应的参数值（β）的乘积，当某个体的某个特征值刚好等于其均值时，其中心化值刚好等于0。所有特征的中心化值都为0时，S(t|X)=S₀(t)。如果在summary指定了times=t，返回值是唯一的一个值，就是时间t时的基线生存函数；如果不指定times，默认输出所有结局事件发生时的生存时间的基线生存函数，如果也想显示删失时间的基线生存函数，需要censored=TRUE。

所以，根据S(t|X)=S0(t)^exp(Xβ)，模型BC.cox预测的1年生存率用以下两种命令是等价的：

summary(survfit(BC.cox,newdata=trainset),times=365.25)$surv

summary(survfit(BC.cox),times=365.25)$surv^exp(BC.cox$linear.predictors)

#训练集线性预测值获得方法：①BC.cox$linear.predictors；②fitted(BC.cox)；③predict(BC.cox,newdata=trainset,type="lp")

#测试集每个个体在1年、3年和5年的结局事件的发生率

testset$Pr.e1Y<-1-summary(survfit(BC.cox,newdata=testset),times=365.25)$surv[1,]

testset$Pr.e3Y<-1-summary(survfit(BC.cox,newdata=testset),times=365.25*3)$surv[1,]

testset$Pr.e5Y<-1-summary(survfit(BC.cox,newdata=testset),times=365.25*5)$surv[1,]

#在整个数据集上计算每个个体在1年、3年和5年的结局事件的发生率，再从中提取训练集和测试集

BC$Pr.e1Y<-1-summary(survfit(BC.cox,newdata=BC),times=365.25)$surv[1,]

BC$Pr.e3Y<-1-summary(survfit(BC.cox,newdata=BC),times=365.25*3)$surv[1,]

BC$Pr.e5Y<-1-summary(survfit(BC.cox,newdata=BC),times=365.25*5)$surv[1,]

trainset<-BC[BC$Subset=="trainset",]

testset<-BC[BC$Subset=="testset",]

3.1.1.3 broom::augment计算事件发生率。思路同3.1，只是换用了不同的函数。在dcurves包的Vignettes文件（https://cran./web/packages/dcurves/vignettes/dca.html）中使用的是这个方法，同时借助了管道函数。

但这里（蓝色背景）的代码实际上存在问题，这个命令会导致新生成的数据集df_surv_updated的生存时间ttcancer全部等于1.5，进而后面绘制的DCA曲线也不对。正确的方法应该是在函数mutate中指定数据集df_surv_updated，这样才可以在数据集df_surv_updated中增加新变量pr_failure18。

在训练集中增加1年的时间发生率命令如下：

library(dplyr)

library(broom)

trainset<-augment(BC.cox,newdata=trainset %>% mutate(rfstime=365.25),type.predict="expected") %>% mutate(trainset,Pr.e1Y=1-exp(-.fitted))

如不使用管道函数，命令可拆解为：

library(dplyr)

trainset.1Y<-mutate(trainset,rfstime=365.25) #trainset.1Y<-trainset$rfstime=365.25

library(broom)

epd<-augment(BC.cox,newdata=trainset.1Y,type.predict="expected")

epd<-mutate(epd,Pr.e1Y=1-exp(-.fitted))

trainset$Pr.e1Y<-epd$Pr.e1Y #mutate(trainset,Pr.e1Y=epd$Pr.e1Y)

其实使用参数type.predict直接指定为"survival"命令更简洁一些：

library(dplyr)

BC.1Y<-mutate(BC,rfstime=365.25) #BC.1Y<-BC$rfstime=365.25

library(broom)

epd<-augment(BC.cox,newdata=BC.1Y,type.predict="survival")

BC$Pr.e1Y<-1-epd$.fitted

以上命令可以用管道函数进行简化：

BC<-augment(BC.cox,newdata=BC %>% mutate(rfstime=365.25),type.predict="survival") %>% mutate(BC,Pr.e1Y=1-.fitted)

trainset<-BC[BC$Subset=="trainset",]

testset<-BC[BC$Subset=="testset",]

3.1.1.4 riskRegression::predictRisk计算事件发生率

library(riskRegression)

trainset$Pr.e1Y<-predictRisk(BC.cox,newdata=trainset,times=365.25)

trainset$Pr.e3Y<-predictRisk(BC.cox,newdata=trainset,times=365.25*3)

trainset$Pr.e5Y<-predictRisk(BC.cox,newdata=trainset,times=365.25*5)

testset$Pr.e5Y<-predictRisk(BC.cox,newdata=testset,times=365.25*5)

3.1.2 绘制决策曲线

dcurves::dca参数介绍见《二分类logistic回归模型决策曲线的绘制与解读》中的5.2，此不赘述。

library(dcurves)

#模型预测训练集1年结局事件发生率的决策曲线

BCM.1Y<-dca(Surv(rfstime,status)~Pr.e1Y,data=trainset,thresholds=seq(0,1,by=0.01),label=list(Pr.e1Y="BC.Model.1Y"),time=365.25)

plot(BCM.1Y,type='net_benefit',smooth=TRUE,show_ggplot_code=FALSE)

#使用通道函数简化命令

dca(Surv(rfstime,status)~Pr.e1Y,data=trainset,thresholds=seq(0,1,by=0.01),label=list(Pr.e1Y="BC.Model.1Y"),time=365.25) %>% plot(smooth=TRUE)

#两个不同模型预测训练集3年结局事件发生率的决策曲线

dca(Surv(rfstime,status)~Pr.e3Y+Pr2.e3Y,data=trainset,thresholds=seq(0,1,by=0.01),label=list(Pr.e3Y="BC.Model1.3Y",Pr2.e3Y="BC.Model2.3Y"),time=365.25*3) %>% plot(smooth=TRUE)

#模型预测测试集5年结局事件发生率的决策曲线

dca(Surv(rfstime,status)~Pr.e5Y,data=testset,thresholds=seq(0,1,by=0.01),label=list(Pr.e5Y="BC.Model.5Y"),time=365.25*5) %>% plot(smooth=TRUE)

3.2 ggDCA：Calculate and Plot Decision Curve

dac::ggDCA绘制DCA非常简单容易理解，而且形式多样，但是笔者认为目前仍存在Bug，希望能在后续的程序包更新中有修正。

library(ggDCA)

library(ggplot2)

#训练集DCA

BCM.trainset<-dca(BC.cox,model.names="BC.Model.1Y.trainset",times=365.25)

#不指定new.data默认拟合模型的数据集；不指定times则默认中位生存时间

ggplot(BCM.trainset)

#测试集DCA

BCM.testset<-dca(BC.cox,model.names="BC.Model.5Y.testset",new.data=testset,times=365.25*5)

ggplot(BCM.testset)

只验证1个模型时，函数会同时绘制训练集模型和验证模型（validate）。同时绘制多个验证模型（2个及以上），结果图中就只包含多个验证模型的DCA结果，不再包含相应的训练模型。这我们在《二分类logistic回归模型决策曲线的绘制与解读》一文中介绍过，同时也指出采用dac::ggDCA绘制的logistic模型测试集决策曲线是有问题的，同样在绘制Cox模型在测试集上的DCA也存在问题。如上图所示，按照教程橘红色应为训练集模型的DCA，而黄绿色的虚线则是测试集的DCA，但是很显然这两条曲线弄反了。去掉上述命令中的new.data参数即可很容易地绘制5年时训练集的DCA：

ggplot(dca(BC.cox,times=365.25*5))

用dca::dcurves函数得到测试集5年非平滑DCA：

dca(Surv(rfstime,status)~Pr.e5Y,data=testset,thresholds=seq(0,1,by=0.01),time=365.25*5) %>% plot

对比一下三幅图，很显然那一条黄绿色的虚线才是训练集模型5年的DCA，而橘红色实线才是测试集的DCA。通过直接比对dac::ggDCA的结果数据也可以知道，橘红色实线数据同相应训练集DCA绘制数据是一致的。除此之外，还有另外一处问题：与单独绘制的训练集DCA相比，测试集DCA中All与黄绿色的DCA曲线（注：如上所言这条线应该是训练集的DCA）之间的距离在阈值概率（x轴）较低的时候距离要远得多（图中蓝色椭圆处）。