Evaluate trtswitch::rpsftm()

This markdown is to validate trtswitch::rpsftm() with trtswitch::shilong and trtswitch::immdef data.

1 Executive Summary

• 1. One-way, with recensoring: Consistent results between trtswitch::rpsftm() and rpsftm::rpsftm() for both shilong and immdef datasets, with very similar psi estimates, hazard ratios, and p-values.
• 2. One-way, without recensoring: Consistent results between trtswitch::rpsftm() and rpsftm::rpsftm() for both shilong and immdef datasets, with very similar psi estimates, hazard ratios, and p-values.
• 3. Two-way, with recensoring: Similar psi estimates but notable differences in hazard ratios and p-values between implementations for shilong dataset, with some discrepancies in counterfactual times for treatment arm patients.
• 4. Two-way, without recensoring: Similar psi estimates but notable differences in hazard ratios and p-values between implementations for shilong dataset, with some discrepancies in counterfactual times for treatment arm patients.

Validation Summary: trtswitch::rpsftm() vs rpsftm::rpsftm()

Switching Type	Recensoring	Dataset	Psi	HR	P-value	Counterfactual Time
1-way	YES	shilong	✅	✅	✅	✅
1-way	NO	shilong	✅	✅	✅	✅
1-way	YES	immdef	✅	✅	✅	✅
1-way	NO	immdef	✅	✅	✅	✅
1-way	YES	anonymised_trial_data	✅	✅	✅	✅
1-way	NO	anonymised_trial_data	✅	✅	✅	✅
2-way	YES	shilong	✅	⚠️	❌	Control: ✅; Treatment: ❌/⚠️
2-way	NO	shilong	✅	⚠️	❌	Control: ✅; Treatment: ❌/⚠️

Legend: - ✅ = Good agreement between implementations - ❌ = Notable differences between implementations
- ⚠️ = Partial agreement or minor differences

2 Data Processing

2.1 Prepare shilong and immdef datasets

library(tidyverse)
library(dplyr)
library(survival)
library(rpsftm)
library(trtswitch) 
library(boot)
library(ggplot2) 
shilong1 <- trtswitch::shilong %>%
  arrange(bras.f, id, ady) %>%
  group_by(bras.f, id) %>%
  slice(n()) %>%
  select(-c("ps", "ttc", "tran")) %>%
  ungroup() %>%
  # two-way switch
    mutate(rx = ifelse(co, ifelse(bras.f == "MTA", dco/ady, 1 - dco/ady),
                     ifelse(bras.f == "MTA", 1, 0))) %>%
  # one-way switch
  mutate(rx1 = ifelse(bras.f == "MTA", 1, 
                     ifelse(co, 1 - dco/ady, 0)))  %>%
  mutate(xo1 = ifelse(bras.f == "MTA",  FALSE, co)) %>%
  mutate(arm = case_when(bras.f == "MTA" ~ 1, TRUE ~ 0)) %>% 
  mutate(survtime = ady, status = event, xoyrs = dco)

immdef1 <- trtswitch::immdef %>% 
  mutate(rx1 = 1-xoyrs/progyrs) %>%
  mutate(xo1 = xo) %>%
  mutate(survtime = progyrs, status = prog, arm = imm)


# Prepare internal_data dataset
anonymised_trial_data <- read_csv("processed_data.csv") %>%
   mutate(event = as.numeric(event)) %>%
  
  # one-way switch  
  mutate(rx1 = ifelse(arm == 1, 1, 
                     ifelse(xo == TRUE, 1 - dco/ady, 0))) %>%
  mutate(xo1 = ifelse(arm == 1, FALSE, xo)) %>%
 mutate(survtime = ady, status = as.integer(event), xoyrs = dco) %>%
mutate(dcut = (ady +30)) 

3 Validation Codes

3.1 Validation code/function for trtswitch::rpsftm()

# input data requires columns:
# ady: survival time
# event: status
# arm: 1/0
# rx: Proportion of duration of treatment switch  
# dcut: the relative day of administrative cutoff (if recensor = TRUE)
validate_rpsftm <- function(data, rx = "rx1", recensor = TRUE){
  fml <- as.formula(paste0("Surv(ady, event)~rand(arm,", rx, ")")) 
  if (recensor){
    Result.RPSFT = rpsftm::rpsftm(fml, data = data, censor_time = dcut, low_psi=-3,hi_psi=3)
  } else if (!recensor){
        Result.RPSFT = rpsftm::rpsftm(fml, data = data, low_psi=-3, hi_psi=3)
  }
  tstar2 = Result.RPSFT$Sstar[,1]
  event2 = Result.RPSFT$Sstar[,2]
  d_add = data.frame(data,tstar2, event2)
  Result.Cox=coxph(Surv(tstar2, event2)~arm , data = d_add, ties="exact")
  Hazard.PointEstimate = exp(Result.Cox$coefficients)
  Hazard.CI = exp(confint(Result.Cox))
  Logrank = survdiff(Surv(tstar2, event2)~arm, data = d_add) 
  Pvalue = 1 - pchisq(Logrank$chisq, 1)
  return(list("psi" = Result.RPSFT$psi, 
              "hazard"= as.numeric(Hazard.PointEstimate),
              "pvalue"= Pvalue,  
              data_outcome = d_add %>% select(id, tstar2, event2)))
}

4 Validation Results

4.1 One-way, with recensoring

4.1.1 shilong Dataset

4.1.1.1 fit1 (trtswitch::rpsftm()) with shilong dataset, recensor is TRUE, one-way switch scenario

fit1 <- trtswitch::rpsftm(
  shilong1, id = "id", time = "ady", event = "event",
  treat = "arm", rx = "rx1", censor_time = "dcut", 
  low_psi = -3, hi_psi = 3, boot = FALSE, recensor = TRUE)

# censor_time needs to be specified even if recensor = FALSE 
fit1_v <- validate_rpsftm(shilong1, rx = "rx1", recensor = TRUE)

• Compare fit1 and fit1_v: psi, hr, p-value

c(fit1$psi, fit1$hr, fit1$cox_pvalue) 

## [1] 0.48211548 1.61653092 0.00871751

c(fit1_v$psi, fit1_v$hazard, fit1_v$pvalue)

## [1] 0.4820509 1.0000598 0.9997310

The estimated causal parameter psi, the estimated hazard ratio and p-value from the Cox model applied to counterfactual unswitched survival times are very close between the result from trtswitch::rpsftm() and validation result using rpsftm::rpsftm() in one-way switch scenario.

