Direct, total and separable effects, using inverse probability weighting (IPW)

Note: Most of the code in this tutorial is written in R base, to mimic the original code from the published papers (see References). Dr. Elena Dudukina has a version of this code in R tidyverse here.

Conditional pooled logistic models

  • First, we run a set of conditional pooled logistic models, both as a first step to calculate the risk of prostate cancer death, and to calculate weights for death due to other causes, and weights for censoring due to loss to follow-up
### Pooled logistic regression for prostate cancer death (this will be needed in the next section)
plrFitP <-
  glm(
    prostateDeath ~ rx * (dtime + I(dtime ^ 2) + I(dtime ^ 3)) + normalAct + ageCat + hx + hgBinary ,
    data = longProstRed,
    family = binomial())

### Pooled logistic regression for death due to other causes

plrFitO <-
  glm(
    otherDeath ~ dtime + I(dtime ^ 2) + normalAct + ageCat + hx + hgBinary + rx ,
    data = longProstRed,
    family = binomial())


### Pooled logistic regression for censoring after time 50

temp50 <- longProstRed[longProstRed$dtime > 50,]

plrFitC <-
  glm(eventCens ~  normalAct + ageCat + hx + rx,
      data = temp50,
      family = binomial())
  • Fit models for the numerators in the stabilized weights corresponding to death due to other causes, and for censoring
plrFitCnum <-
  glm(eventCens ~ rx , data = temp50, family = binomial())

plrFitOnum <-
  glm(otherDeath ~ dtime + I(dtime ^ 2) + rx,
      data = longProstRed,
      family = binomial())

Weights calculation

Calculate weights for censoring due to loss to follow-up

### Denominator

longProstRed$predC <-
  rep(1, dim(longProstRed)[1]) # The probability of *not* being censored is 1 for all t <= 50

longProstRed$predC[longProstRed$dtime > 50] <- # update this probability from time 50 and onwards
  1 - predict(plrFitC, newdata = temp50, type = 'response')

### Numerator

longProstRed$predCnum <-
  rep(1, dim(longProstRed)[1])

longProstRed$predCnum [longProstRed$dtime > 50] <-
  1 - predict(plrFitCnum, newdata = temp50, type = 'response') 

predC <- longProstRed$predC

predCnum <- longProstRed$predCnum
  • Glimpse of data:
patno dtime predC predCnum
8 44 1.000000 1.0000000
8 45 1.000000 1.0000000
8 46 1.000000 1.0000000
8 47 1.000000 1.0000000
8 48 1.000000 1.0000000
8 49 1.000000 1.0000000
8 50 1.000000 1.0000000
8 51 0.963553 0.9693252
8 52 0.963553 0.9693252
8 53 0.963553 0.9693252
8 54 0.963553 0.9693252
8 55 0.963553 0.9693252
8 56 0.963553 0.9693252
8 57 0.963553 0.9693252
8 58 0.963553 0.9693252
8 59 0.963553 0.9693252
  • Calculate the cumulative probability over time:
cumPredC <-
  unlist(aggregate(predC ~ longProstRed$patno, FUN = cumprod)$predC,
         use.names = F)

cumPredCnum <-
  unlist(aggregate(predCnum ~ longProstRed$patno, FUN = cumprod)$predC,
         use.names = F)
  • Glimpse of data:
patno dtime predC predCnum cumPredC
8 44 1.000000 1.0000000 1.0000000
8 45 1.000000 1.0000000 1.0000000
8 46 1.000000 1.0000000 1.0000000
8 47 1.000000 1.0000000 1.0000000
8 48 1.000000 1.0000000 1.0000000
8 49 1.000000 1.0000000 1.0000000
8 50 1.000000 1.0000000 1.0000000
8 51 0.963553 0.9693252 0.9635530
8 52 0.963553 0.9693252 0.9284343
8 53 0.963553 0.9693252 0.8945956
8 54 0.963553 0.9693252 0.8619902
8 55 0.963553 0.9693252 0.8305733
8 56 0.963553 0.9693252 0.8003013
8 57 0.963553 0.9693252 0.7711327
8 58 0.963553 0.9693252 0.7430272
8 59 0.963553 0.9693252 0.7159460
  • Calculate unstabilized and stabilized weights
### Unstabilized weights
weights_cens <- 1 / cumPredC

### Stabilized weights

