Fine-Gray检验、竞争风险模型、列线图绘制

关注

转载最新推荐文章于 2026-08-30 16:22:44 发布 5.8k 阅读 AI 写同款· GEO 优化 ›

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

竞争风险模型(Competing Risk Model)适用于多个终点的生存数据，传统的生存分析(survival analysis) 一般只关心一个终点事件(即研究者感兴趣的结局)。将其他事件均按删失数据(Censored Data)处理，要求个体删失情况与个体终点事件相互独立，结局不存在竞争风险。

竞争风险模型(Competing Risk Model) : 指的是在观察队列中，存在某种已知事件可能会影响另一种事件发生的概率或者是完全阻碍其发生,则可认为前者与后者存在竞争风险。

比如我们研究某药物和肠癌复发的关系，感兴趣事件是肠癌复发，但是研究过程中病人因为心梗死了，这样就观察不到感兴趣事件了，那心梗死亡就可以被叫做竞争风险事件。

本文目录：

加载数据和R包
Fine-Gray检验（单因素分析）
- 图形展示结果
- ggplot2
竞争风险模型（多因素分析）
列线图
参考资料

加载数据和R包

探讨骨髓移植和血液移植治疗白血病的疗效，结局事件定义为复发，某些患者因为移植不良反应死亡，定义为竞争风险事件。

rm(list = ls())
data("bmtcrr",package = "casebase")
str(bmtcrr)
## 'data.frame': 177 obs. of  7 variables:
##  $ Sex   : Factor w/ 2 levels "F","M": 2 1 2 1 1 2 2 1 2 1 ...
##  $ D     : Factor w/ 2 levels "ALL","AML": 1 2 1 1 1 1 1 1 1 1 ...
##  $ Phase : Factor w/ 4 levels "CR1","CR2","CR3",..: 4 2 3 2 2 4 1 1 1 4 ...
##  $ Age   : int  48 23 7 26 36 17 7 17 26 8 ...
##  $ Status: int  2 1 0 2 2 2 0 2 0 1 ...
##  $ Source: Factor w/ 2 levels "BM+PB","PB": 1 1 1 1 1 1 1 1 1 1 ...
##  $ ftime : num  0.67 9.5 131.77 24.03 1.47 ...

这个数据一共7个变量，177行。

Sex: 性别，F是女，M是男
D: 疾病类型，ALL是急性淋巴细胞白血病，AML是急性髓系细胞白血病。
Phase: 不同阶段，4个水平，CR1，CR2，CR3，Relapse。
Age: 年龄。
Status: 结局变量，0=删失，1=复发，2=竞争风险事件。
Source: 因子变量，2个水平：BM+PB(骨髓移植+血液移植)，PB(血液移植)。
ftime: 生存时间。

# 竞争风险分析需要用的R包
library(cmprsk)
## Loading required package: survival

Fine-Gray检验（单因素分析）

在普通的生存分析中，可以用log-rank检验做单因素分析，在竞争风险模型中，使用Fine-Gray检验进行单因素分析。

比如现在我们想要比较不同疾病类型（D）有没有差异，可以进行Fine-Gray检验：

bmtcrr$Status <- factor(bmtcrr$Status)
f <- cuminc(bmtcrr$ftime, bmtcrr$Status, bmtcrr$D)
f
## Tests:
##        stat         pv df
## 1 2.8623325 0.09067592  1
## 2 0.4481279 0.50322531  1
## Estimates and Variances:
## $est
##              20        40        60        80       100       120
## ALL 1 0.3713851 0.3875571 0.3875571 0.3875571 0.3875571 0.3875571
## AML 1 0.2414530 0.2663827 0.2810390 0.2810390 0.2810390        NA
## ALL 2 0.3698630 0.3860350 0.3860350 0.3860350 0.3860350 0.3860350
## AML 2 0.4439103 0.4551473 0.4551473 0.4551473 0.4551473        NA
## 
## $var
##                20          40          60          80         100         120
## ALL 1 0.003307032 0.003405375 0.003405375 0.003405375 0.003405375 0.003405375
## AML 1 0.001801156 0.001995487 0.002130835 0.002130835 0.002130835          NA
## ALL 2 0.003268852 0.003373130 0.003373130 0.003373130 0.003373130 0.003373130
## AML 2 0.002430406 0.002460425 0.002460425 0.002460425 0.002460425          NA