• Compare fit1 and fit1_v: counterfactual survival time, event

compare_fit1 <- fit1$data_outcome %>%
  left_join(fit1_v$data_outcome, by = "id")
ggplot(compare_fit1, aes(x = t_star, y = tstar2)) +
  geom_point() +
  geom_abline(intercept = 0, slope = 1) + 
  labs(title = "Counterfactual survival time comparison", x = "From trtswitch::rpsftm", y = "From validation")

table(compare_fit1$event2, compare_fit1$d_star) 

##    
##       0   1
##   0  64   0
##   1   0 129

4.1.2 immdef Dataset

4.1.2.1 (trtswitch::rpsftm()) with immdef dataset, recensor = TRUE, one-way switch

# RUN TRTSWITCH RPSFTM ON IMMDEF DATASET
fit5 <- trtswitch::rpsftm(
  immdef1, id = "id", time = "progyrs", event = "prog", 
  treat = "imm", rx = "rx1", censor_time = "censyrs", 
  low_psi = -3, hi_psi = 3, boot = FALSE, recensor = TRUE)

# RUN REFERENCE IMPLEMENTATION ON IMMDEF DATASET
# First need to prepare data with correct column names for validation function
immdef_rpsftm <- immdef1 %>%
  transmute(
    id = id,
    ady = progyrs, 
    event = prog, 
    arm = imm, 
    dcut = censyrs,
    rx1 = rx1
  )

fit5_v <- validate_rpsftm(immdef_rpsftm, rx = "rx1", recensor = TRUE)

Compare RPSFTM results on immdef dataset: psi, hazard ratio, p-value

c(fit5$psi, fit5$hr, fit5$cox_pvalue) 

## [1] -0.18117698  0.76109922  0.02408344

c(fit5_v$psi, fit5_v$hazard, fit5_v$pvalue)

## [1] -0.1814538  1.0035655  0.9760373

Compare counterfactual survival times and events

compare_fit5 <- fit5$data_outcome %>%
  left_join(fit5_v$data_outcome, by = "id")

# PLOT COMPARISON OF COUNTERFACTUAL SURVIVAL TIMES
ggplot(compare_fit5, aes(x = t_star, y = tstar2)) +
  geom_point() +
  geom_abline(intercept = 0, slope = 1) + 
  labs(title = "Counterfactual survival time comparison - immdef dataset", 
       x = "From trtswitch::rpsftm", 
       y = "From validation")

# COMPARE EVENT INDICATORS
table(compare_fit5$event2, compare_fit5$d_star) 

##    
##       0   1
##   0 714   1
##   1   0 285

Check results for patients who DIDN’T switch (immdef dataset)

compare_fit5_nonswt <- compare_fit5 %>% 
  left_join(immdef1 %>% select(id, rx1, progyrs), by = "id") %>% 
  filter((rx1 == 1 & imm == 1) | (rx1 == 0 & imm == 0))
head(compare_fit5_nonswt, n = 5)

##   id    t_star d_star treated imm    tstar2 event2 rx1   progyrs
## 1  1 3.0000000      0       1   1 2.5021704      0   1 3.0000000
## 2  3 1.7378377      1       1   1 1.4494553      1   1 1.7378377
## 3  4 2.1662905      1       1   1 1.8068093      1   1 2.1662905
## 4  7 2.1894703      1       0   0 2.1894703      1   0 2.1894703
## 5  8 0.9226239      1       1   1 0.7695208      1   1 0.9226239

Check results for patients who DID switch (immdef dataset)

compare_fit5_swt <- compare_fit5 %>% 
  left_join(immdef1 %>% select(id, rx1, progyrs, xoyrs), by = "id") %>% 
  filter(rx1 > 0 & imm == 0)
head(compare_fit5_swt, n = 5)

##   id   t_star d_star treated imm   tstar2 event2       rx1  progyrs     xoyrs
## 1  2 2.502863      0       0   0 2.502170      0 0.1157343 3.000000 2.6527972
## 2  5 2.502863      0       0   0 2.502170      0 0.2643466 2.884646 2.1220999
## 3  6 2.502863      0       0   0 2.502170      0 0.8142027 3.000000 0.5573920
## 4 13 2.502863      0       0   0 2.502170      0 0.4089745 3.000000 1.7730765
## 5 19 2.023826      1       0   0 2.023395      1 0.8006016 2.333397 0.4652756

Mathematical Check for immdef dataset: Verify the relationship between counterfactual time and psi parameter.

if(nrow(compare_fit5_swt) > 0) {
  # Take first patient who switched as an example
  compare_fit5_example <- compare_fit5_swt %>% slice(1)
  
  # Calculate the scaling factor applied to the post-switch period
  scaling_factor <- (compare_fit5_example$t_star - compare_fit5_example$xoyrs) / 
                   (compare_fit5_example$progyrs - compare_fit5_example$xoyrs)
  
  cat("Scaling factor for patient", compare_fit5_example$id, ":", scaling_factor, "\n")
  cat("Expected value exp(psi):", exp(fit5$psi), "\n")
  cat("Difference:", abs(scaling_factor - exp(fit5$psi)), "\n")
} else {
  cat("No patients switched in the immdef dataset")
}

## Scaling factor for patient 2 : -0.4318344 
## Expected value exp(psi): 0.8342877 
## Difference: 1.266122

4.1.3 Anonymised Trial Dataset

4.1.3.1 (trtswitch::rpsftm()) with internal_data dataset, recensor = TRUE, one-way switch

# RUN TRTSWITCH RPSFTM ON INTERNAL_DATA DATASET
fit7 <- trtswitch::rpsftm(
  anonymised_trial_data, id = "id", time = "survtime", event = "status", 
  treat = "arm", rx = "rx1", censor_time = "dcut", 
  low_psi = -3, hi_psi = 3, boot = FALSE, recensor = TRUE)

# RUN REFERENCE IMPLEMENTATION ON INTERNAL_DATA DATASET
#First need to prepare data with correct column names for validation function
internal_data_rpsftm <- anonymised_trial_data %>%
  transmute(
   id = id,
   ady = ady, 
   event = event, 
  arm = arm, 
   dcut = dcut,
   rx1 = rx1
  )

fit7_v <- validate_rpsftm(internal_data_rpsftm, rx = "rx1", recensor = TRUE)

Compare RPSFTM results on internal_data dataset: psi, hazard ratio, p-value

c(fit7$psi, fit7$hr, fit7$cox_pvalue) 

## [1] -0.46988825  0.69457888  0.05584737

