admin 管理员组文章数量: 887021
Fine
竞争风险模型(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()函数进行竞争风险分析，不然最后画列线图会报错。
# 参考资料
#
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)
在这个列线图中，将数据集中id=25的患者各协变量的取值映射到相应的得分，并计算总得分,并分别计算其在36个月和60个月的累计复发概率，此概率即为控制了竞争风险的累计复发概率，分别为：0.134和0.146。
参考资料
本文标签： Fine
版权声明：本文标题：Fine 内容由网友自发贡献，该文观点仅代表作者本人， 转载请联系作者并注明出处：http://www.freenas.com.cn/jishu/1686703025h26976.html， 本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌抄袭侵权/违法违规的内容，一经查实，本站将立刻删除。
发表评论