weights_cens_stab <- cumPredCnum / cumPredC
patno dtime predC predCnum cumPredC weights_cens weights_cens_stab
8 44 1.000000 1.0000000 1.0000000 1.000000 1.000000
8 45 1.000000 1.0000000 1.0000000 1.000000 1.000000
8 46 1.000000 1.0000000 1.0000000 1.000000 1.000000
8 47 1.000000 1.0000000 1.0000000 1.000000 1.000000
8 48 1.000000 1.0000000 1.0000000 1.000000 1.000000
8 49 1.000000 1.0000000 1.0000000 1.000000 1.000000
8 50 1.000000 1.0000000 1.0000000 1.000000 1.000000
8 51 0.963553 0.9693252 0.9635530 1.037826 1.005991
8 52 0.963553 0.9693252 0.9284343 1.077082 1.012017
8 53 0.963553 0.9693252 0.8945956 1.117824 1.018080
8 54 0.963553 0.9693252 0.8619902 1.160106 1.024178
8 55 0.963553 0.9693252 0.8305733 1.203988 1.030314
8 56 0.963553 0.9693252 0.8003013 1.249529 1.036486
8 57 0.963553 0.9693252 0.7711327 1.296794 1.042695
8 58 0.963553 0.9693252 0.7430272 1.345846 1.048941
8 59 0.963553 0.9693252 0.7159460 1.396753 1.055225

Calculate unstabilized and stabilized weights for death from other causes

### Denominator

predO <-
  1 - predict(plrFitO, newdata = longProstRed, type = 'response')

### Numerator
predOnum <-
  1 - predict(plrFitOnum, newdata = longProstRed, type = 'response')

### Calculate cumulative probability

cumPredO <-
  unlist(aggregate(predO ~ longProstRed$patno, FUN = cumprod)$predO,
         use.names = F)

cumPredOnum <-
  unlist(aggregate(predOnum ~ longProstRed$patno, FUN = cumprod)$predO,
         use.names = F) 

### Weight calculation

weights_compev <- 1 / cumPredO

weights_compev_stab <-
  cumPredOnum / cumPredO
  • Glimpse of data
patno dtime cumPredO cumPredOnum weights_compev weights_compev_stab
487 0 0.9892248 0.9816935 1.010893 0.9923866
487 1 0.9787493 0.9641208 1.021712 0.9850538
487 2 0.9685558 0.9472325 1.032465 0.9779844
487 3 0.9586275 0.9309829 1.043158 0.9711624
487 4 0.9489487 0.9153301 1.053798 0.9645728
487 5 0.9395048 0.9002352 1.064391 0.9582018
487 6 0.9302818 0.8856621 1.074943 0.9520364
487 7 0.9212666 0.8715776 1.085462 0.9460645

  • Multiply weights for censoring and for death due to other causes
weights_de <- weights_cens * weights_compev

weights_de_stab <- weights_cens_stab * weights_compev_stab
  • Assess weights
summary(weights_de)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.003   1.112   1.266   1.484   1.542  74.950
summary(weights_de_stab)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.4028  0.8615  0.9733  1.0156  1.0681 31.7854

Direct effect estimation

In this section and the following ones, we will apply the nonParametricCumHaz and nonParametricCumInc functions from the utility_functions script available in this repository. These functions are necessary to calculate hazards and cumulative incidence. Link

The functions look like this:

nonParametricCumHaz <-
  function(weightVector,
           inputdata,
           grp,
           outcomeProstate = TRUE) {
    outputHazards <- rep(NA, length.out = length(cutTimes))
    counter <- 1
    for (i in cutTimes) {
      if (outcomeProstate) {
        indices <-
          inputdata$dtime == i &
          inputdata$rx == grp &
          inputdata$eventCens == 0 & inputdata$otherDeath == 0
        eventIndicator <- indices & inputdata$prostateDeath == 1
      } else{
        indices <-
          inputdata$dtime == i & inputdata$rx == grp & inputdata$eventCens == 0
        eventIndicator <- indices & inputdata$otherDeath == 1
      }
      outputHazards[counter] <-
        sum(weightVector[eventIndicator]) / sum(weightVector[indices])
      counter <- counter + 1
    }
    return(outputHazards)
  }