c(fit7_v$psi, fit7_v$hazard, fit7_v$pvalue)

## [1] -0.4700738  1.0000633  0.9997328

Compare counterfactual survival times and events

compare_fit7 <- fit7$data_outcome %>%
  left_join(fit7_v$data_outcome, by = "id")

# PLOT COMPARISON OF COUNTERFACTUAL SURVIVAL TIMES
ggplot(compare_fit7, aes(x = t_star, y = tstar2)) +
  geom_point() +
  geom_abline(intercept = 0, slope = 1) + 
  labs(title = "Counterfactual survival time comparison - internal_data dataset", 
       x = "From trtswitch::rpsftm", 
      y = "From validation")

# COMPARE EVENT INDICATORS
table(compare_fit7$event2, compare_fit7$d_star) 

##    
##       0   1
##   0  67   0
##   1   0 113

Check results for patients who DIDN’T switch (internal_data dataset)

compare_fit7_nonswt <- compare_fit7 %>% 
  left_join(anonymised_trial_data %>% select(id, rx1, ady), by = "id") %>% 
  filter((rx1 == 1 & arm == 1) | (rx1 == 0 & arm == 0))
head(compare_fit7_nonswt, n = 5)

##      id    t_star d_star treated arm    tstar2 event2 rx1       ady
## 1 10067 26.670376      1       1   1 16.667816      1   1 26.670376
## 2 10075 23.769947      0       1   1 14.855175      0   1 23.769947
## 3 10109 40.926910      0       1   1 25.577525      0   1 40.926910
## 4 10072  3.373891      1       1   1  2.108534      1   1  3.373891
## 5 10031 24.297373      1       1   1 15.184793      1   1 24.297373

Check results for patients who DID switch (internal_data dataset)

compare_fit7_swt <- compare_fit7 %>% 
  left_join(anonymised_trial_data %>% select(id, rx1, ady, xoyrs), by = "id") %>% 
 filter(rx1 > 0 & arm == 0)
head(compare_fit7_swt, n = 5)

##      id   t_star d_star treated arm   tstar2 event2       rx1      ady   xoyrs
## 1 10191 21.72202      0       0   0 21.72043      0 0.5108907 26.86863 13.1417
## 2 10032 18.16075      1       0   0 18.15894      1 0.6496853 24.00899  8.4107
## 3 10015 13.30396      1       0   0 13.30310      1 0.4606494 16.08138  8.6735
## 4 10022 33.07687      0       0   0 33.07414      0 0.5616927 41.90097 18.3655
## 5 10126 13.14653      1       0   0 13.14583      1 0.3926068 15.41572  9.3634

Mathematical Check for internal_data dataset: Verify the relationship between counterfactual time and psi parameter.

if(nrow(compare_fit7_swt) > 0) {
  # Take first patient who switched as an example
  compare_fit7_example <- compare_fit7_swt %>% slice(1)
  
  # Calculate the scaling factor applied to the post-switch period
  scaling_factor <- (compare_fit7_example$t_star - compare_fit7_example$xoyrs) / 
                  (compare_fit7_example$ady - compare_fit7_example$xoyrs)
  
 cat("Scaling factor for patient", compare_fit7_example$id, ":", scaling_factor, "\n")
cat("Expected value exp(psi):", exp(fit7$psi), "\n")
 cat("Difference:", abs(scaling_factor - exp(fit7$psi)), "\n")
} else {
 cat("No patients switched in the anonymised_trial_data dataset")
}

## Scaling factor for patient 10191 : 0.6250721 
## Expected value exp(psi): 0.6250721 
## Difference: 1.110223e-16

4.2 One-way, without recensoring

4.2.1 shilong Dataset

4.2.1.1 fit2 (trtswitch::rpsftm()) with shilong dataset, recensor = FALSE, one-way switch scenario

fit2 <- trtswitch::rpsftm(
  shilong1, id = "id", time = "ady", event = "event",
  treat = "arm", rx = "rx1", censor_time = "dcut", 
  low_psi = -3, hi_psi = 3, boot = FALSE, recensor = FALSE)
fit2_v<- validate_rpsftm(shilong1, rx = "rx1", recensor = FALSE) 

## Warning in rpsftm::rpsftm(fml, data = data, low_psi = -3, hi_psi = 3):
## Evaluation of a limit of the Confidence Interval failed.  It is set to NA

• Compare fit2 and fit2_v: psi, hr, pvalue

c(fit2$psi, fit2$hr, fit2$cox_pvalue) 

## [1] 0.578078205 1.700890346 0.003967654

c(fit2_v$psi, fit2_v$hazard, fit2_v$pvalue)

## [1] 0.5780950 0.9999659 0.9998461

Similar as model with recensor, results of model without recensoring are also very comparable between trtswitch::rpsftm() and rpsftm::rpsftm() for one-way switch. Improvement suggestion for trtswitch::rpsftm(): parameter censor_time is required even recensor is FALSE.

• Compare fit2 and fit2_v: counterfactual survival time, event

compare_fit2 <- fit2$data_outcome %>%
  left_join(fit2_v$data_outcome, by = "id")
ggplot(compare_fit2, aes(x = t_star, y = tstar2)) +
  geom_point() +
  geom_abline(intercept = 0, slope = 1) + 
  labs(title = "Counterfactual survival time comparison", x = "From trtswitch::rpsftm", y = "From validation")

table(compare_fit2$event2, compare_fit2$d_star) 

##    
##       0   1
##   0  63   0
##   1   0 130

4.2.1.2 fit2 vs. fit2_v: Compare counterfactual survival time to original survival time for subjects who didn’t switch.

compare_fit2_nonswt <- compare_fit2 %>% left_join(shilong1 %>% select(id, rx1, ady), by = "id") %>% filter(rx1 == 1)
head(compare_fit2_nonswt, n = 5)

##   id t_star d_star treated arm   tstar2 event2 rx1 ady
## 1  2     64      1       1   1 114.0889      1   1  64
## 2  4    156      1       1   1 278.0917      1   1 156
## 3  5    436      0       1   1 777.2307      0   1 436
## 4  6    500      1       1   1 891.3197      1   1 500
## 5 12    139      1       1   1 247.7869      1   1 139

In one-way switch scenario, counterfactual survival time (t_star from trtswitch::rpsftm() and tstar2 from rpsftm::rpsftm()) for subjects who didn’t switch (treatment arm in the shilong dataset) is the same as the original survival time ady, which is expected.