结果中1代表复发,2代表竞争风险事件。

第一行统计量=2.8623325, P=0.09067592,表示在控制了竞争风险事件（即第二行计算的统计量和P值）后，两种疾病类型ALL和AML的累计复发风险无统计学差异P=0.09067592。

第2行说明ALL和AML的累计竞争风险无统计学差异。

$est表示估计的各时间点ALL和AML组的累计复发率与与累计竞争风险事件发生率（分别用1和2来区分，与第一行第二行一致）。

$var表示估计的各时间点ALL和AML组的累计复发率与与累计竞争风险事件发生率的方差（分别用1和2来区分，与第一行第二行一致）。

图形展示结果

对于上述结果可以使用图形展示：

plot(f,xlab = 'Month', ylab = 'CIF',lwd=2,lty=1,
     col = c('red','blue','black','forestgreen'))

unnamed-chunk-4-160276184

图形解读：

纵坐标表示累计发生率CIF，横坐标是时间。我们从ALL1对应的红色曲线和AML1对应的蓝色曲线可以得出，ALL组的复发风险较AML 组高，但无统计学意义，P=0.09067592。同理，ALL2对应的黑色曲线在AML2对应的草绿色曲线下方，我们可以得出，ALL组的竞争风险事件发生率较AML组低，同样无统计学意义，P=0.50322531。

简单来讲，这个图可以用一句话来概括：在控制了竞争风险事件后，ALL和AML累计复发风险无统计学差异P=0.09067592。

ggplot2

这个图不好看，非常的不ggplot，所以我们要用ggplot2重新画它！所以首先要提取数据，因为数就是图，图就是数。但是万能的broom包竟然没有不能提取这个对象的数据，只能手动来，太不优雅了！

# 提取数据
ALL1 <- data.frame(ALL1_t = f[[1]][[1]], ALL1_C = f[[1]][[2]])
AML1 <- data.frame(AML1_t = f[[2]][[1]], AML1_C = f[[2]][[2]])
ALL2 <- data.frame(ALL2_t = f[[3]][[1]], ALL2_C = f[[3]][[2]])
AML2 <- data.frame(AML2_t = f[[4]][[1]], AML2_C = f[[4]][[2]])

library(ggplot2)

ggplot()+
  geom_line(data = ALL1, aes(ALL1_t,ALL1_C))+
  geom_line(data = ALL2, aes(ALL2_t,ALL2_C))+
  geom_line(data = AML1, aes(AML1_t,AML1_C))+
  geom_line(data = AML2, aes(AML2_t,AML2_C))+
  labs(x="month",y="cif")+
  theme_bw()

unnamed-chunk-5-160276184

但是这种不好上色，所以我们美化一下，变成长数据再画图即可。

tmp <- data.frame(month = c(ALL1$ALL1_t,AML1$AML1_t,ALL2$ALL2_t,AML2$AML2_t),
                  cif = c(ALL1$ALL1_C,AML1$AML1_C,ALL2$ALL2_C,AML2$AML2_C),
                  type = rep(c("ALL1","AML1","ALL2","AML2"), c(58,58,58,88))
                  )

ggplot(tmp, aes(month, cif))+
  geom_line(aes(color=type, group=type),size=1.2)+
  theme_bw()+
  theme(legend.position = "top")

unnamed-chunk-6-160276184

但是也不太好看，大家自己修改即可。

竞争风险模型（多因素分析）

