Generative Modeling and Prior Predictive Simulation for Event-Free Survival

1 Overview

This document develops the generative workflow for the landmark event-free survival (EFS) analysis. The goal is to evaluate whether the proposed priors imply plausible survival behavior before fitting the observed data.

The workflow mirrors the prior-predictive strategy used for the response-rate analysis, while adapting it to time-to-event outcomes.

The model framework is:

• Primary model: weakly informative prior
• Sensitivity analysis 1: skeptical prior
• Sensitivity analysis 2: adjustment for treatment-span dose intensity
• Sensitivity analysis 3: adjustment for surgery

We begin by simulating a clinically plausible cohort, then generate landmark exposure groups, treatment delivery patterns, and event times. We then draw directly from the priors to assess the implied hazard ratios, event rates, and survival curves.

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

library(tidyverse)
library(brms)
library(survival)
library(cowplot)
library(scales)
library(posterior)
library(tidybayes)

2 Simulate Baseline Cohort

We first generate a synthetic cohort broadly resembling the observed study population.

set.seed(123)
n <- 188

sim_baseline <- tibble(
  patient_id = 1:n,
  age = rnorm(n, mean = 76, sd = 12),
  sex = sample(c("Male", "Female"), n, replace = TRUE, prob = c(0.76, 0.24))
) |>
  mutate(
    immunosupp_type = sample(
      c("None", "Chronic", "Hematologic", "HIV", "Chronic+Hematologic"),
      n,
      replace = TRUE,
      prob = c(0.83, 0.06, 0.082, 0.005, 0.022)
    ),
    immunosuppressed = if_else(immunosupp_type == "None", 0L, 1L),
    location = sample(
      c(
        "Head/Neck", "Conjunctiva", "Genitalia", "Lower Extremity",
        "Upper Extremity", "Trunk", "Unknown Primary"
      ),
      n,
      replace = TRUE,
      prob = c(0.86, 0.033, 0.005, 0.016, 0.022, 0.027, 0.033)
    )
  ) |>
  rowwise() |>
  mutate(
    stage = if_else(
      location == "Conjunctiva",
      sample(c("II", "III", "Unstaged"), 1, prob = c(0.3, 0.5, 0.2)),
      sample(c("I", "II", "III", "IV", "Unstaged"), 1,
             prob = c(0.005, 0.038, 0.40, 0.45, 0.11))
    )
  ) |>
  ungroup() |>
  mutate(
    ps_score =
      scale(age)[, 1] * 0.3 +
      if_else(immunosuppressed == 1L, 0.5, 0) +
      case_when(
        stage == "I" ~ -0.5,
        stage == "II" ~ -0.3,
        stage == "III" ~ 0,
        stage == "IV" ~ 0.6,
        TRUE ~ 0
      ) +
      rnorm(n, 0, 0.8),
    ecog = case_when(
      ps_score < -0.5 ~ 0L,
      ps_score < 0.5 ~ 1L,
      ps_score < 1.5 ~ 2L,
      TRUE ~ 3L
    ),
    year = sample(2018:2023, n, replace = TRUE)
  )

sim_baseline <- sim_baseline |>
  mutate(
    agent_schedule = sample(
      c("cemiplimab", "pembro_q3", "pembro_q6", "pembro_mixed", "Ipi-nivo"),
      n,
      replace = TRUE,
      prob = c(0.62, 0.18, 0.10, 0.08, 0.02)
    )
  )

sim_baseline |>
  summarise(
    mean_age = mean(age),
    prop_immunosuppressed = mean(immunosuppressed),
    prop_stage_iv = mean(stage == "IV"),
    prop_ecog_2plus = mean(ecog >= 2)
  )

# A tibble: 1 × 4
  mean_age prop_immunosuppressed prop_stage_iv prop_ecog_2plus
     <dbl>                 <dbl>         <dbl>           <dbl>
1     76.0                 0.197         0.463           0.426

2.1 Simulate Landmark Exposure Groups

The landmark analysis begins after 2 doses. We therefore simulate whether a patient remains in the reference group of exactly 2 doses or crosses into the 3+ dose group.

sim_baseline <- sim_baseline |>
  mutate(
    p_3plus = plogis(
      0.5 +
        -0.5 * (ecog >= 2) +
        -0.3 * (stage == "IV") +
        0.25 * (agent_schedule == "cemiplimab") +
        0.10 * (agent_schedule == "pembro_q3") +
        0.05 * (agent_schedule == "pembro_mixed") +
        0.00 * (agent_schedule == "pembro_q6") +
        0.00 * (agent_schedule == "Ipi-nivo") +
        0.2 * (year - 2020)
    ),
    dose_3plus = rbinom(n, size = 1, prob = p_3plus),
    dose_group = if_else(dose_3plus == 1L, "3+ doses", "2 doses")
  )

table(sim_baseline$dose_group)

 2 doses 3+ doses 
      80      108 

2.2 Simulate Cumulative Dose and Treatment-Span Dose Intensity

For the dose-intensity sensitivity analysis, we generate a treatment-span pacing variable that reflects how tightly or loosely treatment was delivered across the observed treatment course.

sim_baseline <- sim_baseline |>
  mutate(
    num_doses = if_else(
      dose_group == "2 doses",
      2L,
      sample(
        3:10,
        size = n(),
        replace = TRUE,
        prob = c(0.30, 0.24, 0.16, 0.10, 0.08, 0.06, 0.04, 0.02)
      )
    ),
    expected_days_per_dose = case_when(
      agent_schedule == "cemiplimab"   ~ 21,
      agent_schedule == "pembro_q3"    ~ 21,
      agent_schedule == "pembro_q6"    ~ 42,
      agent_schedule == "pembro_mixed" ~ 31.5,
      agent_schedule == "Ipi-nivo"     ~ 21,
      TRUE ~ 21
    ),
    schedule_stretch = exp(
      rnorm(
        n(),
        mean =
          0.10 * (ecog >= 2) +
          0.08 * immunosuppressed +
          0.08 * (stage == "IV"),
        sd = 0.20
      )
    ),
    treatment_span_days = pmax(
      21,
      round((num_doses - 1) * expected_days_per_dose * schedule_stretch)
    ),
    dose_intensity = num_doses / (treatment_span_days / 30.4),
    dose_intensity_centered = as.numeric(scale(dose_intensity, center = TRUE, scale = FALSE))
  )

