R语言生存分析：Cox回归

最新推荐文章于 2024-10-31 23:06:23 发布
最新推荐文章于 2024-10-31 23:06:23 发布
阅读量1.1w 收藏 171
点赞数 18
分类专栏： 医学统计学 文章标签： r语言 回归 开发语言
原文链接：https://mp.weixin.qq.com/s/eeQ8PJQZunYzQmssGTimNA
版权

上次介绍了生存分析中的寿命表、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)

 

可以看到sex和age的回归系数随着时间变化基本没啥变化，稳定在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

医学和生信笔记
关注
专栏目录
R语言 所有分类数据回归
08-30
#glmnet包含有线性回归，逻辑回归，泊松计数模型，cox回归模型，多分类逻辑回归多响应线性回归 #阿法系数=0是岭回归，阿法系数=1，是lasso回归 ############################################################ data(QuickStartExample) #20 Independent variables and 1 dependent variable fit<-glmnet(x,y) #查看计算路径 print(fit)#Df是非零系数的个数，第三列是惩罚力度，第二列是相应的惩罚力度对应的解释偏差百分比 #利用coef来提取模型的系数 coef(fit,s=0.1) #可视化 plot(fit) #预测 #参数newx用来设置输入数据，s用来设置朗姆塔（惩罚力度）的值 nx<-matrix(rnorm(10*20),10,20) predict(fit,newx=nx,s=c(0.1,0.05)) #交叉验证，筛选最优的朗姆塔的值，cv.glmnet函数实现这个功能 cvfit<-cv.glmnet(x,y) plot(cvfit)#红色散点为交叉验证的散点图，两条垂直的虚线为交叉验证后选定的朗姆塔 #提取最优的朗姆塔的值 cvfit$lambda.min#均方误差最小 cvfit$lambda.1se#一倍标准差 #根据朗姆塔来取得系数 coef(cvfit,s='lambda.min') #将稀疏矩阵转化为传统矩阵 as.matrix(coef(cvfit,s='lambda.min'))
Cox回归 R语言
qq_36481674的博客
10-29 2584
R数据读取，除重，进行合并，进行Cox回归
如何用R语言做cox先单后多回归分析！（附全套代码）
qq_37610365的博客
08-13 1069
根据文献中的纳入排除标准，提取涉及的相关数据，最终共纳入1，574名患者（原文献1，641）。案例文献是沈阳医学院公共卫生学院学者基于SEER数据库的一项回顾性研究，旨在建立一个列线图来预测老年恶性骨肿瘤（MBT）患者的总生存期（OS）。我们建立了一个新的列线图来预测老年MBT患者的1年，3年，5年的OS。在影响因素研究中，大家较为熟知的筛选方式是先单后多，但是在预测模型文章中，更为常用的筛选方法是逐步回归法。R包，请注意模型变量的筛选是根据训练集进行的，在进行数据拆分后，需要用训练集来进行Cox回归。
Cox风险回归分析
robustness博客
04-24 9340
（event = 1,指的是确实事件；0是死亡，1是缺失 1.用的是ovarian自带数据 ,首先画出生存曲线survfit(),但是曲线上不显示确实值的符号+号，怎么在函数中设置合适的缺失值？？？ 2接下来对rx进行log-rank检验 p取的是0.05，所以拒绝原假设，两组实验没有显著性区别， 若是有显著性，可以通过中位生存期，或是RR值比较两种情况对生存状况的
r语言回归分析_R语言之cox回归分析
weixin_39626237的博客
12-08 1万+
Cox比例风险模型(Cox proportional-hazards model，也称为Cox回归)，主要用于带有时间的生存结局的影响因素研究，或评价某个临床治疗措施对患者生存的影响。 Cox模型可以由hazard function表示，h(t)；简单的说就是t时刻死亡的风险，公式如下：h(t)=h0(t) × exp(b1x1 + b2x2 +…+ bpxp)t...
R语言：多因素Cox回归森林图 (基于forestplot包) 森林图 cox可视化
生信小博士的博客
10-15 1万+
计划是介绍一下单因素、多因素、亚组分析、其他类型森林图绘制；列线图部分介绍下基本的变量筛选方法、绘制过程、验证方法（ROC/DCA/校准图/分组比较等）、网页版列线图。其中，1.2是对表的修改，改完表，图就改了；3-5步是通过改字体、线宽等优化的。三-4、字体大小增粗、线条粗细、线型（包括置信区间）、行间距、列间距；三-3、加水平线、垂直辅助线、x周标签、大标题、森林图占比；三-5、颜色、HR方块形状等。三-2、显示所有亚组的患者数。三-1、显示所有变量；三-2、加入患者数目；三-1、显示所有变量。
R进行各种变量下的cox回归分析，重要的是还将输出结果以表格形式转出
宽嘴鱼汤的博客
02-05 3385
关于单变量cox回归、批量单变量cox回归、有序变量下的cox回归、将cox回归结果以表格形式输出的系列内容
（生物信息学）R语言与统计学入门（十）—— 多因素Cox回归分析
weixin_46500027的博客
05-07 1万+
上次提到单因素cox回归分析： （生物信息学）R语言与统计学入门（九）—— 单因素cox回归分析_Lijingxian教你学生信的博客-CSDN博客_r语言单因素回归分析​COX回归模型，又称“比例风险回归模型(proportional hazards model，简称Cox模型)”，是由英国统计学家D.R.Cox(1972)年提出的一种半参数回归模型。该模型以生存结局和生存时间为因变量，可同时分析众多因素对生存期的影响，能分析带有截尾生存时间的资料，且不要求估计资料的生存分布类型。由于上述优良性质，该模
R语言中的生存分析：从理论到实践
10-25
通过本文的介绍，读者应该能够理解生存分析的基本概念，并能够在R语言中实现基本的生存分析流程，包括创建生存对象、估计生存曲线、使用Cox比例风险模型进行回归分析，以及进行模型诊断。 通过不断的实践和学习，...
Cox生存分析代码
09-02
Cox生存分析，分析疾病的风险因素，得出疾病的生存曲线。
利用R对癌症数据进行生存分析
07-22
本教程将详细介绍如何利用R语言对癌症数据集进行生存分析。 生存分析的核心概念是生存时间（Survival Time），通常表示为从研究开始到某个事件发生的时间，比如癌症患者的存活时间。在R中，`survival`包是进行生存...
(52条消息) R语言生存分析COX回归分析实战：以乳腺癌数据为例_Data+Science+Insight的博客-CSDN博客_利用r对癌症数据进行生存分析.mhtml
01-12
(52条消息) R语言生存分析COX回归分析实战：以乳腺癌数据为例_Data+Science+Insight的博客-CSDN博客_利用r对癌症数据进行生存分析.mhtml
R语言在统计中的应用课件
03-01
- 生存分析：寿命表法、Cox比例风险模型等生存数据分析方法。 - 贝叶斯统计：贝叶斯定理及其在参数估计、假设检验中的应用。 ### 四、特色亮点 - **实用性强**：所有内容均围绕实际案例展开，注重解决实际问题的...
R语言进行COX时变系数模型（含时间依存协变量的Cox回归模型）
dege857的博客
09-13 1万+
0 = 无症状，1 = 有症状但完全不卧床，2 = 卧床时间 < 50%，3 = 卧床时间 > 50% 但不卧床，4 = 卧床，ph.karno: 由医师评定的 Karnofsky 表现评分（差=0-好=100），pat.karno: 由患者评定的 Karnofsky 性能评分，meal.cal: 用餐时消耗的卡路里，wt.loss: 过去六个月的体重减轻（磅）第一种我们可以看到时间t和X有关系，对系数β没什么影响，第二种可以看到时间t对β有影响，随着时间变化β是不同的，即系数是不同的。
用R语言进行Cox回归生存分析
热门推荐
庐州月光的博客
06-29 3万+
欢迎关注”生信修炼手册”!在生存分析中，探究生存时间的影响因素是一个重要的研究内容，通过KM和log-rank test检验的方法，只能够处理单个二分类因素的生存数据。当想探究多个因素或...
r语言 svycoxph_生存分析的Cox回归模型(比例风险模型)R语言实现及结果解读
weixin_39927993的博客
12-21 5018
R语言中的Cox模型分析单变量Cox回归library("survival")library("survminer")res.cox res.cox summary(res.cox)summary的结果：Call:coxph(formula = Surv(time, status) ~ sex, data = lung)n= 228, number of events= 165coef exp(c...
AIGC 浪潮下，R 语言奏响数据分析新乐
ya766829的博客
10-31 332
在当今的 AIGC（人工智能生成内容）时代，数据如同璀璨的星辰，散落在各个领域的浩瀚宇宙中。而在这场数据的盛宴中，R 语言以其独特的魅力和强大的功能，引领着数据分析的新风尚。一、AIGC 时代的数据洪流随着人工智能技术的飞速发展，AIGC 如同一股汹涌的浪潮，席卷了各个行业。从自然语言处理到图像识别，从金融领域到医疗健康，大量的数据被不断地生成、收集和存储。这些数据蕴含着丰富的信息和价值，等待着被挖掘和利用。在这个时代，数据的规模呈指数级增长。传统的数据处理方法已经难以满足需求，迫切需要更高效、更强大的数据
[java][高级]Request&Response
最新发布
一名不太专业的程序员
10-31 738
=Request是请求对象，Response是响应对象。==这两个对象在我们使用Servlet的时候有看到：此时，我们就需要思考一个问题request和response这两个参数的作用是什么?request:==获取==请求数据浏览器会发送HTTP请求到后台服务器[Tomcat]HTTP的请求中会包含很多请求数据[请求行+请求头+请求体]后台服务器[Tomcat]会对HTTP请求中的数据进行解析并把解析结果存入到一个对象中。
CST软件如何理解Axial Ratio轴比
a18665053298的博客
10-31 645
所以：AR= Sqrt[(|E1|^2+|E2| ^2+|E1| ^2-|E2| ^2)/(|E1| ^2+|E2| ^2-|E1| ^2+|E2| ^2)]所以： |EL|=1/sqrt(2)*sqrt[(|E1|+|E2|)^2+(|E1|+|E2|)^2]= |E1|+|E2|同理：|ER|=1/sqrt(2)*sqrt[(|E1|-|E2|)^2+(|E1|-|E2|)^2] = |E1|-|E2|答案当然是可以等效。原因很简单，由定义式可知，AR是远场最基础的两个正交电场分量E1和E2算出来的。
schoenfeld 残差图
01-12
Schoenfeld残差图是生存分析中用来检验Cox比例风险模型中的共线性假设的一种工具。该图用于评估模型中的自变量是否与风险残差存在系统性的关联。 Schoenfeld残差是指在Cox模型中计算得到的观测值的风险残差，它们是模型中的每个个体所产生的有序列的残差。而Schoenfeld残差图则是以时间为横轴，Schoenfeld残差为纵轴进行绘制的。 在Schoenfeld残差图中，如果模型中的自变量与风险残差呈现明显的非零关系，那么可以说明模型存在共线性问题。如果残差连续分布在水平线附近，则表示模型中的自变量对残差没有系统性的影响。 Schoenfeld残差图通常会绘制每个自变量的残差与时间的关系，并通过平滑的曲线来观察残差的变化趋势。如果曲线呈现水平分布，即残差与时间无关，则证明了共线性假设的合理性。而如果曲线呈现明显的趋势或波动，那么可以怀疑模型中的自变量存在共线性问题。 通过分析Schoenfeld残差图，研究者可以检查模型中自变量的线性关系及时间变化性质，以提供决策者选择最佳模型的依据。此外，Schoenfeld残差图还可以用于发现风险比率随时间变化的不规则情况，从而更全面地解释数据。 总之，Schoenfeld残差图在生存分析中具有重要的意义，可以帮助评估模型中的共线性假设，同时也为选择和优化模型提供了可靠的参考。
评论 1
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值