做完了单因素分析，再看看竞争风险模型的多因素分析。

首先要把自变量单独放在一个数据框里，使用中发现一个问题，这里如果把分类变量变为因子型不会自动进行哑变量编码，所以需要手动进行哑变量编码！

但是我这里偷懒了，并没有进行哑变量设置！实际中是需要的哦！！

covs <- subset(bmtcrr, select = - c(ftime,Status))
covs[,c(1:3,5)] <- lapply(covs[,c(1:3,5)],as.integer)

str(covs)
## 'data.frame': 177 obs. of  5 variables:
##  $ Sex   : int  2 1 2 1 1 2 2 1 2 1 ...
##  $ D     : int  1 2 1 1 1 1 1 1 1 1 ...
##  $ Phase : int  4 2 3 2 2 4 1 1 1 4 ...
##  $ Age   : int  48 23 7 26 36 17 7 17 26 8 ...
##  $ Source: int  1 1 1 1 1 1 1 1 1 1 ...

指定failcode=1, cencode=0, 分别代表结局事件1与截尾0，其他默认为竞争风险事件2。

# 构建竞争风险模型
f2 <- crr(bmtcrr$ftime, bmtcrr$Status, covs, failcode=1, cencode=0)
summary(f2)
## Competing Risks Regression
## 
## Call:
## crr(ftime = bmtcrr$ftime, fstatus = bmtcrr$Status, cov1 = covs, 
##     failcode = 1, cencode = 0)
## 
##           coef exp(coef) se(coef)      z p-value
## Sex     0.0494     1.051   0.2867  0.172 0.86000
## D      -0.4860     0.615   0.3040 -1.599 0.11000
## Phase   0.4144     1.514   0.1194  3.470 0.00052
## Age    -0.0174     0.983   0.0118 -1.465 0.14000
## Source  0.9526     2.592   0.5469  1.742 0.08200
## 
##        exp(coef) exp(-coef)  2.5% 97.5%
## Sex        1.051      0.952 0.599  1.84
## D          0.615      1.626 0.339  1.12
## Phase      1.514      0.661 1.198  1.91
## Age        0.983      1.018 0.960  1.01
## Source     2.592      0.386 0.888  7.57
## 
## Num. cases = 177
## Pseudo Log-likelihood = -267 
## Pseudo likelihood ratio test = 23.6  on 5 df,

结果解读：在控制了竞争分险事件后，phase变量，即疾病所处阶段是患者复发的独立影响因素(p =0.00052)。

列线图

regplot包绘制列线图。但是它目前只适用coxph()、lm()和glm()返回的对象。

因此我们需要对原数据集加权创建一个新数据集用于为竞争风险模型分析，使用mstate包中的crprep()创建加权数据集,然后使用coxph()对加权数据集进行竞争风险模型拟合，这样就可以画列线图了。

首先是加载数据和R包：

rm(list = ls())
data("bmtcrr",package = "casebase") # 还是这个数据

library(mstate) # 加权用到的R包

bmtcrr$id <- 1:nrow(bmtcrr) # 创建id

# phase变为2分类，不然列线图不好解释
bmtcrr$Phase <- factor(ifelse(bmtcrr$Phase=="Relapse",1,0)) 
str(bmtcrr)
## 'data.frame': 177 obs. of  8 variables:
##  $ Sex   : Factor w/ 2 levels "F","M": 2 1 2 1 1 2 2 1 2 1 ...
##  $ D     : Factor w/ 2 levels "ALL","AML": 1 2 1 1 1 1 1 1 1 1 ...
##  $ Phase : Factor w/ 2 levels "0","1": 2 1 1 1 1 2 1 1 1 2 ...
##  $ Age   : int  48 23 7 26 36 17 7 17 26 8 ...
##  $ Status: int  2 1 0 2 2 2 0 2 0 1 ...
##  $ Source: Factor w/ 2 levels "BM+PB","PB": 1 1 1 1 1 1 1 1 1 1 ...
##  $ ftime : num  0.67 9.5 131.77 24.03 1.47 ...
##  $ id    : int  1 2 3 4 5 6 7 8 9 10 ...

