之前分别介绍了生存分析中的寿命表法、K-M曲线、logrank检验：R语言生存分析的实现
以及Cox回归的构建、可视化以及比例风险检验的内容：R语言生存分析：Cox回归
本次主要介绍如果数据不符合PH假设时采取的方法。
时间依存协变量的Cox回归和时间依存系数Cox回归
关于时依协变量、时依系数的基础知识，大家可以参考这几篇文章：
- survival包的案例介绍：Using Time Dependent Covariates and Time Dependent
Coefcients in the Cox Model[1]
如果不能满足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 ...
这个数据集中的变量解释如下图：
veteran首先构建普通的Cox回归，进行等比例风险假设，这里只选择了trt/prior/karno3个变量，而且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假设的方法，实际还有很多种方法，比较常用的是对时间进行分层，其他方法有机会继续介绍！
参考资料
[1]survival包: https://cran./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
