COX时变系数模型（含时间协变量的Cox回归模型）的R语言实现

Cox回归模型有效地解决了对生存资料进行多因素分析的问题，但是应用Cox回归模型有一个非常重要的前提条件，即比例风险（Proportional hazards）假设，简称PH假设，其基本假设为：协变量对生存率的影响不随时间的改变而改变。只有当PH假定得到满足时，Cox回归模型的结果才有意义。

比例风险假设可以通过检验模型的残差来检验。零假设的拒绝导致使用时变系数来描述数据。时变系数可以用阶跃函数或参数时间函数来描述。

当数据不险合比例风险假设时，可以引入时间变化分析，这是一种基于cox回归模型的拓展，通过加入时间协变量以估计对生存的影响。

时间协变量可以分为内部和外部两种，内部协变量通常是所研究个体生成的随机过程的输出，并且仅在受试者生存且未经审查的情况下才能观察到，外部协变量X（X）可能会影响随时间推移的故障率，但其路径直到时间t >v不受时间v发生故障时间的影响。它也是一个衍生或预定的协变量。一般外部协变量比较常见，即系数β随着时间变化。

这一时变系数可以通过两种方式来描述。使用阶跃函数就是将分析时间分成几个区间，进行Cox比例模型分层。固定基线协变量的影响随着时间的推移变得更强或更弱，这可以通过按时间分层来探索。另一种方法是使用由用户指定的参数连续函数构造一个交互项，r语言的在coxph（）函数中，有一个tt参数来指定时间的特定转换。

本文以survival包自带的肺癌数据集距离，我们先导入R包和数据

library(survival)
lung<-lung
View(lung)

在数据集中，inst为机构代码，time表示以天为单位的生存时间，status为结局变量，age为 年龄，sex: 男=1 女=2，ph.ecog: 由医师评定的 ECOG 表现评分。 0 = 无症状，1 = 有症状但完全不卧床，2 = 卧床时间 < 50%，3 = 卧床时间 > 50% 但不卧床，4 = 卧床，ph.karno: 由医师评定的 Karnofsky 表现评分（差=0-好=100），pat.karno: 由患者评定的 Karnofsky 性能评分，meal.cal: 用餐时消耗的卡路里，wt.loss: 过去六个月的体重减轻（磅）

接下来，可以拟合cox模型并使用survival包自带的cox.zph函数检验PH假设

#拟合cox模型
fit2 <- coxph(Surv(time, status) ~age +ph.karno+sex,data=lung)
#检查PH假设
zph <- cox.zph(fit2)
zph

 

 在输出结果中可以看到ph.karno的P值小于0.05，表明ph.karno没有通过PH假设，全局显着性检验给出的 P 值 (0.0157) 表明模型缺乏拟合。我们可以通过plot函数对ph.karno这个变量进行可视化

plot(zph[2],lwd=2)
abline(0,0, col=1,lty=3,lwd=2)
abline(h= fit2$coef[2], col=3, lwd=2, lty=2)
legend("bottomright",
       legend=c('Reference line for null effect',
                "Average hazard over time",
                "Time-varying hazard"),
       lty=c(3,2,1), col=c(1,3,1), lwd=2)

其中黑色虚线表示系数为0，就是ph.karno无效的情况，绿色虚线表示的是ph.karno的平均效应值，我们可以看到均线和实际背离是相当大的，因此这个情况使用cox.zph 函数来处理是非常恰当的。
对时变系数建模的一种方法是使用阶跃函数，例如 (g(t) = I(t≥to))，其中to是指定值。

该方法的思想是将分析时间划分为多个区间，并对这些时间区间进行Cox比例模型分层。固定基线协变量的影响随着时间的推移而增强或减弱，这可以通过时间分层来探索，基线风险因素ph_karno的影响随着时间的推移而变化，从而导致一系列hr。使用survSplit()函数，可以按照我们看到的ph.karno大约在180和350左右的转折，将每个记录分割成子记录。

lung.split <- survSplit(Surv(time, status) ~ .,
                        data= lung, cut=c(180, 350),
                        episode= "tgroup", id="id")
head(lung.split[-c(1,4,6:8)])

新增指标tgroup对time的时区进行了划分，分割成 (0, 180), (180, 350) 和 (350, 455) 的3个区间，tgroup=1 标识第一个时间间隔(0, 180)，tgroup=2 标识第二个时间间隔(180, 350)，3就是(350, 455)。
现在我们把lung.split的数据拿来重新进行模型拟合，建立一个ph.karno和tgroup的交互分层并重新进行PH假设的检验。

#重新拟合，建立交互分层
fit.split <- coxph(Surv(tstart, time, status) ~
                     age + ph.karno:strata(tgroup)+
                     sex,
                   data=lung.split)
fit.split
#重新进行PH假设评估
cox.zph(fit.split)

可以看到tgroup=1的时候ph.karno的P值是有意义的，HR=0.965594 

这次全部通过PH假设了，GLOBAL=0.2（大于0.05）表明模型是拟合的。

 另一种方式是使用由用户指定的参数连续函数构造一个交互项，在 coxph() 函数中，有一个tt参数来指定时间的具体变换。在本例中，tt函数定义为“tt = function(x, t, …) x * log(t+20)”，其中 x 是具有时变效应的固定协变量，t是分析时间。tt() 函数应用于变量ph.karno。

fit.tt <- coxph(Surv(time, status) ~
                  age + ph.karno + tt(ph.karno)+ sex,
                data=lung,tt = function(x, t, ...) x * log(t+20))
fit.tt

ph.karno 和 tt(ph.karno) 的系数都具有统计显着性，这意味着 ph.karno 的影响随时间而变化。ph.karno 的时变效应可以写为 β(t) = -0.098+0.015×log( t + 20)。
我们可以使用 abline() 函数在 ph.karno 对生存的时变影响的 cox.zph 图中添加一条线。

zph.tt <- cox.zph(fit2,
                  transform=function(t) log(t+20))
plot(zph.tt[2])
abline(coef(fit.tt)[2:3], col=2)

 

一个参数时间函数被分配给ph.karno。如果用函数log(t+20)变换时间轴，则效果为线性，斜率为0.015(红线)。 与水平线(斜率=0)有显著差异。黑线表示变量ph.karno的时变系数。