题记:本文成形于2015年并于2017年收录于由笔者与胡志德医生共同主编,中南大学出版社出版的《疯狂统计学》,此次属旧文重发,较原文有所改动。
1. 引言
对于结果变量为二分类资料的数据,连续型自变量截断值的确定一般通过ROC分析,我们通常选用约登指数(敏感度+特异度-1)最大的点为最佳cut-off值点,这些都是常用的统计学方法。有时我们面对的问题要更复杂,假定我们的结果已经不是单纯的二分类资料,而是包含有时间因素的分类资料(Time to event data),即我们常说的生存资料。
举个简单例子,假定在某研究中我们定义生存资料的结局是死亡,那作为研究者来说不仅关心研究对象是否死亡,而且关心研究对象从入组开始到死亡的时间长度。比如某研究中试验组共入组100人,假定在入组后的第1年、第2年、第3年死亡人数分别为:0、0、90人;对照组同样入组100人,假定入组后第1年、第二年、第3年死亡人数分别为:90、0、0。在这样一个例子中,如果我们只看重死亡的人数,那么试验组与对照组结果没有差别,如果我们同时关注死亡与发生死亡的时间,那显然试验组的结局要优于对照组。从结局变量的维度上讲,生存资料在二分类的资料的基础上又增加了时间的维度。那么对于生存资料中的连续型自变量是否还可以直接采用常规ROC分析来确定截断值呢?在既往已经发表的文献中我们有时会看到有些作者确实直接采用常规ROC分析方法确定生存资料中连续型自变量的截断值,那么这样的做法是否妥当?笔者认为,这种做法目前不好判断正确与否,但至少是不妥当的,因为我们有更科学的方法。本文中,我们将介绍用R语言的survivalROC包来确定生存资料的截断值。
2. 案例与操作
[案例1] 笔者在The Cancer Genome Atlas(TCGA)数据中下载了1215例乳腺癌患者的临床资料与相对基因表达信息。下载网址:https://genome-cancer.ucsc.edu/。临床数据与基因表达信息按照芯片编号进行严格匹配后,我们从中提取了这1215例患者的生存时间,结局状态、X基因的相对表达水平。我们的研究目的是把X基因通过合适的截断值二分为X基因低表达和X基因高表达组,然后通过生存分析判断这两个组的预后是否有差异。本例中生存时间单位:天,NA表示数据缺失;我们定义的结局为死亡,1:死亡,0:截尾;X基因表达水平为经过标准化之后的基因含量,是连续型数据。读者朋友请注意此处X基因可以替换为与预后可能有关的任何连续型自变量,比如任何分子标志物含量、实验室检查指标的数值、年龄、体重、血压等等。数据整理如表1所示。
表1. 1215例患者X基因表达水平与生存资料
SurvivalROC法确定截断值
有些读者除了需要获得生存资料中连续型自变量的最佳截断值点,可能还需要像普通ROC分析一样画出ROC曲线并计算出曲线下的面积,这时我们也可以通过R软件(version 3.5.1)的survivalROC程序包进行计算(version 1.0.3)。我们继续使用[案例1]的数据,R软件操作代码如下:
library(survivalROC) # 载入程序包
data<-read.csv(file.choose()) # 导入外部数据,首先把整理后的数据另存为.csv格式
nobs<- NROW(data) # 定义数据集的行数
cutoff<- 3650 # 设定为10年生存时间,读者可根据需求设定为3年、5年等。
#Delete “NA”
data<-data[which(data$Status!=”NA”),]
#将结果变量中缺失数据删除,读者根据自己数据特点决定是否需要此命令。
# Fit survival ROC model with method of “KM”
SROC= survivalROC(Stime = data$Time,
Status = data$Status,
marker = data$Xgene,
predict.time = cutoff, method= “KM” ) #构建生存函数
cut.op= SROC$cut.values[which.max(SROC$TP-SROC$FP)] # 计算最佳截断值
cut.op # 输出最佳截断值
# plot survival ROC #绘制survivalROC曲线,并设置图形参数。
plot(SROC$FP,SROC$TP, type=”l”, xlim=c(0,1), ylim=c(0,1),
xlab = paste( “FP”,”n”, “AUC = “,round(SROC$AUC,3)),
ylab = “TP”,main = “10-yearsurvival ROC”, col=”red”)
abline(0,1)
legend(“bottomright”,c(“Xgene expression”),col=”red”,lty=c(1,1))
按照10年生存时间,最后计算的截断值为11.0686,后续可以按照前述方法整理数据进行生存分析等。同时我们也使用R软件画出了survivalROC曲线,计算了曲线下面积(AUC),如图1所示。
图1. survivalROC曲线,AUC=0.601。
3. 总结
survivalROC是一种较好的确定生存资料的方法,但需要有一定R语言的基础,对于临床医生来说可能有一定难度。我们设定不同的时间节点,比如把本例中的10-year survivalROC更换为5-year survivalROC(调整predict.time= cutoff设置),计算的截断值也不同,读者可以自行尝试。关于截断值的确定,从既往发表文献看,大体有以下几种方法:1. 按照连续变量的均数或者中位数;2. 按照连续变量的四分位数间距,比如把第25%分数作为截断值,或者把75%分数作为截断值;3. 把连续变量按照固定比例等分;4. 利用x-tile软件寻找所谓最佳cut-off;5.用普通ROC分析寻找最佳cut-off。当然还有其他方法,关于如何确定截断值,只要言之有据,能够经得住审稿人质询即可。当然在论文写作时我们也可以同时采用多种方法确定截断值,如果多种可以得到类似的结论,那说明结论比较稳健。
4. 参考文献
[1]. 胡志德, 周支瑞. 傻瓜统计学. 长沙:中南大学出版社, 2015.
[2]. 周支瑞, 胡志德. 聪明统计学. 长沙:中南大学出版社, 2016.
[3]. 周支瑞, 胡志德. 疯狂统计学. 长沙:中南大学出版社, 2018.
老师您好,新年快乐,可以提供联系邮箱吗,我的[email protected],谢谢。