KM生存分析
如何取最佳cutoff
- KM生存分析中通常取中位值作为cutoff，但也并不一定是这样
- 必要时可以选择最佳的cutoff值，如何实现？ ### survival package
library(survival)
library(survminer)
## Loading required package: ggplot2
## Loading required package: ggpubr
## Loading required package: magrittr
surv_cutpoint函数
- data：包含生存数据和连续变量的的数据框
- time, event：column names containing time and event data, respectively. Event values sould be 0 or 1.
官方示例
# 0. Load some data
data(myeloma)
head(myeloma)
## molecular_group chr1q21_status treatment event time CCND1
## GSM50986 Cyclin D-1 3 copies TT2 0 69.24 9908.4
## GSM50988 Cyclin D-2 2 copies TT2 0 66.43 16698.8
## GSM50989 MMSET 2 copies TT2 0 66.50 294.5
## GSM50990 MMSET 3 copies TT2 1 42.67 241.9
## GSM50991 MAF <NA> TT2 0 65.00 472.6
## GSM50992 Hyperdiploid 2 copies TT2 0 65.20 664.1
## CRIM1 DEPDC1 IRF4 TP53 WHSC1
## GSM50986 420.9 523.5 16156.5 10.0 261.9
## GSM50988 52.0 21.1 16946.2 1056.9 363.8
## GSM50989 617.9 192.9 8903.9 1762.8 10042.9
## GSM50990 11.9 184.7 11894.7 946.8 4931.0
## GSM50991 38.8 212.0 7563.1 361.4 165.0
## GSM50992 16.9 341.6 16023.4 2096.3 569.2
1. Determine the optimal cutpoint of variables
res.cut <- surv_cutpoint(myeloma, time = "time", event = "event",
variables = c("DEPDC1", "WHSC1", "CRIM1"))
summary(res.cut)
## cutpoint statistic
## DEPDC1 279.8 4.275452
## WHSC1 3205.6 3.361330
## CRIM1 82.3 1.968317
2. Plot cutpoint for DEPDC1
# palette = "npg" (nature publishing group), see ?ggpubr::ggpar
plot(res.cut, "DEPDC1", palette = "npg")
## $DEPDC1
image.png
3. Categorize variables
res.cat <- surv_categorize(res.cut)
head(res.cat)
## time event DEPDC1 WHSC1 CRIM1
## GSM50986 69.24 0 high low high
## GSM50988 66.43 0 low low low
## GSM50989 66.50 0 low high high
## GSM50990 42.67 1 low high low
## GSM50991 65.00 0 low low low
## GSM50992 65.20 0 high low low
4. Fit survival curves and visualize
library("survival")
fit <- survfit(Surv(time, event) ~DEPDC1, data = res.cat)
ggsurvplot(fit, data = res.cat, risk.table = TRUE, conf.int = TRUE)
image.png
参考资料
官方文档