4.2.1.3 fit2 vs. fit2_v: Compare counterfactual survival time to original survival time for subjects who switched.

compare_fit2_swt <- compare_fit2 %>% left_join(shilong1 %>% select(id, rx1, ady, dco), by = "id") %>% filter(rx1 < 1)
head(compare_fit2_swt, n = 5)

##   id   t_star d_star treated arm   tstar2 event2       rx1 ady dco
## 1  1 234.2175      1       0   0 234.2209      1 0.7862069 145  31
## 2  3 412.2175      1       0   0 412.2223      1 0.5574913 287 127
## 3  8 815.2611      0       0   0 815.2738      0 0.8701031 485  63
## 4  9 560.3480      1       0   0 560.3564      1 0.8157895 342  63
## 5 10  37.0000      1       0   0  37.0000      1 0.0000000  37  NA

In one-way switch scenario, counterfactual survival time t_star from trtswitch::rpsftm() and tstar2 from rpsftm::rpsftm() are similar. Both of the counterfactual survival times for subjects who switched (control arm in the shilong dataset) are larger than the original survival time ady. Take subject id = 1 as an example to check how counterfactual survival time is related to the estimated psi from trtswitch::rpsftm().

compare_fit2_1 <- compare_fit2_swt %>% filter(id == "1") 
(compare_fit2_1$t_star- compare_fit2_1$dco)/(compare_fit2_1$ady - compare_fit2_1$dco)

## [1] 1.782609

equals to exp(-fit2$psi)

exp(fit2$psi) 

## [1] 1.782609

4.2.2 immdef Dataset

4.2.2.1 fit6 (trtswitch::rpsftm()) with immdef dataset, recensor = FALSE, one-way switch

# RUN TRTSWITCH RPSFTM ON IMMDEF DATASET
fit6 <- trtswitch::rpsftm(
  immdef1, id = "id", time = "progyrs", event = "prog", 
  treat = "imm", rx = "rx1", censor_time = "censyrs", 
  low_psi = -3, hi_psi = 3, boot = FALSE, recensor = FALSE)

# RUN REFERENCE IMPLEMENTATION ON IMMDEF DATASET
fit6_v <- validate_rpsftm(immdef_rpsftm, rx = "rx1", recensor = FALSE)

Compare RPSFTM results on immdef dataset: psi, hazard ratio, p-value

c(fit6$psi, fit6$hr, fit6$cox_pvalue) 

## [1] -0.18505941  0.76677041  0.01972649

c(fit6_v$psi, fit6_v$hazard, fit6_v$pvalue)

## [1] -0.1848298  0.9999392  0.9995780

Compare counterfactual survival times and events

compare_fit6 <- fit6$data_outcome %>%
  left_join(fit6_v$data_outcome, by = "id")

# PLOT COMPARISON OF COUNTERFACTUAL SURVIVAL TIMES
ggplot(compare_fit6, aes(x = t_star, y = tstar2)) +
  geom_point() +
  geom_abline(intercept = 0, slope = 1) + 
  labs(title = "Counterfactual survival time comparison - immdef dataset", 
       x = "From trtswitch::rpsftm", 
       y = "From validation")

# COMPARE EVENT INDICATORS
table(compare_fit6$event2, compare_fit6$d_star) 

##    
##       0   1
##   0 688   0
##   1   0 312

Check results for patients who DIDN’T switch (immdef dataset)

compare_fit6_nonswt <- compare_fit6 %>% 
  left_join(immdef1 %>% select(id, rx1, progyrs), by = "id") %>% 
  filter((rx1 == 1 & imm == 1) | (rx1 == 0 & imm == 0))
head(compare_fit6_nonswt, n = 5)

##   id    t_star d_star treated imm    tstar2 event2 rx1   progyrs
## 1  1 3.0000000      0       1   1 2.4937373      0   1 3.0000000
## 2  3 1.7378377      1       1   1 1.4445702      1   1 1.7378377
## 3  4 2.1662905      1       1   1 1.8007198      1   1 2.1662905
## 4  7 2.1894703      1       0   0 2.1894703      1   0 2.1894703
## 5  8 0.9226239      1       1   1 0.7669272      1   1 0.9226239

Check results for patients who DID switch (immdef dataset)

compare_fit6_swt <- compare_fit6 %>% 
  left_join(immdef1 %>% select(id, rx1, progyrs, xoyrs), by = "id") %>% 
  filter(rx1 > 0 & imm == 0)
head(compare_fit6_swt, n = 5)

##   id   t_star d_star treated imm   tstar2 event2       rx1  progyrs     xoyrs
## 1  2 2.941342      0       0   0 2.941408      0 0.1157343 3.000000 2.6527972
## 2  5 2.755818      1       0   0 2.755963      1 0.2643466 2.884646 2.1220999
## 3  6 2.587333      0       0   0 2.587800      0 0.8142027 3.000000 0.5573920
## 4 13 2.792717      0       0   0 2.792951      0 0.4089745 3.000000 1.7730765
## 5 19 2.017787      1       0   0 2.018144      1 0.8006016 2.333397 0.4652756

Mathematical Check for immdef dataset: Verify the relationship between counterfactual time and psi parameter.

if(nrow(compare_fit6_swt) > 0) {
  # Take first patient who switched as an example
  compare_fit6_example <- compare_fit6_swt %>% slice(1)
  
  # Calculate the scaling factor applied to the post-switch period
  scaling_factor <- (compare_fit6_example$t_star - compare_fit6_example$xoyrs) / 
                   (compare_fit6_example$progyrs - compare_fit6_example$xoyrs)
  
  cat("Scaling factor for patient", compare_fit6_example$id, ":", scaling_factor, "\n")
  cat("Expected value exp(psi):", exp(fit6$psi), "\n")
  cat("Difference:", abs(scaling_factor - exp(fit6$psi)), "\n")
} else {
  cat("No patients switched in the immdef dataset")
}

## Scaling factor for patient 2 : 0.8310549 
## Expected value exp(psi): 0.8310549 
## Difference: 0

4.2.3 anonymised_trial_data Dataset

4.2.3.1 fit8 (trtswitch::rpsftm()) with internal_data dataset, recensor = FALSE, one-way switch

# RUN TRTSWITCH RPSFTM ON INTERNAL_DATA DATASET
fit8 <- trtswitch::rpsftm(
  anonymised_trial_data, id = "id", time = "survtime", event = "status", 
  treat = "arm", rx = "rx1", censor_time = "ady", 
  low_psi = -3, hi_psi = 3, boot = FALSE, recensor = FALSE)