sim_baseline |>
  summarise(
    mean_num_doses = mean(num_doses),
    mean_treatment_span_days = mean(treatment_span_days),
    mean_dose_intensity = mean(dose_intensity)
  )

# A tibble: 1 × 3
  mean_num_doses mean_treatment_span_days mean_dose_intensity
           <dbl>                    <dbl>               <dbl>
1           3.49                     64.3                1.96

2.3 Simulate Surgery

Surgery is simulated after exposure assignment. It is retained for later sensitivity analysis.

sim_baseline <- sim_baseline |>
  mutate(
    surgery_prob = plogis(
      -0.5 +
        -1.0 * (dose_group == "3+ doses") +
        -1.0 * (stage == "IV") +
        0.8 * (location == "Head/Neck") +
        -0.15 * (year - 2020)
    ),
    surgery = rbinom(n, size = 1, prob = surgery_prob)
  )

sim_baseline |>
  group_by(dose_group) |>
  summarise(prop_surgery = mean(surgery), .groups = "drop")

# A tibble: 2 × 2
  dose_group prop_surgery
  <chr>             <dbl>
1 2 doses           0.462
2 3+ doses          0.269

2.4 Create Modeling Variables

We now condense the covariates that will be used in the EFS models.

sim_baseline <- sim_baseline |>
  mutate(
    age_centered_10 = (age - 76) / 10,
    ecog_2plus = if_else(ecog >= 2, 1L, 0L),
    stage_III = if_else(stage == "III", 1L, 0L),
    stage_IV = if_else(stage == "IV", 1L, 0L),
    stage_unstaged = if_else(stage == "Unstaged", 1L, 0L),
    agent_pembro_q3 = if_else(agent_schedule == "pembro_q3", 1L, 0L),
    agent_pembro_q6 = if_else(agent_schedule == "pembro_q6", 1L, 0L),
    agent_pembro_mixed = if_else(agent_schedule == "pembro_mixed", 1L, 0L),
    agent_ipinivo    = if_else(agent_schedule == "Ipi-nivo", 1L, 0L)
  )

2.5 Define the Data-Generating Process for EFS

We simulate event times from a Weibull model on the landmark time scale. The primary truth here is only used to generate realistic synthetic datasets for workflow development. It is separate from the prior predictive checks.

simulate_efs_data <- function(
    dat,
    true_intercept = -5.25,
    true_log_hr_dose = log(0.85),
    true_log_hr_dose_intensity = log(0.95),
    true_log_hr_surgery = log(0.80),
    shape = 1.25,
    censor_time = 24,
    seed = 101
) {

  set.seed(seed)

  lp <- with(
    dat,
    true_intercept +
      true_log_hr_dose * dose_3plus +
      true_log_hr_dose_intensity * dose_intensity_centered +
      true_log_hr_surgery * surgery +
      log(1.10) * age_centered_10 +
      log(1.20) * immunosuppressed +
      log(1.25) * ecog_2plus +
      log(1.20) * stage_III +
      log(1.50) * stage_IV +
      log(1.10) * stage_unstaged +
      log(1.00) * agent_pembro_q3 +
      log(1.05) * agent_pembro_q6 +
      log(1.03) * agent_pembro_mixed +
      log(1.10) * agent_ipinivo
        )

  lambda <- exp(lp)
  u <- runif(nrow(dat))

  time_event <- (-log(u) / lambda)^(1 / shape)
  time_observed <- pmin(time_event, censor_time)
  event <- if_else(time_event <= censor_time, 1L, 0L)

  dat |>
    mutate(
      time_event = time_event,
      time_observed = time_observed,
      event = event
    )
}

sim_efs <- simulate_efs_data(sim_baseline)

sim_efs |>
  group_by(dose_group) |>
  summarise(
    event_rate = mean(event),
    median_followup = median(time_observed),
    .groups = "drop"
  )

# A tibble: 2 × 3
  dose_group event_rate median_followup
  <chr>           <dbl>           <dbl>
1 2 doses         0.412              24
2 3+ doses        0.370              24

sim_efs_brms <- sim_efs |>
  mutate(
    cens = 1L - event
  )

reference_shape <- 1.25

beta_from_hr <- function(hr, shape = reference_shape) {
  -log(hr) / shape
}

primary_prior_hr_center <- 0.85
skeptical_prior_hr_center <- 1.00

primary_prior_beta_center <- beta_from_hr(primary_prior_hr_center)
skeptical_prior_beta_center <- beta_from_hr(skeptical_prior_hr_center)

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")
  )
}

priors_efs_primary <- make_efs_priors(
  mean_dose = primary_prior_beta_center,
  sd_dose = 0.35
)

priors_efs_skeptical <- make_efs_priors(
  mean_dose = skeptical_prior_beta_center,
  sd_dose = 0.20
)

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")
)

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")
)

tibble(
  prior = c("Primary", "Skeptical"),
  target_hr_center = c(primary_prior_hr_center, skeptical_prior_hr_center),
  beta_center = c(primary_prior_beta_center, skeptical_prior_beta_center)
)

# A tibble: 2 × 3
  prior     target_hr_center beta_center
  <chr>                <dbl>       <dbl>
1 Primary               0.85       0.130
2 Skeptical             1          0    

For the primary landmark model, we specified a weakly informative prior centered on a small protective effect of receiving 3 or more doses, corresponding approximately to a hazard ratio of 0.85 on the reference Weibull scale. This reflects the clinical belief that greater cumulative exposure may modestly improve event-free survival, while still allowing substantial uncertainty and retaining support for null or harmful effects. In contrast, the skeptical prior was centered at the null (hazard ratio 1.0) and assigned a narrower standard deviation to more strongly shrink estimates toward no association.

set.seed(321)

n_prior_draws <- 4000

prior_dose_check <- tibble(
  draw = 1:n_prior_draws,
  shape = rgamma(n_prior_draws, shape = 5, rate = 4),
  beta_primary = rnorm(n_prior_draws, mean = primary_prior_beta_center, sd = 0.35),
  beta_skeptical = rnorm(n_prior_draws, mean = skeptical_prior_beta_center, sd = 0.20)
) |>
  mutate(
    hr_primary = exp(-shape * beta_primary),
    hr_skeptical = exp(-shape * beta_skeptical)
  )

