Modeling Event-Free Survival with Bayesian Regression Methods

Landmark-based Bayesian survival modeling with dose-response interpretation

1 Overview

This document presents a Bayesian generative modeling framework for event-free survival (EFS) in patients with cutaneous squamous cell carcinoma treated with immune checkpoint inhibitors.

Rather than treating survival modeling as a purely statistical exercise, this workflow is structured around three core principles:

  • Explicit data-generating assumptions
    Survival times are modeled using parametric distributions (Weibull), allowing direct interpretation of how covariates—including treatment exposure—shape the underlying hazard function.

  • Causal alignment through study design
    A landmark analysis anchored at the second dose is used to reduce immortal time bias and better approximate the causal effect of treatment intensity (2 doses vs ≥3 doses).

  • Full posterior inference
    Model outputs are expressed as posterior distributions, enabling probabilistic statements about clinically meaningful thresholds (for example, the probability that the hazard ratio is less than 0.8), rather than reliance on dichotomous significance testing.

The workflow proceeds from:

  1. Construction of a cleaned and auditable modeling dataset

  2. Definition of a primary landmark Weibull survival model

  3. Estimation of posterior hazard ratios and standardized survival curves

  4. A structured set of sensitivity analyses, including:

    • skeptical priors,
    • dose-intensity adjustment,
    • surgery adjustment,
    • and a no-landmark analysis for comparison

Importantly, this document is designed not only to produce results, but to serve as a transparent and reproducible record of modeling decisions, linking clinical questions, causal assumptions, and statistical implementation in a single framework.

Show R Code 
source(here::here("scripts", "load_packages.R"))

2 Load and Wrangle Modeling Dataset

We begin by loading the processed modeling dataset (df_model_raw), which is generated by the data-wrangling script:

dose_intensity_and_pre_response_exposure_sensitivity.R

Show R Code 
load_latest_rds <- function(folder, pattern) {
  files <- list.files(folder, pattern = pattern, full.names = TRUE)

  if (length(files) == 0) {
    stop(
      paste0(
        "No files found in:\n  ", folder,
        "\nmatching pattern:\n  ", pattern
      )
    )
  }

  file <-
    tibble(file = files) |>
    mutate(
      timestamp = stringr::str_extract(
        basename(file),
        "\\d{4}-\\d{2}-\\d{2}-\\d{2}-\\d{2}-\\d{2}"
      ),
      timestamp = lubridate::ymd_hms(timestamp)
    ) |>
    arrange(desc(timestamp)) |>
    slice(1) |>
    pull(file)

  readRDS(file)
}

df_model_raw <- load_latest_rds(
  here::here("files", "Dose Intensity Sensitivity"),
  "^Dose Intensity Sensitivity Model Data.*\\.rds$"
)
Show R Code 
df_analysis <- df_model_raw |>
  select(-patient, -MRN) |>
  mutate(record_id = row_number())
Show Variable Names 
colnames(df_model_raw)
##   [1] "MRN"                                             
##   [2] "patient"                                         
##   [3] "provider"                                        
##   [4] "primary_oncologist_final"                        
##   [5] "Primary Diagnosis"                               
##   [6] "immune_status"                                   
##   [7] "earliest_immunosupp_date"                        
##   [8] "immunosupp_chronic"                              
##   [9] "immunosupp_sotr"                                 
##  [10] "immunosupp_hematologic"                          
##  [11] "immunosupp_hiv"                                  
##  [12] "immunosupp_other"                                
##  [13] "ecog"                                            
##  [14] "ecog_date"                                       
##  [15] "mgmt_start_date"                                 
##  [16] "mgmt_end_date"                                   
##  [17] "sys_start_date"                                  
##  [18] "last_sys_date"                                   
##  [19] "surg_date"                                       
##  [20] "last_event_date"                                 
##  [21] "last_status_date"                                
##  [22] "last_status"                                     
##  [23] "death_date"                                      
##  [24] "has_sys"                                         
##  [25] "has_surg"                                        
##  [26] "has_rt"                                          
##  [27] "has_il"                                          
##  [28] "has_sx"                                          
##  [29] "has_other"                                       
##  [30] "regimen"                                         
##  [31] "init_trt_intent"                                 
##  [32] "num_doses"                                       
##  [33] "first_resp_date"                                 
##  [34] "diag_final"                                      
##  [35] "sys_med"                                         
##  [36] "sys_med_other"                                   
##  [37] "ned_status"                                      
##  [38] "sys_intent"                                      
##  [39] "sys_resp_path"                                   
##  [40] "sys_resp_clin"                                   
##  [41] "first_resp_clin"                                 
##  [42] "hospitalization_note"                            
##  [43] "toxicity_note"                                   
##  [44] "clinical_stage"                                  
##  [45] "primary_tumor_location"                          
##  [46] "reason_other_text"                               
##  [47] "immunosuppression_details"                       
##  [48] "other_notes"                                     
##  [49] "stage_recurrent"                                 
##  [50] "stage_pni"                                       
##  [51] "stage_locally_advanced"                          
##  [52] "stage_intransit_metastases"                      
##  [53] "stage_regionally_metastatic"                     
##  [54] "stage_distant_metastatic"                        
##  [55] "death_known"                                     
##  [56] "has_recurrence"                                  
##  [57] "has_progression"                                 
##  [58] "has_hospitalization"                             
##  [59] "has_toxicity"                                    
##  [60] "has_death"                                       
##  [61] "disease_specific_death"                          
##  [62] "marjolin_ulcer"                                  
##  [63] "initial_mgmt_complete"                           
##  [64] "reason_no_response"                              
##  [65] "reason_prog"                                     
##  [66] "reason_clin_response"                            
##  [67] "reason_toxicity"                                 
##  [68] "reason_other"                                    
##  [69] "reason_planned_course"                           
##  [70] "rt_start_date"                                   
##  [71] "rt_end_date"                                     
##  [72] "ned_date"                                        
##  [73] "recurrence_date"                                 
##  [74] "death_date_event"                                
##  [75] "progression_date"                                
##  [76] "hospitalization_date"                            
##  [77] "toxicity_date"                                   
##  [78] "sys_resp_clin_date"                              
##  [79] "diag_final_manual"                               
##  [80] "marjolin_ulcer_manual"                           
##  [81] "primary_tumor_location_manual"                   
##  [82] "immune_status_manual"                            
##  [83] "immunosupp_chronic_manual"                       
##  [84] "immunosupp_sotr_manual"                          
##  [85] "immunosupp_hematologic_manual"                   
##  [86] "immunosupp_hiv_manual"                           
##  [87] "immunosupp_other_manual"                         
##  [88] "immunosuppression_details_manual"                
##  [89] "ecog_manual"                                     
##  [90] "ecog_date_manual"                                
##  [91] "clinical_stage_manual"                           
##  [92] "stage_recurrent_manual"                          
##  [93] "stage_pni_manual"                                
##  [94] "stage_locally_advanced_manual"                   
##  [95] "stage_intransit_metastases_manual"               
##  [96] "stage_regionally_metastatic_manual"              
##  [97] "stage_distant_metastatic_manual"                 
##  [98] "mgmt_start_date_manual"                          
##  [99] "initial_mgmt_complete_manual"                    
## [100] "reason_no_response_manual"                       
## [101] "reason_prog_manual"                              
## [102] "reason_clin_response_manual"                     
## [103] "reason_toxicity_manual"                          
## [104] "reason_other_manual"                             
## [105] "reason_other_text_manual"                        
## [106] "reason_planned_course_manual"                    
## [107] "mgmt_end_date_manual"                            
## [108] "surg_date_manual"                                
## [109] "sys_start_date_manual"                           
## [110] "rt_start_date_manual"                            
## [111] "rt_end_date_manual"                              
## [112] "last_status_manual"                              
## [113] "last_status_date_manual"                         
## [114] "death_date_event_manual"                         
## [115] "recurrence_date_manual"                          
## [116] "ned_status_manual"                               
## [117] "ned_date_manual"                                 
## [118] "progression_date_manual"                         
## [119] "hospitalization_date_manual"                     
## [120] "hospitalization_note_manual"                     
## [121] "toxicity_date_manual"                            
## [122] "toxicity_note_manual"                            
## [123] "regimen_manual"                                  
## [124] "sys_med_manual"                                  
## [125] "sys_med_other_manual"                            
## [126] "sys_intent_manual"                               
## [127] "num_doses_manual"                                
## [128] "sys_resp_path_manual"                            
## [129] "sys_resp_clin_manual"                            
## [130] "sys_resp_clin_date_manual"                       
## [131] "first_resp_clin_manual"                          
## [132] "first_resp_date_manual"                          
## [133] "init_trt_intent_manual"                          
## [134] "has_recurrence_manual"                           
## [135] "has_progression_manual"                          
## [136] "has_hospitalization_manual"                      
## [137] "has_toxicity_manual"                             
## [138] "has_death_manual"                                
## [139] "disease_specific_death_manual"                   
## [140] "event_after_mgmt_manual"                         
## [141] "has_sys_manual"                                  
## [142] "has_surg_manual"                                 
## [143] "has_rt_manual"                                   
## [144] "has_sx_manual"                                   
## [145] "other_notes_manual"                              
## [146] "manual_override_flag"                            
## [147] "timestamp"                                       
## [148] "event_after_mgmt"                                
## [149] "n"                                               
## [150] "last_dose_date"                                  
## [151] "sys_before_surg"                                 
## [152] "sys_after_surg"                                  
## [153] "rt_after_surg"                                   
## [154] "dose_1_date"                                     
## [155] "dose_2_date"                                     
## [156] "dose_3_date"                                     
## [157] "dose_4_date"                                     
## [158] "dose_5_date"                                     
## [159] "dose_6_date"                                     
## [160] "dose_7_date"                                     
## [161] "dose_8_date"                                     
## [162] "dose_9_date"                                     
## [163] "dose_10_date"                                    
## [164] "dose_11_date"                                    
## [165] "dose_12_date"                                    
## [166] "dose_13_date"                                    
## [167] "dose_14_date"                                    
## [168] "dose_15_date"                                    
## [169] "dose_16_date"                                    
## [170] "dose_17_date"                                    
## [171] "dose_18_date"                                    
## [172] "dose_19_date"                                    
## [173] "dose_20_date"                                    
## [174] "dose_21_date"                                    
## [175] "dose_22_date"                                    
## [176] "dose_23_date"                                    
## [177] "dose_24_date"                                    
## [178] "dose_25_date"                                    
## [179] "dose_26_date"                                    
## [180] "dose_27_date"                                    
## [181] "dose_28_date"                                    
## [182] "dose_29_date"                                    
## [183] "dose_30_date"                                    
## [184] "dose_31_date"                                    
## [185] "dose_32_date"                                    
## [186] "dose_33_date"                                    
## [187] "dose_34_date"                                    
## [188] "dose_35_date"                                    
## [189] "dose_36_date"                                    
## [190] "dose_37_date"                                    
## [191] "dose_38_date"                                    
## [192] "dose_39_date"                                    
## [193] "dose_40_date"                                    
## [194] "dose_41_date"                                    
## [195] "dose_42_date"                                    
## [196] "dose_43_date"                                    
## [197] "dose_44_date"                                    
## [198] "dose_45_date"                                    
## [199] "dose_46_date"                                    
## [200] "Patient"                                         
## [201] "EMPI"                                            
## [202] "Sex_At_Birth"                                    
## [203] "Date_of_Birth"                                   
## [204] "Race1"                                           
## [205] "Vital_status"                                    
## [206] "Date_Of_Death"                                   
## [207] "sex"                                             
## [208] "age_at_treatment"                                
## [209] "immuno_pairs"                                    
## [210] "immuno_dx_1"                                     
## [211] "immuno_dx_2"                                     
## [212] "immuno_tx_1"                                     
## [213] "immuno_tx_2"                                     
## [214] "immunosupp_category"                             
## [215] "immunosupp_treatment_combined"                   
## [216] "surg_lab"                                        
## [217] "dose_group"                                      
## [218] "surg_dose_group"                                 
## [219] "expected_interval_days"                          
## [220] "interval_source"                                 
## [221] "likely_schedule"                                 
## [222] "any_200mg"                                       
## [223] "any_400mg"                                       
## [224] "has_med_match"                                   
## [225] "med_episode_match_7d"                            
## [226] "response_cutoff_date"                            
## [227] "num_doses_recorded"                              
## [228] "dose_intensity_final"                            
## [229] "mean_interval_days"                              
## [230] "interval_deviation_days"                         
## [231] "num_doses_pre_response"                          
## [232] "dose_intensity_pre_response"                     
## [233] "dose_intensity_pre_response_uncapped"            
## [234] "dose_intensity_pre_response_uncapped_wins"       
## [235] "dose_intensity_pre_response_uncapped_centered"   
## [236] "n_intervals_observed"                            
## [237] "n_intervals_with_expected"                       
## [238] "n_recorded_200mg"                                
## [239] "n_recorded_400mg"                                
## [240] "n_inferred_q3"                                   
## [241] "n_inferred_q6"                                   
## [242] "n_unclassified"                                  
## [243] "expected_days_sum"                               
## [244] "observed_days_sum"                               
## [245] "dose_intensity_pembro_interval"                  
## [246] "dose_intensity_pembro_interval_uncapped"         
## [247] "dose_intensity_pembro_interval_uncapped_wins"    
## [248] "dose_intensity_pembro_interval_uncapped_centered"
## [249] "extra_doses_after_response"                      
## [250] "dose_intensity_source"                           
## [251] "agent_schedule"
Show Code to Build Model Data 
df_analysis <- df_analysis |>
  mutate(
    event_date_raw = pmin(
      recurrence_date,
      progression_date,
      death_date_event,
      na.rm = TRUE
    ),
    event_date_raw = as.Date(event_date_raw, origin = "1970-01-01"),
    event_date = if_else(
      is.na(event_date_raw),
      as.Date(NA),
      event_date_raw
    ),
    event = as.integer(!is.na(event_date)),
    time_observed = as.numeric(
      difftime(
        if_else(event == 1, event_date, last_status_date),
        sys_start_date,
        units = "days"
      )
    ) / 30.44,

    age_centered = age_at_treatment - 76,

    immunosuppressed = factor(
      immune_status == "Immunosuppressed",
      levels = c(FALSE, TRUE)
    ),

    location_condensed = case_when(
      str_detect(primary_tumor_location, "Head|Neck|Face|Scalp|Ear") ~ "Head/Neck",
      str_detect(primary_tumor_location, "Conjunctiva|Eye") ~ "Conjunctiva",
      str_detect(primary_tumor_location, "Unknown") ~ "Unknown Primary",
      TRUE ~ "Other"
    ),
    location_condensed = factor(
      location_condensed,
      levels = c("Head/Neck", "Conjunctiva", "Other", "Unknown Primary")
    ),

    stage_condensed = case_when(
      clinical_stage %in% c("I", "II") ~ "Stage I-II",
      clinical_stage == "III" ~ "Stage III",
      clinical_stage == "IV" ~ "Stage IV",
      TRUE ~ "Unstaged"
    ),
    stage_condensed = factor(
      stage_condensed,
      levels = c("Stage I-II", "Stage III", "Stage IV", "Unstaged")
    ),

    ecog_num = suppressWarnings(as.numeric(ecog)),
    ecog_condensed = case_when(
      ecog_num %in% c(0, 1) ~ "ECOG 0-1",
      ecog_num %in% c(2, 3) ~ "ECOG 2-3",
      TRUE ~ NA_character_
    ),
    ecog_condensed = factor(
      ecog_condensed,
      levels = c("ECOG 0-1", "ECOG 2-3")
    ),

    agent_schedule = case_when(
      agent_schedule %in% c("cemiplimab", "Cemiplimab") ~ "cemiplimab",
      agent_schedule %in% c("pembro_q3", "Pembrolizumab_q3", "pembro q3") ~ "pembro_q3",
      agent_schedule %in% c("pembro_q6", "Pembrolizumab_q6", "pembro q6") ~ "pembro_q6",
      agent_schedule %in% c("pembro_mixed", "pembro mixed") ~ "pembro_mixed",
      agent_schedule %in% c("Ipi-nivo", "ipi-nivo", "Ipi+Nivo", "ipinivo") ~ "Ipi-nivo",
      TRUE ~ NA_character_
    ),
    agent_schedule = factor(
      agent_schedule,
      levels = c("cemiplimab", "pembro_q3", "pembro_q6", "pembro_mixed", "Ipi-nivo")
    ),

    surgery = as.integer(has_surg)
  ) |>
  filter(!is.na(num_doses), num_doses >= 1) |>
  filter(!is.na(agent_schedule))
