A possible bug of `prepare_hr_forestly`

[LittleBeannie]

The following code is a minimal sufficient example for debugging. If I did not factorize TRT01P in ADSL, then the error message disappeared.

adsl <- r2rtf::r2rtf_adsl |> 
  filter(TRT01P %in% c("Placebo", "Xanomeline High Dose")) |>
  mutate(TRT01P = factor(TRT01P,
                         levels = c("Placebo", "Xanomeline High Dose"),
                         labels = c("Control", "Experimental")),
         TRT01A = factor(TRT01A,
                         levels = c("Placebo", "Xanomeline High Dose"),
                         labels = c("Control", "Experimental")),
         RACE = stringr::str_to_sentence(RACE))

set.seed(2025)

adtte <- kmcurvely_adtte |> 
  filter(TRTP %in% c("Placebo", "Xanomeline High Dose")) |>
  rename(TRT01P = TRTP) |>
  left_join(adsl |> select(USUBJID, EFFFL))
  
# manually add a new endpoint, i.e., PFS
adtte_single_endpt <- adtte |> dplyr::filter(PARAMCD == unique(adtte$PARAMCD)[1])
adtte_pfs <- adtte_single_endpt
adtte_pfs$fail_time <- -100
adtte_pfs$fail_time[which(adtte_pfs$TRT01P == "Placebo")] <- rexp(adtte_single_endpt |> filter(TRT01P == "Placebo") |> nrow(), log(2) / 6)
adtte_pfs$fail_time[which(adtte_pfs$TRT01P == "Xanomeline High Dose")] <- rexp(adtte_single_endpt |> dplyr::filter(TRT01P == "Xanomeline High Dose") |> nrow(), log(2) / (6 / 0.6))

adtte_pfs <- adtte_pfs |>
  select(-c(AVAL, CNSR)) |>
  mutate(PARAMCD = "PFS",
         enroll_time = simtrial::rpwexp_enroll(nrow(adtte_single_endpt), 
                                               enroll_rate = data.frame(duration = 6, 
                                                                        rate = nrow(adtte_single_endpt) / 6)),
         dropout_time = rexp(nrow(adtte_single_endpt), -log(0.98) / 12),
         cte = pmin(dropout_time, fail_time) + enroll_time,
         fail = (fail_time <= dropout_time) * 1,
         tte = pmin(cte, 36) - enroll_time,
         event = fail * (cte <= 36)) |>
    rename(AVAL = tte) |>
    mutate(CNSR = 1 - event) |>
    select(-c(enroll_time, dropout_time, cte, fail, event))

# manually add a new endpoint, i.e., OS
adtte_single_endpt <- adtte |> filter(PARAMCD == unique(adtte$PARAMCD)[1])
adtte_os <- adtte_single_endpt
adtte_os$fail_time <- -100
adtte_os$fail_time[which(adtte_os$TRT01P == "Placebo")] <- rexp(adtte_single_endpt |> filter(TRT01P == "Placebo") |> nrow(), log(2) / 20)
adtte_os$fail_time[which(adtte_os$TRT01P == "Xanomeline High Dose")] <- rexp(adtte_single_endpt |> filter(TRT01P == "Xanomeline High Dose") |> nrow(), log(2) / (20 / 0.6))

adtte_os <- adtte_os |>
  select(-c(AVAL, CNSR)) |>
  mutate(PARAMCD = "OS",
         enroll_time = simtrial::rpwexp_enroll(nrow(adtte_single_endpt), 
                                               enroll_rate = data.frame(duration = 6, rate = nrow(adtte_single_endpt) / 6)),
         dropout_time = rexp(nrow(adtte_single_endpt), -log(0.98) / 12),
         cte = pmin(dropout_time, fail_time) + enroll_time,
         fail = (fail_time <= dropout_time) * 1,
         tte = pmin(cte, 36) - enroll_time,
         event = fail * (cte <= 36)
         ) |>
  rename(AVAL = tte) |>
  mutate(CNSR = 1 - event) |>
  select(-c(enroll_time, dropout_time, cte, fail, event))

adtte <- rbind(adtte_pfs, adtte_os)

adtte |> 
  mutate(TRT01P = factor(TRT01P,
                         levels = c("Placebo", "Xanomeline High Dose"),
                         labels = c("Control", "Experimental"))) 

meta_km <- meta_adam(population = adsl,
                     observation = adtte) |>
  define_plan(plan = plan(analysis = "interactive_km_curve", 
                          population = "itt",
                          observation = "eff-pop", 
                          parameter = "pfs;os;male;female") |> 
                     add_plan(analysis = "hr_forestly", 
                              population = "itt",
                              observation = "eff-pop", 
                              parameter = "pfs;os;male;female;age<65;age65-80;age>80")) |>
  define_population(name = "itt",
                    group = "TRT01P",
                    subset = EFFFL == "Y",
                    var = c("USUBJID", "TRT01P", "SEX", "AGEGR1"),
                    label = "ITT") |>
  define_observation(name = "eff-pop",
                     group = "TRT01P",
                     subset = EFFFL == "Y",
                     var = c("USUBJID", "TRT01P", "SEX", "PARAMCD", "AVAL", "CNSR"),
                     label = "Efficacy population") |>
  define_parameter(name = "pfs",
                   subset = PARAMCD == "PFS",
                   label = "PFS") |>
  define_parameter(name = "os",
                   subset = PARAMCD == "OS",
                   label = "Overall Survival") |>
  define_parameter(name = "male",
                   subset = SEX == "M",
                   label = "Male") |>
  define_parameter(name = "female",
                   subset = SEX == "F",
                   label = "Female") |>
  define_parameter(name = "age<65",
                   subset = AGEGR1 == "<65",
                   label = "Age < 65") |>
  define_parameter(name = "age65-80",
                   subset = AGEGR1 == "65-80",
                   label = "Age 65-80") |>
  define_parameter(name = "age>80",
                   subset = AGEGR1 == ">80",
                   label = "Age > 80") |>
  define_analysis(name = "interactive_km_curve",
                  title = "KM curves",
                  label = "km curve") |>
  define_analysis(name = "hr_forestly",
                  title = "HR forest plot",
                  label = "HR forest plot") |>
  meta_build()

prepare_hr_forestly(meta = meta_km,
                    population = "itt",
                    observation = "eff-pop",
                    endpoint = "pfs;os",
                    subgroup = "male;female",
                    arm_levels = c("Placebo", "Xanomeline High Dose")) 

Error in survival::coxph(survival::Surv(time, event) ~ treatment, data = data_sub_trt) :
No (non-missing) observations