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.）为
生存函数survival function定义如下：
风险函数hazard function定义如下：
即：
一般我们设
基于以上式子，
通过cox回归
为了求得此式，我们需要知道
基于上述过程，我们可以得到通过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)
上述代码运行结果如下：
更多阅读
医咖会微信：medieco-ykh
关注医咖会，提高临床研究水平！
快加小咖个人微信（xys2018ykf），拉你进统计讨论群和众多热爱研究的小伙伴们一起交流学习。
点击左下角“阅读原文”，看看医咖会既往推送了哪些统计教程。或者使用电脑打开网址：http://www.mediecogroup.com/，查看70种SPSS教程。