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.7812829 0.6184085 0.5329110
#> 
#> $se.S1
#> [1] 0.02098232 0.02641465 0.02895948
#> 
#> $est.S0
#> [1] 0.7226016 0.5496908 0.4618324
#> 
#> $se.S0
#> [1] 0.01190855 0.01045965 0.01074519
#> 
#> $est.diff
#> [1] 0.05868132 0.06871770 0.07107869
#> 
#> $se.diff
#> [1] 0.02405836 0.02818650 0.03073770

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.02561835 0.03040323 0.03038897
#> 
#> $est.S0
#> [1] 0.7264196 0.5509058 0.4616483
#> 
#> $se.S0
#> [1] 0.011559520 0.009288898 0.008860051
#> 
#> $est.diff
#> [1] 0.05503432 0.07004544 0.07727005
#> 
#> $se.diff
#> [1] 0.02707052 0.02860870 0.03050650