医学大数据|R|竞争风险模型：基础、R操作与结果解读

久菜盒子工作室

已于 2024-02-27 17:24:14 修改

阅读量1.8k

点赞数 37

分类专栏： spss/stata/R/python与医学（大）数据分析文章标签： r语言人工智能开发语言

于 2024-02-26 23:38:26 首次发布

本文链接：https://blog.csdn.net/weixin_68126662/article/details/136303906

版权

spss/stata/R/python与医学（大）数据分析专栏收录该内容

16 篇文章 0 订阅

订阅专栏

责任编辑：医学大数据刘刘老师：头部医疗大数据公司医学科学部研究员

邮箱：897282268@qq.com

久菜盒子工作室我们是：985硕博/美国全奖doctor/计算机7年产品负责人/医学大数据公司医学研究员/SCI一区2篇/Nature子刊一篇/中文二区核心一篇/都是我们

主要领域：医学大数据分析/经管数据分析/金融模型/统计数理基础/统计学/卫生经济学/流行与统计学/

擅长软件：R/python/stata/spss/matlab/mySQL

团队理念：从零开始，让每一个人都得到优质的科研教育

代码整理与总结：

library("survival")
library("cmprsk")
library("mgus2")
data(mgus2)
#预处理
mgus2<-as.data.frame(mgus2)
data<-as.data.frame(mgus2)

mgus2$etime <- with(mgus2, ifelse(pstat==0, futime, ptime))  #当pstat==0时，etime=futime，否则etime=ptime 
#实际上这个地方，etime当发生竞争事件的时候，比如发生死亡，那么etime等于0-死亡时间
#当没有发生竞争事件的时候，etime等于0-跌落时间

mgus2$event <- with(mgus2, ifelse(pstat==0, 2*death, 1)) #当pstat==0时，event=2*death ，否则event=1
# 0 1 2 labels分别对应("censor", "pcm", "death")
event <- factor(event, 0:2)
单因素分析(cuminc)
cumic1<- cuminc(etime,event)


plot(cumic1,xlab = 'Month', ylab = '单因素的Fine-Gray检验',lwd=2,lty=1,
     col = c('red','blue'))
     
legend(0,800,
       c("PCM","death"))
print(cumic1)
cumic2<- cuminc(etime,event,mgus2$sex)
print(cumic2)
plot(cumic2,xlab = 'Month', ylab = '单因素的Fine-Gray检验',lwd=2,lty=1,
     col = c('red','blue','black','forestgreen'))


#多因素分析
dt<-na.omit(mgus2)
dt<-as.data.frame(dt)
cov <- data.frame(age = dt$age,
                  sex = ifelse(dt$sex=='M',1,0), ## 设置哑变量
                  hgb = dt$hgb) 
##构建多因素的竞争风险模型。此处需要指定failcode=1, cencode=0,
#分别代表结局事件赋值1与截尾赋值0，其他赋值默认为竞争风险事件2。

fit<- cuminc(dt$etime,dt$event,dt$sex)
summary(dt)
crr<-crr(dt$etime, dt$event, cov, failcode=1, cencode=0)
print(crr)

1.竞争风险模型基础

参考文献《基于竞争风险模型分析超高龄人群跌倒致长期卧床的危险因素及预测模型研究》

竞争风险模型，指的是在随访过程中，除了结局事件，还有其他事件，这些事件的发生，会“竞争”结局事件的发生，比如与结局事件无关的死亡，失访等，这些事件就称为竞争事件。当存在竞争事件时，传统的做法，通常是将竞争事件的发生视作删失(censor)，来计算结局事件的生存概率，但这会导致结局事件的生存概率估计有偏差，因此在存在竞争风险时，我们优先使用竞争风险模型，而不是cox等其他风险模型。

2.文献工作

2.1文献工作

于 2015 年 3—11 月选取上海市及杭州市 5 个区县每年定期参与体检的超高龄人群作为研究对象，建立前瞻性队列研究并基于中国健康与养老追踪调查（CHARLS）问卷收集研究对象相关信息，随访观察跌倒致长期卧床（目标事件）和死亡（竞争事件）。采用竞争风险模型分析超高龄人群跌倒致长期卧床的影响因素；将竞争风险模型筛选出的独立危险因素构建超高龄人群跌倒致长期卧床风险预测模型列线图并绘制受试者工作特征（ROC）曲线来评价模型的准确性。

运用R 4.2.2 软件统计描述和分析，符合正态分布的计量资料以（x±s）表示，计数资料以构成比表示，单因素分析组间比较采用 Fine-Gray 检验。采用竞争风险模型分析超高龄人群跌倒致长期卧床的影响因素，使用“cmprsk”、“rms”、“mstate” 等包构建列线图预测模型，并绘制受试者工作特征（ROC）曲线并计算 ROC 曲线下面积（AUC）来评价模型的准确性。检验水准α=0.05，以P<0.05 为差异有统计学意义。

