R语言生存分析-Cox比例风险模型诊断
欢迎来到医科研，这里是白介素2的读书笔记，跟我一起聊临床与科研的故事, 生物医学数据挖掘，R语言，TCGA、GEO数据挖掘。
Cox比例风险模型诊断
Cox比例风险模型的建立是基于几个假设之上的，因此一般建好模型后需要进行诊断，评估拟合的模型是否能够用于描述数据
诊断的内容包括:
比例风险假定；
模型影响点（异常值）识别；
比例风险的对数值与协变量之间的非线性关系识别；
对上述三方面的诊断，常见的方法为残差法。
Schoenfeld 残差用于检验比例风险假定；
Deviance 残差用于影响点（异常值）识别；
Martingale残差用于非线性检验；
载入survival and survminer
library("survival")
library("survminer")
计算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)
##
## coef exp(coef) se(coef) z p
## age 0.0200882 1.0202913 0.0096644 2.079 0.0377
## sex -0.5210319 0.5939074 0.1743541 -2.988 0.0028
## wt.loss 0.0007596 1.0007599 0.0061934 0.123 0.9024
##
## Likelihood ratio test=14.67 on 3 df, p=0.002122
## n= 214, number of events= 152
## (14 observations deleted due to missingness)
检验比例风险假设-PH假设
PH假设可通过假设检验和残差图检验。正常情况下，Schoenfeld残差应该与时间无关，如果残差与时间有相关趋势，则违反PH假设的证据。残差图上的横轴代表时间，如果残差均匀的分布则表示残差与时间相互独立。
R语言survival包中的函数cox.zph函数提供简便的实现这一过程的方法
test.ph <- cox.zph(res.cox)
test.ph
## rho chisq p
## age -0.0483 0.378 0.538
## sex 0.1265 2.349 0.125
## wt.loss 0.0126 0.024 0.877
## GLOBAL NA 2.846 0.416
从输出的结果看，三个协变量的P值都大于0.05，说明每个变量均满足PH检验，而模型的整体检验P值0.416也没有统计学意义，因此我们认为模型整体满足PH检验。
图形诊断
survminer包中的ggcoxzph()函数可以绘制每个协变量随时间变化的Schoenfeld残差图
ggcoxzph(test.ph)
上图中实线是与曲线拟合的平滑样条曲线，虚线表示拟合周围的+/- 2标准误差带。
没有与时间相关变化模式，个各个协变量满足风险比例假设
检验异常的的观测
绘制Deviance残差图或者dfbeta值实现,以下选择dfbeta，改为deviance即残差图
survminer中的ggcoxdiagnostics()函数
ggcoxdiagnostics(res.cox, type = "dfbeta",
linear.predictions = FALSE, ggtheme = theme_bw())
上图表示，将最大dfbeta值的大小与回归系数进行比较表明，没有一个观察结果是单独影响的，即使年龄和重量损失的某些dfbeta值与其他值相比较大。
非线性诊断- non linearity
通常，我们假设连续协变量具有线性形式。但是，应该检验这个假设是否成立。使用连续协变量绘制Martingale残差是用于检测非线性的常用方法，或者换句话说，用于评估协变量的函数形式。对于给定的连续协变量，图中的模式可能表明变量不适合。
R语言survminer中的ggcoxfunctional()函数可以绘图
例如我们检验age变量，可使用如下代码
ggcoxfunctional(Surv(time, status) ~ age + log(age) + sqrt(age), data = lung)
结果显示，有一定程度的非线性存在
参考资料
(http://www.sthda.com/english/wiki/cox-model-assumptions)