代码如下:

rm(list=ls())
cat("\014") 
library(survminer)
library(survival)
data1=cancer
str(data1)

datasurv<- Surv(time=data1$time,event=data1$status==2)
data1$datasurv<-with(data1,datasurv)
colnames(data1)

Multicox1<- coxph(datasurv~age+sex+ph.ecog+ph.karno+meal.cal+wt.loss,data=data1)
summary(Multicox1)
testPH1 <- cox.zph(Multicox1) 
print(testPH1)
ggcoxzph(testPH1)
ggforest(model = Multicox1, data = data1, cpositions = c(0.05, 0.15, 0.35),refLabel = 'reference', noDigits = 3,fontsize = 0.8)

残差检验结果如下:

print(testPH1)
            chisq df     p
age       0.21311  1 0.644
sex       1.24215  1 0.265
ph.ecog   2.40880  1 0.121
ph.karno  5.82165  1 0.016
meal.cal  6.44932  1 0.011
wt.loss   0.00982  1 0.921
GLOBAL   13.19744  6 0.040

结果显示,“ph.karno”、“meal.cal ”以及global都是P<0.05,不符合PH假定。

  1. 要检验PH假定,仅靠cox.zph()足够吗?分类变量和连续变量我都可以用这个函数吗?
  2. 残差图没看懂,感觉所有变量的曲线都很平,为什么有的P<0.05, 有的P>0.05呢?
  3. 假如“ph.karno”和“meal.cal ”不是我关注的因子,我可以简单粗暴的删掉它们吗? 假如它们是我关注的因子(或者我不想删掉它们),那我应该怎么办,把它们转换成分类变量再试一次,或者把其中一个strata(我尝试过好像不能把两个都strata),或者做依时协变量?
7 天 后
  1. cox.zph() 只是提供了一个最基本的标准方法,想做的话当然可以做更多其他检验。CRAN 生存分析 task view 上列举了一些,"The proportionality assumption can be checked" 一段。这里从 coxph() 开始好像是都直接作为连续型处理了。
  2. 这种诊断图也就是起个辅助作用。和检验结果数值没有绝对敏感的对应关系,主要是为了在直观上排除明显的随时间变化的趋势,所以其实用处有限,最后以 p 值为准吧。另外这几个 p 值说实话比较边缘,如果你把阈值设为 0.01,那是不是一下就都不显著了,全部符合 PH 假定,哈哈哈哈 ……
  3. 不是关注的变量在一开始就不应该放进去建模,不过这样做是否合理要看你建模的目的。另外要删变量也不需要在这一步删,在模型 summary 里 Wald test 的 age、meal.cal、wt.loss 就已经不显著了 —— 我觉得先确定一个大致的模型,把那些明显无关的变量去掉,再检验所得模型是否符合各种详细的假定,这样顺序比较合理,而不是反过来 —— 毕竟检验 PH 假定是需要先确定一个模型的。在 cox.zph() 这一步如果一定要把违反 PH 假定的变量留在模型里,就像你说的,要么是增加变量 * 时间的交互作用,要么是做 stratification。