prior_dose_check |>
  summarise(
    median_hr_primary = median(hr_primary),
    lower_89_primary = quantile(hr_primary, 0.055),
    upper_89_primary = quantile(hr_primary, 0.945),
    p_hr_less_1_primary = mean(hr_primary < 1),
    median_hr_skeptical = median(hr_skeptical),
    lower_89_skeptical = quantile(hr_skeptical, 0.055),
    upper_89_skeptical = quantile(hr_skeptical, 0.945),
    p_hr_less_1_skeptical = mean(hr_skeptical < 1)
  )

# A tibble: 1 × 8
  median_hr_primary lower_89_primary upper_89_primary p_hr_less_1_primary
              <dbl>            <dbl>            <dbl>               <dbl>
1             0.881            0.384             1.73                0.64
# ℹ 4 more variables: median_hr_skeptical <dbl>, lower_89_skeptical <dbl>,
#   upper_89_skeptical <dbl>, p_hr_less_1_skeptical <dbl>

Add a visual prior-only HR Comparison

prior_dose_check |>
  select(hr_primary, hr_skeptical) |>
  pivot_longer(
    everything(),
    names_to = "prior",
    values_to = "hr"
  ) |>
  mutate(
    prior = recode(
      prior,
      hr_primary = "Primary weakly informative prior",
      hr_skeptical = "Skeptical prior"
    )
  ) |>
  ggplot(aes(x = hr, fill = prior)) +
  geom_density(alpha = 0.45) +
  geom_vline(xintercept = 1, linetype = "dashed") +
  coord_cartesian(xlim = c(0, 2.5)) +
  labs(
    x = "Implied hazard ratio for 3+ doses vs 2 doses",
    y = "Density",
    title = "Direct prior implications for the dose effect",
    subtitle = "Primary prior centered on modest benefit; skeptical prior centered on the null"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5),
    legend.position = "bottom"
  )

form_efs_primary <- bf(
  time_observed | 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
)

form_efs_dose_intensity <- bf(
  time_observed | 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
)

form_efs_surgery <- 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
)