Show Raw EFS Variable Checks 
df_analysis %>%
  summarise(
    n = n(),
    missing_event = sum(is.na(event)),
    missing_dose_2_date = sum(is.na(dose_2_date)),
    missing_last_status_date = sum(is.na(last_status_date)),
    missing_recurrence_date = sum(is.na(recurrence_date)),
    missing_progression_date = sum(is.na(progression_date)),
    missing_death_date_event = sum(is.na(death_date_event))
  )
## # A tibble: 1 × 7
##       n missing_event missing_dose_2_date missing_last_status_date
##   <int>         <int>               <int>                    <int>
## 1   189             0                  11                        0
## # ℹ 3 more variables: missing_recurrence_date <int>,
## #   missing_progression_date <int>, missing_death_date_event <int>
Show Code to Transform Data for Landmark Analysis 
df_landmark_pre <- df_analysis %>%
  filter(num_doses >= 2, !is.na(dose_2_date)) %>%
  mutate(
    raw_time_from_dose2 = case_when(
      event == 1 ~ as.numeric(difftime(event_date, dose_2_date, units = "days")) / 30.44,
      TRUE ~ as.numeric(difftime(last_status_date, dose_2_date, units = "days")) / 30.44
    )
  )

df_landmark <- df_landmark_pre %>%
  mutate(
    time_from_dose2 = case_when(
      event == 1 & raw_time_from_dose2 == 0 ~ 0.016,
      TRUE ~ raw_time_from_dose2
    ),
    event_landmark = if_else(event == 1 & time_from_dose2 > 0, 1L, 0L),
    time_observed_landmark = pmax(time_from_dose2, 0.01),
    dose_group_landmark = case_when(
      num_doses == 2 ~ "2 doses",
      num_doses >= 3 ~ "3+ doses",
      TRUE ~ NA_character_
    ),
    dose_group_landmark = factor(
      dose_group_landmark,
      levels = c("2 doses", "3+ doses")
    ),
    dose_3plus = case_when(
      num_doses == 2 ~ 0L,
      num_doses >= 3 ~ 1L,
      TRUE ~ NA_integer_
    )
  ) %>%
  filter(time_from_dose2 >= 0)
Evaluate Number of Patients Per Dose Level 
df_landmark %>%
  count(num_doses, dose_group_landmark)
## # A tibble: 17 × 3
##    num_doses dose_group_landmark     n
##        <int> <fct>               <int>
##  1         2 2 doses                81
##  2         3 3+ doses               19
##  3         4 3+ doses               19
##  4         5 3+ doses               10
##  5         6 3+ doses                9
##  6         7 3+ doses                8
##  7         8 3+ doses                8
##  8         9 3+ doses                6
##  9        10 3+ doses                1
## 10        11 3+ doses                3
## 11        13 3+ doses                2
## 12        14 3+ doses                3
## 13        15 3+ doses                1
## 14        16 3+ doses                4
## 15        20 3+ doses                1
## 16        28 3+ doses                1
## 17        46 3+ doses                1
Show landmark dose-group summary
Landmark cohort dose-group distribution
Total patients 2 doses 3+ doses Missing dose group
177 81 96 0
Show landmark timing audit
Landmark Timing Audit
Metric Value
Patients with >=2 doses and dose 2 date 178
Events before dose 2 1
Events on dose 2 date 2
Missing last status date 0
Missing event date among events 0

2.1 Summarize Exclusions from the Landmark Cohort

Show landmark exclusion summary code
Exclusions from the Landmark Cohort
Exclusion reason Count
<2 doses 11
Event or censoring before landmark / other 1

The table below summarizes why patients from the full analytic dataset were not included in the landmark cohort.

2.2 Descriptive Checks

Show overall descriptive summary
Overall Landmark Cohort Summary
Metric Value
Total patients 177.00
Number of events 72.00
Event rate 0.41
Median follow-up (months) 9.72
Show outcomes by dose group
Outcomes by Dose Group
Dose group Patients Events Event rate Median follow-up (months)
2 doses 81 31 0.38 5.98
3+ doses 96 41 0.43 15.08

2.3 Verify Landmark Exposure Group Coding

Before fitting survival models, we confirm that the landmark exposure groups align with the observed number of doses.

Show dose-group definition audit
Dose Group by Observed Number of Doses
Dose group Observed number of doses Count
2 doses 2 81
3+ doses 3 19
3+ doses 4 19
3+ doses 5 10
3+ doses 6 9
3+ doses 7 8
3+ doses 8 8
3+ doses 9 6
3+ doses 10 1
3+ doses 11 3
3+ doses 13 2
3+ doses 14 3
3+ doses 15 1
3+ doses 16 4
3+ doses 20 1
3+ doses 28 1
3+ doses 46 1
Show dose-group count code 
df_landmark %>%
  count(dose_group_landmark, sort = TRUE, name = "Count") %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Landmark Cohort by Dose Group"
  ) %>%
  gt::cols_label(
    dose_group_landmark = "Dose group"
  ) %>%
  gt::fmt_number(
    columns = Count,
    decimals = 0
  ) %>%
  gt::cols_align(
    align = "center",
    columns = Count
  )
Landmark Cohort by Dose Group
Dose group Count
3+ doses 96
2 doses 81
Show dose-group level check code 
tibble::tibble(
  `Defined factor levels` = paste(levels(df_landmark$dose_group_landmark), collapse = ", "),
  `Observed values` = paste(unique(df_landmark$dose_group_landmark), collapse = ", ")
) %>%
  gt::gt() %>%
  gt::tab_header(
    title = "Dose Group Level Check"
  ) %>%
  gt::cols_align(
    align = "left",
    columns = everything()
  )
Dose Group Level Check
Defined factor levels Observed values
2 doses, 3+ doses 3+ doses, 2 doses

2.3.1 Unadjusted Kaplan-Meier Sanity Check

Before fitting the Bayesian Weibull model, we generate an unadjusted Kaplan-Meier curve as a basic sanity check on the landmark dataset and survival coding.

