How can I calculate the measures of effect modification on the "additives scale" and "multiplicative scale" using a coxph?

Question

I read the paper entitled "Recommendations for presenting analyses of effect modification and interaction" published by Knol MJ and VanderWeele TJ Link to download the paper.

I would like to follow their recommendations. I think I get everything to follow their recommendation (1), (2), and (4). However, I do not know how to calculate the measure of effect modification on "additives scale" and "multiplicative scale" using a coxph?

I wrote a simple example.

Can anybody help me out here?

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?

Update #1

The answer of AdamO works fine Thank you a lot! However, it does not work with a multivariable model adjusted for several covariables. May you please show how one have to adapt your code in this situation?

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

Answer 1

RERI: Relative excess risk due to interaction is defined in Rothman's Epidemiology an Introduction as:

A multiplicative scale measure of interaction would be defined as:

I'm using RR to refer to relative risk, but it could be risk difference (RD) instead. If the model using a log-link to transform the outcome like with Poisson regression or Cox models, then the multiplicative scale difference is just the product-term parameter. Similarly, if the identity-link is used like in an additive risk model, the additive scale of interaction is just the product-term parameter. That doesn't mean that multiplicative interactions are not interesting in identity link models, or additive interactions are not interesting in log link models.

With the Cox model, the multiplicative interaction is just the coefficient for the product term. You'd calculate the RERI as:

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

And multiplicative interaction:

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