nonParametricCumInc <- function(hazard1, hazard2, competing = FALSE) {
  inc <- rep(NA, length.out = length(cutTimes))
  cumulativeSurvival <- c(1, cumprod((1 - hazard1) * (1 - hazard2)))
  counter <- 1
  for (i in 1:length(cutTimes)) {
    if (!competing) {
      inc[i] <- hazard1[i] * (1 - hazard2[i]) * cumulativeSurvival[i]
    } else{
      inc[i] <- hazard1[i] * cumulativeSurvival[i]
    }
  }
  cumInc <- cumsum(inc)
  return(cumInc)
}

Following equation (37) from Young et al. Statistics in Medicine. 2020.

### Calculate cumulative hazards and cumulative incidence for the treated arm

#### Unstabilized
Treated_de_hazardP <-
  nonParametricCumHaz(
    weights_de,
    inputdata = longProstRed,
    grp = 1,
    outcomeProstate = TRUE)

Treated_de_hazardO <- rep(0, length.out = (length(cutTimes)))

cumIncTreated_de <-
  nonParametricCumInc(Treated_de_hazardP, Treated_de_hazardO)

#### Stabilized
Treated_de_hazardPs <-
  nonParametricCumHaz(
    weights_de_stab,
    inputdata = longProstRed,
    grp = 1,
    outcomeProstate = TRUE)

Treated_de_hazardOs <- rep(0, length.out = (length(cutTimes)))

cumIncTreated_des <-
  nonParametricCumInc(Treated_de_hazardPs, Treated_de_hazardOs)

### Calculate cumulative hazards and cumulative incidence for the placebo arm

#### Unstabilized

Placebo_de_hazardP <-
  nonParametricCumHaz(
    weights_de,
    inputdata = longProstRed,
    grp = 0,
    outcomeProstate = TRUE
  )

Placebo_de_hazardO <- rep(0, length.out = (length(cutTimes)))

cumIncPlacebo_de <-
  nonParametricCumInc(Placebo_de_hazardP, Placebo_de_hazardO)

#### Stabilized

Placebo_de_hazardPs <-
  nonParametricCumHaz(
    weights_de_stab,
    inputdata = longProstRed,
    grp = 0,
    outcomeProstate = TRUE
  )

Placebo_de_hazardOs <- rep(0, length.out = (length(cutTimes)))

cumIncPlacebo_des <-
  nonParametricCumInc(Placebo_de_hazardPs, Placebo_de_hazardOs)
  • Risks under each treatment arm at 60 months of follow-up, risk difference and risk ratio
## Unstabilized
cumIncTreated_de_60 <- cumIncTreated_de[length(cutTimes)] 

cumIncPlacebo_de_60 <- cumIncPlacebo_de[length(cutTimes)]

de_rr <- cumIncTreated_de_60 / cumIncPlacebo_de_60

de_rd <- cumIncTreated_de_60 - cumIncPlacebo_de_60

results_de <- data.frame(cumIncTreated_de_60, cumIncPlacebo_de_60, de_rr, de_rd)

knitr::kable(round(results_de, 2))
cumIncTreated_de_60 cumIncPlacebo_de_60 de_rr de_rd
0.36 0.37 0.97 -0.01 
## Stabilized

cumIncTreated_des_60 <- cumIncTreated_des[length(cutTimes)] 

cumIncPlacebo_des_60 <- cumIncPlacebo_des[length(cutTimes)]

des_rr <- cumIncTreated_des_60 / cumIncPlacebo_des_60

des_rd <- cumIncTreated_des_60 - cumIncPlacebo_des_60

results_des <- data.frame(cumIncTreated_des_60, cumIncPlacebo_des_60, des_rr, des_rd)

knitr::kable(round(results_des, 2))
cumIncTreated_des_60 cumIncPlacebo_des_60 des_rr des_rd
0.36 0.37 0.97 -0.01
  • Cumulative incidence curves
plot(
  cutTimes, cumIncTreated_de, type = "s",
  ylim = c(0, 1), ylab = "Risk", xlab = "Month",
  lty = 1, lwd = 1, xlim = c(0, 60),
  cex.lab = 1.2, cex.main = 1.2,
  main = c(
    paste("Risk of prostate cancer death"),
    paste("under elimination of competing events")
  )
)
lines(
  cutTimes,
  cumIncPlacebo_de,
  type = "s", col = 2, ylim = c(0, 1), lwd = 1)

