基于R语言实现COX模型诊断

一般在建立好Cox模型之后，需要对模型进行诊断。诊断内容包括模型的前提条件，诸如Cox模型的PH假定（比例风险假定），共线性假定等。本篇我们通过合实际例子讲解Cox模型诊断过程，实现软件R语言。

1.1 COX模型的诊断内容

Cox模型的诊断一般包括三方面的内容：

比例风险假定；

模型影响点（异常值）识别；

比例风险的对数值与协变量之间的非线性关系识别；

对上述三方面的诊断，常见的方法为残差法。

Schoenfeld残差用于检验比例风险假定；

Deviance残差用于影响点（异常值）识别；

Martingale残差用于非线性检验；

1.2 R中用于评估Cox模型的包

我们将会用到以下两个包：

survival #用于cox模型建立

survminer #用于cox模型诊断结果的可视化

安装包

install.packages(c("survival","survminer"))

加载包

library("survival")

library("survminer")

1.3 建立Cox模型

我们利用survial包中自带的肺癌数据“data(lung)”建立cox模型。

library("survival")

res.cox <- coxph(Surv(time, status) ~ age + sex +wt.loss, data =lung)#模型中有三个变量；

res.cox#显示模型结果

Call:

coxph(formula = Surv(time, status) ~ age + sex + wt.loss,data = lung)

 

            coefexp(coef) se(coef)     z      p

age      0.02009   1.02029 0.00966  2.08 0.0377

sex     -0.52103   0.59391 0.17435 -2.99 0.0028

wt.loss  0.00076   1.00076 0.00619  0.12 0.9024

 

Likelihood ratio test=14.7 on 3 df, p=0.00212

n= 214, number of events= 152

   (14 observationsdeleted due to missingness)

1.4 模型诊断——PH假定

PH假定可通过假设检验和残差图检验。正常情况下，Schoenfeld残差应该与时间无关，如果残差与时间有相关趋势，则违反PH假设的证据。残差图上，横轴代表时间，如果残差均匀的分布则，表示残差与时间相互独立。

R语言survival包中的函数cox.zph（）可以实现这一个检验过程。

test.ph <- cox.zph(res.cox)

test.ph

            rhochisq     p

age     -0.0483 0.3780.538

sex      0.1265 2.3490.125

wt.loss  0.0126 0.0240.877

GLOBAL       NA 2.8460.416

从上面的结果可以看出，三个变量的P值都大于0.05，说明每个变量均满足PH检验，而模型的整体检验P值0.416，模型整体满足PH检验。

在R语言 survminer中ggcoxzph( )函数可以画出Schoenfeld残差图。

ggcoxzph(test.ph)

上图中实线是拟合的样条平滑曲线，虚线表示拟合曲线上下2个单位的标准差。如果曲线偏离2个单位的标准差则表示不满足比例风险假定。从上图中可见，各协变量满足PH风险假设。

另一种检查比例风险假定的图形方法是绘制log（-log（S（t）））与t或log（t）是非平行，这个方法只能用于协变量是分类变量的情形。

如果违反比例风险假设可以通过以下方式解决：

模型中添加协变量与时间的交互相应；

分层分析；

至于如何实现，我们后期再做介绍。

1.5 模型诊断——模型影响点（异常值）识别

我们可以通过绘制Deviance残差图或者dfbeta值实现上述诊断。在R语言survminer中ggcoxdiagnostics()函数可以画出Deviance残差图。

ggcoxdiagnostics(res.cox,type = "deviance",

                 linear.predictions = FALSE,ggtheme = theme_bw())

残差值均匀的分布在0上下，表明满足上述假定。

ggcoxdiagnostics(res.cox,type = "dfbeta",

                 linear.predictions = FALSE,ggtheme = theme_bw())

影响点的可能来源于数据录入错误，样本中的极值点、协变量不均衡，数据不足等。对本例，上图显示，将dfbeta值大小与回归系数比较表明，即使某些dfbeta值非常大，但它们不足以对模型系数的估计值产生影响。

1.6 模型诊断——非线性诊断

一般情况下，我们假设协变量与-log(s(t))之间是线性关系。通过绘制Martingale残差图可以实现模型协变量的非线性诊断。非线性诊断一般是针对模型中的连续型变量。

在R语言survminer中ggcoxfunctional()函数可以画出Martingale残差图。

ggcoxfunctional(Surv(time, status) ~ age + log(age) + sqrt(age),data = lung)

图中显示年龄局部有非线性趋势，但整体表现出线性趋势。