Show R Code 
fit_km_landmark <- survfit(
  Surv(time_observed_landmark, event_landmark) ~ dose_group_landmark,
  data = df_landmark
)

ggsurvplot(
  fit_km_landmark,
  data = df_landmark,
  risk.table = TRUE,
  conf.int = TRUE,
  xlab = "Months from dose 2 landmark",
  ylab = "Event-free survival probability",
  legend.title = NULL
)

3 Prepare Modeling Variables for EFS

We next create the variables needed for the landmark Weibull models. These transformations mirror the coding used in the response analyses so that coefficient interpretation remains consistent across endpoints.

Show EFS modeling variable construction
Show R Code 
df_landmark_model <- df_landmark %>%
  mutate(
    cens = 1L - event_landmark,

    age_centered_10 = age_centered / 10,

    immunosuppressed = as.integer(immunosuppressed),

    ecog_2plus = case_when(
      ecog_condensed == "ECOG 2-3" ~ 1L,
      ecog_condensed == "ECOG 0-1" ~ 0L,
      TRUE ~ NA_integer_
    ),

    stage_III = if_else(stage_condensed == "Stage III", 1L, 0L, missing = 0L),
    stage_IV = if_else(stage_condensed == "Stage IV", 1L, 0L, missing = 0L),
    stage_unstaged = if_else(stage_condensed == "Unstaged", 1L, 0L, missing = 0L),

    agent_pembro_q3 = if_else(agent_schedule == "pembro_q3", 1L, 0L, missing = 0L),
    agent_pembro_q6 = if_else(agent_schedule == "pembro_q6", 1L, 0L, missing = 0L),
    agent_pembro_mixed = if_else(agent_schedule == "pembro_mixed", 1L, 0L, missing = 0L),
    agent_ipinivo = if_else(agent_schedule == "Ipi-nivo", 1L, 0L, missing = 0L),

    dose_3plus = case_when(
      num_doses == 2 ~ 0L,
      num_doses >= 3 ~ 1L,
      TRUE ~ NA_integer_
    ),

    dose_group_landmark = factor(
      if_else(dose_3plus == 1L, "3+ doses", "2 doses"),
      levels = c("2 doses", "3+ doses")
    ),

    dose_intensity_centered = as.numeric(
      scale(dose_intensity_final, center = TRUE, scale = FALSE)
    ),

    surgery = as.integer(surgery)
  )
Show dose intensity variable check
Dose Intensity (Centered) — Distribution Check
Metric Value
Total patients 177.00
Missing dose intensity 0.00
Mean (centered) 0.00
SD (centered) 0.15
Min −0.86
Max 0.08
Show modeling dataset completeness
Modeling Dataset Completeness
Metric Value
Total patients 177
2 doses 81
3+ doses 96
Missing dose group 0
Show outcomes by dose group (modeling dataset)
Outcomes by Dose Group (Modeling Dataset)
Dose group Patients Events Event rate Median follow-up (months)
2 doses 81 31 0.38 5.98
3+ doses 96 41 0.43 15.08

3.1 Restrict to Complete Cases for the Primary Landmark Model

For the initial EFS model, we remove observations with missing covariates required by the model. This keeps the fitting dataset explicit and auditable.

Show R Code 
df_efs_primary <- df_landmark_model %>%
  filter(
    !is.na(time_observed_landmark),
    !is.na(event_landmark),
    !is.na(cens),
    !is.na(dose_3plus),
    !is.na(age_centered_10),
    !is.na(immunosuppressed),
    !is.na(ecog_2plus)
  )

4 Summarize the Primary Modeling Dataset

4.1 Overall primary modeling dataset

Show primary modeling dataset summary
Primary Modeling Dataset Summary
Metric Value
Total patients 177.00
Number of events 72.00
Event rate 0.41
Median follow-up (months) 9.72

4.2 By dose group

Show primary dataset outcomes by dose group
Primary Dataset Outcomes by Dose Group
Dose group Patients Events Event rate Median follow-up (months)
2 doses 81 31 0.38 5.98
3+ doses 96 41 0.43 15.08

5 Define the Primary Weibull Landmark Model

The primary analysis models event-free survival using a parametric Weibull regression, which directly specifies the hazard function as a function of treatment exposure and baseline covariates.

This formulation allows estimation of the association between receiving ≥3 doses of immunotherapy (versus 2 doses) and the hazard of recurrence, progression, or death, while adjusting for key clinical factors including age, immunosuppression status, performance status, stage, and treatment regimen.

The model is implemented using the brms package. The function bf() (“Bayesian formula”) defines the model structure, including: - the outcome variable, - censoring information, - and the linear predictor relating covariates to the hazard.

In this specification, the outcome is time from the landmark (dose 2) to event or censoring, with right-censoring incorporated via the cens() modifier.

Show primary Weibull model specification
Show R Code 
# Define the Bayesian model formula using brms
# bf() specifies the likelihood structure and linear predictor

form_efs_primary <- bf(
  # Survival outcome with censoring indicator
  # time_observed_landmark: time from landmark to event/censoring
  # cens(cens): 1 = censored, 0 = event
  time_observed_landmark | cens(cens) ~

    # Primary exposure: treatment intensity
    # 1 = ≥3 doses, 0 = 2 doses
    dose_3plus +

    # Demographic covariate (scaled per 10 years)
    age_centered_10 +

    # Clinical covariates
    immunosuppressed +
    ecog_2plus +

    # Disease stage (indicator variables)
    stage_III +
    stage_IV +
    stage_unstaged +

    # Treatment regimen indicators
    agent_pembro_q3 +
    agent_pembro_q6 +
    agent_pembro_mixed +
    agent_ipinivo
)

5.1 Inspect Prior Slots

Inspect Priors 
get_prior(
  formula = form_efs_primary,
  data = df_efs_primary,
  family = weibull()
)
##                   prior     class               coef group resp dpar nlpar lb
##                  (flat)         b                                            
##                  (flat)         b    age_centered_10                         
##                  (flat)         b      agent_ipinivo                         
##                  (flat)         b agent_pembro_mixed                         
##                  (flat)         b    agent_pembro_q3                         
##                  (flat)         b    agent_pembro_q6                         
##                  (flat)         b         dose_3plus                         
##                  (flat)         b         ecog_2plus                         
##                  (flat)         b   immunosuppressed                         
##                  (flat)         b          stage_III                         
##                  (flat)         b           stage_IV                         
##                  (flat)         b     stage_unstaged                         
##  student_t(3, 2.3, 2.5) Intercept                                            
##       gamma(0.01, 0.01)     shape                                           0
##  ub       source
##          default
##     (vectorized)
##     (vectorized)
##     (vectorized)
##     (vectorized)
##     (vectorized)
##     (vectorized)
##     (vectorized)
##     (vectorized)
##     (vectorized)
##     (vectorized)
##     (vectorized)
##          default
##          default

5.2 Specify Priors

We specify weakly informative priors for the Weibull survival model, with particular attention to the prior on the primary exposure (receipt of ≥3 doses).

Because the Weibull model operates on the log-hazard scale, we define priors in terms of hazard ratios (HRs) and then transform them to the corresponding regression coefficients. This improves interpretability and allows direct encoding of clinically meaningful assumptions.

Specifically:

  • A primary prior is centered at a modest protective effect (HR = 0.85)
  • A skeptical prior is centered at no effect (HR = 1.00)

These are mapped to the linear predictor scale using the Weibull shape parameter, ensuring coherence between the prior and the model parameterization.

All other coefficients are assigned mean-zero regularizing priors, and the Weibull shape parameter is given a weakly informative gamma prior.

Show prior specification code
Show R Code 
# Reference Weibull shape parameter used for HR → beta transformation
reference_shape <- 1.25

# Convert hazard ratio to regression coefficient (log-hazard scale)
beta_from_hr <- function(hr, shape = reference_shape) {
  -log(hr) / shape
}

# Define target prior centers on the hazard ratio scale
primary_prior_hr_center <- 0.85
skeptical_prior_hr_center <- 1.00

# Convert to beta scale for use in the model
primary_prior_beta_center <- beta_from_hr(primary_prior_hr_center)
skeptical_prior_beta_center <- beta_from_hr(skeptical_prior_hr_center)

# Helper function to construct prior sets
make_efs_priors <- function(mean_dose = 0, sd_dose = 0.35) {
  c(
    prior_string("normal(3.0, 0.35)", class = "Intercept"),
    prior_string(
      paste0("normal(", round(mean_dose, 4), ", ", sd_dose, ")"),
      class = "b",
      coef = "dose_3plus"
    ),
    prior_string("normal(0, 0.30)", class = "b"),
    prior_string("gamma(5, 4)", class = "shape")
  )
}

# Primary prior specification
priors_efs_primary <- make_efs_priors(
  mean_dose = primary_prior_beta_center,
  sd_dose = 0.35
)

# Skeptical prior specification
priors_efs_skeptical <- make_efs_priors(
  mean_dose = skeptical_prior_beta_center,
  sd_dose = 0.20
)

# Dose-intensity model priors
priors_efs_dose_intensity <- c(
  prior_string("normal(3.0, 0.35)", class = "Intercept"),
  prior_string(
    paste0("normal(", round(primary_prior_beta_center, 4), ", 0.35)"),
    class = "b",
    coef = "dose_3plus"
  ),
  prior_string("normal(0, 0.30)", class = "b", coef = "dose_intensity_centered"),
  prior_string("normal(0, 0.30)", class = "b"),
  prior_string("gamma(5, 4)", class = "shape")
)

# Surgery-adjusted model priors
priors_efs_surgery <- c(
  prior_string("normal(3.0, 0.35)", class = "Intercept"),
  prior_string(
    paste0("normal(", round(primary_prior_beta_center, 4), ", 0.35)"),
    class = "b",
    coef = "dose_3plus"
  ),
  prior_string("normal(0, 0.30)", class = "b", coef = "surgery"),
  prior_string("normal(0, 0.30)", class = "b"),
  prior_string("gamma(5, 4)", class = "shape")
)

# Summarize prior centers for interpretability
tibble(
  Prior = c("Primary", "Skeptical"),
  `Target HR (center)` = c(primary_prior_hr_center, skeptical_prior_hr_center),
  `Beta (log-hazard scale)` = c(primary_prior_beta_center, skeptical_prior_beta_center)
)
# A tibble: 2 × 3
  Prior     `Target HR (center)` `Beta (log-hazard scale)`
  <chr>                    <dbl>                     <dbl>
1 Primary                   0.85                     0.130
2 Skeptical                 1                        0

5.3 Fit the Primary Weibull Model

We next fit the primary Bayesian Weibull survival model using the landmark cohort, the prespecified covariate set, and the primary prior distribution for the dose effect.

This model represents the main inferential analysis for event-free survival in the document. It combines:

  • the landmark-based cohort definition,
  • the Weibull survival likelihood,
  • adjustment for baseline clinical and treatment-related covariates,
  • and a prior centered on a modest protective association for receipt of ≥3 doses.

The model is fit using brm() from the brms package, which serves as the main interface for fitting Bayesian regression models in Stan. In this context, brm() compiles the model, performs Hamiltonian Monte Carlo sampling, and returns posterior draws for all model parameters.

The settings below use 4 Markov chains with 4,000 iterations per chain, including 2,000 warmup iterations, providing a stable posterior sample for downstream summaries, hazard-ratio calculations, and posterior standardized survival curves.

