4.临床预测模型——竞争风险模型

基本概念

竞争风险事件：指出现研究对象感兴趣事件的同时，其他终点事件也有可能出现，这些终点事件将阻止感兴趣事件的出现，或使其发生概率降低，各终点事件之间形成所谓的竞争关系。

竞争风险模型仅仅关心研究对象发生的第一个终点事件，而后发生的其他终点事件称为删失事件（censoring）。

如基线未发生心血管疾病的研究对象在观察期内死于癌症、车祸等其他原因前并未发生心血管疾病，就不能为CVD的发病做出贡献，传统生存分析将其他原因死亡的个体，失访个体和存活个体记为删失数据，会高估CVD的累积发病率。

累积发生函数（CIF）：

CIF（t）表示在时间t及其他事件之前第k类事件的概率。

在竞争风险中，结局不再仅仅是生存或者死亡，此时CIF（t）不等于F（t）。

1.Gray检验：

Gray检验属于单因素分析方法，可以用来检验存在竞争风险事件情况下，两组或两组以上感兴趣事件的累积发生率是否存在统计学差异。

2. 风险函数回归：

多因素分析需要使用风险函数回归。原因别风险函数（CS），部分分布风险函数（SD）后者又称为CIF回归模型、Fine-Gray模型。CS适合病因学研究，回归系数反映协变量对尚未发生任何事件人群中主要终点事件发生率增加的相对作用，SD适合风险预测研究，建立临床预测模型及风险评分，仅对终点事件绝对发生率感兴趣。

案例1：理解竞争风险模型

某研究人员收集了本市2007年确诊为轻度认知损害（MCI）的518例老年患者资料，包括基本人口学特征、生活方式、体格检查和合并疾病信息等，并于2010-2013年完成6次随访调查，主要观察结局为发生阿尔兹海默病（AD）。

随访期间，共发生AD 78例，失访84例，其中28例搬迁、31例退出、25例死亡。试问影响MCI向AD转归的因素都有哪些？

思考：生存分析是预后研究中比较常见的统计分析方法，但是经典的生存分析一般只关心一个终点事件（即研究者感兴趣的结局），而医学研究中观察的终点往往并不唯一（即出现不感兴趣的结局）。比如MCI患者在观察期间死于癌症、心血管疾病、车祸等原因而未发生AD，就不能为AD的发病做出贡献，即死亡“竞争”了AD的发生。传统统计方法将发生AD前死亡的个体、失访个体和未发生AD个体均按删失数据（censored data）处理，可能会导致估计偏差。

——本例中可以将发生AD前死亡作为AD的竞争风险事件，采用竞争风险模型进行分析。竞争风险的单因素分析常用来估计关心终点事件的发生率，多因素分析常用来探索预后影响因素及效应值。

对于死亡率较高的老年人群，当有竞争风险事件存在时，采用传统生存分析方法（K-M法、Cox比例风险回归模型）会高估所研究疾病的发生风险，产生竞争风险偏倚，有人专门研究发现约46%的文献可能存在这种偏倚。

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

竞争风险模型需要使用cmprsk包中的cuminc（）函数

案例：研究骨髓移植相较血液移植治疗白血病的疗效，结局事件定义为“复发”，某些患者移植后不幸因为移植的不良反应死亡，那对于这些发生移植相关死亡的患者，就无法观察到其复发的终点。“移植相关死亡”和“复发”之间存在竞争关系。

bmt <-read.csv('bmtcrr.csv',stringsAsFactors = TRUE)
install.packages("cmprsk")
library(cmprsk) 
str(bmt)

stringsAsFactors = TRUE表示若数据集中存在字符型变量，则将其处理为因子型变量。

数据框中性别sex，2个水平，F表示女性，M表示男性

D表示疾病类型，2个水平，ALL表示急性淋巴细胞白血病，AML表示急性髓系细胞白血病

Age表示年龄

Phrase表示疾病阶段，CR1表示初始发病治疗后缓解，CR2表示复发后再次缓解，CR3表示子啊此复发后缓解，Relapse表示复发

Source表示治疗方式，BM+PB表示骨髓移植+血液移植，PB表示血液移植

Status表示结局，3个水平，0表示删失，1表示复发，2表示竞争风险事件。

ftime表示时间

