R语言生存分析:Cox回归

article/2025/8/28 10:48:04

上次介绍了生存分析中的寿命表、K-M曲线、logrank检验、最佳切点的寻找等,本次主要介绍Cox回归。

本推文不涉及理论,只有实操,想要了解生存分析的理论的请自行学习。

Cox回归

使用survival包中的lung数据集用于演示,这是一份关于肺癌患者的生存数据。time是生存时间,以天为单位,status是生存状态,1代表删失,2代表死亡。

rm(list = ls())
library(survival)
library(survminer)str(lung)
## 'data.frame': 228 obs. of  10 variables:
##  $ inst     : num  3 3 3 5 1 12 7 11 1 7 ...
##  $ time     : num  306 455 1010 210 883 ...
##  $ status   : num  2 2 1 2 2 1 2 2 2 2 ...
##  $ age      : num  74 68 56 57 60 74 68 71 53 61 ...
##  $ sex      : num  1 1 1 1 1 1 2 2 1 1 ...
##  $ ph.ecog  : num  1 0 0 1 0 1 2 2 1 2 ...
##  $ ph.karno : num  90 90 90 90 100 50 70 60 70 70 ...
##  $ pat.karno: num  100 90 90 60 90 80 60 80 80 70 ...
##  $ meal.cal : num  1175 1225 NA 1150 NA ...
##  $ wt.loss  : num  NA 15 15 11 0 0 10 1 16 34 ...

可以使用cox回归探索危险因素。分类变量需要变为因子型,这样在进行回归时会自动进行哑变量设置。

lung$sex <- factor(lung$sex, labels = c("female","male"))
lung$ph.ecog <- factor(lung$ph.ecog, labels = c("asymptomatic", "symptomatic",'in bed <50%','in bed >50%'))str(lung)
## 'data.frame': 228 obs. of  10 variables:
##  $ inst     : num  3 3 3 5 1 12 7 11 1 7 ...
##  $ time     : num  306 455 1010 210 883 ...
##  $ status   : num  2 2 1 2 2 1 2 2 2 2 ...
##  $ age      : num  74 68 56 57 60 74 68 71 53 61 ...
##  $ sex      : Factor w/ 2 levels "female","male": 1 1 1 1 1 1 2 2 1 1 ...
##  $ ph.ecog  : Factor w/ 4 levels "asymptomatic",..: 2 1 1 2 1 2 3 3 2 3 ...
##  $ ph.karno : num  90 90 90 90 100 50 70 60 70 70 ...
##  $ pat.karno: num  100 90 90 60 90 80 60 80 80 70 ...
##  $ meal.cal : num  1175 1225 NA 1150 NA ...
##  $ wt.loss  : num  NA 15 15 11 0 0 10 1 16 34 ...

拟合多因素Cox回归模型,这里我们只用sex/age/ph.karno3个变量做演示:

fit.cox <- coxph(Surv(time, status) ~ sex + age + ph.karno, data = lung)# 查看结果
summary(fit.cox)
## Call:
## coxph(formula = Surv(time, status) ~ sex + age + ph.karno, data = lung)
## 
##   n= 227, number of events= 164 
##    (1 observation deleted due to missingness)
## 
##               coef exp(coef)  se(coef)      z Pr(>|z|)   
## sexmale  -0.497170  0.608249  0.167713 -2.964  0.00303 **
## age       0.012375  1.012452  0.009405  1.316  0.18821   
## ph.karno -0.013322  0.986767  0.005880 -2.266  0.02348 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##          exp(coef) exp(-coef) lower .95 upper .95
## sexmale     0.6082     1.6441    0.4378    0.8450
## age         1.0125     0.9877    0.9940    1.0313
## ph.karno    0.9868     1.0134    0.9755    0.9982
## 
## Concordance= 0.637  (se = 0.025 )
## Likelihood ratio test= 18.81  on 3 df,   p=3e-04
## Wald test            = 18.73  on 3 df,   p=3e-04
## Score (logrank) test = 19.05  on 3 df,   p=3e-04