legend(
  "topleft", c("Treatment", "Placebo"),
  col = c(1, 2, 1, 2),
  lty = c(1, 1, 2, 2),
  cex = 1.2, pt.cex = 1.2, lwd = 2)

Total effect estimation.

For this effect, we will only use the unstabilized weights we calculated for censoring due to loss to follow-up weights_cens, but using weights_cens_stab is possible as well.

Following equation (40) from Young et al. Statistics in Medicine. 2020.

### Calculate hazards and cumulative incidence for the treated arm

Treated_te_hazardP <-
  nonParametricCumHaz(
    weights_cens,
    inputdata = longProstRed,
    grp = 1,
    outcomeProstate = TRUE)

Treated_te_hazardO <-
  nonParametricCumHaz(
    weights_cens,
    inputdata = longProstRed,
    grp = 1,
    outcomeProstate = FALSE)

cumIncTreated_te <-
  nonParametricCumInc(Treated_te_hazardP, Treated_te_hazardO)

### Calculate cumulative hazards and cumulative incidence for the placebo arm

Placebo_te_hazardP <-
  nonParametricCumHaz(
    weights_cens,
    inputdata = longProstRed,
    grp = 0,
    outcomeProstate = TRUE)

Placebo_te_hazardO <-
  nonParametricCumHaz(
    weights_cens,
    inputdata = longProstRed,
    grp = 0,
    outcomeProstate = FALSE)

cumIncPlacebo_te <-
  nonParametricCumInc(Placebo_te_hazardP, Placebo_te_hazardO)
  • Risks under each treatment arm at 60 months of follow-up, risk difference and risk ratio
cumIncTreated_te_60 <- cumIncTreated_te[length(cutTimes)] 

cumIncPlacebo_te_60 <- cumIncPlacebo_te[length(cutTimes)]

te_rr <- cumIncTreated_te_60 / cumIncPlacebo_te_60

te_rd <- cumIncTreated_te_60 - cumIncPlacebo_te_60

results_te <- data.frame(cumIncTreated_te_60, cumIncPlacebo_te_60, te_rr, te_rd)

knitr::kable(round(results_te, 2))
cumIncTreated_te_60 cumIncPlacebo_te_60 te_rr te_rd
0.22 0.28 0.78 -0.06
  • Cumulative incidence curves
plot(
  cutTimes, cumIncTreated_te, type = "s",
  ylim = c(0, 1), ylab = "Risk", xlab = "Month",
  lty = 1, lwd = 1, xlim = c(0, 60),
  cex.lab = 1.2, cex.main = 1.2,
  main = c(
    paste("Risk of prostate cancer death"),
    paste("without elimination of competing events")
  )
)
lines(
  cutTimes,
  cumIncPlacebo_te,
  type = "s", col = 2, ylim = c(0, 1), lwd = 1)

legend(
  "topleft", c("Treatment", "Placebo"),
  col = c(1, 2, 1, 2), lty = c(1, 1, 2, 2),
  cex = 1.2, pt.cex = 1.2, lwd = 2)

Total effect on death for other causes

We will use the hazards obtained for the total effect, from the previous section.

cumIncTreated_te_comp <-
  nonParametricCumInc(Treated_te_hazardO, Treated_te_hazardP, competing = TRUE)

cumIncPlacebo_te_comp <-
  nonParametricCumInc(Placebo_te_hazardO, Placebo_te_hazardP, competing = TRUE)
  • Risks under each treatment arm at 60 months of follow-up, risk difference and risk ratio
cumIncTreated_te_comp_60 <- cumIncTreated_te_comp[length(cutTimes)] 

cumIncPlacebo_te_comp_60 <- cumIncPlacebo_te_comp[length(cutTimes)]

te_comp_rr <- cumIncTreated_te_comp_60 / cumIncPlacebo_te_comp_60

te_comp_rd <- cumIncTreated_te_comp_60 - cumIncPlacebo_te_comp_60

results_te_comp <- data.frame(cumIncTreated_te_comp_60, cumIncPlacebo_te_comp_60, te_comp_rr, te_comp_rd)

knitr::kable(round(results_te_comp, 2))
cumIncTreated_te_comp_60 cumIncPlacebo_te_comp_60 te_comp_rr te_comp_rd
0.51 0.43 1.19 0.08
  • Cumulative incidence curves
