不满足proportional hazards假定的生存分析
不满足proportional hazards假定的生存分析
2019中秋节我们的全国巡讲的优秀学员开启了她的临床统计笔记分享，前面的两个是：
下面仍然是杜云的投稿
近期做生存分析的数据，遇到了不满足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.000473 不满足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最后的处理方法部分，是浅显的总结部分。
实际操作，尤其是方法二是需要花一些功夫的。
用到的朋友们，若遇到问题或心得，欢迎留言交流。