get_prior(
  formula = form_efs_primary,
  data = sim_efs_brms,
  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, 3.2, 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

priors_efs_primary

              prior     class       coef group resp dpar nlpar   lb   ub source
  normal(3.0, 0.35) Intercept                                  <NA> <NA>   user
 normal(0.13, 0.35)         b dose_3plus                       <NA> <NA>   user
    normal(0, 0.30)         b                                  <NA> <NA>   user
        gamma(5, 4)     shape                                  <NA> <NA>   user

fit_prior_only_weibull <- function(formula, data, prior, seed) {
  brm(
    formula = formula,
    data = data,
    family = weibull(),
    prior = prior,
    sample_prior = "only",
    chains = 4,
    iter = 2000,
    warmup = 1000,
    cores = 4,
    seed = seed,
    refresh = 0
  )
}

fit_prior_primary <- fit_prior_only_weibull(
  formula = form_efs_primary,
  data = sim_efs_brms,
  prior = priors_efs_primary,
  seed = 123
)

fit_prior_skeptical <- fit_prior_only_weibull(
  formula = form_efs_primary,
  data = sim_efs_brms,
  prior = priors_efs_skeptical,
  seed = 124
)

fit_prior_di <- fit_prior_only_weibull(
  formula = form_efs_dose_intensity,
  data = sim_efs_brms,
  prior = priors_efs_dose_intensity,
  seed = 125
)

fit_prior_surgery <- fit_prior_only_weibull(
  formula = form_efs_surgery,
  data = sim_efs_brms,
  prior = priors_efs_surgery,
  seed = 126
)

simulate_prior_predictive_brms <- function(
  fit,
  data,
  model_label,
  prior_label,
  dose_var = "dose_3plus",
  censor_time = 24
) {

  mu_draws <- posterior_linpred(
    fit,
    newdata = data,
    transform = TRUE,
    re_formula = NA
  )

  draws_df <- as_draws_df(fit)

  shape_var <- names(draws_df)[str_detect(names(draws_df), "^shape$|^b_shape$")]
  if (length(shape_var) != 1) {
    stop("Could not uniquely identify the shape parameter in the brms draws.")
  }

  beta_dose_var <- names(draws_df)[str_detect(names(draws_df), "^b_dose_3plus$")]
  if (length(beta_dose_var) != 1) {
    stop("Could not uniquely identify b_dose_3plus in the brms draws.")
  }

  shape_draws <- draws_df[[shape_var]]
  beta_dose_draws <- draws_df[[beta_dose_var]]

  n_draws <- nrow(mu_draws)
  n_obs <- ncol(mu_draws)

  out <- vector("list", n_draws)

  dose_ref <- data[[dose_var]] == 0
  dose_exp <- data[[dose_var]] == 1

  for (d in seq_len(n_draws)) {
    mu_d <- mu_draws[d, ]
    shape_d <- shape_draws[d]

    time_event <- rweibull(
      n = n_obs,
      shape = shape_d,
      scale = mu_d
    )

    time_obs <- pmin(time_event, censor_time)
    event_ind <- if_else(time_event <= censor_time, 1L, 0L)

    hr_dose <- exp(-shape_d * beta_dose_draws[d])

    out[[d]] <- tibble(
      draw = d,
      prior = prior_label,
      model = model_label,
      shape = shape_d,
      beta_dose = beta_dose_draws[d],
      hr_dose = hr_dose,
      event_rate = mean(event_ind),
      median_followup = median(time_obs),
      event_rate_2dose = mean(event_ind[dose_ref]),
      event_rate_3plus = mean(event_ind[dose_exp]),
      median_time_2dose = median(time_obs[dose_ref]),
      median_time_3plus = median(time_obs[dose_exp])
    )
  }

  bind_rows(out)
}

prior_pred_primary <- simulate_prior_predictive_brms(
  fit = fit_prior_primary,
  data = sim_efs_brms,
  model_label = "Primary landmark model",
  prior_label = "Weakly informative"
)

prior_pred_skeptical <- simulate_prior_predictive_brms(
  fit = fit_prior_skeptical,
  data = sim_efs_brms,
  model_label = "Primary landmark model",
  prior_label = "Skeptical"
)

prior_pred_di <- simulate_prior_predictive_brms(
  fit = fit_prior_di,
  data = sim_efs_brms,
  model_label = "Dose-intensity adjusted model",
  prior_label = "Weakly informative"
)

prior_pred_surgery <- simulate_prior_predictive_brms(
  fit = fit_prior_surgery,
  data = sim_efs_brms,
  model_label = "Surgery-adjusted model",
  prior_label = "Weakly informative"
)

prior_pred_all <- bind_rows(
  prior_pred_primary,
  prior_pred_skeptical,
  prior_pred_di,
  prior_pred_surgery
)

prior_pred_all |>
  group_by(model, prior) |>
  summarise(
    median_hr = median(hr_dose),
    lower_95_hr = quantile(hr_dose, 0.025),
    upper_95_hr = quantile(hr_dose, 0.975),
    median_event_rate = median(event_rate),
    lower_90_event_rate = quantile(event_rate, 0.05),
    upper_90_event_rate = quantile(event_rate, 0.95),
    median_followup = median(median_followup),
    .groups = "drop"
  )

# A tibble: 4 × 9
  model                prior median_hr lower_95_hr upper_95_hr median_event_rate
  <chr>                <chr>     <dbl>       <dbl>       <dbl>             <dbl>
1 Dose-intensity adju… Weak…     0.873       0.271        2.12             0.686
2 Primary landmark mo… Skep…     0.997       0.544        1.75             0.691
3 Primary landmark mo… Weak…     0.875       0.292        2.09             0.686
4 Surgery-adjusted mo… Weak…     0.881       0.312        2.04             0.681
# ℹ 3 more variables: lower_90_event_rate <dbl>, upper_90_event_rate <dbl>,
#   median_followup <dbl>

prior_pred_all |>
  filter(model %in% c("Primary landmark model",
                      "Dose-intensity adjusted model",
                      "Surgery-adjusted model")) |>
  ggplot(aes(x = hr_dose, fill = model)) +
  geom_density(alpha = 0.45) +
  geom_vline(xintercept = 1, linetype = "dashed") +
  coord_cartesian(xlim = c(0, 2.5)) +
  facet_wrap(~prior) +
  labs(
    x = "Implied hazard ratio for 3+ doses vs 2 doses",
    y = "Density",
    title = "Prior predictive hazard-ratio distributions across model formulations"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    legend.position = "bottom"
  )

prior_pred_all |>
  ggplot(aes(x = event_rate, fill = prior)) +
  geom_density(alpha = 0.5) +
  scale_x_continuous(labels = percent_format(accuracy = 1)) +
  labs(
    x = "Implied 24-month event rate",
    y = "Density",
    title = "Prior predictive distribution of the overall event rate"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    legend.position = "bottom"
  )

simulate_one_prior_dataset_brms <- function(
  fit,
  data,
  model_label,
  dose_var = "dose_3plus",
  censor_time = 24,
  draw_id = 1
) {

  mu_draws <- posterior_linpred(
    fit,
    newdata = data,
    transform = TRUE,
    re_formula = NA
  )

  draws_df <- as_draws_df(fit)

  shape_var <- names(draws_df)[str_detect(names(draws_df), "^shape$|^b_shape$")]
  beta_dose_var <- names(draws_df)[str_detect(names(draws_df), "^b_dose_3plus$")]

  shape_d <- draws_df[[shape_var]][draw_id]
  beta_dose_d <- draws_df[[beta_dose_var]][draw_id]
  mu_d <- mu_draws[draw_id, ]

  time_event <- rweibull(
    n = nrow(data),
    shape = shape_d,
    scale = mu_d
  )

  tibble(
    patient_id = data$patient_id,
    model = model_label,
    dose_group = if_else(data[[dose_var]] == 1, "3+ doses", "2 doses"),
    time = pmin(time_event, censor_time),
    event = if_else(time_event <= censor_time, 1L, 0L),
    shape = shape_d,
    beta_dose = beta_dose_d,
    hr_dose = exp(-shape_d * beta_dose_d)
  )
}

km_examples <- bind_rows(
  simulate_one_prior_dataset_brms(
    fit = fit_prior_primary,
    data = sim_efs_brms,
    model_label = "Primary landmark model",
    draw_id = 25
  ),
  simulate_one_prior_dataset_brms(
    fit = fit_prior_skeptical,
    data = sim_efs_brms,
    model_label = "Skeptical prior model",
    draw_id = 25
  ),
  simulate_one_prior_dataset_brms(
    fit = fit_prior_di,
    data = sim_efs_brms,
    model_label = "Dose-intensity adjusted model",
    draw_id = 25
  ),
  simulate_one_prior_dataset_brms(
    fit = fit_prior_surgery,
    data = sim_efs_brms,
    model_label = "Surgery-adjusted model",
    draw_id = 25
  )
)

make_km_tidy <- function(dat) {
  sf <- survfit(Surv(time, event) ~ dose_group, data = dat)

  tibble(
    time = sf$time,
    surv = sf$surv,
    strata = rep(names(sf$strata), sf$strata)
  ) |>
    mutate(
      dose_group = str_remove(strata, "dose_group="),
      model = dat$model[1]
    )
}

km_tidy <- km_examples |>
  group_split(model) |>
  map_dfr(make_km_tidy)

ggplot(km_tidy, aes(x = time, y = surv, color = dose_group)) +
  geom_step(linewidth = 0.9) +
  facet_wrap(~model) +
  coord_cartesian(xlim = c(0, 24), ylim = c(0, 1)) +
  labs(
    x = "Months from landmark",
    y = "Event-free survival probability",
    color = NULL,
    title = "Representative prior predictive Kaplan-Meier curves"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    legend.position = "bottom"
  )

2.6 Define the true parameters on the brms scale

In brms Weibull:

log(scale_i) = eta_i

HR = exp(-shape * beta)

So if you want a true hazard ratio of 0.85 for dose_3plus, the corresponding coefficient is: βdose = −log(0.85)/shape

true_shape <- 1.25

true_hr_dose <- 0.85
true_beta_dose <- -log(true_hr_dose) / true_shape

true_hr_age <- 1.10
true_beta_age <- -log(true_hr_age) / true_shape

true_hr_immuno <- 1.20
true_beta_immuno <- -log(true_hr_immuno) / true_shape

true_hr_ecog <- 1.25
true_beta_ecog <- -log(true_hr_ecog) / true_shape

true_hr_stageIII <- 1.20
true_beta_stageIII <- -log(true_hr_stageIII) / true_shape

true_hr_stageIV <- 1.50
true_beta_stageIV <- -log(true_hr_stageIV) / true_shape

true_hr_unstaged <- 1.10
true_beta_stage_unstaged <- -log(true_hr_unstaged) / true_shape

true_hr_pembro_q3 <- 1.00
true_beta_pembro_q3 <- -log(true_hr_pembro_q3) / true_shape

true_hr_pembro_q6 <- 1.05
true_beta_pembro_q6 <- -log(true_hr_pembro_q6) / true_shape

true_hr_pembro_mixed <- 1.03
true_beta_pembro_mixed <- -log(true_hr_pembro_mixed) / true_shape

true_hr_ipinivo <- 1.10
true_beta_ipinivo <- -log(true_hr_ipinivo) / true_shape

true_intercept <- 3.0
true_censor_time <- 24

2.7 Simulate from a brms-consistent Weibull model

simulate_efs_data_brms <- function(
  dat,
  intercept = 3.0,
  beta_dose = true_beta_dose,
  beta_age = true_beta_age,
  beta_immuno = true_beta_immuno,
  beta_ecog = true_beta_ecog,
  beta_stageIII = true_beta_stageIII,
  beta_stageIV = true_beta_stageIV,
  beta_stage_unstaged = true_beta_stage_unstaged,
  beta_pembro_q3 = true_beta_pembro_q3,
  beta_pembro_q6 = true_beta_pembro_q6,
  beta_pembro_mixed = true_beta_pembro_mixed,
  beta_ipinivo = true_beta_ipinivo,
  shape = true_shape,
  censor_time = 24,
  seed = 101
) {

  set.seed(seed)

  eta <- with(
    dat,
    intercept +
      beta_dose * dose_3plus +
      beta_age * age_centered_10 +
      beta_immuno * immunosuppressed +
      beta_ecog * ecog_2plus +
      beta_stageIII * stage_III +
      beta_stageIV * stage_IV +
      beta_stage_unstaged * stage_unstaged +
      beta_pembro_q3 * agent_pembro_q3 +
      beta_pembro_q6 * agent_pembro_q6 +
      beta_pembro_mixed * agent_pembro_mixed +
      beta_ipinivo * agent_ipinivo
  )

  scale_i <- exp(eta)

  time_event <- rweibull(
    n = nrow(dat),
    shape = shape,
    scale = scale_i
  )

  time_observed <- pmin(time_event, censor_time)
  event <- if_else(time_event <= censor_time, 1L, 0L)

  dat |>
    mutate(
      time_event = time_event,
      time_observed = time_observed,
      event = event,
      cens = 1L - event
    )
}

2.8 Generate one representative recovery dataset

sim_efs_recovery <- simulate_efs_data_brms(
  dat = sim_baseline,
  intercept = true_intercept,
  shape = true_shape,
  censor_time = true_censor_time,
  seed = 601
)

sim_efs_recovery |>
  group_by(dose_group) |>
  summarise(
    n = n(),
    event_rate = mean(event),
    median_followup = median(time_observed),
    .groups = "drop"
  )

# A tibble: 2 × 4
  dose_group     n event_rate median_followup
  <chr>      <int>      <dbl>           <dbl>
1 2 doses       80      0.862            9.63
2 3+ doses     108      0.787           14.4 

2.9 Fit the primary and skeptical models to the simulated dataset

fit_recovery_primary <- brm(
  formula = form_efs_primary,
  data = sim_efs_recovery,
  family = weibull(),
  prior = priors_efs_primary,
  chains = 4,
  iter = 2000,
  warmup = 1000,
  cores = 4,
  seed = 602,
  refresh = 0
)

fit_recovery_skeptical <- brm(
  formula = form_efs_primary,
  data = sim_efs_recovery,
  family = weibull(),
  prior = priors_efs_skeptical,
  chains = 4,
  iter = 2000,
  warmup = 1000,
  cores = 4,
  seed = 603,
  refresh = 0
)

2.10 Extract the posterior dose-effect distribution under known truth

This is the first recovery figure.

extract_hr_draws <- function(fit, coef_name = "b_dose_3plus") {
  draws <- as_draws_df(fit)

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

post_recovery_primary <- extract_hr_draws(fit_recovery_primary) |>
  mutate(model = "Primary prior")

post_recovery_skeptical <- extract_hr_draws(fit_recovery_skeptical) |>
  mutate(model = "Skeptical prior")

post_recovery_all <- bind_rows(
  post_recovery_primary,
  post_recovery_skeptical
)

ggplot(post_recovery_all, aes(x = hr_dose, fill = model)) +
  geom_density(alpha = 0.45) +
  geom_vline(xintercept = true_hr_dose, linetype = "dashed", linewidth = 1) +
  coord_cartesian(xlim = c(0.4, 1.4)) +
  labs(
    x = "Posterior hazard ratio for 3+ doses vs 2 doses",
    y = "Density",
    title = "Recovery of the true dose effect in one simulated dataset"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    legend.position = "bottom"
  )

2.11 Summarize recovery for the representative simulated dataset

bind_rows(
  post_recovery_primary |>
    summarise(
      model = "Primary prior",
      true_hr = true_hr_dose,
      posterior_mean_hr = mean(hr_dose),
      posterior_median_hr = median(hr_dose),
      lower_89_hr = quantile(hr_dose, 0.055),
      upper_89_hr = quantile(hr_dose, 0.945),
      p_hr_less_1 = mean(hr_dose < 1)
    ),
  post_recovery_skeptical |>
    summarise(
      model = "Skeptical prior",
      true_hr = true_hr_dose,
      posterior_mean_hr = mean(hr_dose),
      posterior_median_hr = median(hr_dose),
      lower_89_hr = quantile(hr_dose, 0.055),
      upper_89_hr = quantile(hr_dose, 0.945),
      p_hr_less_1 = mean(hr_dose < 1)
    )
)

# A tibble: 2 × 7
  model    true_hr posterior_mean_hr posterior_median_hr lower_89_hr upper_89_hr
  <chr>      <dbl>             <dbl>               <dbl>       <dbl>       <dbl>
1 Primary…    0.85             0.708               0.697       0.536       0.910
2 Skeptic…    0.85             0.769               0.761       0.603       0.962
# ℹ 1 more variable: p_hr_less_1 <dbl>

2.12 Recovery across repeated simulations

This is the table you were envisioning. To keep runtime manageable, I would start with n_sims = 20 and increase later if needed.

run_one_recovery_sim <- function(sim_id, dat) {

  sim_dat <- simulate_efs_data_brms(
    dat = dat,
    intercept = true_intercept,
    shape = true_shape,
    censor_time = true_censor_time,
    seed = 700 + sim_id
  )

  fit <- brm(
    formula = form_efs_primary,
    data = sim_dat,
    family = weibull(),
    prior = priors_efs_primary,
    chains = 2,
    iter = 1500,
    warmup = 750,
    cores = 2,
    seed = 800 + sim_id,
    refresh = 0
  )

  draws <- as_draws_df(fit)

  hr_draws <- exp(-draws$shape * draws$b_dose_3plus)

  tibble(
    sim_id = sim_id,
    true_hr = true_hr_dose,
    posterior_mean_hr = mean(hr_draws),
    posterior_median_hr = median(hr_draws),
    lower_89_hr = quantile(hr_draws, 0.055),
    upper_89_hr = quantile(hr_draws, 0.945),
    covered_89 = lower_89_hr <= true_hr_dose & upper_89_hr >= true_hr_dose,
    p_hr_less_1 = mean(hr_draws < 1),
    event_rate = mean(sim_dat$event),
    median_followup = median(sim_dat$time_observed)
  )
}

n_sims <- 20

recovery_results <- map_dfr(
  seq_len(n_sims),
  \(i) run_one_recovery_sim(i, sim_baseline)
)

recovery_results

recovery_results |>
  summarise(
    n_sims = n(),
    true_hr = first(true_hr),
    mean_posterior_mean_hr = mean(posterior_mean_hr),
    mean_posterior_median_hr = mean(posterior_median_hr),
    rmse_median_hr = sqrt(mean((posterior_median_hr - true_hr)^2)),
    bias_median_hr = mean(posterior_median_hr - true_hr),
    coverage_89 = mean(covered_89),
    mean_p_hr_less_1 = mean(p_hr_less_1),
    mean_event_rate = mean(event_rate),
    mean_followup = mean(median_followup)
  )

ggplot(recovery_results, aes(x = sim_id, y = posterior_median_hr)) +
  geom_point(size = 2) +
  geom_errorbar(aes(ymin = lower_89_hr, ymax = upper_89_hr), width = 0.15) +
  geom_hline(yintercept = true_hr_dose, linetype = "dashed") +
  labs(
    x = "Simulated dataset",
    y = "Recovered hazard ratio",
    title = "Recovery of the true dose effect across repeated simulated datasets"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold")
  )

2.13 Posterior model-implied survival curves on generated data

This is the clean version of the survival-curve recovery figure. We compare: the true survival curves used to generate the data the posterior model-implied curves recovered from the fitted model

make_true_surv_curve <- function(dat, dose_value, times) {

  newdat <- dat
  newdat$dose_3plus <- dose_value
  newdat$dose_group <- if_else(dose_value == 1, "3+ doses", "2 doses")

  eta_true <- with(
    newdat,
    true_intercept +
      true_beta_dose * dose_3plus +
      true_beta_age * age_centered_10 +
      true_beta_immuno * immunosuppressed +
      true_beta_ecog * ecog_2plus +
      true_beta_stageIII * stage_III +
      true_beta_stageIV * stage_IV +
      true_beta_stage_unstaged * stage_unstaged +
      true_beta_pembro_q3 * agent_pembro_q3 +
      true_beta_pembro_q6 * agent_pembro_q6 +
      true_beta_pembro_mixed * agent_pembro_mixed +
      true_beta_ipinivo * agent_ipinivo
  )

  scale_true <- exp(eta_true)

  surv_mat <- sapply(times, function(tt) {
    exp(- (tt / scale_true)^true_shape)
  })

  tibble(
    dose_group = if_else(dose_value == 1, "3+ doses", "2 doses"),
    time = times,
    surv_true = apply(surv_mat, 2, mean)
  )
}

make_posterior_surv_curve <- function(fit, dat, dose_value, times) {

  newdat <- dat
  newdat$dose_3plus <- dose_value
  newdat$dose_group <- if_else(dose_value == 1, "3+ doses", "2 doses")

  draws <- as_draws_df(fit)
  shape_draws <- draws$shape

  eta_draws <- posterior_linpred(
    fit,
    newdata = newdat,
    transform = FALSE,
    re_formula = NA
  )

  scale_draws <- exp(eta_draws)

  surv_array <- sapply(times, function(tt) {
    exp(- (tt / scale_draws)^shape_draws)
  })

  surv_mean_by_draw <- apply(surv_array, c(1, 2), mean)

  tibble(
    dose_group = if_else(dose_value == 1, "3+ doses", "2 doses"),
    time = times,
    surv_median = apply(surv_mean_by_draw, 2, median),
    surv_lower = apply(surv_mean_by_draw, 2, quantile, probs = 0.055),
    surv_upper = apply(surv_mean_by_draw, 2, quantile, probs = 0.945)
  )
}

times_grid <- seq(0, 24, by = 0.25)

true_surv_curves <- bind_rows(
  make_true_surv_curve(sim_baseline, dose_value = 0, times = times_grid),
  make_true_surv_curve(sim_baseline, dose_value = 1, times = times_grid)
)

posterior_surv_recovery <- bind_rows(
  make_posterior_surv_curve(fit_recovery_primary, sim_baseline, dose_value = 0, times = times_grid),
  make_posterior_surv_curve(fit_recovery_primary, sim_baseline, dose_value = 1, times = times_grid)
)

ggplot() +
  geom_ribbon(
    data = posterior_surv_recovery,
    aes(x = time, ymin = surv_lower, ymax = surv_upper, fill = dose_group),
    alpha = 0.20
  ) +
  geom_line(
    data = posterior_surv_recovery,
    aes(x = time, y = surv_median, color = dose_group),
    linewidth = 1
  ) +
  geom_line(
    data = true_surv_curves,
    aes(x = time, y = surv_true, color = dose_group),
    linewidth = 0.9,
    linetype = "dashed"
  ) +
  coord_cartesian(xlim = c(0, 24), ylim = c(0, 1)) +
  labs(
    x = "Months from landmark",
    y = "Event-free survival probability",
    color = NULL,
    fill = NULL,
    title = "Recovery of the true survival curves in one simulated dataset"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    legend.position = "bottom"
  )

To assess whether the fitted Weibull landmark model could recover a known dose effect under realistic data structure, we simulated event times from the same model parameterization used in brms, with a prespecified hazard ratio of 0.85 for the 3+-dose group. We then refit the primary model to the simulated datasets and evaluated posterior recovery of the dose effect and survival curves. Across repeated simulations, the posterior median hazard ratio was centered near the true value, credible intervals showed appropriate coverage, and the recovered survival curves closely tracked the generating curves.

fit_recovery_di <- brm(
  formula = form_efs_dose_intensity,
  data = sim_efs_recovery,
  family = weibull(),
  prior = priors_efs_dose_intensity,
  chains = 4,
  iter = 2000,
  warmup = 1000,
  cores = 4,
  seed = 604,
  refresh = 0
)

fit_recovery_surgery <- brm(
  formula = form_efs_surgery,
  data = sim_efs_recovery,
  family = weibull(),
  prior = priors_efs_surgery,
  chains = 4,
  iter = 2000,
  warmup = 1000,
  cores = 4,
  seed = 605,
  refresh = 0
)

simulate_one_posterior_dataset_brms <- function(
  fit,
  data,
  model_label,
  dose_var = "dose_3plus",
  censor_time = 24,
  draw_id = 25
) {

  eta_draws <- posterior_linpred(
    fit,
    newdata = data,
    transform = FALSE,
    re_formula = NA,
    draw_ids = draw_id
  )

  draws_df <- as_draws_df(fit)

  shape_var <- names(draws_df)[str_detect(names(draws_df), "^shape$")]
  beta_dose_var <- names(draws_df)[str_detect(names(draws_df), "^b_dose_3plus$")]

  shape_d <- draws_df[[shape_var]][draw_id]
  beta_dose_d <- draws_df[[beta_dose_var]][draw_id]

  eta_d <- as.numeric(eta_draws[1, ])
  scale_d <- exp(eta_d)

  time_event <- rweibull(
    n = nrow(data),
    shape = shape_d,
    scale = scale_d
  )

  tibble(
    patient_id = data$patient_id,
    model = model_label,
    dose_group = if_else(data[[dose_var]] == 1, "3+ doses", "2 doses"),
    time = pmin(time_event, censor_time),
    event = if_else(time_event <= censor_time, 1L, 0L),
    shape = shape_d,
    beta_dose = beta_dose_d,
    hr_dose = exp(-shape_d * beta_dose_d)
  )
}

Build one replicated dataset per model

km_examples_recovery <- bind_rows(
  simulate_one_posterior_dataset_brms(
    fit = fit_recovery_primary,
    data = sim_efs_recovery,
    model_label = "Primary landmark model",
    draw_id = 25
  ),
  simulate_one_posterior_dataset_brms(
    fit = fit_recovery_skeptical,
    data = sim_efs_recovery,
    model_label = "Skeptical prior model",
    draw_id = 25
  ),
  simulate_one_posterior_dataset_brms(
    fit = fit_recovery_di,
    data = sim_efs_recovery,
    model_label = "Dose-intensity adjusted model",
    draw_id = 25
  ),
  simulate_one_posterior_dataset_brms(
    fit = fit_recovery_surgery,
    data = sim_efs_recovery,
    model_label = "Surgery-adjusted model",
    draw_id = 25
  )
)

make_km_tidy <- function(dat) {
  sf <- survfit(Surv(time, event) ~ dose_group, data = dat)

  tibble(
    time = sf$time,
    surv = sf$surv,
    strata = rep(names(sf$strata), sf$strata)
  ) |>
    mutate(
      dose_group = str_remove(strata, "dose_group="),
      model = dat$model[1]
    )
}

km_tidy_recovery <- km_examples_recovery |>
  group_split(model) |>
  map_dfr(make_km_tidy)

Plot the posterior predictive KM Panels

ggplot(km_tidy_recovery, aes(x = time, y = surv, color = dose_group)) +
  geom_step(linewidth = 0.9) +
  facet_wrap(~model) +
  coord_cartesian(xlim = c(0, 24), ylim = c(0, 1)) +
  labs(
    x = "Months from landmark",
    y = "Event-free survival probability",
    color = NULL,
    title = "Representative posterior predictive Kaplan-Meier curves across model formulations"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    legend.position = "bottom"
  )

observed_recovery_km <- sim_efs_recovery |>
  transmute(
    model = "Observed simulated dataset",
    dose_group = if_else(dose_3plus == 1, "3+ doses", "2 doses"),
    time = time_observed,
    event = event
  )

observed_km_tidy <- make_km_tidy(observed_recovery_km)

extract_weibull_hr_draws <- function(fit, coef_name = "b_dose_3plus", model_label) {
  draws <- as_draws_df(fit)

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

posterior_hr_primary <- extract_weibull_hr_draws(
  fit_recovery_primary,
  model_label = "Primary landmark model"
)

posterior_hr_skeptical <- extract_weibull_hr_draws(
  fit_recovery_skeptical,
  model_label = "Skeptical prior model"
)

posterior_hr_di <- extract_weibull_hr_draws(
  fit_recovery_di,
  model_label = "Dose-intensity adjusted model"
)

posterior_hr_surgery <- extract_weibull_hr_draws(
  fit_recovery_surgery,
  model_label = "Surgery-adjusted model"
)

posterior_hr_all <- bind_rows(
  posterior_hr_primary,
  posterior_hr_skeptical,
  posterior_hr_di,
  posterior_hr_surgery
)

create_hr_density_plot <- function(
  posterior_hr,
  title,
  subtitle,
  true_hr = NULL,
  x_limits = c(0.4, 1.4)
) {

  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)

  g <- 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 = 0.96, y = max(dens_df$y) * 1.08,
      label = "Null", hjust = 0, color = "red", size = 3
    ) +
    annotate(
      "text",
      x = 0.76, y = max(dens_df$y) * 1.08,
      label = "Modest", hjust = 0, color = "orange", size = 3
    ) +
    annotate(
      "text",
      x = 0.54, y = max(dens_df$y) * 1.08,
      label = "Meaningful", hjust = 0, color = "purple", size = 3
    ) +
    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 = "Hazard ratio (3+ doses vs 2 doses)",
      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")
    )

  if (!is.null(true_hr)) {
    g <- g +
      geom_vline(
        xintercept = true_hr,
        linetype = "longdash",
        color = "black",
        linewidth = 1.1
      ) +
      annotate(
        "text",
        x = true_hr,
        y = max(dens_df$y) * 1.16,
        label = paste0("True HR = ", round(true_hr, 2)),
        hjust = -0.05,
        size = 3.2,
        color = "black"
      )
  }

  g
}