crrmod<-cuminc(ftime=bmt$ftime,fstatus=bmt$Status,group=bmt$D,cencode=0)
crrmod

ftime指定数据集中的生存时间，fstatus指定结局，group指定数据集中的分类变量，选项cencode指定数据集bmt中的结局变量status的删失值为0。

cuminc函数只能进行分类变量的单因素分析，无法进行连续性自变量的单因素分析，也无法进行任何形式自变量的多因素分析。

 这种时候无法呈现特定时间点的结果，若要呈现特定时间点的结果，可以利用timepoints（）函数

timepoints(w=crrmod,times=c(12,24,36,48,60))

w指定需要看的模型，times指定观测特定的时间点。

$est表示估计各时间点ALL和AML组的累积复发率（后缀1表示），后缀为2的表示累积竞争风险事件的发生率。例如ALL组12个月的累积复发率为34.2%，AML组中的累积12个月复发率为22.1% ， ALL组中12个月累积竞争风险事件发生率约为35.6%，AML组中12个月累积竞争风险事件发生率为43.4%。

$var表示估计各时间点ALL和AML组的累积复发率、累积竞争风险事件发生率的方差。

方差开二次方为标准差

利用点估计值及其标准差可计算其95%参考范围。

如果想要检验两个水平ALL和AML组是否存在统计学差异。

crrmod$Tests

第一行统计量为2.86，p值为0.09，表示ALL组和AML组的累积复发风险无统计学差异。第二行则表示累积竞争风险无统计学差异。

计算置信区间

x<-timepoints(w=crrmod,times=c(12,24,36,48,60))
x$est-1.96*sqrt(x$var)
x$est+1.96*sqrt(x$var)

 

生存曲线

1.生存曲线

par(mar=c(5,5,1,1))
plot(crrmod,
     xlab="Month",
     ylab = "CIF",
     lty=c(1,3,1,3),
     lwd = 2,
     col=c("#E7B800", "#E7B800","#2E9FDF","#2E9FDF"),
     cex=1.5,
     cex.axis=1.5,
     cex.lab=1.5)

设置图形边界，下、左、上、右依次为5，5，1，1

从图中看到ALL组的复发风险较AML组高，但未达到统计学差异。ALL组的竞争风险事件较AML组低，同样也未达到统计学意义。

2.复发曲线/竞争风险曲线

将复发风险和竞争风险的曲线分别绘制

par(mar=c(5,5,1,1))
plot(x=crrmod$`AML 1`$time,
     y=crrmod$`AML 1`$est,
     type = "s",
     col="#E7B800",
     lty=1,
     lwd=3,
     xlab="Month",
     ylab = "CIF", 
     ylim = c(0,1),
     cex=1.5,
     cex.axis=1.5,
     cex.lab=1.5)
lines(x=crrmod$`ALL 1`$time,
      y=crrmod$`ALL 1`$est,
      type = "s",
      col="#2E9FDF",
      lty=3,
      lwd=3)
legend(0,0.98,
       legend =c("AML","ALL"),
       col=c("#E7B800","#2E9FDF"),
       lty = c(1,3),
       lwd = 3,
       bty = "n")

bty=“n“表示不显示图例边框

par(mar=c(5,5,1,1))
plot(x=crrmod$`AML 2`$time,
     y=crrmod$`AML 2`$est,
     type = "s",
     col="#E7B800",
     lty=1,
     lwd=3,
     xlab="Month",
     ylab = "CIF", 
     ylim = c(0,1),
     cex=1.5,
     cex.axis=1.5,
     cex.lab=1.5)
lines(x=crrmod$`ALL 2`$time,
      y=crrmod$`ALL 2`$est,
      type = "s",
      col="#2E9FDF",
      lty=3,
      lwd=3)
legend(0,0.98,
       legend =c("AML","ALL"),
       col=c("#E7B800","#2E9FDF"),
       lty = c(1,3),
       lwd = 3,
       bty = "n")

3.添加置信区间

par(mar=c(5,5,1,1))
plot(x=crrmod[[1]][[1]],
     y=crrmod[[1]][[2]],
     type="s",
     ylim=c(0,1), 
     ylab="Cumulative Incidence",
     xlab="Time (Months)",
     cex=1.5,
     cex.axis=1.5,
     cex.lab=1.5)