结果解读和logistic回归的结果解读类似:超级详细的logistic细节解读

  • coef是回归系数,

  • exp(coef)是HR值,

  • se(coef)是回归系数的标准误,

  • z是Wald检验的z值,

  • Pr(>|z|)是回归系数的P值,

  • lower .95/upper .95是HR值的95%可信区间。

Concordance= 0.645是Cox回归的C-index,最后给出了Likelihood ratio test似然比检验的统计量、自由度、P值;Wald test的统计量、自由度、P值;Score (logrank) test的统计量、自由度、P值。

想获得整洁的结果不需要自己提取,只要用神包broom即可:

broom::tidy(fit.cox, exponentiate = T, conf.int = T)
## # A tibble: 3 × 7
##   term     estimate std.error statistic p.value conf.low conf.high
##   <chr>       <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
## 1 sexmale     0.608   0.168       -2.96 0.00303    0.438     0.845
## 2 age         1.01    0.00940      1.32 0.188      0.994     1.03 
## 3 ph.karno    0.987   0.00588     -2.27 0.0235     0.975     0.998
  • estimate:HR值(exp(coef))

  • std.error:回归系数的标准误(se(coef))

  • statistic:Wald检验的z值

  • p.value:回归系数的P值

  • conf.low/conf.high:HR的95%的可信区间

构建好Cox回归后,也可以用函数单独提取想要的结果,以下图片展示了可用于提取模型信息的函数,和logistic回归差不多:

进行Cox回归必须要符合等比例风险假设,关于什么是等比例风险假设,可以参考这篇文章:https://mp.weixin.qq.com/s/l1Sd_nB9XY11FZMuda7xBg

等比例风险的检验可以通过很多方法进行,比如K-M曲线,一般如果有交叉,那么可能不符合等比例风险假设,还可以通过各种残差分布来检验。

下面是Cox回归的等比例风险假设检验,检验方法是基于Schoenfeld残差:

ftest <- cox.zph(fit.cox)
ftest
##           chisq df      p
## sex       3.085  1 0.0790
## age       0.478  1 0.4892
## ph.karno  8.017  1 0.0046
## GLOBAL   10.359  3 0.0157

可以看到ph.karno的P值是小于0.05的,其实是不满足等比例风险假设的,下一篇推文会说到不符合等比例风险假设时该怎么办。

这种方法是基于Schoenfeld残差,检验结果可以通过图示画出来:

library(survminer)ggcoxzph(ftest)

 

可以看到sexage的回归系数随着时间变化基本没啥变化,稳定在0水平线上,和上面的检验结果一样。

还可以通过以下方式查看残差的变化:

ggcoxdiagnostics(fit.cox, type = "schoenfeld")
## `geom_smooth()` using formula 'y ~ x'

这张图反映的也是回归系数随时间的变化趋势,和上面的图意思一样,如果符合比例风险假设,那么结果应该是一条水平线,从图示来看,这3个变量都是有点问题的,但是真实数据往往不可能是完美的,很少有完全符合要求的数据。

除了Schoenfeld残差外,ggcoxdiagnostics()还支持其他类型,比如:"martingale", "deviance", "score","dfbetas" and "scaledsch"在,只需要在type参数中提供合适的类型即可。

cox回归也是回归分析的一种,可以计算出回归系数和95%的可信区间,因此结果可以通过森林图展示:

# 为了森林图好看点,多选几个变量
fit.cox <- coxph(Surv(time, status) ~ . , data = lung)ggforest(fit.cox, data = lung,main = "Hazard ratio",cpositions = c(0.01, 0.15, 0.35), # 更改前三列的相对位置fontsize = 0.7,refLabel = "reference",noDigits = 2)

这个结果如果你觉得不好看,或者你还有其他的森林图想做到统一的样式,可以考虑之前的介绍的画森林图的方法进行个性化定制:

画一个好看的森林图

用更简单的方式画森林图

R语言画森林图系列3!

R语言画森林图系列4!

以上是Cox回归的主要内容,大家有问题可以加群或者评论区留言,下次继续介绍时依协变量Cox回归和时依系数Cox回归。

参考资料

  1. http://www.sthda.com/english/wiki/survival-analysis-basics

  2. https://www.emilyzabor.com/tutorials/survival_analysis_in_r_tutorial.html

  3. survival包帮助文档

  4. https://mp.weixin.qq.com/s/2rwxeaF_M0UnqPi2F9JNxA


http://chatgpt.dhexx.cn/article/cj3YZbno.shtml

相关文章

SPSS教学—Cox回归模型探索多因素对生存期的影响

Cox回归模型又称为比例风险回归模型&#xff0c;该模型以生存结局和生存时间作为因变量&#xff0c;进而分析众多因素对生存期的影响&#xff0c;是一个典型的多因素分析方法。 SPSS中就带有Cox回归模型方法&#xff0c;本节将带大家进行深入的了解与探索&#xff0c;话不多说…

oracle dump enq hw,等待事件enq:HW–contention说明及解决方法

一、今天在查看awr报告中&#xff0c;发现Top 5 Timed Foreground Events发现enq: HW - contention的等待事件&#xff1b; 二、enq: HW - contention的官方说明&#xff1a; The HW enqueue is used to serialize the allocation of space beyond the high water mark of a se…

关于enq: TX - allocate ITL entry的问题分析

今天发现系统在下午1点左右的时候负载比较高,就抓取了一个最新的awr报告. Snap Id Snap Time Sessions Cursors/Session Begin Snap: 20892 26-Nov-14 13:20:17 3623 5.4 End Snap: 20893 26-Nov-14 13:30:17 3602 5.4 Elapsed: 10.01 (mins) DB Time…

enq: TX - index contention

解决方案&#xff1a;alter system set “_lm_drm_disable”5 sid’*’;&#xff08;重启库&#xff09;

oracle enq: tm,Tuning enq: TM – contention with foreign key (外键引起的队列)

TM – Enqueue contention 与Table Manipulation相关的入队争用&#xff0c;可以在使用需要锁定表的重组活动显式锁定表时看到。 ID1 ID2的含义 ID1 : 0(LGWR) or object_number&#xff0c; 即DBA_OBJECTS.OBJECT_ID ID2: 0 for a normal table / partition lock ; 1 for…

Oracle-enq: TX - row lock contention 等待事件分析

什么是enq:TX - row lock contention等待: 等待事件enq:TX - row lock contention 是Oracle常见的几大等待事件之一,在开启的事务中&#xff0c;为了维护事务数据的一致性&#xff0c;会在事务所涉及的修改行中添加TX锁以防止其他会话同时修改数据&#xff0c;当其他会话等待该…

Oracle死锁问题: enq: TX - row lock contention

前言 这篇文章也是记录近期遇到的问题以及从中学到的知识 &#xff0c;近期一直在救火&#xff0c;有些问题自认为还是挺有代表性的&#xff0c;有兴趣的话再继续向下看 问题现象 线上反馈&#xff0c;执行批量处理EXCEL数据时&#xff0c;系统一直卡在进度滚动条界面。处理任务…

oracle dump enq hw,enq:HW–contention 故障处理

enq: HW - contention 说明&#xff1a; 硬件队列用于序列化段的高水位线以外的空间分配。 可以用手动分配范围即可解决问题。 V$SESSION_WAIT,P2 / V$LOCK.ID1 is the tablespace number. V$SESSION_WAIT.P3 / V$LOCK.ID2 is the relative data block address (dba) of segmen…

oracle enq ta,Oracle 的 enq: TT - contention 等待事件

Oracle 的 enq: TT - contention 等待事件 在遇到 TT 锁等待时,你可能会被Oracle的文档所困扰。即便是在Oracle Database 12c的文档中,关于TT锁的描述也是:Temporary table enqueue。 这其实已经不准确了,从数据库中可以得到更详细和准确的描述,V$LOCK_TYPE中有着准确记录…

oracle enq ta,enq: TA – contention 等待事件

enq: TA – contention 等待事件 This enqueue is used when undo tablespace operations are being performed. Some examples of such operations are: When dropping an undo tablespace we acquire the enqueue in exclsuive mode to stop other sessions using the undo t…

关于AQS中的enq方法的理解

自己太笨了&#xff0c;总感觉有点绕&#xff0c;就整理下吧~ private Node enq(final Node node) {//自旋锁for (;;) {//tail默认就是nullNode t tail;if (t null) { // Must initialize//因为tail默认是null&#xff0c;所以首次一定会进来//compareAndSetHead在下面//也就…

队列等待之enq: TX - row lock contention

【性能优化】队列等待之enq: TX - row lock contention 问题背景&#xff1a; 客户反映某条sql DELETE SHAREINNERDOC WHERE SOURCEID:B1<br/>这个执行时间太长 问题解决 1> 查看awr报告&#xff1a; 有队列等待之enq: TX - row lock contention&#xff0c;对应的sq…

等待事件 enq:TX - row lock contention分析与解决

6月30日&#xff0c;数据库发生了大量锁表。大概持续1小时&#xff0c;并且越锁越多。后来通过业务人员停掉程序&#xff0c;并kill掉会话后解决。 几天后再EM上查看CPU占用&#xff1a; CPU发生了明显等待。 主要是由于enq:TX - row lock contention等待事件造成。 等待事…

java -- 随机获取字母或者数字

java只有涉及到随机的&#xff0c;最经常用到的方法就是Math.random()&#xff0c;这个方法会返回一个大于0小于1的随机数( 能取0不能取1 )&#xff0c;如果我们要随机0-9&#xff0c;就可以用&#xff08;Math.random()*10&#xff09;来表示&#xff0c;随机0-99也类似如此操…

JavaScript生成随机字母数字字符串

如何使用javascript生成随机字母数字字符串&#xff1f;下面本篇文章就来给大家介绍一下使用JavaScript生成随机字母数字字符串的方法&#xff0c;希望对大家有所帮助。 方法一&#xff1a;Math.random()方法和Math.floor()方法 ● 创建一个函数&#xff0c;该函数有两个参数…

Python - 怎么将一个数字拆分成多个随机数字

前情提要 使用numpy.random.choice()的时候&#xff0c;通过参数p&#xff08;一个列表&#xff09;来指定所给选择元素的选择概率。但参数p&#xff08;选择概率&#xff09;要保证和为1&#xff0c;这时我又想随机生成选择概率&#xff0c;所以现在的问题就是怎么将1拆分成多…

python随机生成一个数字_如何实现python随机生成数字?

今天小编就生成随机数,整理了多个方式,方便大家在项目时,根据自己的需求,直接拿来套用即可,以下内容相当详细,具体来看看吧~ 说明:python中生成随机数主要用到random模块,方法主要包括:randint、uniform、random、sample、choice等几种常用方法; 环境:Mac OS 10.1…

python随机生成三位数字_python3 随机生成数字

原博文 2019-11-25 10:07 − random模块 random.randint(1,10)--随机生成0-10之间的随机整数 random.uniform(1,10)--随机生成0-10之间的实数 random.randrange(9,100,10)--从9-100之间随机选取一个实数,差为10,也就是说从9,19,29,39,49... 0 3530 相关推荐 2019-12-0…

C#生成含数字字母的随机字符串

C#生成含数字字母的随机字符串 要求生成的字符串是随机的&#xff0c;也就是字母和数字都需要随机&#xff0c;既可能只包含数字&#xff0c;也可能只包含字母&#xff0c;也可能两者都有。 实现方式如下: 首先定义一个包含所有字母和数字的字符串和一个空串&#xff0c;随后…

css 随机 数,纯CSS实现随机效果

最近在Codepen上看到了Adir写的随机翻牌和找蛋蛋(可以想象是砸金蛋)效果&#xff0c;让我再次刷新了对CSS的认知。看到这两个效果之后我才知道&#xff0c;在CSS中除了可以实现 来自其他语言的随机性 众所周知&#xff0c;CSS是一种声明式的标记语言。在很多同学的认知中&#…