不满足proportional hazards假定的生存分析
近期做生存分析的数据,遇到了不满足proportional hazards
assumption (简称PH 假定)的问题,折腾了一圈,总结下这部分内容。
- 什么是proportional hazards?
- 如何判断是否满足PH 假定?
- 遇到不满足PH 假定的如何处理。
我们一步步来看,希望对大家有用。若有错误,欢迎批评指正。
注:有代码的部分,用的是r语言。
1 什么是proportional hazards (PH) 假定。
1.1 hazard曲线。
从字面上理解,成比例的hazards,也就是两组的hazards的比值是恒定的。
举个例子,如下图。随访10年以后两组的hazard曲线有交叉,显然不是恒定的PH,不满足PH假定。
而在曲线交叉前,肉眼看是基本满足PH假定的 (随访刚刚开始的那部分可以不计)。
1.2 与随访时间的关系。
从上面可以知道,PH假定也就是两组hazard的比值不随着随访时间的变化而变化,是一个恒定的值。
再深入一下,Drug组与Placebo组的HRs,不随随访时间变化而变化,也就是Drug的作用与随访时间没有交互作用 (interaction)。
1.3 何时需要满足PH假定?
Kaplan Meier 曲线,比较组间的long-rank test。需满足PH假定。
Poisson regression,估计HRs,需满足PH假定。
Cox regression.需满足PH假定。
可以看出,我们常用的生存分析模型,都有PH的假定。
所以,在我们用这些模型时,要判断,我们感兴趣的两组的HRs,是否满足PH假定。
2 如何判断是否满足PH 假定?
2.1 曲线直观判断
2.1.1 第一种,hazard curve。
上面的例子图形就是一种图形判断方法。但是这种方法比较粗糙,不相交的曲线也可能不满足PH假定。
2.1.2 第二种, log (cumulative hazard)图形
log (cumulative hazard)作为y轴, 随访时间作为x轴。若两组的曲线基本平行,则满足PH假定。如下图。
从上图可以看出,两组曲线基本平行的,满足PH假定。这种方法是比较直观靠谱的图形判断方法。
图形判断方法虽然直观,但是给不了我们钟爱的p 值。下面介绍可以给出p 值的方法。
2.2 检验interaction term
上面我们提到,PH假定等同于某个因素对hazard的作用不与随访时间发生交互作用。所以,我们可以检验加入interaction term和不加入interaction的两个model是否有统计学差异。
举个例子
anova(fit1,fit2)
## Analysis of Deviance Table
## Cox model: response is Surv(start, trunc_yy, death_cancer)
## Model 1: ~ agegrp:fu + agegrp + fu
## Model 2: ~ agegrp + fu
## loglik Chisq Df P(>|Chi|)
## 1 -7830.8
## 2 -7836.0 10.319 3 0.01604 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
两个model比较,得到的p值为0.016.也就是说,年龄对预后的作用,随着随访时间变化而变化。不满足PH假定。
2.3 Schoenfeld’s residuals 检验
这种方法适用于cox 回归。看结果中的p值。该方法的H0假设是满足PH假定,p值小于0.05,就是拒绝H0,不满足PH假定了。
分享一个链接http://www.sthda.com/english/wiki/cox-model-assumptions,理解更多。
coxph(Surv(trunc_yy,death_cancer)~year8594+sex+agegrp+stage, data=localised) %>% cox.zph()
## chisq df p
## year8594 2.632 1 0.10471
## sex 0.823 1 0.36431
## agegrp 16.858 3 0.00076
## GLOBAL 22.247 5 0.00047
3 不满足PH假定如何处理。
3.1 方法一:随访时间截断
把PH假定不满足随访时间去掉,比如Figure 1中,随访10年后的时间,PH假定肯定不满足。那么我们可以把随访的最长时间定为10年,得到10-year 的HRs。这种方法用于随访时间足够,而去掉的随访时间不影响我们对研究问题的回答。
3.2 方法二:加入interaction term
把随访时间分成几段,加入interaction term,得到不同时间段的HRs。比如上面的model1.
3.3 方法三:Stratified cox 回归
由于Cox 回归比较常用,介绍下这种方法。把不满足PH假定的变量放到strata()
参数中。这种方法适用于不满足PH假定的是exposure之外的变量。因为放到strata()
当中的变量,我们得到该变量对结局的coefficient了。
比如, 下面的这个例子,agegrp
不满足PH假定,如果我们感兴趣的是stage
对hazards的作用,完全可以把agegrp
放进strata()
. Stratified之后,就没有不满足PH假定的变量了。
coxph(Surv(trunc_yy,death_cancer)~year8594+sex+agegrp+stage, data=localised) %>% cox.zph()
## chisq df p
## year8594 2.632 1 0.10471
## sex 0.823 1 0.36431
## agegrp 16.858 3 0.00076
## GLOBAL 22.247 5 0.00047coxph(Surv(trunc_yy,death_cancer)~year8594+sex+strata(agegrp)+stage, data=localised) %>% cox.zph()
## chisq df p
## year8594 3.71 1 0.054
## sex 1.48 1 0.224
## GLOBAL 5.48 2 0.064
最后的处理方法部分,是浅显的总结部分。
实际操作,尤其是方法二是需要花一些功夫的。
用到的朋友们,若遇到问题或心得,欢迎留言交流。