Show primary Weibull model fitting code
Show R Code 
fit_efs_primary <- brm(
  formula = form_efs_primary,
  data = df_efs_primary,
  family = weibull(),
  prior = priors_efs_primary,
  chains = 4,
  iter = 4000,
  warmup = 2000,
  cores = 4,
  seed = 2201,
  refresh = 0
)
Show R Code 
fit_efs_primary_landmark_weibull <- fit_efs_primary
Show R Code 
summary(fit_efs_primary)
 Family: weibull 
  Links: mu = log; shape = identity 
Formula: time_observed_landmark | cens(cens) ~ dose_3plus + age_centered_10 + immunosuppressed + ecog_2plus + stage_III + stage_IV + stage_unstaged + agent_pembro_q3 + agent_pembro_q6 + agent_pembro_mixed + agent_ipinivo 
   Data: df_efs_primary (Number of observations: 177) 
  Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
         total post-warmup draws = 8000

Regression Coefficients:
                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept              4.25      0.42     3.43     5.08 1.00    13752     6520
dose_3plus             0.39      0.25    -0.11     0.88 1.00    16644     5961
age_centered_10       -0.19      0.14    -0.47     0.09 1.00    15066     6805
immunosuppressed      -0.20      0.26    -0.69     0.30 1.00    16304     5816
ecog_2plus            -0.11      0.27    -0.64     0.41 1.00    16767     6027
stage_III              0.20      0.24    -0.28     0.68 1.00    13146     6161
stage_IV              -0.38      0.24    -0.85     0.09 1.00    14209     6877
stage_unstaged         0.03      0.26    -0.49     0.54 1.00    15899     5510
agent_pembro_q3       -0.33      0.24    -0.80     0.16 1.00    15991     5877
agent_pembro_q6       -0.01      0.29    -0.57     0.54 1.00    16383     6039
agent_pembro_mixed     0.09      0.27    -0.46     0.61 1.00    14974     6105
agent_ipinivo          0.02      0.30    -0.55     0.60 1.00    16232     6237

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape     0.62      0.05     0.53     0.72 1.00    10656     6414

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

5.4 Extract Posterior Hazard-Ratio Draws

Under the brms Weibull parameterization, the hazard ratio for 3+ doses versus 2 doses is:

HR=exp(−shape×β dose_3plus)

Show R Code 
extract_weibull_hr_draws <- function(fit, coef_name = "b_dose_3plus") {
  draws <- as_draws_df(fit)

  tibble(
    beta_dose = draws[[coef_name]],
    shape = draws[["shape"]],
    hr = exp(-draws[["shape"]] * draws[[coef_name]])
  )
}

posterior_hr_primary <- extract_weibull_hr_draws(fit_efs_primary)
Show R Code 
posterior_hr_efs_primary_landmark <- posterior_hr_primary
Show R Code 
posterior_hr_primary %>%
  summarise(
    mean_hr = mean(hr),
    median_hr = median(hr),
    lower_89 = quantile(hr, 0.055),
    upper_89 = quantile(hr, 0.945),
    p_hr_less_1 = mean(hr < 1),
    p_hr_less_0_8 = mean(hr < 0.8),
    p_hr_less_0_6 = mean(hr < 0.6)
  )
# A tibble: 1 × 7
  mean_hr median_hr lower_89 upper_89 p_hr_less_1 p_hr_less_0_8 p_hr_less_0_6
    <dbl>     <dbl>    <dbl>    <dbl>       <dbl>         <dbl>         <dbl>
1   0.795     0.786    0.606     1.01       0.934         0.543        0.0508

5.5 Plot the Posterior Dose-Effect Distribution

create_hr_density_plot <- function(
  posterior_hr,
  title,
  subtitle,
  x_limits = c(0.4, 1.4),
  x_label = "Hazard ratio (3+ doses vs 2 doses)"
) {

  dens <- density(posterior_hr)
  dens_df <- tibble(x = dens$x, y = dens$y)

  ci_lower <- quantile(posterior_hr, 0.055)
  ci_upper <- quantile(posterior_hr, 0.945)

  prob_protective <- mean(posterior_hr < 1.0)
  prob_modest <- mean(posterior_hr < 0.8)
  prob_meaningful <- mean(posterior_hr < 0.6)

  ggplot(dens_df, aes(x = x, y = y)) +
    geom_area(fill = "steelblue", alpha = 0.7) +
    geom_area(data = subset(dens_df, x < 1.0), fill = "green", alpha = 0.22) +
    geom_area(data = subset(dens_df, x < 0.8), fill = "#CFEFE5", alpha = 0.35) +
    geom_area(data = subset(dens_df, x < 0.6), fill = "purple", alpha = 0.18) +
    geom_vline(xintercept = 1.0, linetype = "dashed", color = "red", linewidth = 1) +
    geom_vline(xintercept = 0.8, linetype = "dashed", color = "orange", linewidth = 0.8) +
    geom_vline(xintercept = 0.6, linetype = "dashed", color = "purple", linewidth = 0.8) +
    geom_vline(xintercept = c(ci_lower, ci_upper), linetype = "dotted", color = "darkblue") +
    annotate(
      "text",
      x = x_limits[1], y = max(dens_df$y) * 0.95,
      label = paste0("P(HR < 1.0) = ", round(prob_protective * 100, 1), "%"),
      hjust = 0, size = 3.3, color = "darkgreen"
    ) +
    annotate(
      "text",
      x = x_limits[1], y = max(dens_df$y) * 0.85,
      label = paste0("P(HR < 0.8) = ", round(prob_modest * 100, 1), "%"),
      hjust = 0, size = 3.3, color = "#2E7C71"
    ) +
    annotate(
      "text",
      x = x_limits[1], y = max(dens_df$y) * 0.75,
      label = paste0("P(HR < 0.6) = ", round(prob_meaningful * 100, 1), "%"),
      hjust = 0, size = 3.3, color = "purple"
    ) +
    annotate(
      "text",
      x = mean(x_limits), y = max(dens_df$y) * 0.48,
      label = paste0(
        "89% CrI: [", round(ci_lower, 2), ", ", round(ci_upper, 2), "]"
      ),
      size = 3.5
    ) +
    labs(
      title = title,
      subtitle = subtitle,
      x = x_label,
      y = "Density"
    ) +
    scale_x_continuous(limits = x_limits) +
    theme_minimal(base_size = 12) +
    theme(
      plot.title = element_text(hjust = 0.5, face = "bold"),
      plot.subtitle = element_text(hjust = 0.5, face = "bold")
    )
}

create_hr_density_plot(
  posterior_hr_primary$hr,
  title = "EFS dose effect",
  subtitle = "Primary landmark model"
)

5.6 Code for posterior survival curves

Show R Code 
library(posterior)

draws <- as_draws_df(fit_efs_primary)

df_pred_2 <- df_efs_primary %>%
  mutate(
    dose_3plus = 0L,
    dose_group_landmark = factor("2 doses", levels = c("2 doses", "3+ doses"))
  )

df_pred_3plus <- df_efs_primary %>%
  mutate(
    dose_3plus = 1L,
    dose_group_landmark = factor("3+ doses", levels = c("2 doses", "3+ doses"))
  )
Show R Code 
lp_2 <- posterior_linpred(
  fit_efs_primary,
  newdata = df_pred_2,
  re_formula = NA
)

lp_3plus <- posterior_linpred(
  fit_efs_primary,
  newdata = df_pred_3plus,
  re_formula = NA
)

mu_2 <- exp(lp_2)
mu_3plus <- exp(lp_3plus)

shape_draws <- draws$shape

5.7 Compute Survival Curves

Show R Code 
time_grid <- seq(0, 60, by = 1)

compute_survival <- function(mu_mat, shape_vec, times) {

  n_draws <- nrow(mu_mat)
  n_times <- length(times)

  surv_mat <- matrix(NA, nrow = n_draws, ncol = n_times)

  for (i in seq_len(n_times)) {
    t <- times[i]

    surv_mat[, i] <- exp(
      -((t / mu_mat)^shape_vec)
    ) |> rowMeans()
  }

  surv_mat
}

surv_2 <- compute_survival(mu_2, shape_draws, time_grid)
surv_3plus <- compute_survival(mu_3plus, shape_draws, time_grid)

5.8 Summarize Curves

Show R Code 
summarize_surv <- function(surv_mat, label) {
  tibble(
    time = time_grid,
    median = apply(surv_mat, 2, median),
    lower = apply(surv_mat, 2, quantile, 0.055),
    upper = apply(surv_mat, 2, quantile, 0.945),
    group = label
  )
}

surv_df_efs_primary_landmark <- bind_rows(
  summarize_surv(surv_2, "2 doses"),
  summarize_surv(surv_3plus, "3+ doses")
)
Show R Code 
hr_probs_primary <- posterior_hr_primary %>%
  summarise(
    p_hr_lt_1   = mean(hr < 1.0),
    p_hr_lt_09  = mean(hr < 0.9),
    p_hr_lt_08  = mean(hr < 0.8),
    p_hr_lt_06  = mean(hr < 0.6),
    median_hr   = median(hr),
    lower_89_hr = quantile(hr, 0.055),
    upper_89_hr = quantile(hr, 0.945)
  )

hr_probs_primary
# A tibble: 1 × 7
  p_hr_lt_1 p_hr_lt_09 p_hr_lt_08 p_hr_lt_06 median_hr lower_89_hr upper_89_hr
      <dbl>      <dbl>      <dbl>      <dbl>     <dbl>       <dbl>       <dbl>
1     0.934      0.805      0.543     0.0508     0.786       0.606        1.01
Show R Code 
library(scales)
ann_text_primary <- paste0(
  "Median HR = ", round(hr_probs_primary$median_hr, 2),
  " (89% CrI ", round(hr_probs_primary$lower_89_hr, 2),
  "\u2013", round(hr_probs_primary$upper_89_hr, 2), ")\n",
  "P(HR < 1.0) = ", percent(hr_probs_primary$p_hr_lt_1, accuracy = 0.1), "\n",
  "P(HR < 0.9) = ", percent(hr_probs_primary$p_hr_lt_09, accuracy = 0.1), "\n",
  "P(HR < 0.8) = ", percent(hr_probs_primary$p_hr_lt_08, accuracy = 0.1), "\n",
  "P(HR < 0.6) = ", percent(hr_probs_primary$p_hr_lt_06, accuracy = 0.1)
)

5.9 Plot

Show R Code 
ggplot(
  surv_df_efs_primary_landmark,
  aes(x = time, y = median, color = group, fill = group)
) +
  geom_ribbon(
    aes(ymin = lower, ymax = upper),
    alpha = 0.15,
    color = NA
  ) +
  geom_line(linewidth = 1.3) +
  annotate(
    "label",
    x = 36,
    y = 1.0,
    label = ann_text_primary,
    hjust = 0,
    vjust = 1,
    size = 3.7,
    label.size = 0.25,
    fill = "white",
    color = "black"
  ) +
  scale_x_continuous(
    breaks = seq(0, 60, by = 12),
    limits = c(0, 60)
  ) +
  coord_cartesian(ylim = c(0, 1)) +
  labs(
    x = "Months from dose 2 landmark",
    y = "Predicted event-free survival",
    color = NULL,
    fill = NULL,
    title = "Posterior standardized EFS curves",
    subtitle = "Primary landmark Weibull model"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, face = "bold"),
    legend.position = "right"
  )

