科研星球

R语言统计与绘图:Cox回归中,不满足PH假定时该怎么处理?

作为一个临床研究工作者,在日常分析数据过程中,我们大量地使用Cox比例风险模型,却往往忽略Cox比例风险模型的一个重要假设-PH假设。这就导致我们在投文章的时候,审稿人经常会要求文章中补充如何判断PH假设的部分。


笔者就曾经遇到过这样的问题,被审稿人要求增加PH假设的判断。审稿人的comments如下:


The assumption for the proportional hazards model might have been violated, as seen in Figure 4 where the KM curves crossed each other, So a single summary measure using a proportional hazards model might not be appropriate.

 

所以本期我们将从三个方面展开,力图一次性地解决PH假设的相关问题,首先介绍什么是PH假设,然后介绍PH假设的判断方法,最后介绍PH假设不满足时该如何处理


本文会尽量使用通俗的语言来说明PH假设的问题,但是还是会保留一定的数学公式,供想深入了解的读者学习,只是想在自己文章中解决PH假设问题的读者可以跳过这些公式。




什么是PH假设


PH假设,即 proportional hazard assumption,指的是在使用Cox比例风险模型处理数据时,干预组与对照组的风险比(hazard ratio)在整个随访期间为常数,不随随访时间的变化而变化。(这里进行了简化处理,因为在投文章的时候reviewer主要关注的就是对于treatment来说PH假设是否成立,其实应该是所有纳入Cox模型的协变量都要满足PH假设)


Cox回归模型如下:

 

0 (2).png


从上面的式子可以看出,当X1,...,Xk一定时,treatment等于1时的hazard为

 

0.png


treatment等于0时的hazard为

 

0 (1).png

 

同理,对于其他协变量X1,...,Xk也一样。




PH假设的判断方法


PH假设的判断方法有两种:一种是基于回归方程的,一种是基于残差的。


1)基于回归方程的判断,这种方法的思路说起来非常简单,假设我们不确定对于treatment来说PH假定是否成立,那么我们不妨假设它是不成立的,这时我们就在Cox回归模型中引入一个treatment和时间的交互项,得到如下回归模型:

 

0 (3).png


0 (4).png一般可以选择t或者log(t)。


接下来就是常规的假设检验了,原假设为0 (5).png,当p值小于0.05时,拒绝原假设,即对于treatment来说PH假设不成立;当p值大于等于0.05,接受原假设,即对于treatment来说PH假设成立。


R代码如下:

# cox model

res <- coxph(Surv(tstop, status) ~ treatment + X1 + X2 + X3, dta)

summary(res)

 

# cox model with treatment * time interaction

## f(t) = t

res.t <- coxph(Surv(tstop, status) ~ treatment + X1 + X2 + X3 + tt(treatment), dta, 

               tt = function(x, t, ...) x * t)

summary(res.t)

## f(t) = log(t)

res.logt <- coxph(Surv(tstop, status) ~ treatment + X1 + X2 + X3 + tt(treatment), dta, 

                  tt = function(x, t, ...) x * log(t))

summary(res.logt)


注意,res.t <- coxph(Surv(tstop, status) ~ treatment + X1 + X2 + X3 + treatment * t, dta)是错误的,R也会给出如下错误信息


Error in model.frame.default(formula = Surv(tstop, status) ~ treatment +  : 

  参数't'的种类(closure)不对


2)基于残差的,Cox模型中用来判断PH假设的残差为Schoenfeld residual。

为了展开说明的方便,这里有一个重要结论需要记住,假设对于treatment来说,hazard ratio随时间变化的函数形式为0 (8).png,那么treatment在各个时间点k上scaled Schoenfeld residuals的期望0 (7).png(r代表残差,S代表Schoenfeld,*代表scaled)随时间变化的函数形式近似为0 (9).png


因此,验证PH假设是否成立,这里只需要拟合0 (7).png0 (4).png的回归模型,判断原假设0 (5).png即可,原理同上。同样,0 (4).png一般可以选择t或者log(t)。


R代码如下


# Schoenfeld residual test

## f(t) = t

cox.zph(res, transform = 'identity')

## f(t) = log(t)

cox.zph(res, transform = log)


但是在实际数据分析过程中,仅仅通过基于回归方程或基于残差的假设检验是不够的,有时甚至得到误导性的结果,原因主要有二,一是假设检验的p值受样本量的影响较大,二是从上面描述中我们可以知道,我们必须先指定的函数形式之后,才能进行假设检验,但问题是我们如何保证指定的的函数形式是正确的呢,这里存在较大的model misspecification的风险。


所以,在实际数据分析过程中还要结合图示法对PH假设进行判断。


1) Cumulative incidence plot,,其中


结论:图上两条cumulative incidence curve相交-->PH假设不满足,PH假设不满足-\->图上两条cumulative incidence curve相交,即图上两条cumulative incidence curve相交是PH假设不满足的充分非必要条件。


说明:这里为什么不是0 (13).png,而是0 (6).png,这主要是与我们分析临床研究数据时的习惯相匹配,我们分析临床研究数据时,一般会在unadjusted analysis部分绘制cumulative incidence plot,即我们常说的KM曲线。但是无论是0 (13).png,还是0 (6).png上述结论均成立。


R代码如下:

# cumulative incidence plot

