R生存分析|关心的变量KM曲线不显著,还有救吗?

简介: R生存分析|关心的变量KM曲线不显著,还有救吗?


如果想查看某些因素,如年龄,性别,分期,肿瘤数目,大小,实验室指标 或者 通过生信手(tao)段(lu)构建的模型和评分是否对预后有影响时候,经常会把连续变量变为分类变量,然后绘制KM曲线或者列线图等。

这时候会有一些常用的方法:

(1)实验室指标:根据正常范围进行分类

(2)临床指标:根据临床意义进行分类

(3)生信模型评分:根据中位数,平均值等进行分类

(4)生信模型评分:根据统计上的最优cutoff来分类


本次主要介绍基于统计上的最优cutoff分类的方法,并与常见的中位数进行简单的比较。


一 载入数据,R包

为了复现方便,使用内置myeloma数据集

#载入所需的R包
library("survival")
library("survminer")
#查看myeloma数据集
data(myeloma)
head(myeloma)


二 KM-中位数分类

以TP53基因为例,按照常用的中位数分为表达量高,低两组

#按照中位数进行分类
myeloma <- myeloma %>%
 mutate(TP53_cat = ifelse(TP53 > median(TP53) ,"High" , "Low"))
head(myeloma)

构建模型,并绘制KM曲线

#构建模型
fit <- survfit(Surv(time, event) ~ TP53_cat, data = myeloma)
#绘制生存曲线并显示P值
ggsurvplot(fit,
          data = myeloma,
          palette=c("red", "blue"), #自定义颜色
          legend.labs=c("TP53_High","TP53_Low"), #自定义标签
          risk.table = TRUE, 
          break.x.by = 6, #横坐标刻度间隔
          pval = T) #是否显示P值

如图显示P值不显著,这时候可以试一下最优cutoff。

更多调整可参考R|生存分析 - KM曲线 ,必须拥有姓名和颜值


三 KM-最优cutoff分类

3.1 计算最优cutoff

使用surv_cutpoint函数找到最优cutoff

res.cut <- surv_cutpoint(myeloma, 
                        time = "time",
                        event = "event",
                        variables = c("TP53", "WHSC1")) #可以添加多列
summary(res.cut)#查看最佳cutoff
      cutpoint statistic
TP53     748.3  2.928906
WHSC1   3205.6  3.361330

可以看到TP53 和 WHSC1 基因统计得到的最优cutoff。


3.2 根据最优cutoff分类

A:根据得到的最优cutoff 自行分类

myeloma <- myeloma %>%
 mutate(TP53_cutoff = ifelse(TP53 > 748.3 ,"High" , "Low"))
head(myeloma)

构建模型,并绘制KM曲线

#构建模型
fit <- survfit(Surv(time, event) ~ TP53_cutoff, data = myeloma)
#绘制生存曲线
ggsurvplot(fit,
          data = myeloma,
          palette=c("red", "blue"), #自定义颜色
          legend.labs=c("TP53_High","TP53_Low"), #自定义标签
          risk.table = TRUE,  
          break.x.by = 6,##横坐标间隔
          pval = T) #是否展示P值

可以看到P值 < 0.05了,变化还是很明显的。


B:根据surv_categorize函数获取重新构建的矩阵

此处推荐这种方法,能比较简单的获取重新构建的矩阵

##重新构建的矩阵
res.cat <- surv_categorize(res.cut)
head(res.cat)
          time event TP53 WHSC1
GSM50986 69.24     0  low   low
GSM50988 66.43     0 high   low
GSM50989 66.50     0 high  high
GSM50990 42.67     1 high  high
GSM50991 65.00     0  low   low
GSM50992 65.20     0 high   low

构建模型,并绘制KM曲线

fit <- survfit(Surv(time, event) ~TP53, data = res.cat) 
#绘制生存曲线
ggsurvplot(fit,
          data = res.cat,
          palette=c("red", "blue"),
          legend.labs=c("TP53_High","TP53_Low"), #标签
          risk.table = TRUE,
          break.x.by = 6, ##横坐标间隔
          pval = T)

结果和自行根据最优cutoff,使用ifelse进行分类得到的结果一致,此处不展示了。

相关文章
|
6月前
|
图形学
【unity实战】时间控制 昼夜交替 四季变化 天气变化效果
【unity实战】时间控制 昼夜交替 四季变化 天气变化效果
141 0
|
7月前
|
算法 测试技术 C#
【树上倍增】【内向基环树】【 图论 】2836. 在传球游戏中最大化函数值
【树上倍增】【内向基环树】【 图论 】2836. 在传球游戏中最大化函数值
|
7月前
R语言用GAM广义相加模型研究公交专用道对行程时间变异度数据的影响
R语言用GAM广义相加模型研究公交专用道对行程时间变异度数据的影响
|
7月前
R语言使用限制平均生存时间RMST比较两条生存曲线分析肝硬化患者
R语言使用限制平均生存时间RMST比较两条生存曲线分析肝硬化患者
|
7月前
|
数据挖掘
R语言绘制生存曲线估计|生存分析|如何R作生存曲线图
R语言绘制生存曲线估计|生存分析|如何R作生存曲线图
R|生存分析 - KM曲线 ,必须拥有姓名和颜值
R|生存分析 - KM曲线 ,必须拥有姓名和颜值
|
编解码 Java BI
10米精度NPP净初级生产力数据/NDVI数据/植被类型数据/土地利用数据/降雨气温分布数据/太阳辐射分布数据
10米精度NPP净初级生产力数据/NDVI数据/植被类型数据/土地利用数据/降雨气温分布数据/太阳辐射分布数据
10米精度NPP净初级生产力数据/NDVI数据/植被类型数据/土地利用数据/降雨气温分布数据/太阳辐射分布数据
|
API
DJI开发之航线重叠率的计算
DJI开发之航线重叠率的计算
990 0
DJI开发之航线重叠率的计算
|
数据处理 ice
Google Earth Engine——2004-2010年时间平均基线的月度引力异常值,该数据集所包含的数据是以 “等水厚度 “为单位,以厘米为单位表示水的垂直范围的质量偏差
Google Earth Engine——2004-2010年时间平均基线的月度引力异常值,该数据集所包含的数据是以 “等水厚度 “为单位,以厘米为单位表示水的垂直范围的质量偏差
185 0
Google Earth Engine——2004-2010年时间平均基线的月度引力异常值,该数据集所包含的数据是以 “等水厚度 “为单位,以厘米为单位表示水的垂直范围的质量偏差