Posterior standardized event-free survival curves from the primary landmark Weibull model. Survival curves represent posterior predictions from the fitted Bayesian Weibull regression model after standardizing over the observed covariate distribution in the study cohort. For each posterior draw, predicted survival was computed after assigning all patients to either two doses or three or more doses of immunotherapy while leaving all other covariates unchanged. The resulting curves therefore represent model-estimated survival under each dosing scenario averaged across the empirical patient population. Shaded regions denote 89% credible intervals. Because the model uses a Weibull survival specification rather than a Cox proportional hazards model, the hazard functions are not constrained to be proportional over time. The annotation box summarizes the posterior distribution of the hazard ratio comparing ≥3 doses with 2 doses, including the median hazard ratio, its 89% credible interval, and posterior probabilities that the hazard ratio is below several clinically interpretable thresholds.
Show R Code 
# -----------------------------
# Primary EFS model object
# -----------------------------
fit_efs_primary_landmark_weibull <- fit_efs_primary
Show R Code 
# -----------------------------
# Posterior HR draws
# -----------------------------
posterior_hr_efs_primary_landmark <- posterior_hr_primary
Show R Code 
# -----------------------------
# HR summary table
# -----------------------------
hr_summary_efs_primary_landmark <- posterior_hr_efs_primary_landmark %>%
  summarise(
    mean_hr = mean(hr),
    median_hr = median(hr),
    lower_89 = quantile(hr, 0.055),
    upper_89 = quantile(hr, 0.945),
    p_hr_less_1 = mean(hr < 1),
    p_hr_less_0_9 = mean(hr < 0.9),
    p_hr_less_0_8 = mean(hr < 0.8),
    p_hr_less_0_6 = mean(hr < 0.6),
    p_hr_less_0_5 = mean(hr < 0.5)
  )
Show R Code 
# -----------------------------
# Annotation text for figure
# -----------------------------
ann_text_efs_primary_landmark <- paste0(
  "Median HR = ", round(hr_summary_efs_primary_landmark$median_hr, 2), "\n",
  "89% CrI: ", round(hr_summary_efs_primary_landmark$lower_89, 2),
  "\u2013", round(hr_summary_efs_primary_landmark$upper_89, 2), "\n\n",
  "P(HR < 1.0) = ", scales::percent(hr_summary_efs_primary_landmark$p_hr_less_1, accuracy = 0.1), "\n",
  "P(HR < 0.9) = ", scales::percent(hr_summary_efs_primary_landmark$p_hr_less_0_9, accuracy = 0.1), "\n",
  "P(HR < 0.8) = ", scales::percent(hr_summary_efs_primary_landmark$p_hr_less_0_8, accuracy = 0.1), "\n",
  "P(HR < 0.6) = ", scales::percent(hr_summary_efs_primary_landmark$p_hr_less_0_6, accuracy = 0.1)
)
Show R Code 
surv_2_efs_primary_landmark <- surv_2
surv_3plus_efs_primary_landmark <- surv_3plus
time_grid_efs_primary_landmark <- time_grid

6 Sensitivity Analysis: Fit EFS Model with Skeptical prior

Show R Code 
fit_efs_skeptical <- brm(
  formula = form_efs_primary,
  data = df_efs_primary,
  family = weibull(),
  prior = priors_efs_skeptical,
  chains = 4,
  iter = 4000,
  warmup = 2000,
  cores = 4,
  seed = 2202,
  refresh = 0
)
Show R Code 
fit_efs_skeptical_landmark_weibull <- fit_efs_skeptical
Show R Code 
posterior_hr_skeptical <- extract_weibull_hr_draws(fit_efs_skeptical)

posterior_hr_efs_skeptical_landmark <- posterior_hr_skeptical
Show R Code 
hr_summary_efs_skeptical_landmark <- posterior_hr_efs_skeptical_landmark %>%
  summarise(
    mean_hr = mean(hr),
    median_hr = median(hr),
    lower_89 = quantile(hr, 0.055),
    upper_89 = quantile(hr, 0.945),
    p_hr_less_1 = mean(hr < 1),
    p_hr_less_0_9 = mean(hr < 0.9),
    p_hr_less_0_8 = mean(hr < 0.8),
    p_hr_less_0_6 = mean(hr < 0.6),
    p_hr_less_0_5 = mean(hr < 0.5)
  )

hr_summary_efs_skeptical_landmark
# A tibble: 1 × 9
  mean_hr median_hr lower_89 upper_89 p_hr_less_1 p_hr_less_0_9 p_hr_less_0_8
    <dbl>     <dbl>    <dbl>    <dbl>       <dbl>         <dbl>         <dbl>
1   0.914     0.908    0.763     1.08       0.819         0.468         0.116
# ℹ 2 more variables: p_hr_less_0_6 <dbl>, p_hr_less_0_5 <dbl>
Show R Code 
create_hr_density_plot(
  posterior_hr_skeptical$hr,
  title = "EFS dose effect",
  subtitle = "Skeptical prior landmark model"
)

Show R Code 
ann_text_efs_skeptical_landmark <- paste0(
  "Median HR = ", round(hr_summary_efs_skeptical_landmark$median_hr, 2), "\n",
  "89% CrI: ", round(hr_summary_efs_skeptical_landmark$lower_89, 2),
  "\u2013", round(hr_summary_efs_skeptical_landmark$upper_89, 2), "\n\n",
  "P(HR < 1.0) = ", scales::percent(hr_summary_efs_skeptical_landmark$p_hr_less_1, accuracy = 0.1), "\n",
  "P(HR < 0.9) = ", scales::percent(hr_summary_efs_skeptical_landmark$p_hr_less_0_9, accuracy = 0.1), "\n",
  "P(HR < 0.8) = ", scales::percent(hr_summary_efs_skeptical_landmark$p_hr_less_0_8, accuracy = 0.1), "\n",
  "P(HR < 0.6) = ", scales::percent(hr_summary_efs_skeptical_landmark$p_hr_less_0_6, accuracy = 0.1)
)
Show R Code 
df_pred_2_skeptical <- df_efs_primary %>%
  mutate(
    dose_3plus = 0L,
    dose_group_landmark = factor("2 doses", levels = c("2 doses", "3+ doses"))
  )

df_pred_3plus_skeptical <- df_efs_primary %>%
  mutate(
    dose_3plus = 1L,
    dose_group_landmark = factor("3+ doses", levels = c("2 doses", "3+ doses"))
  )
Show R Code 
draws_skeptical <- as_draws_df(fit_efs_skeptical)

lp_2_skeptical <- posterior_linpred(
  fit_efs_skeptical,
  newdata = df_pred_2_skeptical,
  re_formula = NA
)

lp_3plus_skeptical <- posterior_linpred(
  fit_efs_skeptical,
  newdata = df_pred_3plus_skeptical,
  re_formula = NA
)

mu_2_skeptical <- exp(lp_2_skeptical)
mu_3plus_skeptical <- exp(lp_3plus_skeptical)

shape_draws_skeptical <- draws_skeptical$shape
time_grid_skeptical <- seq(0, 60, by = 1)
Show R Code 
surv_2_skeptical <- compute_survival(
  mu_mat = mu_2_skeptical,
  shape_vec = shape_draws_skeptical,
  times = time_grid_skeptical
)

surv_3plus_skeptical <- compute_survival(
  mu_mat = mu_3plus_skeptical,
  shape_vec = shape_draws_skeptical,
  times = time_grid_skeptical
)
Show R Code 
summarize_surv_with_time <- function(surv_mat, label, time_grid) {
  tibble(
    time = time_grid,
    median = apply(surv_mat, 2, median),
    lower = apply(surv_mat, 2, quantile, 0.055),
    upper = apply(surv_mat, 2, quantile, 0.945),
    group = label
  )
}

surv_df_efs_skeptical_landmark <- bind_rows(
  summarize_surv_with_time(surv_2_skeptical, "2 doses", time_grid_skeptical),
  summarize_surv_with_time(surv_3plus_skeptical, "3+ doses", time_grid_skeptical)
)
Show R Code 
ggplot(
  surv_df_efs_skeptical_landmark,
  aes(x = time, y = median, color = group, fill = group)
) +
  geom_ribbon(
    aes(ymin = lower, ymax = upper),
    alpha = 0.15,
    color = NA
  ) +
  geom_line(linewidth = 1.3) +
  annotate(
    "label",
    x = 36,
    y = 1.0,
    label = ann_text_efs_skeptical_landmark,
    hjust = 0,
    vjust = 1,
    size = 3.7,
    label.size = 0.25,
    fill = "white",
    color = "black"
  ) +
  scale_x_continuous(
    breaks = seq(0, 60, by = 12),
    limits = c(0, 60)
  ) +
  coord_cartesian(ylim = c(0, 1)) +
  labs(
    x = "Months from dose 2 landmark",
    y = "Predicted event-free survival",
    color = NULL,
    fill = NULL,
    title = "Posterior standardized EFS curves",
    subtitle = "Skeptical prior landmark Weibull model"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, face = "bold"),
    legend.position = "right"
  )

7 Sensitivity Analysis: Dose-Intensity–Adjusted Event-Free Survival Model

This sensitivity analysis extends the primary landmark Weibull model by adding a dose-intensity parameter. The goal is to assess whether the association between receipt of ≥3 doses and event-free survival remains after accounting for how consistently treatment was delivered over time.

Conceptually, this model separates two related but distinct ideas:

  • Cumulative exposure: whether a patient received 2 doses versus ≥3 doses
  • Delivered treatment intensity: whether treatment was administered on the expected schedule, summarized here as a centered dose-intensity measure

This distinction matters because a patient may receive multiple doses but still experience delays, interruptions, or schedule modifications that alter the biological intensity of treatment. By including both variables in the same model, this sensitivity analysis asks whether the primary dose-group association is robust to additional adjustment for treatment delivery patterns.

As in the primary analysis, the model is fit within the landmark cohort and uses a Weibull survival specification with right-censoring. The interpretation remains probabilistic and model-based, but the estimand is now conditioned not only on dose group and baseline covariates, but also on relative treatment intensity.

7.1 Construct the Dose-Intensity Modeling Dataset

We first restrict the landmark dataset to observations with complete information for the outcome, censoring indicator, primary exposure, dose-intensity variable, and required baseline covariates.

Show dose-intensity modeling dataset code
Show R Code 
df_efs_dose_intensity <- df_landmark_model %>%
  filter(
    !is.na(time_observed_landmark),
    !is.na(event_landmark),
    !is.na(cens),
    !is.na(dose_3plus),
    !is.na(dose_intensity_centered),
    !is.na(age_centered_10),
    !is.na(immunosuppressed),
    !is.na(ecog_2plus)
  )

7.2 Summarize the Dose-Intensity Modeling Dataset

The table below summarizes the resulting analytic dataset used for the dose-intensity sensitivity model.
Show dose-intensity modeling dataset summary
Show R Code 
df_efs_dose_intensity %>%
  summarise(
    n = n(),
    n_events = sum(event_landmark),
    event_rate = mean(event_landmark),
    median_followup = median(time_observed_landmark)
  )
# A tibble: 1 × 4
      n n_events event_rate median_followup
  <int>    <int>      <dbl>           <dbl>
1   177       72      0.407            9.72

7.3 Define the Dose-Intensity Weibull Model

We next define the Weibull regression formula for this sensitivity analysis. As in the primary model, bf() from brms specifies the Bayesian model formula, including the survival outcome, censoring structure, and covariates.