# RUN REFERENCE IMPLEMENTATION ON INTERNAL_DATA DATASET
fit8_v <- validate_rpsftm(internal_data_rpsftm, rx = "rx1", recensor = FALSE)

Compare RPSFTM results on internal_data dataset: psi, hazard ratio, p-value

c(fit8$psi, fit8$hr, fit8$cox_pvalue) 

## [1] -0.47224404  0.69359797  0.05494223

c(fit8_v$psi, fit8_v$hazard, fit8_v$pvalue)

## [1] -0.4700738  1.0000633  0.9997328

Compare counterfactual survival times and events

compare_fit8 <- fit8$data_outcome %>%
 left_join(fit8_v$data_outcome, by = "id")

# PLOT COMPARISON OF COUNTERFACTUAL SURVIVAL TIMES
ggplot(compare_fit8, aes(x = t_star, y = tstar2)) +
 geom_point() +
 geom_abline(intercept = 0, slope = 1) + 
 labs(title = "Counterfactual survival time comparison - internal_data dataset", 
 x = "From trtswitch::rpsftm", 
 y = "From validation")

# COMPARE EVENT INDICATORS
table(compare_fit8$event2, compare_fit8$d_star) 

##    
##       0   1
##   0  67   0
##   1   0 113

Check results for patients who DIDN’T switch (internal_data dataset)

compare_fit8_nonswt <- compare_fit8 %>% 
left_join(anonymised_trial_data %>% select(id, rx1, ady), by = "id") %>% 
filter((rx1 == 1 & arm == 1) | (rx1 == 0 & arm == 0))
head(compare_fit8_nonswt, n = 5)

##      id    t_star d_star treated arm    tstar2 event2 rx1       ady
## 1 10067 26.670376      1       1   1 16.667816      1   1 26.670376
## 2 10075 23.769947      0       1   1 14.855175      0   1 23.769947
## 3 10109 40.926910      0       1   1 25.577525      0   1 40.926910
## 4 10072  3.373891      1       1   1  2.108534      1   1  3.373891
## 5 10031 24.297373      1       1   1 15.184793      1   1 24.297373

Check results for patients who DID switch (internal_data dataset)

compare_fit8_swt <- compare_fit8 %>% 
left_join(anonymised_trial_data %>% select(id, rx1, ady, xoyrs), by = "id") %>% 
filter(rx1 > 0 & arm == 0)
head(compare_fit8_swt, n = 5)

##      id   t_star d_star treated arm   tstar2 event2       rx1      ady   xoyrs
## 1 10191 21.70183      0       0   0 21.72043      0 0.5108907 26.86863 13.1417
## 2 10032 18.13781      1       0   0 18.15894      1 0.6496853 24.00899  8.4107
## 3 10015 13.29306      1       0   0 13.30310      1 0.4606494 16.08138  8.6735
## 4 10022 33.04225      0       0   0 33.07414      0 0.5616927 41.90097 18.3655
## 5 10126 13.13763      1       0   0 13.14583      1 0.3926068 15.41572  9.3634

Mathematical Check for anonymised_trial_data: Verify the relationship between counterfactual time and psi parameter.

if(nrow(compare_fit8_swt) > 0) {
  #Take first patient who switched as an example
compare_fit8_example <- compare_fit8_swt %>% slice(1)
  
  # Calculate the scaling factor applied to the post-switch period
scaling_factor <- (compare_fit8_example$t_star - compare_fit8_example$xoyrs) / 
                  (compare_fit8_example$ady - compare_fit8_example$xoyrs)
  
cat("Scaling factor for patient", compare_fit8_example$id, ":", scaling_factor, "\n")
cat("Expected value exp(psi):", exp(fit8$psi), "\n")
cat("Difference:", abs(scaling_factor - exp(fit8$psi)), "\n")
} else {
  cat("No patients switched in the internal_data dataset")
}

## Scaling factor for patient 10191 : 0.6236013 
## Expected value exp(psi): 0.6236013 
## Difference: 1.110223e-16

4.3 Two-way, with recensoring

4.3.1 Shilong Dataset

4.3.1.1 fit3 (trtswitch::rpsftm()) with shilong dataset, recensor = TRUE, two-way switch

fit3 <- trtswitch::rpsftm(
  shilong1, id = "id", time = "tstop", event = "event",
  treat = "arm", rx = "rx", censor_time = "dcut",
  low_psi = -3, hi_psi = 3, boot = FALSE, recensor = TRUE)
fit3_v <- validate_rpsftm(shilong1, rx = "rx", recensor = TRUE)

• Compare fit3 and fit3_v: psi, hr, pvalue

c(fit3$psi, fit3$hr, fit3$cox_pvalue) 

## [1] 1.007842e+00 2.721078e+00 3.648658e-06

c(fit3_v$psi, fit3_v$hazard, fit3_v$pvalue)

## [1] 1.0079860 0.9953629 0.9798349

In the two-way switch scenario, the estimated causal parameter psi = 1.0078423 from trtswitch::rpsftm() is close to psi = 1.007986 from rpsftm::rpsftm(). However, the estimated hazard ratio and p-value from the Cox model applied to counterfactual unswitched survival times are very different between fit3 and fit3_v.

• Compare fit3 and fit3_v: counterfactual survival time

compare_fit3 <- fit3$data_outcome %>%
  left_join(fit3_v$data_outcome, by = "id") %>%
  select(id, t_star, arm, tstar2, d_star, event2) %>% left_join(shilong1 %>% select(id, ady, dco, rx, status))

## Joining with `by = join_by(id)`

ggplot(compare_fit3, aes(x = t_star, y = tstar2, color = factor(arm))) +
  geom_point() +
  geom_abline(intercept = 0, slope = 1) + 
  labs(title = "Counterfactual survival time comparison", x = "From trtswitch::rpsftm", y = "From validation")

4.3.1.2 Treatment arm patients: Counterfactual survival time for switchers vs non-switchers

# Filter for treatment arm patients only and add switching status
compare_fit3_treatment <- compare_fit3 %>% 
  filter(arm == 1) %>%
  mutate(switching_status = ifelse(rx == 1, "Non-switcher", "Switcher"))

