cause_specific_vs_subdist.R

library("survival")
library("cmprsk")
library("RColorBrewer")

set.seed(66)

# Standard function for simulating survival with a Weibull baseline hazard
# https://stats.stackexchange.com/questions/135124/how-to-create-a-toy-survival-time-to-event-data-with-right-censoring

# baseline hazard: Weibull
# x = feature/covariate to simulate from
# lambda = scale parameter in h0()
# rho = shape parameter in h0()
# beta = fixed effect parameter

simulWeib <- function(x, lambda, rho, beta )
{

  N = length(x)
  
  # Weibull latent event times
  v = runif(n=N)
  tstop = (- log(v) / (lambda * exp(x * beta)))^(1 / rho)
  event = rep(1,N)
  
  # data set
  df = data.frame(id=1:N,
                  tstop=tstop,
                  event=event,
                  x=x)
  return(df)
}

# simulate feature with direct effect on the first outcome
N_tot = 5e3

# causal hazard ratio for the primary event
hr1 = 1.5
all_hr2 = c(0.5,1,2)

# offsets the competing event
# set this to 10 to minimize competing effects
event2_offset = 0

# for plotting
clr = brewer.pal(3,"Set1")
clr = c(clr[2],clr[1])
par(mfrow=c(1,3))

for ( i in 1:length(all_hr2) ) {
  cur_hr2 = all_hr2[i]
  # simulate feature with direct effect on primary outcome
  x = rbinom(N_tot,1,0.5)
  event1 = simulWeib( x , lambda=0.0001, rho=3.5, beta = log(hr1) )
  
  # simulate feature with direct effect on a competing outcome
  event2 = simulWeib( x , lambda=0.0001, rho=3.5, beta = log(cur_hr2) )
  event2$tstop = event2$tstop + event2_offset
  
  # censor on second event and estimate using cause-specific model
  cs_event = event1
  compete = cs_event$tstop > event2$tstop
  cs_event$tstop[ compete ] = event2$tstop[ compete ]
  cs_event$event[ compete ] = 0
  cs_fit <- coxph(Surv(tstop, event) ~ x , data=cs_event )
  cshr = summary(cs_fit)$coef[1,2]
  cat( cshr , '\n' )
  
  # estimate using Fine-Gray model
  cr_event = event1
  compete = cr_event$tstop > event2$tstop
  cr_event$tstop[ compete ] = event2$tstop[ compete ]
  cr_event$event[ compete ] = 2
  
  crr_fit = crr(cr_event$tstop, cr_event$event,cr_event$x)
  sdhr = summary(crr_fit)$coef[1,2]
  cat( sdhr , '\n' )
  cat('\n')
  
  # plot
  kmfit = survfit(Surv(tstop, event) ~ x, data=cs_event)

  ttl = paste("Competing HR=",round(cur_hr2,2),"\nCSHR=",round(cshr,2)," SDHR=",round(sdhr,2),sep='')
  plot( kmfit , fun="event" , col=paste(clr,50,sep=''),lwd=4,las=1 , xlab="Time", ylab="Cumulative Incidence",main=ttl)
  
  crr_ci = cuminc(cr_event$tstop, cr_event$event,cr_event$x)
  lines( crr_ci[["0 1"]]$time , crr_ci[["0 1"]]$est , col=clr[1] , lty=3 , lwd=2 )
  lines( crr_ci[["1 1"]]$time , crr_ci[["1 1"]]$est , col=clr[2] , lty=3 , lwd=2 )
  
  legend("topleft",legend=c("Cause-specific","Subdistribution"),bty="n",lty=c(1,3),lwd=c(3,2))
  
}