💡专注R语言在🩺生物医学中的使用
之前分别介绍了生存分析中的寿命表法、K-M曲线、logrank检验:R语言生存分析的实现
以及Cox回归的构建、可视化以及比例风险检验的内容:R语言生存分析:Cox回归
本次主要介绍如果数据不符合PH假设时采取的方法。
关于时依协变量、时依系数的基础知识,大家可以参考这几篇文章:
如果不能满足PH假设,可以考虑使用时依协变量或者时依系数Cox回归,时依协变量和时依系数是两个概念,简单来说就是如果一个协变量本身会随着时间而改变,这种叫时依协变量,如果是协变量的系数随着时间改变,这种叫时依系数。
这里以survival包的veteran数据集为例,演示如何处理此类不符合PH检验的情况。
rm(list = ls())
library(survival)
str(veteran)
## 'data.frame': 137 obs. of 8 variables:
## $ trt : num 1 1 1 1 1 1 1 1 1 1 ...
## $ celltype: Factor w/ 4 levels "squamous","smallcell",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ time : num 72 411 228 126 118 10 82 110 314 100 ...
## $ status : num 1 1 1 1 1 1 1 1 1 0 ...
## $ karno : num 60 70 60 60 70 20 40 80 50 70 ...
## $ diagtime: num 7 5 3 9 11 5 10 29 18 6 ...
## $ age : num 69 64 38 63 65 49 69 68 43 70 ...
## $ prior : num 0 10 0 10 10 0 10 0 0 0 ...
这个数据集中的变量解释如下图:
首先构建普通的Cox回归,进行等比例风险假设,这里只选择了trt/prior/karno
3个变量,而且trt/prior
作为分类变量并没有转换为因子型,因为二分类变量数值型和因子型的结果是一样的,转不转换没啥影响!如果你还不懂分类变量在r语言中的编码方案,一定要看这篇:分类变量进行回归分析时的编码方案
fit <- coxph(Surv(time, status) ~ trt + prior + karno, data = veteran)
# 进行PH检验
zp <- cox.zph(fit)
zp
## chisq df p
## trt 0.288 1 0.59125
## prior 2.168 1 0.14087
## karno 12.138 1 0.00049
## GLOBAL 18.073 3 0.00042
可以看到变量karno
的P值小于0.05,是不满足PH假设的。
通过图形化方法查看PH检验的结果:
#op <- par(mfrow=c(1,3))
#plot(zp)
#par(op)
#ggcoxdiagnostics(fit, type = "schoenfeld")
plot(zp[3])
abline(0,0, col="red") # 0水平线
abline(h=fit$coef[3], col="green", lwd=2, lty=2)
黑色实线以及两侧的虚线是karno
的系数随着时间变化的曲线,绿色虚线是假设karno
符合PH检验时的总体估计线,红色实线是参考线。
这张图反映了karno
变量的系数随着时间的改变,karno
偏离的比较厉害(上面注释掉的代码可以都运行看看其他变量的情况),系数最开始接近-0.05,然后逐渐趋于0,最后又开始趋向-0.05,所以它的系数是一直在随着时间改变的,不符合比例风险假设。
这种情况下一个比较简单的解决方式是对时间使用分层函数。根据上面的图示我们知道karno
的系数大概分为3层(3段),可以根据两个拐点进行分层,通过survival
中的survSplit()
实现。
vet2 <- survSplit(Surv(time, status) ~ ., data= veteran,
cut=c(90, 180), # 两个拐点把时间分为3层(3段)
episode= "tgroup",
id="id")
vet2[1:7, c("id", "tstart", "time", "status", "tgroup", "age", "karno")]
## id tstart time status tgroup age karno
## 1 1 0 72 1 1 69 60
## 2 2 0 90 0 1 64 70
## 3 2 90 180 0 2 64 70
## 4 2 180 411 1 3 64 70
## 5 3 0 90 0 1 38 60
## 6 3 90 180 0 2 38 60
## 7 3 180 228 1 3 38 60
结果多了两列:tstart/tgroup
。
受试者1(id编号为1)在第72天的时候死了,所以数据和之前一样。受试者2和3(id为2和3)虽然时间在变,但是直到第3层才死去,karno
的值没有变化。
重新拟合Cox模型,此时tgroup
是分好的层,所以要用strata()
,另外karno
会随着时间变化,和时间有交互,所以用karno:strata(tgroup)
。
# 注意此时Surv()的用法!
fit2 <- coxph(Surv(tstart, time, status) ~ trt + prior + karno:strata(tgroup), data = vet2)
fit2
## Call:
## coxph(formula = Surv(tstart, time, status) ~ trt + prior + karno:strata(tgroup),
## data = vet2)
##
## coef exp(coef) se(coef) z p
## trt -0.011025 0.989035 0.189062 -0.058 0.953
## prior -0.006107 0.993912 0.020355 -0.300 0.764
## karno:strata(tgroup)tgroup=1 -0.048755 0.952414 0.006222 -7.836 4.64e-15
## karno:strata(tgroup)tgroup=2 0.008050 1.008083 0.012823 0.628 0.530
## karno:strata(tgroup)tgroup=3 -0.008349 0.991686 0.014620 -0.571 0.568
##
## Likelihood ratio test=63.04 on 5 df, p=2.857e-12
## n= 225, number of events= 128
结果表明karno
这个变量只有在tgroup=1
(第1层,前3个月)才有意义,后面两层是没有意义的。
再次进行PH检验:
cox.zph(fit2)
## chisq df p
## trt 1.72 1 0.189
## prior 3.81 1 0.051
## karno:strata(tgroup) 3.04 3 0.385
## GLOBAL 8.03 5 0.154
这时karno:strata(tgroup)
就满足了等比例风险假设。
除了对时间进行分层外,还有一种解决方法。
上面的图中我们可以看出karno
系数随时间变化的曲线明显不是线性的,我们可以通过数据变换把它变成类似线性的,比如取log
,这种变换通过tt(time transform)
函数实现。
这种方法实际上是通过tt()
函数构建了一个时依协变量,但是这样做是为了解决系数随着时间改变的问题(也就是为了解决时依系数的问题)。
fit3 <- coxph(Surv(time, status) ~ trt + prior + karno + tt(karno), # 对karno进行变换
data = veteran,
tt = function(x, t, ...) x * log(t+20) # 具体变换方式
)
fit3
## Call:
## coxph(formula = Surv(time, status) ~ trt + prior + karno + tt(karno),
## data = veteran, tt = function(x, t, ...) x * log(t + 20))
##
## coef exp(coef) se(coef) z p
## trt 0.016478 1.016614 0.190707 0.086 0.93115
## prior -0.009317 0.990726 0.020296 -0.459 0.64619
## karno -0.124662 0.882795 0.028785 -4.331 1.49e-05
## tt(karno) 0.021310 1.021538 0.006607 3.225 0.00126
##
## Likelihood ratio test=53.84 on 4 df, p=5.698e-11
## n= 137, number of events= 128
此时karno
的时依系数估计为:-0.124662 * log(t + 20)。
在构建时依协变量时,可以选择x * t、x * log(t)、x * log(t + 20)、x * log(t + 200)等等,没有明确的规定,要结合结果和图示进行选择,可以参考冯国双老师的文章:一文详解时依协变量[4]。
我们可以把现在的时依系数估计和经过变换后的的PH检验画在一起,看看变换后的效果:
# 变换后的PH检验
zp <- cox.zph(fit, transform = function(time) log(time + 20))
# 画图
plot(zp[3])
abline(0,0, col="red") # 0水平线
abline(h=fit$coef[3], col="green", lwd=2, lty=2) # 整体估计
abline(coef(fit3)[3:4],lwd=2,lty=3,col="blue") # 现在的估计
可以看到变换后结果好多了(蓝色虚线,和黑色曲线相比较),虽然还是有一点倾斜。
以上是两种处理不满足PH假设的方法,实际还有很多种方法,比较常用的是对时间进行分层,其他方法有机会继续介绍!
survival包: https://cran.r-project.org/web/packages/survival/vignettes/timedep.pdf
[2]医咖会时依协变量: https://mp.weixin.qq.com/s/CSwO42ucfeKipHl32NoOCA
[3]7code时依协变量: https://mp.weixin.qq.com/s/sBJ8Xg_e3uBfUnDOIPoxug
[4]医咖会时依协变量: https://mp.weixin.qq.com/s/CSwO42ucfeKipHl32NoOCA
🔖精选合集
联系客服