cutoff：最佳截断值之cutoff包

关于连续性变量最佳截断值的选择，之前介绍了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':## ## myelomadata(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 DEPDC1plot(res.cut, "DEPDC1", palette = "npg")## $DEPDC1

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

对于DEPDC1来说，它的最佳截断值是279.8：

# 3. Categorize variablesres.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 visualizelibrary("survival")fit <- survfit(Surv(time, event) ~DEPDC1, data = res.cat)ggsurvplot(fit, data = res.cat, pval = T)

cutoff

cutoff使用非常简单，主要是6个函数，并且使用逻辑一致。

• cox:寻找COX回归中P值最小的cutoff
• linear:寻找线性回归中P值最小的cutoff
• logit：寻找logistic回归中P值最小的cutoff
• logrank：寻找logrank检验中P值最小的cutoff
• cutit:根据cutoff对连续性变量分组
• roc:寻找ROC曲线下面积最大的cutoff，只能是分类不能是生存

前4个函数是确定最佳截点用的，支持多个截点。

library(cutoff)

使用此函数找到的最佳截点再去做cox回归得到的P值是最小的。

# 参数很简单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: 58library(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, unioncut_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值哈)

group <- cutoff::cutit(myeloma$DEPDC1,279.8)

再做cox回归：

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包功能强大，使用简答，非常值得大家学习。