The key difference from the primary model is the addition of dose_intensity_centered, which captures variation in regimen delivery beyond the simple 2-dose versus ≥3-dose contrast.

Show dose-intensity Weibull model specification
Show R Code 
# Define the Bayesian survival model formula using brms
# This sensitivity model adds dose intensity to the primary landmark specification

form_efs_dose_intensity <- bf(
  # Survival outcome with right-censoring
  time_observed_landmark | cens(cens) ~

    # Primary exposure: 1 = ≥3 doses, 0 = 2 doses
    dose_3plus +

    # Additional treatment-delivery parameter:
    # centered measure of relative dose intensity
    dose_intensity_centered +

    # Demographic covariate
    age_centered_10 +

    # Clinical covariates
    immunosuppressed +
    ecog_2plus +

    # Disease stage indicators
    stage_III +
    stage_IV +
    stage_unstaged +

    # Treatment regimen indicators
    agent_pembro_q3 +
    agent_pembro_q6 +
    agent_pembro_mixed +
    agent_ipinivo
)

7.4 Inspect Prior Slots

Before fitting the model, we inspect the available prior slots to confirm the parameterization used by brms for the Weibull family and ensure that priors are mapped to the intended coefficients.

Show dose-intensity prior inspection code
Show R Code 
get_prior(
  formula = form_efs_dose_intensity,
  data = df_efs_dose_intensity,
  family = weibull()
)
                  prior     class                    coef group resp dpar nlpar
                 (flat)         b                                              
                 (flat)         b         age_centered_10                      
                 (flat)         b           agent_ipinivo                      
                 (flat)         b      agent_pembro_mixed                      
                 (flat)         b         agent_pembro_q3                      
                 (flat)         b         agent_pembro_q6                      
                 (flat)         b              dose_3plus                      
                 (flat)         b dose_intensity_centered                      
                 (flat)         b              ecog_2plus                      
                 (flat)         b        immunosuppressed                      
                 (flat)         b               stage_III                      
                 (flat)         b                stage_IV                      
                 (flat)         b          stage_unstaged                      
 student_t(3, 2.3, 2.5) Intercept                                              
      gamma(0.01, 0.01)     shape                                              
 lb ub       source
            default
       (vectorized)
       (vectorized)
       (vectorized)
       (vectorized)
       (vectorized)
       (vectorized)
       (vectorized)
       (vectorized)
       (vectorized)
       (vectorized)
       (vectorized)
       (vectorized)
            default
  0         default

7.5 Fit the Dose-Intensity Sensitivity Model

We then fit the Bayesian Weibull survival model using the dose-intensity-adjusted formula and the corresponding prior specification.

Show dose-intensity model fitting code
Show R Code 
fit_efs_dose_intensity <- brm(
  formula = form_efs_dose_intensity,
  data = df_efs_dose_intensity,
  family = weibull(),
  prior = priors_efs_dose_intensity,
  chains = 4,
  iter = 4000,
  warmup = 2000,
  cores = 4,
  seed = 2203,
  refresh = 0
)
Show R Code 
summary(fit_efs_dose_intensity)
 Family: weibull 
  Links: mu = log; shape = identity 
Formula: time_observed_landmark | cens(cens) ~ dose_3plus + dose_intensity_centered + age_centered_10 + immunosuppressed + ecog_2plus + stage_III + stage_IV + stage_unstaged + agent_pembro_q3 + agent_pembro_q6 + agent_pembro_mixed + agent_ipinivo 
   Data: df_efs_dose_intensity (Number of observations: 177) 
  Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
         total post-warmup draws = 8000

Regression Coefficients:
                        Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
Intercept                   4.25      0.43     3.42     5.12 1.00    14359
dose_3plus                  0.39      0.25    -0.10     0.87 1.00    17614
dose_intensity_centered    -0.00      0.29    -0.57     0.57 1.00    16459
age_centered_10            -0.18      0.14    -0.47     0.09 1.00    16682
immunosuppressed           -0.20      0.26    -0.70     0.31 1.00    16061
ecog_2plus                 -0.12      0.27    -0.65     0.42 1.00    18189
stage_III                   0.20      0.25    -0.30     0.68 1.00    14688
stage_IV                   -0.38      0.24    -0.86     0.12 1.00    14410
stage_unstaged              0.03      0.27    -0.48     0.56 1.00    17677
agent_pembro_q3            -0.33      0.24    -0.81     0.15 1.00    17101
agent_pembro_q6            -0.02      0.29    -0.58     0.54 1.00    16741
agent_pembro_mixed          0.08      0.28    -0.47     0.65 1.00    15678
agent_ipinivo               0.02      0.30    -0.56     0.60 1.00    17368
                        Tail_ESS
Intercept                   6186
dose_3plus                  6528
dose_intensity_centered     5741
age_centered_10             6166
immunosuppressed            6206
ecog_2plus                  5940
stage_III                   6282
stage_IV                    6477
stage_unstaged              6920
agent_pembro_q3             5841
agent_pembro_q6             5796
agent_pembro_mixed          5830
agent_ipinivo               5333

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape     0.62      0.05     0.54     0.72 1.00    11619     6831

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
Show R Code 
posterior_hr_dose_intensity <- extract_weibull_hr_draws(fit_efs_dose_intensity)

posterior_hr_efs_dose_intensity_landmark <- posterior_hr_dose_intensity
Show R Code 
hr_summary_efs_dose_intensity_landmark <- posterior_hr_efs_dose_intensity_landmark %>%
  summarise(
    mean_hr = mean(hr),
    median_hr = median(hr),
    lower_89 = quantile(hr, 0.055),
    upper_89 = quantile(hr, 0.945),
    p_hr_less_1 = mean(hr < 1),
    p_hr_less_0_9 = mean(hr < 0.9),
    p_hr_less_0_8 = mean(hr < 0.8),
    p_hr_less_0_6 = mean(hr < 0.6),
    p_hr_less_0_5 = mean(hr < 0.5)
  )

hr_summary_efs_dose_intensity_landmark
# A tibble: 1 × 9
  mean_hr median_hr lower_89 upper_89 p_hr_less_1 p_hr_less_0_9 p_hr_less_0_8
    <dbl>     <dbl>    <dbl>    <dbl>       <dbl>         <dbl>         <dbl>
1   0.794     0.785    0.609     1.01       0.940         0.809         0.551
# ℹ 2 more variables: p_hr_less_0_6 <dbl>, p_hr_less_0_5 <dbl>
Show R Code 
create_hr_density_plot(
  posterior_hr_dose_intensity$hr,
  title = "EFS dose effect",
  subtitle = "Dose-intensity adjusted landmark model"
)

Show R Code 
ann_text_efs_dose_intensity_landmark <- paste0(
  "Median HR = ", round(hr_summary_efs_dose_intensity_landmark$median_hr, 2), "\n",
  "89% CrI: ", round(hr_summary_efs_dose_intensity_landmark$lower_89, 2),
  "\u2013", round(hr_summary_efs_dose_intensity_landmark$upper_89, 2), "\n\n",
  "P(HR < 1.0) = ", scales::percent(hr_summary_efs_dose_intensity_landmark$p_hr_less_1, accuracy = 0.1), "\n",
  "P(HR < 0.9) = ", scales::percent(hr_summary_efs_dose_intensity_landmark$p_hr_less_0_9, accuracy = 0.1), "\n",
  "P(HR < 0.8) = ", scales::percent(hr_summary_efs_dose_intensity_landmark$p_hr_less_0_8, accuracy = 0.1), "\n",
  "P(HR < 0.6) = ", scales::percent(hr_summary_efs_dose_intensity_landmark$p_hr_less_0_6, accuracy = 0.1)
)
Show R Code 
df_pred_2_dose_intensity <- df_efs_dose_intensity %>%
  mutate(
    dose_3plus = 0L,
    dose_group_landmark = factor("2 doses", levels = c("2 doses", "3+ doses"))
  )

df_pred_3plus_dose_intensity <- df_efs_dose_intensity %>%
  mutate(
    dose_3plus = 1L,
    dose_group_landmark = factor("3+ doses", levels = c("2 doses", "3+ doses"))
  )
Show R Code 
draws_dose_intensity <- as_draws_df(fit_efs_dose_intensity)

lp_2_dose_intensity <- posterior_linpred(
  fit_efs_dose_intensity,
  newdata = df_pred_2_dose_intensity,
  re_formula = NA
)

lp_3plus_dose_intensity <- posterior_linpred(
  fit_efs_dose_intensity,
  newdata = df_pred_3plus_dose_intensity,
  re_formula = NA
)

mu_2_dose_intensity <- exp(lp_2_dose_intensity)
mu_3plus_dose_intensity <- exp(lp_3plus_dose_intensity)

shape_draws_dose_intensity <- draws_dose_intensity$shape
time_grid_dose_intensity <- seq(0, 60, by = 1)
Show R Code 
surv_2_dose_intensity <- compute_survival(
  mu_mat = mu_2_dose_intensity,
  shape_vec = shape_draws_dose_intensity,
  times = time_grid_dose_intensity
)

surv_3plus_dose_intensity <- compute_survival(
  mu_mat = mu_3plus_dose_intensity,
  shape_vec = shape_draws_dose_intensity,
  times = time_grid_dose_intensity
)
Show R Code 
surv_df_efs_dose_intensity_landmark <- bind_rows(
  summarize_surv_with_time(surv_2_dose_intensity, "2 doses", time_grid_dose_intensity),
  summarize_surv_with_time(surv_3plus_dose_intensity, "3+ doses", time_grid_dose_intensity)
)
Show R Code 
ggplot(
  surv_df_efs_dose_intensity_landmark,
  aes(x = time, y = median, color = group, fill = group)
) +
  geom_ribbon(
    aes(ymin = lower, ymax = upper),
    alpha = 0.15,
    color = NA
  ) +
  geom_line(linewidth = 1.3) +
  annotate(
    "label",
    x = 36,
    y = 1.0,
    label = ann_text_efs_dose_intensity_landmark,
    hjust = 0,
    vjust = 1,
    size = 3.7,
    label.size = 0.25,
    fill = "white",
    color = "black"
  ) +
  scale_x_continuous(
    breaks = seq(0, 60, by = 12),
    limits = c(0, 60)
  ) +
  coord_cartesian(ylim = c(0, 1)) +
  labs(
    x = "Months from dose 2 landmark",
    y = "Predicted event-free survival",
    color = NULL,
    fill = NULL,
    title = "Posterior standardized EFS curves",
    subtitle = "Dose-intensity adjusted landmark Weibull model"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, face = "bold"),
    legend.position = "right"
  )

8 Sensitivity Analysis: Surgery-Adjusted Event-Free Survival Model

Show R Code 
df_efs_surgery <- df_landmark_model %>%
  filter(
    !is.na(time_observed_landmark),
    !is.na(event_landmark),
    !is.na(cens),
    !is.na(dose_3plus),
    !is.na(surgery),
    !is.na(age_centered_10),
    !is.na(immunosuppressed),
    !is.na(ecog_2plus)
  )
Show R Code 
df_efs_surgery %>%
  summarise(
    n = n(),
    n_events = sum(event_landmark),
    event_rate = mean(event_landmark),
    median_followup = median(time_observed_landmark),
    prop_surgery = mean(surgery == 1)
  )
