【原】不满足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./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

最后的处理方法部分，是浅显的总结部分。

实际操作，尤其是方法二是需要花一些功夫的。

用到的朋友们，若遇到问题或心得，欢迎留言交流。