【原】KM生存分析如何取最佳的cutoff

  

欢迎来到医科研，这里是白介素2的读书笔记，跟我一起聊临床与科研的故事, 生物医学数据挖掘，R语言，TCGA、GEO, SEER数据挖掘。

---

---

KM生存分析

如何取最佳cutoff

• KM生存分析中通常取中位值作为cutoff，但也并不一定是这样

• 必要时可以选择最佳的cutoff值，如何实现？### survival package

1library(survival)
2library(survminer)
3## Loading required package: ggplot2
4## Loading required package: ggpubr
5## Loading required package: magrittr


surv_cutpoint函数

• data：包含生存数据和连续变量的的数据框

• time, event：column names containing time and event data, respectively. Event values sould be 0 or 1.

官方示例

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


1. Determine the optimal cutpoint of variables

1res.cut <- surv_cutpoint(myeloma, time = "time", event = "event",
2   variables = c("DEPDC1", "WHSC1", "CRIM1"))
3
4summary(res.cut)
5##        cutpoint statistic
6## DEPDC1    279.8  4.275452
7## WHSC1    3205.6  3.361330
8## CRIM1      82.3  1.968317


2. Plot cutpoint for DEPDC1

1# palette = "npg" (nature publishing group), see ?ggpubr::ggpar
2plot(res.cut, "DEPDC1", palette = "npg")
3## $DEPDC1


image.png

3. Categorize variables

1res.cat <- surv_categorize(res.cut)
2head(res.cat)
3##           time event DEPDC1 WHSC1 CRIM1
4## GSM50986 69.24     0   high   low  high
5## GSM50988 66.43     0    low   low   low
6## GSM50989 66.50     0    low  high  high
7## GSM50990 42.67     1    low  high   low
8## GSM50991 65.00     0    low   low   low
9## GSM50992 65.20     0   high   low   low


4. Fit survival curves and visualize

1library("survival")
2fit <- survfit(Surv(time, event) ~DEPDC1, data = res.cat)
3ggsurvplot(fit, data = res.cat, risk.table = TRUE, conf.int = TRUE)


image.png

参考资料

官方文档

---

本期内容就到这里，我是白介素2，下期再见，点击下方框框留言。