ggplot(compare_fit3_treatment, aes(x = t_star, y = tstar2, color = switching_status)) +
  geom_point(size = 3, alpha = 0.7) +
  geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "black") + 
  labs(title = "Counterfactual survival time comparison - Treatment arm patients only", 
       subtitle = "Two-way switching scenario with recensoring",
       x = "From trtswitch::rpsftm", 
       y = "From rpsftm::rpsftm (validation)",
       color = "Switching Status") +
  scale_color_manual(values = c("Switcher" = "#E31A1C", "Non-switcher" = "#1F78B4")) +
  theme_minimal() +
  theme(legend.position = "bottom",
        plot.title = element_text(size = 14, face = "bold"),
        plot.subtitle = element_text(size = 12))

4.3.1.3 fit3 vs. fit3_v: Counterfactual survival time for subjects who didn’t switch.

compare_fit3_nonswt <- compare_fit3 %>% filter((arm == 0 & rx == 0) | (arm == 1 & rx == 1))
ggplot(compare_fit3_nonswt, aes(x = t_star, y = tstar2, color = factor(arm))) +
  geom_point() +
  geom_abline(intercept = 0, slope = 1) + 
  labs(title = "Counterfactual survival time comparison for subjects who didn't switch", x = "From trtswitch::rpsftm", y = "From validation")

compare_fit3_nonswt %>% mutate(delta = t_star - tstar2) %>% filter(delta != 0)                                                                                            

##     id   t_star arm     tstar2 d_star event2 ady dco rx status      delta
## 1    2  64.0000   1  175.36492      1      1  64  NA  1      1 -111.36492
## 2   19  73.0000   1  200.02562      0      0  73  NA  1      0 -127.02562
## 3   25  62.0000   1  169.88477      1      1  62  NA  1      1 -107.88477
## 4   26 329.0000   1  901.48531      1      1 329  NA  1      1 -572.48531
## 5   33 167.0000   1  457.59285      1      1 167  NA  1      1 -290.59285
## 6   37  76.0000   1  208.24585      1      1  76  NA  1      1 -132.24585
## 7   51 209.0000   1  572.67608      0      0 209  NA  1      0 -363.67608
## 8   52  65.0000   1  178.10500      1      1  65  NA  1      1 -113.10500
## 9   53 246.0000   1  674.05892      0      0 246  NA  1      0 -428.05892
## 10  55 146.0000   1  400.05123      0      0 146  NA  1      0 -254.05123
## 11  56 210.6083   1  577.00000      0      0 286  NA  1      0 -366.39171
## 12  59 177.0000   1  484.99361      0      0 177  NA  1      0 -307.99361
## 13  60 210.2433   1  576.00000      0      0 280  NA  1      1 -365.75671
## 14  62 109.0000   1  298.66838      1      1 109  NA  1      1 -189.66838
## 15  64  43.0000   1  117.82331      1      1  43  NA  1      1  -74.82331
## 16  66 146.0000   1  400.05123      0      0 146  NA  1      0 -254.05123
## 17  69  76.0000   1  208.24585      1      1  76  NA  1      1 -132.24585
## 18  70 120.0000   1  328.80923      1      1 120  NA  1      1 -208.80923
## 19  72 392.0161   1 1074.00000      0      0 510  NA  1      1 -681.98387
## 20  73  90.0000   1  246.60692      1      1  90  NA  1      1 -156.60692
## 21  77  28.0000   1   76.72215      1      1  28  NA  1      1  -48.72215
## 22  79 162.0000   1  443.89246      0      0 162  NA  1      0 -281.89246
## 23  83 198.0000   1  542.53523      1      1 198  NA  1      1 -344.53523
## 24  87 252.0000   1  690.49938      0      0 252  NA  1      0 -438.49938
## 25  89  27.0000   1   73.98208      1      1  27  NA  1      1  -46.98208
## 26  92  30.0000   1   82.20231      1      1  30  NA  1      1  -52.20231
## 27  93 150.0000   1  411.01154      1      1 150  NA  1      1 -261.01154
## 28  94 105.0000   1  287.70808      1      1 105  NA  1      1 -182.70808
## 29  95 139.0000   1  380.87069      1      1 139  NA  1      1 -241.87069
## 30  96 114.0000   1  312.36877      1      1 114  NA  1      1 -198.36877
## 31  97  98.0000   1  268.52754      1      1  98  NA  1      1 -170.52754
## 32  99 329.9652   1  904.00000      0      0 386  NA  1      0 -574.03484
## 33 101 180.0000   1  493.21385      1      1 180  NA  1      1 -313.21385
## 34 103  97.0000   1  265.78746      0      0  97  NA  1      0 -168.78746
## 35 104  49.0000   1  134.26377      1      1  49  NA  1      1  -85.26377
## 36 105 223.0000   1  611.03715      1      1 223  NA  1      1 -388.03715
## 37 106 103.0000   1  282.22792      1      1 103  NA  1      1 -179.22792
## 38 108 178.0000   1  487.73369      0      0 178  NA  1      0 -309.73369
## 39 110 202.0000   1  553.49554      0      0 202  NA  1      0 -351.49554
## 40 111 132.0000   1  361.69015      1      1 132  NA  1      1 -229.69015
## 41 113 324.0000   1  887.78492      1      1 324  NA  1      1 -563.78492
## 42 118  92.0000   1  252.08708      1      1  92  NA  1      1 -160.08708
## 43 120 314.9999   1  863.00000      0      0 315  NA  1      1 -548.00007
## 44 125 102.0000   1  279.48785      1      1 102  NA  1      1 -177.48785
## 45 126 278.0000   1  761.74138      1      1 278  NA  1      1 -483.74138
## 46 127  58.0000   1  158.92446      1      1  58  NA  1      1 -100.92446
## 47 128 368.2908   1 1009.00000      0      0 389  NA  1      1 -640.70924
## 48 129 170.0000   1  465.81308      1      1 170  NA  1      1 -295.81308
## 49 130  61.0000   1  167.14469      1      1  61  NA  1      1 -106.14469
## 50 134 196.0000   1  537.05508      1      1 196  NA  1      1 -341.05508
## 51 136  55.0000   1  150.70423      1      1  55  NA  1      1  -95.70423
## 52 139  24.0000   1   65.76185      1      1  24  NA  1      1  -41.76185
## 53 140 305.0000   1  835.72346      0      0 305  NA  1      0 -530.72346
## 54 143  55.0000   1  150.70423      0      0  55  NA  1      0  -95.70423
## 55 144 102.0000   1  279.48785      0      0 102  NA  1      0 -177.48785
## 56 146 130.0000   1  356.21000      1      1 130  NA  1      1 -226.21000
## 57 147  83.0000   1  227.42638      1      1  83  NA  1      1 -144.42638
## 58 148 210.0000   1  575.41615      1      1 210  NA  1      1 -365.41615
## 59 149 218.0000   1  597.33677      0      0 218  NA  1      0 -379.33677
## 60 151 230.6836   1  632.00000      0      0 573  NA  1      0 -401.31639
## 61 152  20.0000   1   54.80154      1      1  20  NA  1      1  -34.80154
## 62 156  93.0000   1  254.82715      1      1  93  NA  1      1 -161.82715
## 63 158  39.0000   1  106.86300      1      1  39  NA  1      1  -67.86300
## 64 159  47.0000   1  128.78362      1      1  47  NA  1      1  -81.78362
## 65 161 296.0000   1  811.06277      1      1 296  NA  1      1 -515.06277
## 66 175  66.0000   1  180.84508      1      1  66  NA  1      1 -114.84508
## 67 183 142.0000   1  389.09092      1      1 142  NA  1      1 -247.09092
## 68 184 149.2873   1  409.00000      0      0 261  NA  1      0 -259.71266
## 69 186 205.1332   1  562.00000      0      0 238  NA  1      0 -356.86679
## 70 187  77.0000   1  210.98592      1      1  77  NA  1      1 -133.98592
## 71 189 105.0000   1  287.70808      1      1 105  NA  1      1 -182.70808
## 72 192 214.6234   1  588.00000      0      0 553  NA  1      0 -373.37664
## 73 193  46.0000   1  126.04354      1      1  46  NA  1      1  -80.04354
## 74 194 193.0000   1  528.83485      1      1 193  NA  1      1 -335.83485
## 75 196 239.0000   1  654.87838      0      0 239  NA  1      0 -415.87838

