关于连续性变量最佳截断值的选择，之前介绍了survminer中的surv_cutpoint以及X-tile软件：
- R语言生存分析的实现
- 生存分析最佳截断值的确定
今天再介绍一个非常好用的R包：cutoff
准备数据
使用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
查看根据最佳切点进行分组后的数据分布情况：
代码语言：javascript复制# 2. Plot cutpoint for DEPDC1
plot(res.cut, "DEPDC1", palette = "npg")
## $DEPDC1
然后使用surv_categorize()根据最佳截断值对数据进行分组，这样数据就根据最佳截断值变成了高表达/低表达组。
对于DEPDC1来说，它的最佳截断值是279.8：
代码语言：javascript复制# 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
根据最佳切点绘制生存曲线：
代码语言：javascript复制# 4. Fit survival curves and visualize
library("survival")
fit <- survfit(Surv(time, event) ~DEPDC1, data = res.cat)
ggsurvplot(fit, data = res.cat, pval = T)
cutoff
cutoff使用非常简单，主要是6个函数，并且使用逻辑一致。
cox:寻找COX回归中P值最小的cutofflinear:寻找线性回归中P值最小的cutofflogit：寻找logistic回归中P值最小的cutofflogrank：寻找logrank检验中P值最小的cutoffcutit:根据cutoff对连续性变量分组roc:寻找ROC曲线下面积最大的cutoff，只能是分类不能是生存
前4个函数是确定最佳截点用的，支持多个截点。
代码语言：javascript复制library(cutoff)
使用此函数找到的最佳截点再去做cox回归得到的P值是最小的。
代码语言：javascript复制# 参数很简单
cut_cox <- cutoff::cox(data = myeloma,
time = "time",
y = "event",
x = "DEPDC1",
cut.numb = 1, #截点个数
n.per = 0.2, #分组后每组样本量占总样本量的比例
y.per = 0.1, #分组后因变量比例
round = 5 #保留几位小数
)
##
## 70 rows were deleted because of missing value.
## 186 rows were analysised
##
##
## 1: all combination: 185
## 2: filter by n.per
## -
## =
## Combination: 111
##
## 3: filter by y.per
## Combination: 111
##
## 4: last combination: 58
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
cut_cox %>% arrange(pvalue)
## cut1 n n.per y y.per dump hr pvalue
## 1 279.8 108/78 0.58065/0.41935 21/30 0.19444/0.38462 b2 2.19404 0.00577
## 2 481.8 146/40 0.78495/0.21505 34/17 0.23288/0.42500 b2 2.25001 0.00665
## 3 273.9 107/79 0.57527/0.42473 21/30 0.19626/0.37975 b2 2.15075 0.00713
## 4 466.1 143/43 0.76882/0.23118 33/18 0.23077/0.41860 b2 2.14269 0.00939
## 5 273.3 106/80 0.56989/0.43011 21/30 0.19811/0.37500 b2 2.09152 0.00953
## 6 284.0 109/77 0.58602/0.41398 22/29 0.20183/0.37662 b2 2.06914 0.01015
##
## p.adjust
## 1 0.64062
## 2 0.73839
## 3 0.79196
## 4 1.04224
## 5 1.05836
## 6 1.12655
可以看到这里找到的最佳截断值也是279.8，下面使用这个值进行分组：(注意上面得到的pvalue并不是COX回归的P值哈)
代码语言：javascript复制group <- cutoff::cutit(myeloma$DEPDC1,279.8)
再做cox回归：
代码语言：javascript复制summary(coxph(Surv(myeloma$time,myeloma$event)~factor(group)))
## Call:
## coxph(formula = Surv(myeloma$time, myeloma$event) ~ factor(group))
##
## n= 256, number of events= 70
##
## coef exp(coef) se(coef) z Pr(>|z|)
## factor(group)2 1.020 2.773 0.242 4.215 2.5e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## factor(group)2 2.773 0.3606 1.726 4.456
##
## Concordance= 0.638 (se = 0.03 )
## Likelihood ratio test= 17.87 on 1 df, p=2e-05
## Wald test = 17.76 on 1 df, p=2e-05
## Score (logrank) test = 19.34 on 1 df, p=1e-05
P值是2.5e-05，你可以换其他的截断值试试，得到的P值都比这个大~
根据这个分组做cox回归得到的P值是最小的。
其他几个函数的用法都是一模一样的，就不多做讲解了，这个R包功能强大，使用简答，非常值得大家学习。
