序言
本章主要介绍如果数据不符合PH假设时采取的方法
1.提取数据集
library(mlr3verse)
library(mlr3proba)
library(survival)
task =tsk("gbcs")
data<-task$data()2.Cox回归，进行等比例风险假设
fit <- coxph(Surv(time, status) ~ nodes + prog_recp, data = data)
# 进行PH检验
zp <- cox.zph(fit)
zp chisq df p
nodes 1.05 1 0.30628
prog_recp 12.68 1 0.00037
GLOBAL 14.44 2 0.00073可以看到变量prog_recp的P值小于0.05，是不满足PH假设的
plot(zp[2])这张图反映了prog_recp变量的系数随着时间的改变，prog_recp偏离的比较厉害，系数大概从810天开始增加，趋近于0；1900天又开始下降，所以它的系数是一直在随着时间改变的，不符合比例风险假设。
3.对时间使用分层函数
vet2 <- survSplit(Surv(time, status) ~ ., data= data,
cut=c(810, 1900), # 两个拐点把时间分为3层（3段）
episode= "tgroup",
id="id")
vet2[1:7, c("id", "tstart", "time", "status", "tgroup", "prog_recp", "nodes")] id tstart time status tgroup prog_recp nodes
1 1 0 810 0 1 141 5
2 1 810 1900 0 2 141 5
3 1 1900 2282 0 3 141 5
4 2 0 810 0 1 78 1
5 2 810 1900 0 2 78 1
6 2 1900 2006 0 3 78 1
7 3 0 810 0 1 422 1数据多了三列：id/tstart/tgroup
fit2 <- coxph(Surv(tstart, time, status) ~nodes
+prog_recp:strata(tgroup), data = vet2)
fit2 Call:
coxph(formula = Surv(tstart, time, status) ~ nodes + prog_recp:strata(tgroup),
data = vet2)
coef exp(coef) se(coef) z p
nodes 0.064944 1.067099 0.008907 7.292 3.06e-13
prog_recp:strata(tgroup)tgroup=1 -0.022822 0.977436 0.004791 -4.764 1.90e-06
prog_recp:strata(tgroup)tgroup=2 -0.003196 0.996809 0.001171 -2.728 0.00637
prog_recp:strata(tgroup)tgroup=3 -0.004299 0.995710 0.003534 -1.217 0.22376
Likelihood ratio test=121.1 on 4 df, p=< 2.2e-16
n= 1333, number of events= 171结果表明prog_recp这个变量只有在tgroup=1和2才有意义，后面一层是没有意义的。
zp1<-cox.zph(fit2)
zp1 chisq df p
nodes 1.40 1 0.24
prog_recp:strata(tgroup) 2.66 3 0.45
GLOBAL 4.02 4 0.40这时prog_recp:strata(tgroup)就满足了等比例风险假设
plot(zp1)4.连续性时依系数变换
上面的图中可以看出prog_recp系数随时间变化的曲线明显不是线性的，我们可以通过数据变换把它变成类似线性的，比如取log，这种变换通过tt(time transform)函数实现。
fit3 <- coxph(Surv(time, status) ~ nodes + prog_recp + tt(prog_recp), # 对prog_recp进行变换
data = data,
tt = function(x, t, ...) x * log(t)
)
fit3 Call:
coxph(formula = Surv(time, status) ~ nodes + prog_recp + tt(prog_recp),
data = data, tt = function(x, t, ...) x * log(t))
coef exp(coef) se(coef) z p
nodes 0.064031 1.066125 0.008728 7.337 2.19e-13
prog_recp -0.076417 0.926430 0.019173 -3.986 6.73e-05
tt(prog_recp) 0.009967 1.010017 0.002659 3.749 0.000178
Likelihood ratio test=110.5 on 3 df, p=< 2.2e-16
n= 686, number of events= 171此时prog_recp的时依系数估计为：0.009967 * log(t)。
在构建时依协变量时，可以选择x * t、x * log(t)、x * log(t + 20)、x * log(t + 200)等等，没有明确的规定，要结合结果和图示进行选择
# 变换后的PH检验
zp <- cox.zph(fit, transform = function(time) log(time ))
plot(zp[2])
abline(0,0, col="red") # 0水平线
abline(h=fit$coef[2], col="green", lwd=2, lty=2) # 整体估计
abline(coef(fit3)[2:3],lwd=2,lty=3,col="blue") # 现在的估计样条技术
样条技术是一种用于平滑和逼近数据的数学方法。它常用于统计学和数据分析中，特别是当数据呈现出复杂的模式或非线性关系时。
使用样条技术pspline(time)进行数据变换
fit4 <- coxph(Surv(time, status) ~ nodes + prog_recp + tt(prog_recp), # 对prog_recp进行变换
data = data,
tt = function(x, t, ...) x * pspline(t)
)
fit4 Call:
coxph(formula = Surv(time, status) ~ nodes + prog_recp + tt(prog_recp),
data = data, tt = function(x, t, ...) x * pspline(t))
coef se(coef) se2 Chisq DF p
nodes 6.48e-02 8.91e-03 8.91e-03 5.28e+01 1.00 3.6e-13
prog_recp -5.46e-02 3.74e-02 2.34e-02 2.13e+00 1.00 0.1441
tt(prog_recp), linear 6.55e-06 3.06e-06 3.06e-06 4.56e+00 1.00 0.0326
tt(prog_recp), nonlin 1.45e+01 3.61 0.0042
Iterations: 10 outer, 40 Newton-Raphson
Theta= 1
Degrees of freedom for terms= 1.0 0.4 4.6
Likelihood ratio test=124 on 6 df, p=<2e-16
n= 686, number of events= 171根据结果，可以看出nodes是一个显著的预测因子，而prog_recp在线性平滑函数下对风险没有显著影响，但在非线性平滑函数下具有显著影响。
还可以使用pspline()函数在R代码里主要是修改了COX回归的基线风险函数，这样，基线风险函数不再是固定的，而是可以随时间变化。
fit5 <- coxph(Surv(time, status) ~ nodes + prog_recp + pspline(time),
data = data
)
fit5 Call:
coxph(formula = Surv(time, status) ~ nodes + prog_recp + pspline(time),
data = data)
coef se(coef) se2 Chisq DF p
nodes 3.62e-02 1.06e-02 1.06e-02 1.17e+01 1.00 0.00063
prog_recp -3.75e-03 8.31e-04 8.31e-04 2.04e+01 1.00 6.3e-06
pspline(time), linear -2.83e-02 1.31e-03 1.31e-03 4.65e+02 1.00 < 2e-16
pspline(time), nonlin 2.88e+02 2.99 < 2e-16
Iterations: 3 outer, 60 Newton-Raphson
Theta= 0.273
Degrees of freedom for terms= 1 1 4
Likelihood ratio test=1108 on 5.99 df, p=<2e-16
n= 686, number of events= 171在这个模型中，时间的多项式样条变换（pspline(time)）在统计上是显著的，这表明生存风险是随时间的变化的。在调整了时间的非线性变化（即基线风险的波动）之后nodes和prog_recp的影响在统计上依然显著。
zp5<-cox.zph(fit5)
zp5 chisq df p
nodes 1.97 1.00 0.16
prog_recp 1.07 1.00 0.30
pspline(time) 430.84 3.99 <2e-16
GLOBAL 435.59 5.99 <2e-16plot(zp5)参考
https://zhuanlan.zhihu.com/p/596831293
统计咨询||做COX回归不满足PH假定的处理方法解答