There are a couple of dots that are not on the line. They are all from treatment arm (arm = 1). Counterfactual survival time (t_star) from trtswitch::rpsftm() for these subjects shrunk, while results (tstar2) from rpsftm::rpsftm() are the same as the original survival time (ady). We notice that d_star and event2 stay consistent, whereas the tstar2 and t_star are different. tstar2 is the same as “ady” (original value), whereas t_star is being updated based on psi.

4.3.1.4 A closer look at counterfactual survival time for subjects who switch

head(compare_fit3 %>% filter(arm == 1 & rx < 1), n = 5)

##   id    t_star arm   tstar2 d_star event2 ady dco        rx status
## 1  4  75.99072   1 208.2023      1      1 156  30 0.1923077      1
## 2  5 213.75200   1 585.6466      0      0 436  86 0.1972477      0
## 3  6 263.14713   1 720.9898      1      1 500 127 0.2540000      1
## 4 12  92.64542   1 253.8451      1      1 139  66 0.4748201      1
## 5 23 215.34644   1 590.0295      1      1 376 123 0.3271277      1

In treatment arm, counterfactual survival time (tstar2) doesn’t get updated from rpsftm::rpsftm(), while t_star from trtswitch::rpsftm() is smaller than ady.

Take subject id = 5 from treatment arm as an example to check how counterfactual survival time is related to the estimated psi from trtswitch::rpsftm().

compare_fit3_5 <- compare_fit3 %>% filter(id == "5")
compare_fit3_5

##   id  t_star arm   tstar2 d_star event2 ady dco        rx status
## 1  5 213.752   1 585.6466      0      0 436  86 0.1972477      0

(compare_fit3_5$t_star- compare_fit3_5$dco)/(compare_fit3_5$ady - compare_fit3_5$dco)

## [1] 0.3650057

which equals to exp(-fit3$psi). This is consistent with our expectations.

exp(-fit3$psi)

## [1] 0.3650057

In control arm, both counterfactual survival time (tstar2) from rpsftm::rpsftm() and t_star from trtswitch::rpsftm() are larger than ady.

head(compare_fit3 %>% filter(arm == 0 & rx > 0), n = 5)

##   id   t_star arm    tstar2 d_star event2 ady dco        rx status
## 1  1 343.3239   0 343.36877      1      1 145  31 0.7862069      1
## 2  3 565.3493   0 565.41231      1      1 287 127 0.5574913      1
## 3  8 598.0000   0 598.00000      0      0 485  63 0.8701031      0
## 4  9 827.3716   0 827.48146      1      1 342  63 0.8157895      1
## 5 11  53.4381   0  53.44046      1      1  43  37 0.1395349      1

Take subject id = 9 from control arm as an example to check how counterfactual survival time is related to the estimated psi from trtswitch::rpsftm().

compare_fit3_9 <- compare_fit3 %>% filter(id == "9")
compare_fit3_9

##   id   t_star arm   tstar2 d_star event2 ady dco        rx status
## 1  9 827.3716   0 827.4815      1      1 342  63 0.8157895      1

(compare_fit3_9$t_star- compare_fit3_9$dco)/(compare_fit3_9$ady - compare_fit3_9$dco)

## [1] 2.739683

equals to exp(fit3$psi)

exp(fit3$psi)

## [1] 2.739683

(compare_fit3_9$tstar2- compare_fit3_9$dco)/(compare_fit3_9$ady - compare_fit3_9$dco)

## [1] 2.740077

equals to exp(fit3_v$psi)

exp(fit3_v$psi)

## [1] 2.740077

4.4 Two-way, without recensoring

4.4.1 shilong Dataset

4.4.1.1 fit4 (trtswitch::rpsftm()) with shilong dataset, recensor = FALSE, two-way switch

fit4 <- trtswitch::rpsftm(
  shilong1, id = "id", time = "ady", event = "event",
  treat = "arm", rx = "rx", censor_time = "dcut",
  low_psi = -3, hi_psi = 3, boot = FALSE, recensor = FALSE)
fit4_v <- validate_rpsftm(shilong1, rx = "rx", recensor = FALSE)

## Warning in rpsftm::rpsftm(fml, data = data, low_psi = -3, hi_psi = 3):
## Evaluation of a limit of the Confidence Interval failed.  It is set to NA

• Compare fit4 and fit4_v: psi, hr, pvalue

c(fit4$psi, fit4$hr, fit4$cox_pvalue) 

## [1] 1.119231e+00 2.717309e+00 7.257647e-07

c(fit4_v$psi, fit4_v$hazard, fit4_v$pvalue)

## [1] 1.1188999 1.0000278 0.9998749

• Compare fit4 and fit4_v: counterfactual survival time

compare_fit4 <- fit4$data_outcome %>%
  left_join(fit4_v$data_outcome, by = "id") %>%
  select(id, t_star, arm, tstar2) %>% left_join(shilong1 %>% select(id, ady, dco, rx))

