利用cutoff包验证surv_cutpoint函数的统计学原理

gg系列的survminer包的函数surv_cutpoint可以帮助大家直接计算生存资料中连续自变量的截断值，在之前的的文章和大家分享过，因为篇幅限制，上次没有分析过原理，这里我们用cutoff包说下笔者对其原理的体会，不对的地方请大家批评指正。surv_cutpoint函数具体用法见链接https://www.jianshu.com/p/538f7fa11f6c

首先说下cutoff包，这也是一个可以帮助大家找截断值的R包，为CRAN的包，今天咱们就用下他在寻找生存资料连续自变量截断值的函数logrank，重要形式参数定义如下

cutoff包logrank函数形式参数
cutoff包logrank函数形式参数续

数据集我们用survival包里面的sphase

library(survival)
library(cutoff)
library(ggpubr)
library(survminer)
data(sphase)
head(sphase)

数据集结构

以下利用longrank函数计算截断值

logresult<- logrank(data=sphase,#数据集
                    time = 'RFS',#生存时间
                    y='event', #生存状态
                    x='SPF',#连续自变量
                    cut.numb=1,#截断值选择1个分界点
                    n.per=0.10,#分组后自变量组最少比例
                    y.per=0.10,#分组后因变量组最少比例
                    p.cut=0.05,#检验水准，结果只显示pvalue小于0.05的截断值
                    round=5)#保留5位有效小数
logresult

longrank结果

可以看到，一共75种情况p值小于0.05，根据自变量及因变量分组后比例过滤后就剩下9个了，以下画图看下其P值分布

ggline(logresult,
       x = "cut1",
       y="pvalue",
       palette = "jco",
       xlab = "截断值取值",
       ylab="P值")

不同截断值时P值大小情况

可以看到，在107时的P是最小的了，理论上这时候在自由度固定的情况下logrank计算的卡方值是最大的（这里的自由度为1）

以下我们用surv_cutpoint函数确定下截断值并画出生存曲线

res.cut <-surv_cutpoint(sphase, time = "RFS", event = "event",variables = c("SPF"))
summary(res.cut)

surv_cutpoint确定的截断值

可以看到，其确定的截断值为107，统计量为2.4

以下画图并及计算logrank的p值

res.cat <- surv_categorize(res.cut)
fit <- survfit(Surv(RFS, event) ~SPF, data = res.cat)
ggsurvplot(fit, 
           data = res.cat, 
           risk.table = TRUE, 
           conf.int = TRUE,
           pval = T)

生存曲线，截断值为107分组时的logrank的p值为0.0051，与上述cutoff包计算结果一致

为了进一步查看surv_cutpoint计算出的统计量是否为logrank的卡方值我们把数据导出后利用SPSS软件进行计算验证

write.csv(res.cat,"log.csv")#导出数据集

导出数据集
首先把event变成数值型变量方便后续计算，事件发生定义为1
logrank参数选择，三种计算P值的方法都用，我们最常用的其实是第一种，既不对个案进行加权
SPSS计算结果

可以看到，我们用以上两个包计算的记过，及时基于logrank的P值的，其对应的卡方值为7.8，与surv_cutpoint计算的统计量是不一样的。

实际上，surv_cutpoint计算的统计量是基于maxstat包的maxstat.text函数计算出的Maximally selected rank statistics ，感兴趣的同学可以查看原包的参考文献

原文截图