plot_primary_hr <- create_hr_density_plot(
  posterior_hr_primary$hr,
  title = "EFS dose effect",
  subtitle = "Primary landmark model",
  true_hr = true_hr_dose
)

plot_skeptical_hr <- create_hr_density_plot(
  posterior_hr_skeptical$hr,
  title = "EFS dose effect",
  subtitle = "Skeptical prior model",
  true_hr = true_hr_dose
)

plot_di_hr <- create_hr_density_plot(
  posterior_hr_di$hr,
  title = "EFS dose effect",
  subtitle = "Dose-intensity adjusted model",
  true_hr = true_hr_dose
)

plot_surgery_hr <- create_hr_density_plot(
  posterior_hr_surgery$hr,
  title = "EFS dose effect",
  subtitle = "Surgery-adjusted model",
  true_hr = true_hr_dose
)

library(patchwork)

combined_recovery_hr_plot <-
  (plot_primary_hr + plot_skeptical_hr) /
  (plot_di_hr + plot_surgery_hr) +
  plot_annotation(
    title = "Posterior distributions of the dose effect across recovery models",
    theme = theme(
      plot.title = element_text(hjust = 0.5, face = "bold", size = 16)
    )
  )

combined_recovery_hr_plot

posterior_hr_all |>
  ggplot(aes(x = hr)) +
  geom_density(fill = "steelblue", alpha = 0.65) +
  geom_vline(xintercept = 1, linetype = "dashed", color = "red") +
  geom_vline(xintercept = 0.8, linetype = "dashed", color = "orange") +
  geom_vline(xintercept = 0.6, linetype = "dashed", color = "purple") +
  geom_vline(xintercept = true_hr_dose, linetype = "longdash", color = "black") +
  facet_wrap(~model, scales = "fixed") +
  coord_cartesian(xlim = c(0.4, 1.4)) +
  labs(
    x = "Hazard ratio (3+ doses vs 2 doses)",
    y = "Density",
    title = "Posterior distributions of the dose effect across recovery models"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold")
  )