## Joining with `by = join_by(id)`

ggplot(compare_fit4, aes(x = t_star, y = tstar2, color = factor(arm))) +
  geom_point() +
  geom_abline(intercept = 0, slope = 1) + 
  labs(title = "Counterfactual survival time comparison", x = "From trtswitch::rpsftm", y = "From validation")

5 Unadjusted CoxPH Model Results

5.1 shilong dataset

unadj_fit <- coxph(Surv(ady, event) ~ arm, data = shilong1)
unadj_fit

## Call:
## coxph(formula = Surv(ady, event) ~ arm, data = shilong1)
## 
##       coef exp(coef) se(coef)     z     p
## arm 0.2349    1.2648   0.1777 1.322 0.186
## 
## Likelihood ratio test=1.75  on 1 df, p=0.1859
## n= 193, number of events= 130

5.2 immdef dataset

unadj_fit_immdef <- coxph(Surv(progyrs, prog) ~ imm, data = immdef1)
unadj_fit_immdef

## Call:
## coxph(formula = Surv(progyrs, prog) ~ imm, data = immdef1)
## 
##        coef exp(coef) se(coef)     z      p
## imm -0.2171    0.8048   0.1137 -1.91 0.0561
## 
## Likelihood ratio test=3.66  on 1 df, p=0.05562
## n= 1000, number of events= 312

5.3 anonymised_trial_data dataset

unadj_fit_internal <- coxph(Surv(ady, event) ~ arm, data = anonymised_trial_data)
unadj_fit_internal

## Call:
## coxph(formula = Surv(ady, event) ~ arm, data = anonymised_trial_data)
## 
##        coef exp(coef) se(coef)      z     p
## arm -0.2489    0.7797   0.1887 -1.319 0.187
## 
## Likelihood ratio test=1.74  on 1 df, p=0.1873
## n= 180, number of events= 113

6 Detailed Summary

6.1 Summarize the findings for comparing rpsftm from trtswitch::rpsftm() and rpsftm::rpsftm()

1. The unadjusted CoxPH model gives a HR > 1, indicating a worse treatment effect. Both One-way switch and Two-way switch RPSFTM model further enlarge the HR. With one-way switch scenario, two functions trtswitch::rpsftm() and rpsftm::rpsftm() give very similar hazard ratio. For two-way switch scenario, the trend is similar, but the numbers are quite different regarding hazard ratio. Here is a summary table of HRs from different models.

##                                 Model Package::function  Dataset        HR
## 1                          Unadjusted   survival::coxph shilong1 1.2647965
## 2     One-way switch, recensor = TRUE trtswitch::rpsftm shilong1 1.6165309
## 3     One-way switch, recensor = TRUE    rpsftm::rpsftm shilong1 1.0000598
## 4    One-way switch, recensor = FALSE trtswitch::rpsftm shilong1 1.7008903
## 5    One-way switch, recensor = FALSE    rpsftm::rpsftm shilong1 0.9999659
## 6     Two-way switch, recensor = TRUE trtswitch::rpsftm shilong1 2.7210780
## 7     Two-way switch, recensor = TRUE    rpsftm::rpsftm shilong1 0.9953629
## 8    Two-way switch, recensor = FALSE trtswitch::rpsftm shilong1 2.7173094
## 9    Two-way switch, recensor = FALSE    rpsftm::rpsftm shilong1 1.0000278
## 10  RPSFTM on immdef, recensor = TRUE trtswitch::rpsftm  immdef1 0.7610992
## 11  RPSFTM on immdef, recensor = TRUE    rpsftm::rpsftm  immdef1 1.0035655
## 12 RPSFTM on immdef, recensor = FALSE trtswitch::rpsftm  immdef1 0.7667704
## 13 RPSFTM on immdef, recensor = FALSE    rpsftm::rpsftm  immdef1 0.9999392

2. In one-way switch scenario, trtswitch::rpsftm() and rpsftm::rpsftm() give comparable numbers for estimated psi, hazard ratio and cox p-value. The counterfactual survival time for subjects who didn’t switch (treatment arm) are consistent with original survival time, which is expected. The counterfactual survival time for subjects who switched (control arm) are larger than original survival time, which is also expected since the unadjusted HR is greater than 1. This consistency is observed in both the shilong and immdef datasets.

3. RPSFTM validation on immdef dataset: The newly added validation using the immdef1 dataset (fit5 and fit6) shows that trtswitch::rpsftm() and rpsftm::rpsftm() produce consistent results across different datasets. This cross-dataset validation strengthens confidence in the implementation.

4. In two-way switch scenario (Recensor = TRUE), trtswitch::rpsftm() and rpsftm::rpsftm() give similar estimated psi, but very different numbers for hazard ratio and cox p-value.

   • For subjects who didn’t switch: counterfactual survival time from rpsftm::rpsftm() is the same as original survival time, however, there are a couple of subjects in treatment arm (arm =1) whose counterfactual survival time from trtswitch::rpsftm() is shorter than the original survival time. This needs to be looked into (maybe in trtswitch::rpsftm, it also recensors the non-switchers).

   • For subjects who switched: counterfactual survival time from both trtswitch::rpsftm() and rpsftm::rpsftm() are larger than original survival time for control arm (arm = 0). This is aligned with a worse treatment effect and a positive estimated psi. counterfactual survival time from trtswitch::rpsftm() for treatment arm (arm = 1) are shorter than original survival time. This is also expected and was double-checked with manual calculations. However, counterfactual survival time from rpsftm::rpsftm() for treatment arm (arm = 1) didn’t get updated. The values are the same as original survival time, which may be a bug.

5. In two-way switch scenario (Recensor = FALSE), similar to the recensoring = TRUE scenario, trtswitch::rpsftm() and rpsftm::rpsftm() give comparable estimated psi, but different numbers for hazard ratio and cox p-value. This pattern suggests that the discrepancies between implementations are not driven by the recensoring mechanism, but rather by fundamental differences in the two-way switching algorithms.

   • The counterfactual survival time comparison plots show similar patterns to the recensoring = TRUE scenario, with points colored by treatment arm revealing systematic differences between the two implementations.

   • The scatter plots demonstrate that while most patients fall near the diagonal line of agreement, there are consistent deviations particularly visible when stratified by treatment arm, indicating the same underlying algorithmic differences persist regardless of recensoring setting.This consistency across recensoring settings confirms that the implementation differences observed in two-way switching scenarios are inherent to how each package handles bidirectional treatment switching, rather than being artifacts of the recensoring process.