# A tibble: 1 × 5
      n n_events event_rate median_followup prop_surgery
  <int>    <int>      <dbl>           <dbl>        <dbl>
1   177       72      0.407            9.72        0.486
Show R Code 
df_efs_surgery %>%
  group_by(dose_group_landmark) %>%
  summarise(
    n = n(),
    n_events = sum(event_landmark),
    event_rate = mean(event_landmark),
    median_followup = median(time_observed_landmark),
    prop_surgery = mean(surgery == 1),
    .groups = "drop"
  )
# A tibble: 2 × 6
  dose_group_landmark     n n_events event_rate median_followup prop_surgery
  <fct>               <int>    <int>      <dbl>           <dbl>        <dbl>
1 2 doses                81       31      0.383            5.98        0.728
2 3+ doses               96       41      0.427           15.1         0.281
Show R Code 
form_efs_surgery <- bf(
  time_observed_landmark | cens(cens) ~
    dose_3plus +
    surgery +
    age_centered_10 +
    immunosuppressed +
    ecog_2plus +
    stage_III +
    stage_IV +
    stage_unstaged +
    agent_pembro_q3 +
    agent_pembro_q6 +
    agent_pembro_mixed +
    agent_ipinivo
)
Show R Code 
get_prior(
  formula = form_efs_surgery,
  data = df_efs_surgery,
  family = weibull()
)
                  prior     class               coef group resp dpar nlpar lb
                 (flat)         b                                            
                 (flat)         b    age_centered_10                         
                 (flat)         b      agent_ipinivo                         
                 (flat)         b agent_pembro_mixed                         
                 (flat)         b    agent_pembro_q3                         
                 (flat)         b    agent_pembro_q6                         
                 (flat)         b         dose_3plus                         
                 (flat)         b         ecog_2plus                         
                 (flat)         b   immunosuppressed                         
                 (flat)         b          stage_III                         
                 (flat)         b           stage_IV                         
                 (flat)         b     stage_unstaged                         
                 (flat)         b            surgery                         
 student_t(3, 2.3, 2.5) Intercept                                            
      gamma(0.01, 0.01)     shape                                           0
 ub       source
         default
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
    (vectorized)
         default
         default
Show R Code 
fit_efs_surgery <- brm(
  formula = form_efs_surgery,
  data = df_efs_surgery,
  family = weibull(),
  prior = priors_efs_surgery,
  chains = 4,
  iter = 4000,
  warmup = 2000,
  cores = 4,
  seed = 2204,
  refresh = 0
)
Show R Code 
summary(fit_efs_surgery)
 Family: weibull 
  Links: mu = log; shape = identity 
Formula: time_observed_landmark | cens(cens) ~ dose_3plus + surgery + age_centered_10 + immunosuppressed + ecog_2plus + stage_III + stage_IV + stage_unstaged + agent_pembro_q3 + agent_pembro_q6 + agent_pembro_mixed + agent_ipinivo 
   Data: df_efs_surgery (Number of observations: 177) 
  Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
         total post-warmup draws = 8000

Regression Coefficients:
                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept              4.32      0.45     3.45     5.22 1.00    14631     6482
dose_3plus             0.37      0.25    -0.14     0.86 1.00    17489     6660
surgery               -0.11      0.24    -0.58     0.36 1.00    15328     6697
age_centered_10       -0.20      0.15    -0.49     0.08 1.00    15449     6320
immunosuppressed      -0.20      0.26    -0.71     0.31 1.00    17210     5772
ecog_2plus            -0.12      0.26    -0.63     0.41 1.00    17759     5684
stage_III              0.20      0.25    -0.28     0.70 1.00    16522     6358
stage_IV              -0.39      0.24    -0.86     0.09 1.00    16005     6662
stage_unstaged         0.03      0.27    -0.52     0.56 1.00    16508     6431
agent_pembro_q3       -0.33      0.24    -0.79     0.15 1.00    18136     6300
agent_pembro_q6       -0.02      0.30    -0.61     0.56 1.00    17180     6222
agent_pembro_mixed     0.08      0.28    -0.46     0.63 1.00    18201     5931
agent_ipinivo          0.02      0.30    -0.56     0.61 1.00    16691     5854

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape     0.62      0.05     0.53     0.72 1.00    11928     7410

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
Show R Code 
posterior_hr_surgery <- extract_weibull_hr_draws(fit_efs_surgery)

posterior_hr_efs_surgery_landmark <- posterior_hr_surgery
Show R Code 
hr_summary_efs_surgery_landmark <- posterior_hr_efs_surgery_landmark %>%
  summarise(
    mean_hr = mean(hr),
    median_hr = median(hr),
    lower_89 = quantile(hr, 0.055),
    upper_89 = quantile(hr, 0.945),
    p_hr_less_1 = mean(hr < 1),
    p_hr_less_0_9 = mean(hr < 0.9),
    p_hr_less_0_8 = mean(hr < 0.8),
    p_hr_less_0_6 = mean(hr < 0.6),
    p_hr_less_0_5 = mean(hr < 0.5)
  )

hr_summary_efs_surgery_landmark
# A tibble: 1 × 9
  mean_hr median_hr lower_89 upper_89 p_hr_less_1 p_hr_less_0_9 p_hr_less_0_8
    <dbl>     <dbl>    <dbl>    <dbl>       <dbl>         <dbl>         <dbl>
1   0.805     0.796    0.616     1.03       0.922         0.783         0.514
# ℹ 2 more variables: p_hr_less_0_6 <dbl>, p_hr_less_0_5 <dbl>
Show R Code 
create_hr_density_plot(
  posterior_hr_surgery$hr,
  title = "EFS dose effect",
  subtitle = "Surgery-adjusted landmark model"
)

Show R Code 
ann_text_efs_surgery_landmark <- paste0(
  "Median HR = ", round(hr_summary_efs_surgery_landmark$median_hr, 2), "\n",
  "89% CrI: ", round(hr_summary_efs_surgery_landmark$lower_89, 2),
  "\u2013", round(hr_summary_efs_surgery_landmark$upper_89, 2), "\n\n",
  "P(HR < 1.0) = ", scales::percent(hr_summary_efs_surgery_landmark$p_hr_less_1, accuracy = 0.1), "\n",
  "P(HR < 0.9) = ", scales::percent(hr_summary_efs_surgery_landmark$p_hr_less_0_9, accuracy = 0.1), "\n",
  "P(HR < 0.8) = ", scales::percent(hr_summary_efs_surgery_landmark$p_hr_less_0_8, accuracy = 0.1), "\n",
  "P(HR < 0.6) = ", scales::percent(hr_summary_efs_surgery_landmark$p_hr_less_0_6, accuracy = 0.1)
)
Show R Code 
df_pred_2_surgery <- df_efs_surgery %>%
  mutate(
    dose_3plus = 0L,
    dose_group_landmark = factor("2 doses", levels = c("2 doses", "3+ doses"))
  )

df_pred_3plus_surgery <- df_efs_surgery %>%
  mutate(
    dose_3plus = 1L,
    dose_group_landmark = factor("3+ doses", levels = c("2 doses", "3+ doses"))
  )
Show R Code 
draws_surgery <- as_draws_df(fit_efs_surgery)

lp_2_surgery <- posterior_linpred(
  fit_efs_surgery,
  newdata = df_pred_2_surgery,
  re_formula = NA
)

lp_3plus_surgery <- posterior_linpred(
  fit_efs_surgery,
  newdata = df_pred_3plus_surgery,
  re_formula = NA
)

mu_2_surgery <- exp(lp_2_surgery)
mu_3plus_surgery <- exp(lp_3plus_surgery)

shape_draws_surgery <- draws_surgery$shape
time_grid_surgery <- seq(0, 60, by = 1)
Show R Code 
surv_2_surgery <- compute_survival(
  mu_mat = mu_2_surgery,
  shape_vec = shape_draws_surgery,
  times = time_grid_surgery
)

surv_3plus_surgery <- compute_survival(
  mu_mat = mu_3plus_surgery,
  shape_vec = shape_draws_surgery,
  times = time_grid_surgery
)
Show R Code 
surv_df_efs_surgery_landmark <- bind_rows(
  summarize_surv_with_time(surv_2_surgery, "2 doses", time_grid_surgery),
  summarize_surv_with_time(surv_3plus_surgery, "3+ doses", time_grid_surgery)
)
Show R Code 
ggplot(
  surv_df_efs_surgery_landmark,
  aes(x = time, y = median, color = group, fill = group)
) +
  geom_ribbon(
    aes(ymin = lower, ymax = upper),
    alpha = 0.15,
    color = NA
  ) +
  geom_line(linewidth = 1.3) +
  annotate(
    "label",
    x = 36,
    y = 1.0,
    label = ann_text_efs_surgery_landmark,
    hjust = 0,
    vjust = 1,
    size = 3.7,
    label.size = 0.25,
    fill = "white",
    color = "black"
  ) +
  scale_x_continuous(
    breaks = seq(0, 60, by = 12),
    limits = c(0, 60)
  ) +
  coord_cartesian(ylim = c(0, 1)) +
  labs(
    x = "Months from dose 2 landmark",
    y = "Predicted event-free survival",
    color = NULL,
    fill = NULL,
    title = "Posterior standardized EFS curves",
    subtitle = "Surgery-adjusted landmark Weibull model"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, face = "bold"),
    legend.position = "right"
  )

Show R Code 
surv_2_efs_surgery_landmark <- surv_2_surgery
surv_3plus_efs_surgery_landmark <- surv_3plus_surgery
time_grid_efs_surgery_landmark <- time_grid_surgery

9 Sensitivity Analysis: Non-Landmark (Treatment-Initiation) Model

As a sensitivity analysis, we fit a model without a landmark restriction, including all patients from treatment initiation. This model evaluates the association between receipt of 3 or more doses and event-free survival without conditioning on survival to a fixed time point. Because treatment exposure is defined post-baseline, this analysis is subject to immortal time bias and is interpreted as a prognostic association rather than a causal estimate.

Show R Code 
df_efs_no_landmark <- df_analysis %>%
  mutate(
    cens = 1L - event,

    age_centered_10 = age_centered / 10,

    immunosuppressed = as.integer(immunosuppressed),

    ecog_2plus = case_when(
      ecog_condensed == "ECOG 2-3" ~ 1L,
      ecog_condensed == "ECOG 0-1" ~ 0L,
      TRUE ~ NA_integer_
    ),

    stage_III = if_else(stage_condensed == "Stage III", 1L, 0L, missing = 0L),
    stage_IV = if_else(stage_condensed == "Stage IV", 1L, 0L, missing = 0L),
    stage_unstaged = if_else(stage_condensed == "Unstaged", 1L, 0L, missing = 0L),

    agent_pembro_q3 = if_else(agent_schedule == "pembro_q3", 1L, 0L, missing = 0L),
    agent_pembro_q6 = if_else(agent_schedule == "pembro_q6", 1L, 0L, missing = 0L),
    agent_pembro_mixed = if_else(agent_schedule == "pembro_mixed", 1L, 0L, missing = 0L),
    agent_ipinivo = if_else(agent_schedule == "Ipi-nivo", 1L, 0L, missing = 0L),

    dose_3plus = if_else(num_doses >= 3, 1L, 0L),

    dose_group = factor(
      if_else(dose_3plus == 1L, "3+ doses", "0–2 doses"),
      levels = c("0–2 doses", "3+ doses")
    ),

    surgery = as.integer(surgery)
  ) %>%
  filter(
    !is.na(time_observed),
    !is.na(event),
    !is.na(cens),
    !is.na(dose_3plus),
    !is.na(age_centered_10),
    !is.na(immunosuppressed),
    !is.na(ecog_2plus)
  )