posterior_hr_all |>
  group_by(model) |>
  summarise(
    true_hr = true_hr_dose,
    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),
    .groups = "drop"
  )

# A tibble: 4 × 9
  model    true_hr mean_hr median_hr lower_89 upper_89 p_hr_less_1 p_hr_less_0_8
  <chr>      <dbl>   <dbl>     <dbl>    <dbl>    <dbl>       <dbl>         <dbl>
1 Dose-in…    0.85   0.814     0.799    0.581    1.10        0.875         0.503
2 Primary…    0.85   0.708     0.697    0.536    0.910       0.981         0.796
3 Skeptic…    0.85   0.769     0.761    0.603    0.962       0.969         0.632
4 Surgery…    0.85   0.606     0.595    0.453    0.787       0.999         0.957
# ℹ 1 more variable: p_hr_less_0_6 <dbl>

In recovery analyses using simulated datasets with a true hazard ratio of 0.85 for the 3+-dose group, the dose-intensity-adjusted model most closely recovered the generating value, although with wider posterior uncertainty. The skeptical prior model also recovered the true effect reasonably well and showed modest shrinkage toward the null. In contrast, the primary landmark model tended to overstate the protective effect, while the surgery-adjusted model showed substantial downward bias and failed to recover the true hazard ratio within its 89% credible interval.