然后是对原数据进行加权：

df.w <- crprep("ftime", "Status",
               data=bmtcrr, 
               trans=c(1,2),# 要加权的变量，1表示结局事件，2表示竞争风险事件
               cens=0, # 删失
               id="id",
               
               # 要保留的协变量
               keep=c("Age","Sex","D","Source","Phase"))

df.w
##       id Tstart  Tstop status weight.cens Age Sex   D Source Phase count
## 1      1   0.00   0.67      2   1.0000000  48   M ALL  BM+PB     1     1
## 2      1   0.67   9.50      2   1.0000000  48   M ALL  BM+PB     1     2
## 3      1   9.50  13.07      2   0.9679938  48   M ALL  BM+PB     1     3
## 4      1  13.07  17.23      2   0.8730924  48   M ALL  BM+PB     1     4
## 5      1  17.23  20.83      2   0.8536904  48   M ALL  BM+PB     1     5
## 6      1  20.83  28.53      2   0.8120469  48   M ALL  BM+PB     1     6
## 省略......

df.w$T<- df.w$Tstop - df.w$Tstart

上述代码已经创建一个加权数据集df.w，此时还需要选择failcode == 1的行，然后我们才可以在此数据集上使用coxph()函数进行竞争风险分析，不然最后画列线图会报错。

# 参考资料
# https://blog.csdn.net/zhongkeyuanchongqing/article/details/124086113
df.w2 <- df.w[df.w$failcode == 1,]

构建cox模型：

m.crr<- coxph(Surv(T,status==1)~Age+Sex+D+Source+Phase,
             data=df.w2,
             weight=weight.cens,
             subset=failcode==1)
summary(m.crr)
## Call:
## coxph(formula = Surv(T, status == 1) ~ Age + Sex + D + Source + 
##     Phase, data = df.w2, weights = weight.cens, subset = failcode == 
##     1)
## 
##   n= 686, number of events= 56 
## 
##              coef exp(coef) se(coef) robust se      z Pr(>|z|)    
## Age      -0.02174   0.97850  0.01172   0.01208 -1.800 0.071914 .  
## SexM      0.10551   1.11128  0.27981   0.29571  0.357 0.721247    
## DAML     -0.53163   0.58764  0.29917   0.30613 -1.737 0.082450 .  
## SourcePB  1.06564   2.90269  0.53453   0.56000  1.903 0.057051 .  
## Phase1    1.06140   2.89040  0.27870   0.28129  3.773 0.000161 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##          exp(coef) exp(-coef) lower .95 upper .95
## Age         0.9785     1.0220    0.9556     1.002
## SexM        1.1113     0.8999    0.6225     1.984
## DAML        0.5876     1.7017    0.3225     1.071
## SourcePB    2.9027     0.3445    0.9686     8.699
## Phase1      2.8904     0.3460    1.6654     5.016
## 
## Concordance= 0.737  (se = 0.037 )
## Likelihood ratio test= 28.33  on 5 df,   p=3e-05
## Wald test            = 27.27  on 5 df,   p=5e-05
## Score (logrank) test = 30.49  on 5 df,   p=1e-05,   Robust = 20.2  p=0.001
## 
##   (Note: the likelihood ratio and score tests assume independence of
##      observations within a cluster, the Wald and robust score tests do not).

接下来，我们可以使用regplot()函数绘制nomogram。其实你可以绘制多种不同的列线图，可以参考之前的推文：生存资料列线图的4种绘制方法

library(regplot)
regplot(m.crr,
        observation=df.w2[df.w2$id==25&df.w2$failcode==1,],
        failtime = c(36, 60), 
        prfail = T, 
        droplines=T)