plot(
  cutTimes,
  cumIncTreated_te_comp, type = "s",
  ylim = c(0, 1), ylab = "Risk", xlab = "Month",
  lty = 1, lwd = 1, xlim = c(0, 60),
  cex.lab = 1.2, cex.main = 1.2,
  main = c(
    paste("Risk of other causes of death (the competing event)")
  ))

lines(
  cutTimes,
  cumIncPlacebo_te_comp,
  type = "s", col = 2, ylim = c(0, 1), lwd = 1)

legend(
  "topleft", c("Treatment", "Placebo"),
  col = c(1, 2, 1, 2), lty = c(1, 1, 2, 2),
  cex = 1.2, pt.cex = 1.2, lwd = 2)

Separable effects

Following equation (11) from Stensrud et al. Journal of the American Statistical Association. 2020

To this matter, we will use the long format data set prepared before: longProstRed.

  • Make one extra copy of the exposure rx, named Orx, to be used in the hazard of death due to other causes.
longProstRed$Orx <- longProstRed$rx
  • Fit conditional pooled logistic regression models
plrFitP <-
    glm(
      prostateDeath ~ rx*(dtime + I(dtime ^ 2)+ I(dtime ^ 3)) + normalAct + ageCat + hx + hgBinary,
      data = longProstRed,
      family = binomial()
    )
    
plrFitO <-
  glm(
    otherDeath ~ Orx * (dtime + I(dtime ^ 2) + I(dtime ^ 3)) + normalAct + ageCat + hx + hgBinary,
    data = longProstRed,
    family = binomial()
  )
  • Create three copies of the data. Data set construction where Ay and Ad = 1
treated <-
    baseline[rep(1:n, each = length(cutTimes)), ] #One row for each time

treated$dtime <- rep(cutTimes, n)
  
treated$rx <- 1

treated$Orx <- 1
  • Data set construction where Ay and Ad = 0
placebo <-
    baseline[rep(1:n, each = length(cutTimes)), ] #One row for each time
  
placebo$dtime <- rep(cutTimes, n)

placebo$rx <- 0

placebo$Orx <- 0
  • Dataset construction where Ay = 1, Ad = 0
treatAy <-
    baseline[rep(1:n, each = length(cutTimes)), ]
  
treatAy$dtime <- rep(cutTimes, n)
  
treatAy$rx <- 1
  
treatAy$Orx <- 0
  • Calculate hazards for each time interval in the three cloned data sets
## For Ay = 1 and Ad = 1

treated$hazardP <-
    predict(plrFitP, newdata = treated, type = 'response')
  
treated$hazardO <-
    predict(plrFitO, newdata = treated, type = 'response')
  
treated$s <- (1 - treated$hazardP) * (1 - treated$hazardO)

## For Ay = 0, Ad = 0

placebo$hazardP <-
    predict(plrFitP, newdata = placebo, type = 'response')

placebo$hazardO <-
    predict(plrFitO, newdata = placebo, type = 'response')

placebo$s <- (1 - placebo$hazardP) * (1 - placebo$hazardO)

## For Ay = 1, Ad = 0

treatAy$hazardP <-
    predict(plrFitP, newdata = treatAy, type = 'response')
  
treatAy$hazardO <-
    predict(plrFitO, newdata = treatAy, type = 'response')
  
treatAy$s <- (1 - treatAy$hazardP) * (1 - treatAy$hazardO)
  • Calculate weights for censoring for lost to follow-up
plrFitC <-
  glm(
    eventCens ~ Orx* (dtime + I(dtime ^ 2) + I(dtime ^ 3)) + normalAct + ageCat +
      hx + hgBinary,
    data = longProstRed,
    family = binomial()
  )
  
predC <- 1 - predict(plrFitC, newdata = longProstRed, type = 'response')

# Create a new data frame with predicted hazard estimates from treated and 
weightFrame <-
  data.frame(
    dtime = treated$dtime,
    patno = treated$patno,
    predOa1 = (1 - placebo$hazardO),
    predOa0 = (1 - treated$hazardO)
  )

compWeightsO <-
  merge(longProstRed,
        weightFrame,
        by = c('patno', 'dtime'),
        all.x = TRUE, sort = FALSE)

predOa0 <- compWeightsO$predOa0

predOa1 <- compWeightsO$predOa1

cum_pred_O_0 <-
  unlist(aggregate(predOa0 ~ compWeightsO$patno, FUN = cumprod)$predOa0,
         use.names = F) #unlist because aggregate creates list