2.2文献结果解读

其结果的解读阐述为：在有竞争事件的影响下

本研究共纳入986 名超高龄老年人，其中男431名（43.7%）、女555 名（56.3%），平均年龄为（89.8±5.2）岁。经过 8 年的随访，失访 96 人，失访率为 9.74%；发生目标事件 165 人，发生率为 16.73%；发生竞争事件 134 人，发生率为 13.59%。竞争风险模型分析结果显示，在有竞争事件的影响下，肌力的增加（HR=1.071，95%CI=1.049~1.091）、年龄 >85 岁（HR=1.954，95%CI=1.255~3.042）、居住地为农村（HR=1.946，95%CI=1.385~2.731）、睡眠质量较差（HR=5.756， 95%CI=3.904~8.491）、白内障（HR=1.832，95%CI=1.201~2.794）、糖尿病（HR=1.549，95%CI=1.121~2.143）、认知功能损害（HR=1.717，95%CI=1.258~2.344）为超高龄人群跌倒致长期卧床发生风险的独立危险因素（P<0.05）。超高龄人群跌倒致长期卧床风险预测模型的 ROC 曲线 AUC 为 0.798（95%CI=0.608~0.988），灵敏度为 0.841，特异度为 0.677。

赋值情况，需要注意的是，有三个结局，1个是0，一个是目标事件，一个是竞争事件（也就是死亡）

结果解读其实和cox回归是一样的

将上述单因素分析差异有统计学意义的11个变量纳入竞争风险模型中(变量赋值情况见表2)分析结果显示，在竞争事件的影响下，肌力增加(HR=1.071，95%C=1.049~1.091)、年龄>85 岁(HR=1.954，95%C=1.255~3.042)、居住地为农村(HR=1.946，95%CI=1.385~2.731)、睡眠质量较差(HR=5.756，95%CI=3.904~8.491)、白内障HR=1.832，95%C1=1.201~2.794)、糖尿病(HR=1.54995%CI=1.121~2.143)、认知功能损害(HR=1.717,95%CI=1.258~2.344)为超高龄人群跌倒致长期卧床发生风险的独立危险因素(P<0.05)。

3.R操作与结果解读

使用survival包（v3.1-12）的mgus2 dataset进行分析

hgb 血红蛋白；

creat 肌酐；

ptime 直至发展为浆细胞恶性肿瘤（PCM）或最后一次访视的时间（以月为单位）；

pstat 出现PCM(浆细胞恶性肿瘤)：0 =否，1 =是；

futime 直到死亡或最后一次接触的时间，以月为单位；

death 发生死亡：0 =否，1 =是；这里把PCM作为结局事件，death作为竞争事件。

#加载包

library("survival")

library("cmprsk")

Log

> library("survival")

Warning message:

程辑包‘survival’是用R版本4.3.2 来建造的

> library("cmprsk")

Warning message:

程辑包‘cmprsk’是用R版本4.3.2 来建造的

导入数据

#预处理

mgus2<-as.data.frame(mgus2)

data<-as.data.frame(mgus2)

单因素分析

3.1单因素分析(cuminc)

Cumic1<- cuminc(etime,event)

plot(cumic1,xlab = 'Month', ylab = '单因素的Fine-Gray检验',lwd=2,lty=1,

col = c('red','blue'))

这个图要画的更好看一点

print(cumic1)

est是对应估计出来的子分布累计风险概率，var是各个估计子分布累计风险概率的估计偏差。

加入分组单因素，比如性别

cumic2<- cuminc(etime,event,mgus2$sex)

Tests:

       stat          pv        df

1 1.194508     0.274422157    1

2 11.651259    0.000641591    1

Estimates and Variances:

$est

           100        200       300       400

F 1 0.06196166 0.09535198 0.1151287        NA

M 1 0.04865458 0.08914153 0.1044602 0.1044602

F 2 0.41670174 0.65422980 0.7218621        NA

M 2 0.52056642 0.70771036 0.7994364 0.7994364

$var

                  100          200          300          400

F 1 9.518160e-05    0.0001611587 0.0003038417           NA

M 1 6.479291e-05   0.0001455188 0.0002579428 0.0002579428

F 2 4.039584e-04    0.0004893207 0.0006741219           NA

M 2 3.527969e-04   0.0003697058 0.0005904656 0.0005904656

Test对应的是性别变量分别对于结局1，2的影响是否显著，pv即P-value。可见，性别对于death影响是显著的，对于PCM是不显著的。

