如何使用coxph计算对“添加剂量表”和“乘法量表”进行效果修改的度量？

Posted 2021-04-29

我阅读了由Knol MJ和VanderWeele TJ Link to download the paper发表的题为“提出效果修改与互动分析的建议”的论文。

我想按照他们的建议。我想我得到的一切都遵循他们的建议（1），（2）和（4）。但是，我不知道如何使用coxph计算“添加剂量表”和“乘法量表”的效果修改度量？

我写了一个简单的例子。

有人可以帮帮我吗？

library(survival)
library(Hmisc)
library(papeR)

#########################################################
# loading and preparing data
#########################################################
data(colon)
d <- colon[, Cs(time, status, rx, surg)]
rm(colon)
names(d) <- Cs(time, status, group, modifier)
d$group <- ifelse(d$group == "Obs", 0, 1)

d$group.4.stata <- NA
d$group.4.stata <- ifelse(d$group==0 & d$modifier==0, 0, d$group.4.stata)
d$group.4.stata <- ifelse(d$group==0 & d$modifier==1, 1, d$group.4.stata)
d$group.4.stata <- ifelse(d$group==1 & d$modifier==0, 2, d$group.4.stata)
d$group.4.stata <- ifelse(d$group==1 & d$modifier==1, 3, d$group.4.stata)

d$group.4.stata <- factor(d$group.4.stata, 0:3, c("control + short", "control + long", "intervention + short", "intervention + long"))

#########################################################
# number of patients in each stratum
#########################################################
table(d$group.4.stata)

#########################################################
# Recommendation (1)
# HR for each statum of exposure (control=0 versus intervention=1) and potential effect modifier (time from surgery to registration: 0=short, 1=long)
#########################################################
prettify(fit <- summary(coxph(Surv(time, status) ~ group.4.stata, data=d)), digits=3)

#########################################################
# Recommendation (2)
# HR for each statum of exposure (control=0 versus intervention=1) within stata of potential effect modifier (time from surgery to registration: 0=short, 1=long)
#########################################################
prettify(fit <- summary(coxph(Surv(time, status) ~ group, data=d[which(d$modifier==0),])), digits=3)
prettify(fit <- summary(coxph(Surv(time, status) ~ group, data=d[which(d$modifier==1),])), digits=3)

#########################################################
# Recommendation (3)
# meassure of effect modification on "additives scale" or "multiplicative scale"?
#########################################################
prettify(fit <- summary(coxph(Surv(time, status) ~ group*modifier, data=d)), digits=3) # What kind of measure is this?

# And how can I calculate the other measure?

更新＃1

AdamO的答案很好，非常感谢你！但是，它不适用于针对多个协变量调整的多变量模型。请您说明在这种情况下如何调整代码？

library(survival)
library(Hmisc)
library(papeR)
library(epiR)
#########################################################
# loading and preparing data
#########################################################
data(colon)
d <- colon[, Cs(time, status, rx, surg, sex, age, node4)]
rm(colon)
names(d) <- Cs(time, status, group, modifier, covariate1, covariate2, covariate3)
d$group <- ifelse(d$group == "Obs", 0, 1)

d$group.4.stata <- NA
d$group.4.stata <- ifelse(d$group==0 & d$modifier==0, 0, d$group.4.stata)
d$group.4.stata <- ifelse(d$group==0 & d$modifier==1, 1, d$group.4.stata)
d$group.4.stata <- ifelse(d$group==1 & d$modifier==0, 2, d$group.4.stata)
d$group.4.stata <- ifelse(d$group==1 & d$modifier==1, 3, d$group.4.stata)

d$group.4.stata <- factor(d$group.4.stata, 0:3, c("control + short", "control + long", "intervention + short", "intervention + long"))

#########################################################
# number of patients in each stratum
#########################################################
table(d$group.4.stata)

#########################################################
# HR for each statum of exposure (control=0 versus intervention=1) and potential effect modifier (time from surgery to registration: 0=short, 1=long)
#########################################################
prettify(fit <- summary(coxph(Surv(time, status) ~ group.4.stata+covariate1+covariate2+covariate3, data=d)), digits=3)

#########################################################
# HR for each statum of exposure (control=0 versus intervention=1) within stata of potential effect modifier (time from surgery to registration: 0=short, 1=long)
#########################################################
prettify(fit <- summary(coxph(Surv(time, status) ~ group+covariate1+covariate2+covariate3, data=d[which(d$modifier==0),])), digits=3)
prettify(fit <- summary(coxph(Surv(time, status) ~ group+covariate1+covariate2+covariate3, data=d[which(d$modifier==1),])), digits=3)

#########################################################
# meassure of effect modification on "additives scale" or "multiplicative scale"?
#########################################################
prettify(fit <- summary(coxph(Surv(time, status) ~ group*modifier+covariate1+covariate2+covariate3, data=d)), digits=3) # What kind of measure is this? 

# And how can I calculate the other measure?

#########################################################
# calculating RERI 
#########################################################
f <- function(parms) exp(sum(parms)) - exp(parms[1]) - exp(parms[2]) + 1
df <- function(parms) c( ## derivative of RERI for use in delta-method
  exp(sum(parms)) - exp(parms[1]), ## or use bootstrap
  exp(sum(parms)) - exp(parms[2]),
  exp(sum(parms))
)

fit <- coxph(Surv(time, status) ~ group*modifier+covariate1+covariate2+covariate3, data=d)
parms <- coef(fit)
reri <- f(parms)
v.reri <- t(df(parms)) %*% vcov(fit) %*% matrix(df(parms))

## RERI and 95% CI
reri + c(0,-1,1)*1.96*c(sqrt(v.reri))

答案

RERI：由于相互作用导致的相对超额风险在Rothman的流行病学和引言中定义为：

乘法量表的相互作用衡量标准定义如下：

我使用RR来指代相对风险，但它可能是风险差异（RD）。如果模型使用对数链接来转换结果，例如Poisson回归或Cox模型，则乘法量表差异只是乘积项参数。类似地，如果像添加风险模型那样使用身份链接，则交互的附加规模只是产品期限参数。这并不意味着乘法交互在身份链接模型中不感兴趣，或者附加交互在日志链接模型中不感兴趣。

使用Cox模型，乘法相互作用只是乘积项的系数。您将RERI计算为：

x <- rnorm(100)
w <- rnorm(100)
y <- rexp(100, exp(-3 + 0.4*x + 1.2*w + 0.4*w*x))
fit <- coxph(Surv(y) ~ x*w)

f <- function(parms) exp(sum(parms)) - exp(parms[1]) - exp(parms[2]) + 1
df <- function(parms) c( ## derivative of RERI for use in delta-method
  exp(sum(parms)) - exp(parms[1]), ## or use bootstrap
  exp(sum(parms)) - exp(parms[2]),
  exp(sum(parms))
)
parms <- coef(fit)
reri <- f(parms)
v.reri <- t(df(parms)) %*% vcov(fit) %*% matrix(df(parms))

## RERI and 95% CI
reri + c(0,-1,1)*1.96*c(sqrt(v.reri))

和乘法互动：

exp(c(coef(fit)[3], confint(fit, 3)))