x<- crrmod[names(crrmod) != 'Tests']
col=c("#2E9FDF","#E7B800")
for (i in 1:2) {
  lines(x[[i]][[1]],x[[i]][[2]],lty=1,lwd=2,col=col[i])
  lines(x[[i]][[1]],x[[i]][[2]]-1.96*sqrt(crrmod[[i]][[3]]),lty=3,col=col[i])
  lines(x[[i]][[1]],x[[i]][[2]]+1.96*sqrt(crrmod[[i]][[3]]),lty=3,col=col[i])
} 
legend(0,0.98,
       legend=c("AML","ALL"),
       lty = 1,
       lwd = 3,
       col=c("#E7B800","#2E9FDF"))

crrmod[[1]][[1]],crrmod[[1]][[2]]表示模型crrmod中AML组的各时间点及其对应的复发风险。将才不是“Test”的结果赋值为x，形成一个x的列表。

x[[1]]表示x中ALL1这个列表，x[[1]][[1]]表示ALL1中列表的time的数据，x[[1]][[2]]表示累积发生率，x[[1]][[3]]表示ALL1列表中的vars方差。

循环中1=ALL，2=AML

竞争风险曲线可信区间绘制

par(mar=c(5,5,1,1))
plot(x=crrmod[[3]][[1]],
     y=crrmod[[3]][[2]],
     type="s",
     ylim=c(0,1), 
     ylab="Cumulative Incidence",
     xlab="Time (Months)",
     cex=1.5,
     cex.axis=1.5,
     cex.lab=1.5)
x<- crrmod[names(crrmod) != 'Tests']
col=c("#2E9FDF","#E7B800")
for (i in 1:2) {
  lines(x[[i+2]][[1]],x[[i+2]][[2]],lty=1,lwd=2,col=col[i])
  lines(x[[i+2]][[1]],x[[i+2]][[2]]-1.96*sqrt(crrmod[[i+2]][[3]]),lty=3,col=col[i])
  lines(x[[i+2]][[1]],x[[i+2]][[2]]+1.96*sqrt(crrmod[[i+2]][[3]]),lty=3,col=col[i])
} 
legend(0,0.98,
       legend=c("AML","ALL"),
       lty = 1,
       lwd = 3,
       col=c("#E7B800","#2E9FDF"))

crr（）函数法（多因素分析）

cov <- data.frame(age = bmt$Age,
                  sex_F = ifelse(bmt$Sex=='F',1,0),
                  dis_AML = ifelse(bmt$D=='AML',1,0),
                  phase_cr1 = ifelse(bmt$Phase=='CR1',1,0),
                  phase_cr2 = ifelse(bmt$Phase=='CR2',1,0),
                  phase_cr3 = ifelse(bmt$Phase=='CR3',1,0),
                  source_PB = ifelse(bmt$Source=='PB',1,0)) 

将需要校正的变量提取，形成新的数据集，对phase进行哑变量变换。

crrfit <- crr(ftime=bmt$ftime, fstatus=bmt$Status, cov1=cov, failcode=1, cencode=0)
summary(crrfit)

其中cov1指定自变量的数据集cov，failcode=1表示1为结局事件，cencode=0表示0为删失事件，其他为竞争风险事件。

 结果解读：

phase_cr1和phase_cr2的p值小于0.05，其中phase_cr1的sHR及其可信区间为0.332（0.159-0.695），表示在控制竞争风险事件后，以复发作为参考，cr1阶段的患者复发的风险降低66.8%或者cr1阶段的患者复发的风险是relapse阶段患者0.332倍。

进一步判断phase整体是否有统计学意义，需要使用aod包中的wald.test()函数

#install.packages("aod")
library(aod)
wald.test(Sigma=crrfit$var,b=crrfit$coef,Terms = 4:6)

其中sigma指定模型中的方差，b指定模型中的偏回归系数，terms指定针对模型的4-6行进行整体的假设检验。

p值为0.0029，说明变量phase整体上存在统计学差异。

重点：

如果研究中采用Fine-Gray检验与竞争风险模型，那么COX比例风险模型必须同时拟合，将两种结果进行比对。

什么时候用：

如果存在竞争风险事件，而且极有可能对结论产生影响，这时采用Fine-Gray检验与竞争风险模型才是合适的。