最佳截断值确定之cutoff

2023-10-23 15:54:01 浏览数 (1)

关于连续性变量最佳截断值的选择，之前介绍了survminer中的surv_cutpoint以及X-tile软件：

  • R语言生存分析的实现
  • 生存分析最佳截断值的确定

今天再介绍一个非常好用的R包：cutoff

准备数据

使用myeloma数据进行演示。

代码语言：javascript复制
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()函数找最佳截断值：

代码语言：javascript复制
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值最小的cutoff
  • linear:寻找线性回归中P值最小的cutoff
  • logit：寻找logistic回归中P值最小的cutoff
  • logrank：寻找logrank检验中P值最小的cutoff
  • cutit:根据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包功能强大，使用简答，非常值得大家学习。

线性回归 time 变量 函数 数据

0 人点赞