R语言计算竞争风险数据中累积发生率置信区间并绘图（1）

竞争风险模型就是指在临床事件中出现和它竞争的结局事件，这是事件会导致原有结局的改变，因此叫做竞争风险模型。比如我们想观察患者肿瘤的复发情况，但是患者在观察期突然车祸死亡，或者因其他疾病死亡，这样我们就观察不到复发情况了，这种情况下不能把缺失数据仅仅当做右删失处理，这样的话会造成数据的估值错误。这是我们应该优先选择竞争风险模型来做数据分析，而不是COX回归。我们既往已经多篇文章介绍了R语言竞争风险风险模型，需要的请看文章合集，粉丝后台问如何计算竞争风险模型累积发生率的可信区间，如下面文章所示


今天我们通过R语言来演示一下，

需要用到的是一个妊娠数据（公众号回复：妊娠数据，可以获得数据），我们导入数据和R包看看

library(etm)
library(survival)
abortion <-read.csv("E:/r/test/renshen.csv",sep=',',header=TRUE)
1
2
3

这是一个关于使用香豆素衍生物对妊娠结果影响的数据。 目的是调查暴露于香豆素衍生物是否会增加自然流产的可能性。 除自然流产外，妊娠可能以人工流产或活产告终，从而导致相互竞争的风险局面。Id：编码，entry：进入队列时间，exit：退出队列时间，group组别：0：对照组，1：暴露于香豆素衍生物组，cause：结局变量1：人工流产，2：活产，3：自然流产。
我们先来建立竞争风险模型，在这里要解释一下Surv(entry, exit, cause != 0)这一句和我们的生存分析的代码是一样的entry, exit相当于时间time，cause != 0相当于不等于0等变量，就是就是相当于status. failcode：表示感兴趣的结局指标，这里默认为1，就是人工流产，等下画图也是以1来画。

cif.abortion <- etmCIF(Surv(entry, exit, cause != 0) ~ group,
                       abortion, etype = cause, failcode = 3)
1
2

这样两个组别的标准误就出来了，我们可以通过标准误计算可信区间，更详细的信息可以

summary(cif.abortion)
1

我们还可以进一步画图，默认就是failcode = 3，表示对照组和暴露组的cif

plot(cif.abortion)
1

假设我们想了解第27周的风险和可信区间

plot(cif.abortion, curvlab = c("Control", "Exposed"), ylim = c(0, 0.6),
     ci.type = "bars", pos.ci = 27, col = c(1, 2), ci.lwd = 6,
     lwd = 2, lty = 1, cex = 1.3)
1
2
3

有些杂志是发黑白的图片，我们也可以换成黑白的图

plot(cif.abortion, curvlab = c("Control", "Exposed"), ylim = c(0, 0.6),
     ci.type = "bars", pos.ci = c(27, 28), col = c(1, 1), ci.lwd = 6,
     lwd = 2, lty = c(2, 1), cex = 1.3)
1
2
3

如果我们想了解整个随访期的可信区间，

plot(cif.abortion, curvlab = c("Control", "Exposed"), ylim = c(0, 0.5),
     ci.type = "pointwise", col = c(1, 2), lwd = 2, lty = 1, cex = 1.3)
1
2

如果我们想了解其他指标的的cif,可以使用which.cif 函数指定，对于本例数据which.cif = c(1, 2)

plot(cif.abortion, which.cif = c(1, 2), ylim = c(0, 0.8), lwd = 2,
       col = c(1, 1, 2, 2), lty = c(1, 2, 1, 2), legend = FALSE)
legend(0, 0.8, c("Control", "Exposed"), col = c(1, 2), lty = 1,
         bty = "n", lwd = 2)
legend(0, 0.7, c("ETOP", "Life Birth"), col = 1, lty = c(1, 2),
         bty = "n", lwd = 2)
1
2
3
4
5
6