Show R Code 
form_efs_no_landmark <- bf(
  time_observed | cens(cens) ~
    dose_3plus +
    surgery +
    age_centered_10 +
    immunosuppressed +
    ecog_2plus +
    stage_III +
    stage_IV +
    stage_unstaged +
    agent_pembro_q3 +
    agent_pembro_q6 +
    agent_pembro_mixed +
    agent_ipinivo
)
Show R Code 
fit_efs_no_landmark <- brm(
  formula = form_efs_no_landmark,
  data = df_efs_no_landmark,
  family = weibull(),
  prior = priors_efs_surgery,
  chains = 4,
  iter = 4000,
  warmup = 2000,
  cores = 4,
  seed = 2205,
  refresh = 0
)
Show R Code 
posterior_hr_no_landmark <- extract_weibull_hr_draws(fit_efs_no_landmark)

hr_summary_no_landmark <- posterior_hr_no_landmark %>%
  summarise(
    mean_hr = mean(hr),
    median_hr = median(hr),
    lower_89 = quantile(hr, 0.055),
    upper_89 = quantile(hr, 0.945),
    p_hr_less_1 = mean(hr < 1),
    p_hr_less_0_8 = mean(hr < 0.8),
    p_hr_less_0_6 = mean(hr < 0.6)
  )

hr_summary_no_landmark
# A tibble: 1 × 7
  mean_hr median_hr lower_89 upper_89 p_hr_less_1 p_hr_less_0_8 p_hr_less_0_6
    <dbl>     <dbl>    <dbl>    <dbl>       <dbl>         <dbl>         <dbl>
1   0.756     0.745    0.561    0.980       0.959         0.656         0.111
Show R Code 
df_pred_0_2 <- df_efs_no_landmark %>%
  mutate(dose_3plus = 0L)

df_pred_3plus <- df_efs_no_landmark %>%
  mutate(dose_3plus = 1L)
Show R Code 
draws_no_landmark <- as_draws_df(fit_efs_no_landmark)

lp_0_2 <- posterior_linpred(fit_efs_no_landmark, newdata = df_pred_0_2, re_formula = NA)
lp_3plus <- posterior_linpred(fit_efs_no_landmark, newdata = df_pred_3plus, re_formula = NA)

mu_0_2 <- exp(lp_0_2)
mu_3plus <- exp(lp_3plus)

shape <- draws_no_landmark$shape
time_grid <- seq(0, 60, by = 1)

surv_0_2 <- compute_survival(mu_0_2, shape, time_grid)
surv_3plus <- compute_survival(mu_3plus, shape, time_grid)

surv_df_no_landmark <- bind_rows(
  summarize_surv_with_time(surv_0_2, "0–2 doses", time_grid),
  summarize_surv_with_time(surv_3plus, "3+ doses", time_grid)
)
Show R Code 
summary(fit_efs_no_landmark)
 Family: weibull 
  Links: mu = log; shape = identity 
Formula: time_observed | cens(cens) ~ dose_3plus + surgery + age_centered_10 + immunosuppressed + ecog_2plus + stage_III + stage_IV + stage_unstaged + agent_pembro_q3 + agent_pembro_q6 + agent_pembro_mixed + agent_ipinivo 
   Data: df_efs_no_landmark (Number of observations: 189) 
  Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
         total post-warmup draws = 8000

Regression Coefficients:
                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept              4.25      0.43     3.43     5.11 1.00    14011     6266
dose_3plus             0.40      0.24    -0.05     0.86 1.00    15745     6697
surgery               -0.08      0.22    -0.53     0.35 1.00    14708     6807
age_centered_10       -0.19      0.13    -0.45     0.06 1.00    14616     5934
immunosuppressed      -0.30      0.24    -0.76     0.18 1.00    18369     6165
ecog_2plus            -0.13      0.25    -0.61     0.36 1.00    16215     6531
stage_III              0.19      0.24    -0.27     0.65 1.00    13992     6691
stage_IV              -0.45      0.23    -0.92     0.00 1.00    13317     6792
stage_unstaged         0.05      0.26    -0.46     0.56 1.00    15990     6322
agent_pembro_q3       -0.32      0.23    -0.77     0.12 1.00    16495     5707
agent_pembro_q6       -0.01      0.29    -0.57     0.55 1.00    16553     5837
agent_pembro_mixed     0.09      0.27    -0.42     0.61 1.00    19568     6164
agent_ipinivo          0.02      0.29    -0.55     0.61 1.00    16951     6022

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape     0.73      0.06     0.62     0.85 1.00    11126     6558

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).
Show R Code 
hr_summary_efs_no_landmark <- posterior_hr_efs_no_landmark %>%
  summarise(
    mean_hr = mean(hr),
    median_hr = median(hr),
    lower_89 = quantile(hr, 0.055),
    upper_89 = quantile(hr, 0.945),
    p_hr_less_1 = mean(hr < 1),
    p_hr_less_0_9 = mean(hr < 0.9),
    p_hr_less_0_8 = mean(hr < 0.8),
    p_hr_less_0_6 = mean(hr < 0.6),
    p_hr_less_0_5 = mean(hr < 0.5)
  )

hr_summary_efs_no_landmark
# A tibble: 1 × 9
  mean_hr median_hr lower_89 upper_89 p_hr_less_1 p_hr_less_0_9 p_hr_less_0_8
    <dbl>     <dbl>    <dbl>    <dbl>       <dbl>         <dbl>         <dbl>
1   0.756     0.745    0.561    0.980       0.959         0.862         0.656
# ℹ 2 more variables: p_hr_less_0_6 <dbl>, p_hr_less_0_5 <dbl>
Show R Code 
create_hr_density_plot(
  posterior_hr_no_landmark$hr,
  title = "EFS dose effect",
  subtitle = "No-landmark Weibull model",
  x_label = "Hazard ratio (3+ doses vs 0–2 doses)"
)

Show R Code 
ann_text_efs_no_landmark <- paste0(
  "Median HR = ", round(hr_summary_efs_no_landmark$median_hr, 2), "\n",
  "89% CrI: ", round(hr_summary_efs_no_landmark$lower_89, 2),
  "\u2013", round(hr_summary_efs_no_landmark$upper_89, 2), "\n\n",
  "P(HR < 1.0) = ", scales::percent(hr_summary_efs_no_landmark$p_hr_less_1, accuracy = 0.1), "\n",
  "P(HR < 0.9) = ", scales::percent(hr_summary_efs_no_landmark$p_hr_less_0_9, accuracy = 0.1), "\n",
  "P(HR < 0.8) = ", scales::percent(hr_summary_efs_no_landmark$p_hr_less_0_8, accuracy = 0.1), "\n",
  "P(HR < 0.6) = ", scales::percent(hr_summary_efs_no_landmark$p_hr_less_0_6, accuracy = 0.1)
)
Show R Code 
ggplot(
  surv_df_efs_no_landmark,
  aes(x = time, y = median, color = group, fill = group)
) +
  geom_ribbon(
    aes(ymin = lower, ymax = upper),
    alpha = 0.15,
    color = NA
  ) +
  geom_line(linewidth = 1.3) +
  annotate(
    "label",
    x = 36,
    y = 1.0,
    label = ann_text_efs_no_landmark,
    hjust = 0,
    vjust = 1,
    size = 3.7,
    label.size = 0.25,
    fill = "white",
    color = "black"
  ) +
  scale_x_continuous(
    breaks = seq(0, 60, by = 12),
    limits = c(0, 60)
  ) +
  coord_cartesian(ylim = c(0, 1)) +
  labs(
    x = "Months from treatment initiation",
    y = "Predicted event-free survival",
    color = NULL,
    fill = NULL,
    title = "Posterior standardized EFS curves",
    subtitle = "No-landmark Weibull model"
  ) +
  theme_minimal(base_size = 14) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, face = "bold"),
    legend.position = "right"
  )

Posterior standardized event-free survival curves from the no-landmark Weibull model. Survival curves represent posterior predictions from the fitted Bayesian Weibull regression model after standardizing over the observed covariate distribution in the full study cohort. For each posterior draw, predicted survival was computed after assigning all patients to either 0–2 doses or 3 or more doses of immunotherapy while leaving all other covariates unchanged. Because cumulative dose is defined after treatment initiation, this no-landmark analysis is vulnerable to immortal time bias and should be interpreted as a prognostic association rather than a causal estimate. Shaded regions denote 89% credible intervals. The annotation box summarizes the posterior hazard-ratio distribution comparing ≥3 doses with 0–2 doses.
Show R Code 
surv_0_2_efs_no_landmark <- surv_0_2
surv_3plus_efs_no_landmark <- surv_3plus
time_grid_efs_no_landmark <- time_grid
Show R Code 
df_efs_no_landmark %>%
  summarise(
    n = n(),
    n_events = sum(event),
    event_rate = mean(event),
    median_followup = median(time_observed)
  )
# A tibble: 1 × 4
      n n_events event_rate median_followup
  <int>    <int>      <dbl>           <dbl>
1   189       79      0.418            10.2
Show R Code 
df_efs_no_landmark %>%
  group_by(dose_group) %>%
  summarise(
    n = n(),
    n_events = sum(event),
    event_rate = mean(event),
    median_followup = median(time_observed),
    .groups = "drop"
  )
# A tibble: 2 × 5
  dose_group     n n_events event_rate median_followup
  <fct>      <int>    <int>      <dbl>           <dbl>
1 0–2 doses     92       37      0.402            6.62
2 3+ doses      97       42      0.433           15.7

10 Conclusion

This modeling framework provides a coherent approach to evaluating the relationship between treatment exposure and event-free survival using a Bayesian generative paradigm.

Across the primary landmark analysis and multiple sensitivity specifications, the results consistently frame the dose–outcome relationship in probabilistic and clinically interpretable terms, including:

  • posterior distributions of hazard ratios,
  • probabilities of benefit across clinically meaningful thresholds,
  • and standardized survival curves that reflect the observed patient population.

Several key insights emerge from this approach:

  • Model-based standardization matters
    By averaging predictions over the empirical covariate distribution, the resulting survival curves reflect expected outcomes at the population level rather than for an arbitrary reference patient.

  • Inference is gradient, not binary
    The posterior framework allows direct evaluation of the probability of benefit, avoiding over-reliance on threshold-based interpretations.

  • Design and modeling are inseparable
    The landmark analysis plays a central role in aligning the statistical model with the clinical question, highlighting that causal interpretability is driven as much by study design as by model choice.

  • Sensitivity analyses contextualize rather than replace the primary model
    Variations in priors, covariate adjustment, and cohort definition provide a structured way to assess robustness while preserving interpretability of the primary estimand.

The no-landmark analysis, in particular, illustrates how estimates can shift when temporal alignment is relaxed, reinforcing the importance of design-based control of bias when interpreting treatment effects.

Taken together, this document demonstrates how Bayesian generative modeling can be used not only to estimate survival outcomes, but to clarify the relationship between data, assumptions, and clinical interpretation. The goal is not simply to produce a single estimate, but to provide a framework through which clinicians and researchers can reason more transparently about treatment effects in real-world data.