The two-stage estimation (TSE) approach first estimates the effect of treatment switching on survival, then uses this estimate to derive counterfactual survival times if switching had not occurred. A key requirement of the TSE approach is that treatment switching must occur at or after a disease-related “secondary baseline” time point, typically defined by disease progression.
The TSE method assumes the no unmeasured confounding condition, meaning treatment switching must be independent of potential outcomes, provided all relevant patient characteristics at the secondary baseline are measured and included in the model. The simple TSE method further assumes that no time-dependent confounding exists between the secondary baseline and the time of switching.
The simple TSE method involves applying an accelerated failure time (AFT) model that compares post-progression survival in control group switchers with post-progression survival in control group non-switchers. Prognostic variables measured at the secondary baseline are included to account for differences between the groups. The treatment effect of switching is estimated as a time ratio and used to adjust the survival times of switchers.
Although \(\psi\) is estimated solely from control group patients who experienced disease progression, it will also be applied to adjust the survival times of patients who switched treatments before disease progression, under the assumption that there are only a limited number of such cases.
Once \(\psi\) has been estimated, we can derive an adjusted data set and fit a (potentially stratified) Cox proportional hazards model to the adjusted data set to obtain an estimate of the hazard ratio. The confidence interval for the hazard ratio can be derived by bootstrapping the entire adjustment and subsequent model-fitting process.
We start by preparing the data.
# the eventual survival time
shilong1 <- shilong %>%
arrange(bras.f, id, tstop) %>%
group_by(bras.f, id) %>%
slice(n()) %>%
select(-c("ps", "ttc", "tran"))
# the last value of time-dependent covariates before pd
shilong2 <- shilong %>%
filter(pd == 0 | tstart <= dpd) %>%
arrange(bras.f, id, tstop) %>%
group_by(bras.f, id) %>%
slice(n()) %>%
select(bras.f, id, ps, ttc, tran)
# combine baseline and time-dependent covariates
shilong3 <- shilong1 %>%
left_join(shilong2, by = c("bras.f", "id"))Next we apply the simple TSE method.
fit1 <- tsesimp(
data = shilong3, time = "tstop", event = "event",
treat = "bras.f", censor_time = "dcut", pd = "pd",
pd_time = "dpd", swtrt = "co", swtrt_time = "dco",
base_cov = c("agerand", "sex.f", "tt_Lnum", "rmh_alea.c",
"pathway.f"),
base2_cov = c("agerand", "sex.f", "tt_Lnum", "rmh_alea.c",
"pathway.f", "ps", "ttc", "tran"),
aft_dist = "weibull", alpha = 0.05,
recensor = TRUE, swtrt_control_only = FALSE, offset = 1,
boot = FALSE)We can examine the Weibull AFT model fits and the corresponding value of \(\hat{\psi}\).
# control group
fit1$fit_aft[[1]]$fit$parest[, c("param", "beta", "sebeta", "z")]
#> param beta sebeta z
#> 1 (Intercept) 5.376908384 0.733428671 7.3311947
#> 2 swtrt 1.136190154 0.234370294 4.8478420
#> 3 agerand 0.001992035 0.007352436 0.2709354
#> 4 sex.fFemale 0.352456115 0.239501797 1.4716220
#> 5 tt_Lnum -0.025630824 0.040554642 -0.6320072
#> 6 rmh_alea.c -0.917219583 0.191340548 -4.7936498
#> 7 pathway.fHR -0.442579892 0.347184802 -1.2747675
#> 8 pathway.fPI3K.AKT.mTOR -0.494699137 0.373006209 -1.3262491
#> 9 ps -0.473179085 0.135589101 -3.4898018
#> 10 ttc 0.062793519 0.262081632 0.2395953
#> 11 tran 0.375669348 0.377943241 0.9939835
#> 12 Log(scale) -0.459858421 0.105871866 -4.3435375
fit1$psi
#> [1] -1.13619
# experimental group
fit1$fit_aft[[2]]$fit$parest[, c("param", "beta", "sebeta", "z")]
#> param beta sebeta z
#> 1 (Intercept) 4.819297077 0.760270306 6.3389258
#> 2 swtrt 0.983747428 0.271723181 3.6204030
#> 3 agerand 0.001253863 0.009513484 0.1317985
#> 4 sex.fFemale 0.395790503 0.237599761 1.6657866
#> 5 tt_Lnum 0.066813220 0.051717213 1.2918952
#> 6 rmh_alea.c -0.714216404 0.209567497 -3.4080495
#> 7 pathway.fHR -0.301948747 0.296989947 -1.0166969
#> 8 pathway.fPI3K.AKT.mTOR -0.076210117 0.281075447 -0.2711376
#> 9 ps -0.401333182 0.114040814 -3.5192066
#> 10 ttc 0.414584865 0.377881585 1.0971291
#> 11 tran -0.615854801 0.390436760 -1.5773484
#> 12 Log(scale) -0.414571273 0.116718016 -3.5519047
fit1$psi_trt
#> [1] -0.9837474Now we fit the outcome Cox model and compare the treatment hazard ratio estimate with the reported.
fit1$fit_outcome$parest[, c("param", "beta", "sebeta", "z")]
#> param beta sebeta z
#> 1 treated -0.130373883 0.185429595 -0.7030910
#> 2 agerand 0.005589227 0.008101401 0.6899087
#> 3 sex.fFemale -0.443153389 0.192503561 -2.3020529
#> 4 tt_Lnum -0.010977575 0.037305714 -0.2942599
#> 5 rmh_alea.c 0.708929840 0.192039389 3.6915856
#> 6 pathway.fHR 0.151164769 0.299114625 0.5053740
#> 7 pathway.fPI3K.AKT.mTOR 0.200302360 0.304414444 0.6579923
c(fit1$hr, fit1$hr_CI)
#> [1] 0.8777672 0.6102972 1.2624590Finally, to ensure the uncertainty is accurately represented, the entire adjustment process and subsequent survival modeling can be bootstrapped.
```{r boot}
fit2 <- tsesimp(
data = shilong3, time = "tstop", event = "event",
treat = "bras.f", censor_time = "dcut", pd = "pd",
pd_time = "dpd", swtrt = "co", swtrt_time = "dco",
base_cov = c("agerand", "sex.f", "tt_Lnum", "rmh_alea.c",
"pathway.f"),
base2_cov = c("agerand", "sex.f", "tt_Lnum", "rmh_alea.c",
"pathway.f", "ps", "ttc", "tran"),
aft_dist = "weibull", alpha = 0.05,
recensor = TRUE, swtrt_control_only = FALSE, offset = 1,
boot = TRUE, n_boot = 1000, seed = 12345)
c(fit2$hr, fit2$hr_CI)
```