plot(cumic2,xlab = 'Month', ylab = '单因素的Fine-Gray检验',lwd=2,lty=1,

     col = c('red','blue','black','forestgreen'))

纵坐标表示累计发生率pcm，横坐标是时间轴。我们从F1对应的红色曲线和M1对应的蓝色曲线可以得出，女性组的发生pcm风险较男性组高，但未达到统计学意义，P=0.274422157。同理，F2对应的黑色曲线在M2对应的草绿色曲线下方，我们可以得出，女性组的竞争风险事件发生率较男性组低，具有统计学差异，P=0.0064。这个图可以用一句话来概括：在控制了竞争风险事件后，女性组和男性组累计PCM发生率无统计学差异P=0.274422157。

3.2多因素分析（crr）

在cmprsk包中crr()函数可以很方便的实现多因素分析。

参考链接：https://mp.weixin.qq.com/s/IKsBZMjDzz4RiUrn7RuV_g

###多因素竞争风险

cov <- data.frame(age = dt$age,

                  sex = ifelse(dt$sex=='M',1,0), ## 设置哑变量

                  hgb = dt$hgb)

covv <- mgus2[, c('age', 'hgb', 'creat')]

设置哑变量其实还有很多方法，我们今天暂时用这种方法

crr<-crr(dt$etime, dt$event, cov, failcode=1, cencode=0)  

#第一次试分析，有删除的数据会报错

> fit<- cuminc(dt$etime,dt$event,dt$sex)

Error in data.frame(time = ftime, cause = fstatus, group = as.factor(if (missing(group)) rep(1, :

  参数值意味着不同的行数: 0, 1384

后面发现是因为dt数据集中不包含，event和etime，把这两个放到同一个数据集里。

> mgus2$etime <- with(mgus2, ifelse(pstat==0, futime, ptime)) #当pstat==0时，etime=futime，否则etime=ptime

> mgus2$event <- with(mgus2, ifelse(pstat==0, 2*death, 1)) #当pstat==0时，event=2*death ，否则event=1

> #多因素分析

> dt<-na.omit(mgus2)

> cov <- data.frame(age = dt$age,

+                   sex = ifelse(dt$sex=='M',1,0), ## 设置哑变量

+                   hgb = dt$hgb)

> fit<- cuminc(dt$etime,dt$event,dt$sex)

> crr<-crr(dt$etime, dt$event, cov, failcode=1, cencode=0)

print(crr)

> print(crr)

convergence: TRUE

coefficients:

     age         sex        hgb

-0.01911         -0.27880    -0.01408

standard errors:

[1] 0.005898      0.190200   0.046240

two-sided p-values:

   age         sex      hgb

0.0012         0.1400   0.7600

coefficients为系数，standard errors为标准误差，two-sided p-values为双边P值。

代码整理与总结

library("survival")
library("cmprsk")
library("mgus2")
data(mgus2)
#预处理
mgus2<-as.data.frame(mgus2)
data<-as.data.frame(mgus2)

mgus2$etime <- with(mgus2, ifelse(pstat==0, futime, ptime))  #当pstat==0时，etime=futime，否则etime=ptime 
#实际上这个地方，etime当发生竞争事件的时候，比如发生死亡，那么etime等于0-死亡时间
#当没有发生竞争事件的时候，etime等于0-跌落时间

mgus2$event <- with(mgus2, ifelse(pstat==0, 2*death, 1)) #当pstat==0时，event=2*death ，否则event=1
# 0 1 2 labels分别对应("censor", "pcm", "death")
event <- factor(event, 0:2)
单因素分析(cuminc)
cumic1<- cuminc(etime,event)


plot(cumic1,xlab = 'Month', ylab = '单因素的Fine-Gray检验',lwd=2,lty=1,
     col = c('red','blue'))
     
legend(0,800,
       c("PCM","death"))
print(cumic1)
cumic2<- cuminc(etime,event,mgus2$sex)
print(cumic2)
plot(cumic2,xlab = 'Month', ylab = '单因素的Fine-Gray检验',lwd=2,lty=1,
     col = c('red','blue','black','forestgreen'))


#多因素分析
dt<-na.omit(mgus2)
dt<-as.data.frame(dt)
cov <- data.frame(age = dt$age,
                  sex = ifelse(dt$sex=='M',1,0), ## 设置哑变量
                  hgb = dt$hgb) 
##构建多因素的竞争风险模型。此处需要指定failcode=1, cencode=0,
#分别代表结局事件赋值1与截尾赋值0，其他赋值默认为竞争风险事件2。

fit<- cuminc(dt$etime,dt$event,dt$sex)
summary(dt)
crr<-crr(dt$etime, dt$event, cov, failcode=1, cencode=0)
print(crr)