Skip to content

Commit 01066b0

Browse files
nanxstatsnanxstats
committed
Stylize discrepancy vignette
1 parent dd234a0 commit 01066b0

File tree

1 file changed

+125
-123
lines changed

1 file changed

+125
-123
lines changed

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

Lines changed: 125 additions & 123 deletions
Original file line numberDiff line numberDiff line change
@@ -15,44 +15,46 @@ library(gsDesign2)
1515
library(dplyr)
1616
library(tibble)
1717
library(gt)
18-
#library(ggplot2)
19-
#library(cowplot)
2018
library(simtrial)
2119
library(tidyr)
22-
#library(future.batchtools)
23-
#library(doFuture)
24-
#library(foreach)
25-
26-
#library(tictoc)
27-
2820
library(survival)
29-
30-
3121
```
3222

33-
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.
34-
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*).
23+
In the survival (base R) package, the log-rank and Cox estimation procedures
24+
apply (by default) a correction to "fix" roundoff errors.
25+
These are implemented with the `timefix` option (by default `timefix = TRUE`)
26+
via the `aeqSurv()` function.
27+
However, in the simtrial package, (and also Hmisc), such a correction is not
28+
implemented; Consequently, there can be discrepancies between simtrial and
29+
base R survival (`survdiff()`, `coxph()`, and `survfit()`).
3530

36-
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)
31+
For details on the `aeqSurv()` function, see [Therneau,
32+
2016](https://cran.r-project.org/package=survival/vignettes/tiedtimes.pdf) and
33+
the `?aeqSurv` function documentation.
3734

38-
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.
35+
In the following, we describe a simulation scenario where a discrepancy is
36+
generated and illustrate how discrepancies can be resolved (if desired) by
37+
pre-processing survival times with `aeqSurv()` and thus replicating `survdiff()`
38+
and `coxph()` default calculations.
3939

40-
In the simulated dataset two observations are generated:
41-
42-
- Observation $i=464$ with survival time $Y=0.306132722582$
43-
- Observation $i=516$ with survival time $Y=0.306132604679$
44-
- Per "aeqSurv" these times are tied and set to $Y=0.306132604679$
45-
- The log-rank and Cox estimates can therefore differ between other approaches without the "timefix" correction
40+
In the simulated dataset, two observations are generated:
4641

42+
- Observation $i=464$ with survival time $Y=0.306132722582$.
43+
- Observation $i=516$ with survival time $Y=0.306132604679$.
44+
- Per `aeqSurv()`, these times are tied and set to $Y=0.306132604679$.
45+
- The log-rank and Cox estimates can therefore differ between other approaches
46+
without the "timefix" correction.
4747

4848
## Scenario definitions
4949

50-
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.
51-
50+
We define various true data generating model scenarios and convert for use in
51+
gsDesign2. Here, we are using a single scenario where discrepancies were found.
52+
This is just for illustration to inform the user of simtrial that discrepancies
53+
can occur and how to resolve via `aeqSurv()`, if desired.
5254

5355
```{r}
5456
survival_at_24_months <- 0.35
55-
hr <- log(.35)/log(.25)
57+
hr <- log(.35) / log(.25)
5658
control_median <- 12
5759
control_rate <- c(log(2) / control_median, (log(.25) - log(.2)) / 12)
5860
@@ -79,32 +81,32 @@ scenarios <- tribble(
7981
5, "Weak null", 0, 0, 1,
8082
5, "Weak null", 1, 24, .25,
8183
5, "Weak null", 2, 12, .2,
82-
6, "Strong null", 0, 0, 1,
84+
6, "Strong null", 0, 0, 1,
8385
6, "Strong null", 1, 3, exp(-3 * control_rate[1] * 1.5),
8486
6, "Strong null", 2, 3, exp(-6 * control_rate[1]),
8587
6, "Strong null", 3, 18, .25,
8688
6, "Strong null", 4, 12, .2,
87-
)
89+
)
8890
# scenarios |> gt()
8991
```
9092

91-
9293
```{r}
93-
fr <-
94-
scenarios |>
95-
group_by(Scenario) |>
94+
fr <- scenarios |>
95+
group_by(Scenario) |>
9696
# filter(Scenario == 2) |>
97-
mutate(Month = cumsum(duration),
98-
x_rate = -(log(Survival) - log(lag(Survival, default = 1))) /
99-
duration,
100-
rate = ifelse(Month > 24, control_rate[2], control_rate[1]),
101-
hr = x_rate / rate) |>
102-
select(-x_rate) |>
103-
filter(Period > 0, Scenario > 0) |> ungroup()
104-
#fr |> gt() |> fmt_number(columns = everything(), decimals = 2)
105-
106-
fr <- fr |> mutate(fail_rate = rate, dropout_rate =0.001, stratum = "All")
107-
97+
mutate(
98+
Month = cumsum(duration),
99+
x_rate = -(log(Survival) - log(lag(Survival, default = 1))) /
100+
duration,
101+
rate = ifelse(Month > 24, control_rate[2], control_rate[1]),
102+
hr = x_rate / rate
103+
) |>
104+
select(-x_rate) |>
105+
filter(Period > 0, Scenario > 0) |>
106+
ungroup()
107+
# fr |> gt() |> fmt_number(columns = everything(), decimals = 2)
108+
109+
fr <- fr |> mutate(fail_rate = rate, dropout_rate = 0.001, stratum = "All")
108110
109111
# MWLR
110112
mwlr <- fixed_design_mb(
@@ -115,86 +117,93 @@ mwlr <- fixed_design_mb(
115117
study_duration = 36
116118
) |> to_integer()
117119
118-
119120
er <- mwlr$enroll_rate
120-
121121
```
122122

123-
# A scenario that generates a discrepancy
124-
123+
## A scenario that generates a discrepancy
125124

126125
```{r}
127-
128126
set.seed(3219)
129127
130-
dgm <- fr[c(14:17),]
128+
dgm <- fr[c(14:17), ]
129+
130+
fail_rate <- data.frame(
131+
stratum = rep("All", 2 * nrow(dgm)),
132+
period = rep(dgm$Period, 2),
133+
treatment = c(
134+
rep("control", nrow(dgm)),
135+
rep("experimental", nrow(dgm))
136+
),
137+
duration = rep(dgm$duration, 2),
138+
rate = c(dgm$rate, dgm$rate * dgm$hr)
139+
)
131140
132-
fail_rate <- data.frame(stratum = rep("All", 2 * nrow(dgm)),
133-
period = rep(dgm$Period, 2),
134-
treatment = c(rep("control", nrow(dgm)),
135-
rep("experimental", nrow(dgm))),
136-
duration = rep(dgm$duration, 2),
137-
rate = c(dgm$rate, dgm$rate * dgm$hr)
138-
)
139-
140141
dgm$stratum <- "All"
141142
# Constant dropout rate for both treatment arms and all scenarios
142-
dropout_rate <- data.frame(stratum = rep("All", 2),
143-
period = rep(1, 2),
144-
treatment = c("control", "experimental"),
145-
duration = rep(100, 2),
146-
rate = rep(.001, 2)
147-
)
143+
dropout_rate <- data.frame(
144+
stratum = rep("All", 2),
145+
period = rep(1, 2),
146+
treatment = c("control", "experimental"),
147+
duration = rep(100, 2),
148+
rate = rep(.001, 2)
149+
)
148150
```
149151

150-
151-
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*])
152+
Simulated dataset with discrepancy between logrank test of `simtrial::wlr()`
153+
and `survdiff()` (also compare to score test of `coxph()` [same as `survdiff()`
154+
with default `timefix = TRUE`]).
152155

153156
```{r}
154-
155157
ss <- 395
156158
157-
set.seed(8316951+ss*1000)
158-
159+
set.seed(8316951 + ss * 1000)
160+
159161
# Generate a dataset
160-
dat <- sim_pw_surv(n = 698, enroll_rate = er,
161-
fail_rate = fail_rate, dropout_rate = dropout_rate)
162+
dat <- sim_pw_surv(
163+
n = 698,
164+
enroll_rate = er,
165+
fail_rate = fail_rate,
166+
dropout_rate = dropout_rate
167+
)
162168
163169
analysis_data <- cut_data_by_date(dat, 36)
164170
165171
dfa <- analysis_data
166172
167-
dfa$treat <- ifelse(dfa$treatment=="experimental",1,0)
173+
dfa$treat <- ifelse(dfa$treatment == "experimental", 1, 0)
168174
169-
z1 <- dfa |> wlr(weight=fh(rho=0,gamma=0))
175+
z1 <- dfa |> wlr(weight = fh(rho = 0, gamma = 0))
170176
171-
check <- survdiff(Surv(tte,event)~ treat, data=dfa)
177+
check <- survdiff(Surv(tte, event) ~ treat, data = dfa)
172178
173-
# Note, for coxph use
174-
#cph.score <- summary(coxph(Surv(tte,event)~ treat, data=dfa, control=coxph.control(timefix=TRUE)))$sctest
179+
# Note, for `coxph()`, use
180+
# cph.score <- summary(coxph(Surv(tte, event) ~ treat, data = dfa, control = coxph.control(timefix = TRUE)))$sctest
175181
176-
cat("Log-rank wlr() vs survdiff()",c(z1$z^2,check$chisq),"\n")
182+
cat("Log-rank wlr() vs survdiff()", c(z1$z^2, check$chisq), "\n")
177183
```
178184

179-
Verify that *timefix=FALSE* in *coxph* agrees with *wlr*
185+
Verify that `timefix = FALSE` in `coxph()` agrees with `wlr()`:
180186

181187
```{r}
182-
cph.score <- summary(coxph(Surv(tte,event)~ treat, data=dfa, control=coxph.control(timefix=FALSE)))$sctest
183-
cat("Log-rank wlr() vs Cox score z^2",c(z1$z^2,cph.score["test"]),"\n")
188+
cph.score <- summary(coxph(
189+
Surv(tte, event) ~ treat,
190+
data = dfa,
191+
control = coxph.control(timefix = FALSE)
192+
))$sctest
193+
cat("Log-rank wlr() vs Cox score z^2", c(z1$z^2, cph.score["test"]), "\n")
184194
```
185195

196+
Pre-processing survival times with `aeqSurv()` to implement `timefix = TRUE` procedure.
186197

187-
Pre-processing survival times with *aeqSurv* to implement *timefix=TRUE* procedure.
188-
189-
Verify *wlr* and *survdiff* now agree.
198+
Verify `wlr()` and `survdiff()` now agree.
190199

191200
```{r}
192-
Y <- dfa[,"tte"]
193-
Delta <- dfa[,"event"]
201+
Y <- dfa[, "tte"]
202+
Delta <- dfa[, "event"]
194203
195-
tfixed <- aeqSurv(Surv(Y,Delta))
196-
Y<- tfixed[,"time"]
197-
Delta <- tfixed[,"status"]
204+
tfixed <- aeqSurv(Surv(Y, Delta))
205+
Y <- tfixed[, "time"]
206+
Delta <- tfixed[, "status"]
198207
# Use aeqSurv version
199208
dfa$tte2 <- Y
200209
dfa$event2 <- Delta
@@ -203,72 +212,65 @@ dfa$event2 <- Delta
203212
dfa2 <- dfa
204213
dfa2$tte <- dfa2$tte2
205214
dfa2$event <- dfa2$event2
206-
z1new <- dfa2 |> wlr(weight=fh(rho=0,gamma=0))
207-
cat("Log-rank wlr() with timefix vs survdiff() z^2",c(z1new$z^2,check$chisq),"\n")
208-
215+
z1new <- dfa2 |> wlr(weight = fh(rho = 0, gamma = 0))
216+
cat("Log-rank wlr() with timefix vs survdiff() z^2", c(z1new$z^2, check$chisq), "\n")
209217
```
210218

211-
212-
Where do they differ (tte2 are times after *aeqSurv*) ?
219+
Where do they differ (`tte2` are times after `aeqSurv()`)?
213220

214221
```{r}
222+
dfa <- dfa[order(dfa$tte2), ]
215223
216-
dfa <- dfa[order(dfa$tte2),]
217-
218-
id <- seq(1,nrow(dfa))
224+
id <- seq(1, nrow(dfa))
219225
220226
diff <- exp(dfa$tte) - exp(dfa$tte2)
221-
id_diff <- which(abs(diff)>0)
227+
id_diff <- which(abs(diff) > 0)
222228
223-
tolook <- seq(id_diff-2,id_diff+2)
229+
tolook <- seq(id_diff - 2, id_diff + 2)
224230
225-
dfcheck <- dfa[tolook,c("tte","tte2","event","event2","treatment")]
226-
print(dfcheck,digits=12)
231+
dfcheck <- dfa[tolook, c("tte", "tte2", "event", "event2", "treatment")]
232+
print(dfcheck, digits = 12)
227233
```
228234

235+
Verify `coxph()` (default) and `coxph()` with `aeqSurv()` pre-processing
236+
(using `tte2` as outcome and setting `timefix = FALSE`) are identical:
229237

230-
231-
232-
Verify *coxph* (default) and *coxph* with aeqSurv pre-processing (using tte2 as outcome and setting *timefix=FALSE*) are identical:
233-
234-
Also note that here ties do not have impact because in separate arms
238+
Also note that here ties do not have impact because in separate arms.
235239

236240
```{r}
237241
# Check Cox with ties
238-
cox_breslow <- summary(coxph(Surv(tte,event)~treatment,data=dfa,ties="breslow"))$conf.int
239-
cox_efron <- summary(coxph(Surv(tte,event)~treatment,data=dfa,ties="efron"))$conf.int
240-
cat("Cox Breslow and Efron hr (tte, timefix=TRUE):",c(cox_breslow[1],cox_efron[1]),"\n")
242+
cox_breslow <- summary(coxph(Surv(tte, event) ~ treatment, data = dfa, ties = "breslow"))$conf.int
243+
cox_efron <- summary(coxph(Surv(tte, event) ~ treatment, data = dfa, ties = "efron"))$conf.int
244+
cat("Cox Breslow and Efron hr (tte, timefix=TRUE):", c(cox_breslow[1], cox_efron[1]), "\n")
241245
242246
# Here ties do not have impact because in separate arms
243-
cox_breslow <- summary(coxph(Surv(tte2,event2)~treatment,data=dfa,ties="breslow", control=coxph.control(timefix=FALSE)))$conf.int
244-
cox_efron <- summary(coxph(Surv(tte2,event2)~treatment,data=dfa,ties="efron", control=coxph.control(timefix=FALSE)))$conf.int
245-
cat("Cox Breslow and Efron hr (tte2, timefix=FALSE):",c(cox_breslow[1],cox_efron[1]),"\n")
247+
cox_breslow <- summary(coxph(Surv(tte2, event2) ~ treatment, data = dfa, ties = "breslow", control = coxph.control(timefix = FALSE)))$conf.int
248+
cox_efron <- summary(coxph(Surv(tte2, event2) ~ treatment, data = dfa, ties = "efron", control = coxph.control(timefix = FALSE)))$conf.int
249+
cat("Cox Breslow and Efron hr (tte2, timefix=FALSE):", c(cox_breslow[1], cox_efron[1]), "\n")
246250
```
247251

248-
**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**.
249-
250-
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:
252+
**So here there is a difference between `tte` and `tte2` times, but there is
253+
not an impact of ties for Cox between `"breslow"` and `"efron"` because the ties
254+
(single tie in `tte2`) are in separate arms**.
251255

256+
Lastly, artificially change treatment so that two observations are tied within
257+
the same treatment arm which generates difference between `"breslow"` and
258+
`"efron"` options for `ties`:
252259

253260
```{r}
254-
# Create tie within treatment arm by changing treatment
261+
# Create tie within treatment arm by changing treatment
255262
dfa3 <- dfa
256-
dfa3[19,"treat"] <- 1.0
257-
258-
cox_breslow <- summary(coxph(Surv(tte,event)~treat, data=dfa3,ties="breslow", control=coxph.control(timefix=TRUE)))$conf.int
259-
cox_efron <- summary(coxph(Surv(tte,event)~treat, data=dfa3,ties="efron", control=coxph.control(timefix=TRUE)))$conf.int
260-
cat("Cox Breslow and Efron hr (tte, timefix=TRUE)=",c(cox_breslow[1],cox_efron[1]),"\n")
263+
dfa3[19, "treat"] <- 1.0
261264
265+
cox_breslow <- summary(coxph(Surv(tte, event) ~ treat, data = dfa3, ties = "breslow", control = coxph.control(timefix = TRUE)))$conf.int
266+
cox_efron <- summary(coxph(Surv(tte, event) ~ treat, data = dfa3, ties = "efron", control = coxph.control(timefix = TRUE)))$conf.int
267+
cat("Cox Breslow and Efron hr (tte, timefix=TRUE)=", c(cox_breslow[1], cox_efron[1]), "\n")
262268
```
263269

264-
265270
Same as
266271

267272
```{r}
268-
269-
cox_breslow <- summary(coxph(Surv(tte2,event2)~treat, data=dfa3,ties="breslow", control=coxph.control(timefix=FALSE)))$conf.int
270-
cox_efron <- summary(coxph(Surv(tte2,event2)~treat, data=dfa3,ties="efron", control=coxph.control(timefix=FALSE)))$conf.int
271-
cat("Cox Breslow and Efron hr (tte2, timefix=FALSE)=",c(cox_breslow[1],cox_efron[1]),"\n")
272-
273+
cox_breslow <- summary(coxph(Surv(tte2, event2) ~ treat, data = dfa3, ties = "breslow", control = coxph.control(timefix = FALSE)))$conf.int
274+
cox_efron <- summary(coxph(Surv(tte2, event2) ~ treat, data = dfa3, ties = "efron", control = coxph.control(timefix = FALSE)))$conf.int
275+
cat("Cox Breslow and Efron hr (tte2, timefix=FALSE)=", c(cox_breslow[1], cox_efron[1]), "\n")
273276
```
274-

0 commit comments

Comments
 (0)