plot(survfit(Surv(tstop, status) ~ treatment, data=dta), fun = 'F', conf.int = T, 

     xlab = 'Days', 

     ylab = 'Cumulative Incidence', col = c('red', 'blue'))


2) Log cumulative hazard plot,0 (13).png,其中0 (14).png


结论:图上两条log cumulative hazard curve平行-->PH假设满足,PH假设满足-->图上两条log cumulative hazard curve平行,即图上两条log cumulative hazard curve平行是PH假设的充分且必要条件。


R代码如下:

# log cumulative hazard plot

plot(survfit(coxph(Surv(tstop, status) ~ X1 + X2 + X3 + strata(treatment), data = dta)), 

     fun = 'cloglog', conf.int = T, 

     xlab = 'Days', 

     ylab = 'Log Cumulative Hazard', col = c('red', 'blue'))


3) Schoenfeld residual plot,即0 (15).png


结论:图上0 (16).png均匀分布在y=0两侧-->PH假设满足,PH假设满足-->图上0 (16).png均匀分布在y=0两侧,即图上0 (16).png均匀分布在y=0两侧是PH假设的充分且必要条件。


该图与上面两个图相比,可以看出hazard ratio随时间变化的趋势。


R代码如下:


# shoenfeld residual plot

scaled_sches <- residuals(res, type = "scaledsch" )

scaled_sch_for_trt <- scaled_sches[ , 1]

time <- as.numeric(rownames(scaled_sches))

plot(time, scaled_sch_for_trt, 

     xlab = "Days", 

     ylab = "Scaled Schoenfeld Residual (treatment)")

lines(smooth.spline(time, scaled_sch_for_trt), col = "red", lwd = 2)


综上,判断PH假设是否成立时需要综合假设检验法和图示法。


如在审稿人给了如下comments的文章中:


The assumption for the proportional hazards model might have been violated, as seen in Figure 4where the KM curves crossed each other, So a single summary measure using a proportional hazards model might not be appropriate.


笔者实际上做了基于残差的假设检验,p值为0.27,从假设检验的结果来看似乎满足PH假设,但是我们下结论之前应该万分小心,还需要结合图示法进行判断。

 

0.jpeg


Cumulative incidence plot如上所示,可以看见两条曲线有交叉,因此PH假设不满足,这也是审稿人给的comments。




PH假设不满足时的处理措施


当PH假定不满足时,我们应该怎么处理呢?在上面“PH假设的判断方法”中我们实际上已经介绍了处理方法,在模型中引入treatment和时间的交互项,即0 (17).png


但是在临床研究类文章中,作者需要平衡所用统计模型的严密性和结果可解释性之间的关系,因此在0 (18).png函数模型选择时,我们往往倾向于选择指示函数0 (19).png,即在0 (20).png时,0 (21).png,在0 (22).png时,0 (23).png


这样的话就可以得到两个hazardratio值,在0 (26).png之前有一个,之后有一个。当然也可以设置多个时间分割点,0 (24).png0 (25).png


R代码如下:

# deal with PH assumption unmet

dta2 <- survSplit(Surv(tstop, status) ~ ., dta, cut = c(0.5), episode = "tgroup")

res2 <- coxph(Surv(tstart, tstop, status) ~ treatment * strata(tgroup) + 

                X1 + X2 + X3, data = dta2)

summary(res2)


切分节点的选择,即图片如何确定,常见的作法是:


观察cumulative incidence plot上两曲线相交的点,选此点对应的时间作为切分节点,至于此种选择方法数学上的意义何在,笔者没有找到合适的证明方法,但是似乎已经成为约定俗成的一种常规处理方式了。


本期生存数据生成的R代码如下

generate survival data

set.seed(20190421)

n = 1000

treatment <- rbinom(n, 1, 0.5)

X1 <- rbinom(n, 1, 0.3)

X2 <- rbinom(n, 1, 0.4)

X3 <- rbinom(n, 1, 0.1)

dta <- data.frame(id = NULL, tstart = NULL, tstop = NULL, 

                  treatment = NULL, X1 = NULL, X2 = NULL, X3 = NULL, 

                  status = NULL, Ttime = NULL)

for (i in 1 : 1000) {

  for (j in seq(0, 1.98, 0.02)) {

    Ttime <- rexp(1, rate = exp(1.5 * treatment[i] - 

                                  1.6 * X1[i] - 2.6 * X2[i] - 3.6 * X3[i] + 

                                  2 * j * treatment[i]))

    tmp.dta = c(id = i, tstart = j, tstop = j + 0.02, 

                treatment = treatment[i], X1 = X1[i], X2 = X2[i], X3 = X3[i], 

                status = ifelse(Ttime <= j + 0.02, 1, 0), Ttime = Ttime)

    dta <- rbind(dta, tmp.dta)

    names(dta) <- c('id', 'tstart', 'tstop', 'treatment', 'X1', 'X2', 'X3', 

                    'status', 'Ttime')

    if (Ttime <= j + 0.02) break

  }

}

dta <- dta %>% 

  group_by(id) %>% 

  summarise(tstart = min(tstart), 

            tstop = max(tstop), 

            status = sum(status), 

            treatment = unique(treatment), 

            X1 = unique(X1), 

            X2 = unique(X2), 

            X3 = unique(X3)) %>% 

  ungroup() %>% 

  data.frame()

rm(treatment, X1, X2, X3)



没有账号?