生存分析最佳截断值的确定

最新推荐文章于 2024-01-30 23:31:45 发布

医学和生信笔记

最新推荐文章于 2024-01-30 23:31:45 发布

阅读量2k

点赞数 20

分类专栏：临床预测模型文章标签： r语言

本文链接：https://blog.csdn.net/Ayue0616/article/details/135018136

版权

临床预测模型专栏收录该内容

41 篇文章 101 订阅

订阅专栏

在做生存分析时，连续型变量最佳截断值的选择，是很火爆的问题。因为有的时候即使你通过单因素或者多因素分析得到了结果，但是根据这个结果进行K-M生存分析发现P值并不显著（P>0.05）。那这个时候你可能需要重新寻找最佳的截断值，使得K-M生存分析有意义。

目前比较常见的方法有：

基于X-tile软件实现；
使用survminer包中的surv_cutpoint()函数实现；

这两种方法都很方便，其中X-tile软件其实是基于穷举法，这个软件可以在以下网址免费下载：https://medicine.yale.edu/lab/rimm/research/software/，目前版本是3.6.1（其实很久没更新了，我估计不会更新了，因为功能已经很完备了）。这个软件我不常用，但是网上已经有很多教程了，大家感兴趣的搜索一下即可。

surv_cutpoint()我在之前的推文中详细介绍过用法：

R语言生存分析的实现

surv_cutpoint()和X-tile软件都可以计算连续型变量的最佳截断值，通过这个最佳截断值把所有数据划分为2个组，此时对两组做K-M生存分析，得到的P值是最小的。

但是这两种方法也略有不同：

surv_cutpoint()可以同时计算多个变量的最佳截断值，X-tile软件一次只能计算1个变量的最佳截断值；
surv_cutpoint()只能找1个最佳截断值，X-tile软件还可以找2个最佳截断值，把所有数据分为3组；
surv_cutpoint()是基于maxstat包计算出maximally selected rank statistics，X-tile软件是基于log-rank卡方值

surv_cutpoint()

使用myeloma数据进行演示。

library(survival)
library(survminer)
## Loading required package: ggplot2
## Loading required package: ggpubr
## 
## Attaching package: 'survminer'
## The following object is masked from 'package:survival':
## 
##     myeloma

data(myeloma)
head(myeloma)
##          molecular_group chr1q21_status treatment event  time   CCND1 CRIM1
## GSM50986      Cyclin D-1       3 copies       TT2     0 69.24  9908.4 420.9
## GSM50988      Cyclin D-2       2 copies       TT2     0 66.43 16698.8  52.0
## GSM50989           MMSET       2 copies       TT2     0 66.50   294.5 617.9
## GSM50990           MMSET       3 copies       TT2     1 42.67   241.9  11.9
## GSM50991             MAF           <NA>       TT2     0 65.00   472.6  38.8
## GSM50992    Hyperdiploid       2 copies       TT2     0 65.20   664.1  16.9
##          DEPDC1    IRF4   TP53   WHSC1
## GSM50986  523.5 16156.5   10.0   261.9
## GSM50988   21.1 16946.2 1056.9   363.8
## GSM50989  192.9  8903.9 1762.8 10042.9
## GSM50990  184.7 11894.7  946.8  4931.0
## GSM50991  212.0  7563.1  361.4   165.0
## GSM50992  341.6 16023.4 2096.3   569.2

先使用surv_cutpoint()函数找最佳截断值：

res.cut <- surv_cutpoint(myeloma, time = "time", event = "event",
                         variables = c("CCND1", "CRIM1", "DEPDC1") # 找这3个变量的最佳切点
                         )

summary(res.cut)
##        cutpoint statistic
## CCND1     450.7  1.976398
## CRIM1      82.3  1.968317
## DEPDC1    279.8  4.275452

查看根据最佳切点进行分组后的数据分布情况：

# 2. Plot cutpoint for DEPDC1
plot(res.cut, "DEPDC1", palette = "npg")
## $DEPDC1

然后使用surv_categorize()根据最佳截断值对数据进行分组，这样数据就根据最佳截断值变成了高表达/低表达组。

# 3. Categorize variables
res.cat <- surv_categorize(res.cut)
head(res.cat)
##           time event CCND1 CRIM1 DEPDC1
## GSM50986 69.24     0  high  high   high
## GSM50988 66.43     0  high   low    low
## GSM50989 66.50     0   low  high    low
## GSM50990 42.67     1   low   low    low
## GSM50991 65.00     0  high   low    low
## GSM50992 65.20     0  high   low   high

根据最佳切点绘制生存曲线：

# 4. Fit survival curves and visualize
library("survival")
fit <- survfit(Surv(time, event) ~DEPDC1, data = res.cat)
ggsurvplot(fit, data = res.cat, pval = T)

X-tile

由于X-tile软件需要使用TXT格式，所以我们把myeloma数据存为TXT格式，重新计算一下结果看看。

write.table(myeloma,file = "myeloma.txt",sep = "\t",quote = F,row.names = F)

点点点即可：

对于DEPDC1这个变量来说，X-tile软件得到的结果也是279.8：

可以看到两种方法得到的生存曲线也是一样的。

有些数据两种方法得到的结果可能并不是完全一样，但是并不会相差很大，所以这两种方法的结果基本是一致的。

survivalROC

我们再说下survivalROC，这个包也是寻找最佳截断值，但是使用场景和以上两种并不一样，它需要指定一个时间点，然后它找的最佳截断值是使得ROC曲线下面积最大的值。

所以这个方法其实是寻找ROC曲线的最佳截点，和上面两种方法并不是一个概念~

关于这个包的使用，之前也写过推文了：生存资料ROC曲线的最佳截点和平滑曲线

其他ROC曲线教程：

医学和生信笔记

关注

20
点赞
踩
35

收藏

觉得还不错? 一键收藏
2
评论
生存分析最佳截断值的确定

在做生存分析时，连续型变量最佳截断值的选择，是很火爆的问题。因为有的时候即使你通过单因素或者多因素分析得到了结果，但是根据这个结果进行K-M生存分析发现P值并不显著（P>0.05）。那这个时候你可能需要重新寻找最佳的截断值，使得K-M生存分析有意义。survminer。
复制链接

扫一扫

专栏目录