# Simple Two-Stage Estimation

    library(trtswitch)
    library(dplyr, warn.conflicts = FALSE)

# Introduction

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.

# Estimation of \(\psi\)

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.

# Estimation of Hazard Ratio

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.

# Example

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, id = "id", 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.819297083 0.760270308  6.3389258
    #> 2                   swtrt  0.983747428 0.271723182  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.301948752 0.296989948 -1.0166969
    #> 8  pathway.fPI3K.AKT.mTOR -0.076210121 0.281075448 -0.2711376
    #> 9                      ps -0.401333183 0.114040814 -3.5192066
    #> 10                    ttc  0.414584864 0.377881586  1.0971290
    #> 11                   tran -0.615854802 0.390436760 -1.5773484
    #> 12             Log(scale) -0.414571271 0.116718016 -3.5519047

    fit1$psi_trt
    #> [1] -0.9837474

Now 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.2624590

Finally, to ensure the uncertainty is accurately represented, the entire adjustment process and subsequent survival modeling can be bootstrapped.

    ```{r boot}
    fit2 <- tsesimp(
      data = shilong3, id = "id", 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)
    ```