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)。

為此我們先回顧一下生存分析中各個函數之間的關係：

生存函數survival function定義如下：

風險函數hazard function定義如下：

即：

一般我們設，叫做累計風險函數cumulative hazard function，因此上式可以寫成

基於以上式子，

基於上述過程，我們可以得到通過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)

上述代碼運行結果如下：

更多閱讀

-ykh

關注醫咖會，提高臨床研究水準！

），拉你進統計討論群和眾多熱愛研究的小夥伴們一起交流學習。