Comparison of Methods in the Rotterdam Data

library(drsurv)
#> Loading required package: survival

df <- rotterdam
df$time <- pmin(df$rtime, df$dtime) / 365.25
df$status <- ifelse(df$recur == 1 | df$death == 1, 1, 0)
df$censor <- 1 - df$status
drFit.sjovan <-
 sjovan_survdiff(
   oformula = Surv(time, status) ~ chemo + year + age + meno +
     size + factor(grade) + nodes + pgr + er + hormon,
   ofunc = "survreg",
   cformula = Surv(time, censor) ~ chemo + year + age,
   cfunc = "survreg",
   eformula = chemo ~ year + age + meno + size +
     factor(grade) + nodes + pgr + er + hormon,
   method = "DR",
   times = c(2.5, 5, 7.5),
   se.type = "sandwich",
   data = df
 )
drFit.sjovan
#> $est.S1
#> [1] 0.7661741 0.6032550 0.5224925
#> 
#> $se.S1
#> [1] 0.03440273 0.02906199 0.14520126
#> 
#> $est.S0
#> [1] 0.7190130 0.5515392 0.4670844
#> 
#> $se.S0
#> [1] 0.07008758 0.07427899 0.10500988
#> 
#> $est.diff
#> [1] 0.04716109 0.05171581 0.05540801
#> 
#> $se.diff
#> [1] 0.06863352 0.07611256 0.14153032

df$chemo.fac <- factor(df$chemo)
drFit.blanche <-
 blanche_survdiff(
   oformula = Event(time, status) ~ chemo.fac + year + age + meno +
     size + factor(grade) + nodes + pgr + er + hormon,
   cformula = ~ chemo.fac + meno + size,
   eformula = chemo.fac ~ year + age + meno + size +
     factor(grade) + nodes + pgr + er + hormon,
   times = c(2.5, 5, 7.5),
   data = df
 )
drFit.blanche
#> $est.S1
#> [1] 0.7812813 0.6184004 0.5328929
#> 
#> $se.S1
#> [1] 0.02098236 0.02641496 0.02896027
#> 
#> $est.S0
#> [1] 0.7225998 0.5496811 0.4618080
#> 
#> $se.S0
#> [1] 0.01190859 0.01045986 0.01074570
#> 
#> $est.diff
#> [1] 0.05868142 0.06871926 0.07108498
#> 
#> $se.diff
#> [1] 0.02405841 0.02818684 0.03073858

drFit.pseudo <-
 pseudo_survdiff(
   oformula = Surv(time, status) ~ chemo + year + age + meno +
     size + factor(grade) + nodes + pgr + er + hormon,
   cformula = ~ chemo + meno + size,
   eformula = chemo ~ year + age + meno + size +
     factor(grade) + nodes + pgr + er + hormon,
   times = c(2.5, 5, 7.5),
   data = df, 
   R = 50
 )
drFit.pseudo
#> $est.S1
#> [1] 0.7814539 0.6209512 0.5389184
#> 
#> $se.S1
#> [1] 0.02405820 0.02644088 0.02639502
#> 
#> $est.S0
#> [1] 0.7264196 0.5509058 0.4616483
#> 
#> $se.S0
#> [1] 0.01231595 0.01186391 0.01214011
#> 
#> $est.diff
#> [1] 0.05503432 0.07004544 0.07727005
#> 
#> $se.diff
#> [1] 0.02717281 0.02809705 0.02835318