cum_pred_O_1 <-
  unlist(aggregate(predOa1 ~ compWeightsO$patno, FUN = cumprod)$predOa1,
         use.names = F) #unlist because aggregate creates list

cum_pred_C <-
  unlist(aggregate(predC ~ compWeightsO$patno, FUN = cumprod)$predC,
         use.names = F) #unlist because aggregate creates list

ipw_d <- rep(NA, length.out = length(cum_pred_C))

# index time zeros

t0 <- longProstRed$prostateDeath == 0

ipw_d <- cum_pred_O_1 / cum_pred_O_0

ipw_d <- cum_pred_O_1 / cum_pred_O_0

ipw_cens <- 1 / cum_pred_C

ipw_sep_eff <- ipw_d / cum_pred_C

Last, we will use the function discrete_cuminc_prost, from the utility_functions script available in this repository.

discrete_cuminc_prost <- function(weight_vector, inputdata, grp=0, outcome_y=TRUE,follow_up=1:max_time){
  event_vec <- rep(NA, length.out=length(follow_up))
  counter <- 1 
  # count number of individuals in grp (that is, we cound those who were present at baseline)
  n_grp <- sum(inputdata$dtime==0 & inputdata$rx==grp)
  for(i in follow_up){
    if(outcome_y){
      indices <- inputdata$dtime==i & inputdata$rx == grp & inputdata$eventCens==0 & inputdata$otherDeath==0 
      eventIndicator <- indices & inputdata$prostateDeath==1 
    }else{
      indices <- inputdata$dtime==i & inputdata$rx == grp & inputdata$eventCens==0  
      eventIndicator <- indices &  inputdata$otherDeath==1
    }
    event_vec[counter] <- sum(weight_vector[eventIndicator]) / n_grp
    counter <- counter+1
  }
  output_cuminc <- cumsum(event_vec)
  return(output_cuminc)
}
cumIncTreatAyIPW <-
  discrete_cuminc_prost(ipw_sep_eff, longProstRed, grp = 1, follow_up = cutTimes)

cumIncTreatedIPW <-
  discrete_cuminc_prost(ipw_cens, longProstRed, grp = 1, follow_up = cutTimes)

cumIncPlaceboIPW <-
  discrete_cuminc_prost(ipw_cens, longProstRed, grp = 0, follow_up = cutTimes)


estimatesFullData <-
  data.frame(cumIncTreatedIPW, cumIncTreatAyIPW, cumIncPlaceboIPW)

knitr::kable(round(estimatesFullData[36,], 2))
cumIncTreatedIPW cumIncTreatAyIPW cumIncPlaceboIPW
36 0.14 0.15 0.2
  • Cumulative incidence curves
plot(cutTimes,cumIncTreatedIPW, type="s",ylim=c(0,1), ylab="Cumulative incidence", 
     xlab="Months",lty=1,lwd=1,xlim=c(0,52))
lines(cutTimes,cumIncPlaceboIPW, type="s",col=2,ylim=c(0,1),lwd=1,lty=1)
lines(cutTimes,cumIncTreatAyIPW, type="s",col=3,ylim=c(0,1),lwd=1,lty=1)
legend(x=0,y=1, c(expression('Pr('*'Y'['k+1']^{"a"*'=1'} *'=1' * ')' ),expression('Pr('*'Y'['k+1']^{"a"['Y']*'=1' * ' , ' * "a"['D']*'=0'} *'=1' * ')' ),expression('Pr('*'Y'['k+1']^{"a"*'=0'} *'=1' * ')' )),bty = "n",
       col=c(1,3,2), lty=c(1,1,1),cex=1,pt.cex=1,lwd=1)

References

  • Young JG, Stensrud MJ, Tchetgen Tchetgen EJ, Hernán MA. A causal framework for classical statistical estimands in failure time settings with competing events. Statistics in Medicine 2020. Paper link Code

  • Stensrud MJ, Young JG, Didelez V, Robins JM, Hernán MA. Separable Effects for Causal Inference in the Presence of Competing Events. Journal of the American Statistical Association 2022. Paper link

  • Rojas-Saunero LP, Young JG, Didelez V, Ikram MA, Swanson SA. Considering questions before methods in dementia research with competing events and causal goals. American Journal of Epidemiology. 2023. Paper link Code