R语言生存分析04-Cox比例风险模型诊断

作者:白介素2

相关阅读:
R语言生存分析04-Cox比例风险模型诊断
R语言生存分析03-Cox比例风险模型
R语言生存分析-02-ggforest
R语言生存分析-01
ggpubr-专为学术绘图而生(二)
ggstatsplot-专为学术绘图而生(一)
生存曲线
R语言GEO数据挖掘01-数据下载及提取表达矩阵
R语言GEO数据挖掘02-解决GEO数据中的多个探针对应一个基因
R语言GEO数据挖掘03-limma分析差异基因
R语言GEO数据挖掘04-功能富集分析

如果没有时间精力学习代码,推荐了解:零代码数据挖掘课程

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)
image.png
  • 上图中实线是与曲线拟合的平滑样条曲线,虚线表示拟合周围的+/- 2标准误差带。
  • 没有与时间相关变化模式,个各个协变量满足风险比例假设

检验异常的的观测

  • 绘制Deviance残差图或者dfbeta值实现,以下选择dfbeta,改为deviance即残差图
  • survminer中的ggcoxdiagnostics()函数
ggcoxdiagnostics(res.cox, type = "dfbeta",
                 linear.predictions = FALSE, ggtheme = theme_bw())
image.png

上图表示,将最大dfbeta值的大小与回归系数进行比较表明,没有一个观察结果是单独影响的,即使年龄和重量损失的某些dfbeta值与其他值相比较大。

非线性诊断- non linearity

  • 通常,我们假设连续协变量具有线性形式。但是,应该检验这个假设是否成立。 使用连续协变量绘制Martingale残差是用于检测非线性的常用方法,或者换句话说,用于评估协变量的函数形式。对于给定的连续协变量,图中的模式可能表明变量不适合。
  • R语言survminer中的ggcoxfunctional()函数可以绘图
  • 例如我们检验age变量,可使用如下代码
ggcoxfunctional(Surv(time, status) ~ age + log(age) + sqrt(age), data = lung)
image.png
  • 结果显示,有一定程度的非线性存在

参考资料

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 212,332评论 6 493
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 90,508评论 3 385
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 157,812评论 0 348
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 56,607评论 1 284
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 65,728评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 49,919评论 1 290
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,071评论 3 410
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 37,802评论 0 268
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,256评论 1 303
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,576评论 2 327
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,712评论 1 341
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,389评论 4 332
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 40,032评论 3 316
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,798评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,026评论 1 266
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 46,473评论 2 360
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 43,606评论 2 350

推荐阅读更多精彩内容