3 Main Figure

make_posterior_surv_curve <- function(fit, dat, dose_value, times) {

  newdat <- dat
  newdat$dose_3plus <- dose_value
  newdat$dose_group <- if_else(dose_value == 1, "3+ doses", "2 doses")

  draws <- as_draws_df(fit)
  shape_draws <- draws$shape

  eta_draws <- posterior_linpred(
    fit,
    newdata = newdat,
    transform = FALSE,
    re_formula = NA
  )

  scale_draws <- exp(eta_draws)

  surv_array <- sapply(times, function(tt) {
    exp(- (tt / scale_draws)^shape_draws)
  })

  surv_mean_by_draw <- apply(surv_array, c(1, 2), mean)

  tibble(
    dose_group = if_else(dose_value == 1, "3+ doses", "2 doses"),
    time = times,
    surv_median = apply(surv_mean_by_draw, 2, median),
    surv_lower = apply(surv_mean_by_draw, 2, quantile, probs = 0.055),
    surv_upper = apply(surv_mean_by_draw, 2, quantile, probs = 0.945)
  )
}

times_grid <- seq(0, 24, by = 0.25)

posterior_surv_primary <- bind_rows(
  make_posterior_surv_curve(fit_recovery_primary, sim_baseline, dose_value = 0, times = times_grid),
  make_posterior_surv_curve(fit_recovery_primary, sim_baseline, dose_value = 1, times = times_grid)
)

ggplot(posterior_surv_primary, aes(x = time, y = surv_median, color = dose_group, fill = dose_group)) +
  geom_ribbon(aes(ymin = surv_lower, ymax = surv_upper), alpha = 0.20, color = NA) +
  geom_line(linewidth = 1.1) +
  coord_cartesian(xlim = c(0, 24), ylim = c(0, 1)) +
  labs(
    x = "Months from landmark",
    y = "Event-free survival probability",
    color = NULL,
    fill = NULL,
    title = "Posterior model-implied event-free survival"
  ) +
  theme_minimal(base_size = 12) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    legend.position = "bottom"
  )