R教程:如何通过cox回归模型计算RD(危险差)和NNT?

专题合集更多教程

作者:胥洋

 

RD:risk difference,危险差,指的是干预组与对照组相比,危险减少或增加的绝对值,即RD=Prtreated - Pruntreated

 

NNT/H:number needed to treat/harm,当我们研究保护性因素时,即Prtreated < Pruntreated时,常常报道NNT,指的是为了避免一例结局事件的发生,平均多少个人需要接受有利治疗;当我们研究危险性因素时,即Prtreated > Pruntreated时,常常报道NNH,指的是为了诱发一例结局事件的发生,平均多少人需要接受有害治疗。

 

综上,NNT/H=1÷abs(RD),当RD<0时,为NNT,当RD>0时,为NNH。

 

在生存分析中为什么要报道 RD 和 NNT/H,仅仅报道HR(hazard ratio, 风险比)不够吗?

 

在回答这个问题之前,我们可以将流行病学中常见的效应值分为两大类:

 

一类是相对效应值,如risk ratio,odds ratio和hazard ratio,表示的是与对照组相比,治疗组效应减少或增加的相对值;

 

另一个是绝对效应值,如risk difference,表示的是与对照组相比,治疗组效应减少或增加的绝对值。与相对效应值相比,绝对效应值更具有临床意义。(详见:总结:那些可以评价干预措施效果的指标们

 

举例来说,某种新药的随机对照试验得到该药与安慰剂相比HR为0.5,显示出该药可以显著降低结局事件的发生率。但是,当结局事件在安慰剂组的发生率为1%和50%时,HR等于0.5对于临床实践而言意义却完全不一样。

 

假设结局事件为皮疹,当结局事件在安慰剂组的发生率为1%,结局事件本身发生率不高且不严重,该药物与安慰剂相比降低了0.5%的结局事件发生率,降低有限,所以该药物临床意义不大;

 

当结局事件在安慰剂组的发生率为50%,结局事件本身发生率很高虽不严重,该药物与安慰剂相比降低了25%的结局事件发生率,降低显著,所以该药物临床意义较大。

 

那么,怎么计算生存数据中的 RD 和 NNT/H 呢?

 

在生存数据中,

 

通过上式,我们可以看到要求 RD 和 NNT/H 的关键在于求得Prtreated(T<t)和Pruntreated(T<t)。

 

为此我们先回顾一下生存分析中各个函数之间的关系

 

假设T是一个大于0的随机变量,表示从随访开始到结局事件发生的时间,T的概率密度函数(probability density function, p.d.f.)为

,累积分布函数(cumulative distribution function, c.d.f.)为

 

生存函数survival function定义如下

风险函数hazard function定义如下

即:

一般我们设

 

,叫做累计风险函数cumulative hazard function,因此上式可以写成

 

基于以上式子,

通过cox回归,我们可以得到α和β的估计值,为了求得上式,我们就需要知道S0(t)和Pr(X),显然Pr(X)是未知的,问题到这里似乎无解了,但是我们不妨变换一下思路,

为了求得此式,我们需要知道,因为所有人在Zi=0和Xi=0时,基线生存函数均一致,因此题就简化为求S0(t)了,通过cox回归模型求S0(t)常见的方法有两种,分别为 Breslow‘s estimator 和 Kalbfleish & Prentice estimator,这里就不详细介绍了,感兴趣的读者可以自行阅读。

 

基于上述过程,我们可以得到通过cox回归模型计算 RD 和 NNT/H 的方法

 

1. 拟合cox回归模型;

 

2. 生成治疗数据集,即所有人均接受治疗,设treatment=1,其他协变量取值与原始数据集一致;

 

3. 生成对照数据集,即所有人均不接受治疗,设treatment=0,其他协变量取值与原始数据集一致;

 

4. 指定 RD 和 NNT/H 的预测时间点,即在那些随访时间点上计算 RD 和NNT/H;

 

5. 在治疗数据集上预测每个患者在设定的时间点上的生存概率;

 

6. 在对照数据集上预测每个患者在设定的时间点上的生存概率;

 

7. 计算治疗数据集上所有人在设定的时间点上的生存概率的均值;

 

8. 计算对照数据集上所有人在设定的时间点上的生存概率的均值;

 

9. 将两组均值做差得到在设定的时间点上的RD值,将RD值取倒数得到在设定的时间点上的NNT/H;

 

10. 用bootstrap的方法求得RD和NNT/H的置信区间。

 

上述过程的R代码如下:

 

## import packages

library(rms)

library(withr)

library(dplyr)

library(doParallel)

 

## set cores

cl <- makeCluster(3)

registerDoParallel(cl)

 

## define functions

inverse.logit <- function(x) {

  return(exp(x) / (1 + exp(x)))

}

RD_NNT.H.cal <- function(data, times, treatment, formula) {

  data.treated <- data

  data.treated[ , treatment] <- 1

  data.untreated <- data

  data.untreated[ , treatment] <- 0

  cox.res <- cph(formula, data, surv = T, x = T, y = T)

  treated.sur <- survest(cox.res, newdata = data.treated, times = time) 

  untreated.sur <- survest(cox.res, newdata = data.untreated, times = time)

  treated.Mean <- colMeans(1 - treated.sur$surv) 

  untreated.Mean <- colMeans(1 - untreated.sur$surv)

  RD <- treated.Mean - untreated.Mean

  type <- ifelse(RD >= 0, 'NNH', 'NNT')

  NNT.H <- 1 / abs(RD)

  return(data.frame(time = time, RD = RD, type = type, NNT.H = NNT.H))

}

 

## generate survival data

set.seed(20190401)

n = 1000

alpha0 = log(1)

alpha1 = log(3)

beta0 = log(1)

beta1 = log(3)

beta2 = log(4)

 

X <- rbinom(n, 1, 0.5)

E <- rbinom(n, 1, inverse.logit(alpha0 + alpha1 * X))

Ytime <- rexp(n, rate = exp(beta0 + beta1 * X + beta2 * E))

Ctime <- runif(n, min = 3, max = 6)

Time <- ifelse(Ytime <= Ctime, Ytime, Ctime)

Status <- ifelse(Ytime <= Ctime, 1, 0)

dta <- data.frame(X, E, Time, Status)

 

## calculate NNH

time <- seq(1, 4, length.out = 20)

res.pe <- RD_NNT.H.cal(dta, time, 'E', Surv(Time, Status) ~ X + E)

 

## caculate the CIs

res.CI <- foreach (i = 1 : 1000, 

                   .packages = c('dplyr', 'rms', 'withr'), 

                   .combine = 'rbind') %dopar% {

  with_seed(i, dta.sample <- dta %>% sample_frac(1, replace = T))

  res <- RD_NNT.H.cal(dta.sample, time, 'E', Surv(Time, Status) ~ X + E)

}

res.CI <- res.CI %>% 

  group_by(time) %>% 

  summarise(RD.p2.5 = quantile(RD, 0.025), 

            RD.p97.5 = quantile(RD, 0.975), 

            NNT.H.p2.5 = quantile(NNT.H, 0.025), 

            NNT.H.p97.5 = quantile(NNT.H, 0.975))

 

## combine results

res <- merge(res.pe, res.CI, by = 'time') %>% 

  select(time, type, RD, RD.p2.5, RD.p97.5, NNT.H, NNT.H.p2.5, NNT.H.p97.5)

 

上述代码运行结果如下:

描述问题
选择一个标签 (请选择一个与您问题最相符的标签)
提交问题
我要提问
描述问题
选择一个标签 (请选择一个与您问题最相符的标签)
提交问题
描述问题
选择一个标签 (请选择一个与您问题最相符的标签)
    提交问题