Merge pull request #308 from Merck/larry-leon-patch-1

Add vignette on potential discrepancies between simtrial and survdiff

Author: LittleBeannie

_pkgdown.yml

@@ -92,6 +92,7 @@ articles:
   - maxcombo
   - rmst
   - parallel
+  - discrepancy-between-simtrial-and-survival
 - title: "NPH distribution approximations"
   contents:
   - arbitrary-hazard

vignettes/discrepancy-between-simtrial-and-survival.Rmd

@@ -0,0 +1,274 @@
+---
+title: "Note on potential discrepancies between simtrial and survdiff"
+author: "Larry Leon"
+output: rmarkdown::html_vignette
+vignette: >
+  %\VignetteIndexEntry{Note on potential discrepancies between simtrial and survdiff}
+  %\VignetteEngine{knitr::rmarkdown}
+---
+
+## Overview
+
+```{r, message=FALSE}
+library(gsDesign)
+library(gsDesign2)
+library(dplyr)
+library(tibble)
+library(gt)
+#library(ggplot2)
+#library(cowplot)
+library(simtrial)
+library(tidyr)
+#library(future.batchtools)
+#library(doFuture)
+#library(foreach)
+
+#library(tictoc)
+
+library(survival)
+
+
+```
+
+In the **survival** (base R) package the log-rank and Cox estimation procedures apply (by default) a correction to "fix" roundoff errors.   These are implemented with the *timefix* option (by default *timefix = TRUE*) via the *aeqSurv* function.   
+However in the **simtrial** package (and also **Hmisc**) such a correction is not implemented;  Consequently, there can be discrepancies between **simtrial** and base R *survival* (*survdiff*, *coxph*, and *survfit*).
+
+For details on the *aeqSurv* function see [Therneau, 2016](https://cran.r-project.org/web/packages/survival/vignettes/tiedtimes.pdf) and [R documentation, version 3.803](https://www.rdocumentation.org/packages/survival/versions/3.8-3/topics/aeqSurv)
+
+In the following we describe a simulation scenario where a discrepancy is generated and illustrate how discrepancies can be resolved (if desired) by pre-processing survival times with *aeqSurv* and thus replicating *survdiff* and *coxph* default calculations.
+
+In the simulated dataset two observations are generated:
+
+-   Observation $i=464$ with survival time $Y=0.306132722582$
+-   Observation $i=516$ with survival time $Y=0.306132604679$
+-   Per "aeqSurv" these times are tied and set to $Y=0.306132604679$
+-   The log-rank and Cox estimates can therefore differ between other approaches without the "timefix" correction
+
+
+## Scenario definitions
+
+We define various true data generating model scenarios and convert for use in **gsDesign2**.  Here we are using a single scenario where discrepancies were found.   This is just for illustration to inform the user of **simtrial** that discrepancies can occur and how to resolve via *aeqSurv*, if desired.
+
+
+```{r}
+survival_at_24_months <- 0.35
+hr <- log(.35)/log(.25)
+control_median <- 12
+control_rate <- c(log(2) / control_median, (log(.25) - log(.2)) / 12)
+
+scenarios <- tribble(
+  ~Scenario, ~Name,           ~Period, ~duration, ~Survival,
+  0,         "Control",       0,       0,         1,
+  0,         "Control",       1,       24,        .25,
+  0,         "Control",       2,       12,        .2,
+  1,         "PH",            0,       0,         1,
+  1,         "PH",            1,       24,        .35,
+  1,         "PH",            2,       12,        .2^hr,
+  2,         "3-month delay", 0,       0,         1,
+  2,         "3-month delay", 1,       3,         exp(-3 * control_rate[1]),
+  2,         "3-month delay", 2,       21,        .35,
+  2,         "3-month delay", 3,       12,        .2^hr,
+  3,         "6-month delay", 0,       0,         1,
+  3,         "6-month delay", 1,       6,         exp(-6 * control_rate[1]),
+  3,         "6-month delay", 2,       18,        .35,
+  3,         "6-month delay", 3,       12,        .2^hr,
+  4,         "Crossing",      0,       0,         1,
+  4,         "Crossing",      1,       3,         exp(-3 * control_rate[1] * 1.3),
+  4,         "Crossing",      2,       21,        .35,
+  4,         "Crossing",      3,       12,        .2^hr,
+  5,         "Weak null",     0,       0,         1,
+  5,         "Weak null",     1,       24,        .25,
+  5,         "Weak null",     2,       12,        .2,
+  6,         "Strong null",   0,       0,         1,         
+  6,         "Strong null",   1,       3,         exp(-3 * control_rate[1] * 1.5),
+  6,         "Strong null",   2,       3,         exp(-6 * control_rate[1]),
+  6,         "Strong null",   3,       18,        .25,
+  6,         "Strong null",   4,       12,        .2,
+  )
+# scenarios |> gt()
+```
+
+
+```{r}
+fr <- 
+  scenarios |>
+    group_by(Scenario) |>
+  #  filter(Scenario == 2) |>
+    mutate(Month = cumsum(duration),
+           x_rate = -(log(Survival) - log(lag(Survival, default = 1))) / 
+           duration, 
+           rate = ifelse(Month > 24, control_rate[2], control_rate[1]),
+           hr = x_rate / rate) |>
+    select(-x_rate) |>
+    filter(Period > 0, Scenario > 0) |> ungroup()
+#fr |> gt() |> fmt_number(columns = everything(), decimals = 2)
+
+fr <- fr |> mutate(fail_rate = rate, dropout_rate =0.001, stratum = "All")
+
+
+# MWLR
+mwlr <- fixed_design_mb(
+  tau = 12,
+  enroll_rate = define_enroll_rate(duration = 12, rate = 1),
+  fail_rate = fr |> filter(Scenario == 2),
+  alpha = 0.025, power = .85, ratio = 1,
+  study_duration = 36
+) |> to_integer()
+
+
+er <- mwlr$enroll_rate
+
+```
+
+# A scenario that generates a discrepancy 
+
+
+```{r}
+
+set.seed(3219)
+
+dgm <- fr[c(14:17),]
+
+fail_rate <- data.frame(stratum = rep("All", 2 * nrow(dgm)),
+                          period = rep(dgm$Period, 2),
+                          treatment = c(rep("control", nrow(dgm)), 
+                                        rep("experimental", nrow(dgm))),
+                          duration = rep(dgm$duration, 2),
+                          rate = c(dgm$rate, dgm$rate * dgm$hr)
+                         )
+    
+dgm$stratum <- "All"
+# Constant dropout rate for both treatment arms and all scenarios
+dropout_rate <- data.frame(stratum = rep("All", 2), 
+                           period = rep(1, 2),
+                           treatment = c("control", "experimental"),
+                           duration = rep(100, 2),
+                           rate = rep(.001, 2)
+                          )
+```
+
+
+Simulated dataset with discrepancy between logrank test of *wlr* (**simtrial**) and *survdiff* (also compare to score test of *coxph* [same as survdiff with default *timefix=TRUE*])
+
+```{r}
+
+ss <- 395
+
+set.seed(8316951+ss*1000)  
+  
+# Generate a dataset
+dat <- sim_pw_surv(n = 698, enroll_rate = er,
+fail_rate = fail_rate, dropout_rate = dropout_rate)
+
+analysis_data <- cut_data_by_date(dat, 36)
+
+dfa <- analysis_data
+
+dfa$treat <- ifelse(dfa$treatment=="experimental",1,0)
+
+z1 <- dfa |>  wlr(weight=fh(rho=0,gamma=0))
+
+check <- survdiff(Surv(tte,event)~ treat, data=dfa)
+
+# Note, for coxph use 
+#cph.score <- summary(coxph(Surv(tte,event)~ treat, data=dfa, control=coxph.control(timefix=TRUE)))$sctest
+
+cat("Log-rank wlr() vs survdiff()",c(z1$z^2,check$chisq),"\n")
+```
+
+Verify that *timefix=FALSE* in *coxph* agrees with *wlr* 
+
+```{r}
+cph.score <- summary(coxph(Surv(tte,event)~ treat, data=dfa, control=coxph.control(timefix=FALSE)))$sctest
+cat("Log-rank wlr() vs Cox score z^2",c(z1$z^2,cph.score["test"]),"\n")
+```
+
+
+Pre-processing survival times with *aeqSurv* to implement *timefix=TRUE* procedure.  
+
+Verify *wlr* and *survdiff* now agree.
+
+```{r}
+Y <- dfa[,"tte"]
+Delta <- dfa[,"event"]
+
+tfixed <- aeqSurv(Surv(Y,Delta))
+Y<- tfixed[,"time"]
+Delta <- tfixed[,"status"]
+# Use aeqSurv version
+dfa$tte2 <- Y
+dfa$event2 <- Delta
+
+# wlr() after "timefix"
+dfa2 <- dfa
+dfa2$tte <- dfa2$tte2
+dfa2$event <- dfa2$event2
+z1new <- dfa2 |>  wlr(weight=fh(rho=0,gamma=0))
+cat("Log-rank wlr() with timefix vs survdiff() z^2",c(z1new$z^2,check$chisq),"\n")
+
+```
+
+
+Where do they differ (tte2 are times after *aeqSurv*) ?
+
+```{r}
+
+dfa <- dfa[order(dfa$tte2),]
+
+id <- seq(1,nrow(dfa))
+
+diff <- exp(dfa$tte) - exp(dfa$tte2)
+id_diff <- which(abs(diff)>0)
+
+tolook <- seq(id_diff-2,id_diff+2)
+
+dfcheck <-  dfa[tolook,c("tte","tte2","event","event2","treatment")]
+print(dfcheck,digits=12)
+```
+
+
+
+
+Verify *coxph* (default) and *coxph* with aeqSurv pre-processing (using tte2 as outcome and setting *timefix=FALSE*) are identical:
+
+Also note that here ties do not have impact because in separate arms
+
+```{r}
+# Check Cox with ties
+cox_breslow <- summary(coxph(Surv(tte,event)~treatment,data=dfa,ties="breslow"))$conf.int
+cox_efron <- summary(coxph(Surv(tte,event)~treatment,data=dfa,ties="efron"))$conf.int
+cat("Cox Breslow and Efron hr (tte, timefix=TRUE):",c(cox_breslow[1],cox_efron[1]),"\n")
+
+# Here ties do not have impact because in separate arms
+cox_breslow <- summary(coxph(Surv(tte2,event2)~treatment,data=dfa,ties="breslow", control=coxph.control(timefix=FALSE)))$conf.int
+cox_efron <- summary(coxph(Surv(tte2,event2)~treatment,data=dfa,ties="efron", control=coxph.control(timefix=FALSE)))$conf.int
+cat("Cox Breslow and Efron hr (tte2, timefix=FALSE):",c(cox_breslow[1],cox_efron[1]),"\n")
+```
+
+**So here there is a difference between tte and tte2 times, but there is not an impact of ties for Cox between *breslow* and *efron* because the ties (single tie in tte2) are in separate arms**.
+
+Lastly, artificially change treatment so that two observations are tied within the same treatment arm which generates difference between *breslow* and *efron* options for ties:
+
+
+```{r}
+# Create tie within treatment arm by changing treatment 
+dfa3 <- dfa
+dfa3[19,"treat"]  <- 1.0
+
+cox_breslow <- summary(coxph(Surv(tte,event)~treat, data=dfa3,ties="breslow", control=coxph.control(timefix=TRUE)))$conf.int
+cox_efron <- summary(coxph(Surv(tte,event)~treat, data=dfa3,ties="efron", control=coxph.control(timefix=TRUE)))$conf.int
+cat("Cox Breslow and Efron hr (tte, timefix=TRUE)=",c(cox_breslow[1],cox_efron[1]),"\n")
+
+```
+
+
+Same as
+
+```{r}
+
+cox_breslow <- summary(coxph(Surv(tte2,event2)~treat, data=dfa3,ties="breslow", control=coxph.control(timefix=FALSE)))$conf.int
+cox_efron <- summary(coxph(Surv(tte2,event2)~treat, data=dfa3,ties="efron", control=coxph.control(timefix=FALSE)))$conf.int
+cat("Cox Breslow and Efron hr (tte2, timefix=FALSE)=",c(cox_breslow[1],cox_efron[1]),"\n")
+
+```
+