登录
首页精彩阅读基于R语言实现COX模型诊断
基于R语言实现COX模型诊断
2017-05-11
收藏

基于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)

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

数据分析咨询请扫描二维码

客服在线
立即咨询