---
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(simtrial)
library(tidyr)
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/package=survival/vignettes/tiedtimes.pdf) and
the `?aeqSurv` function documentation.

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 `simtrial::wlr()`
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")
```
