Purpose

This report compares simulated and raw Study 1 overall survival (OS) and progression-free survival (PFS) results for all simulated patient-reported outcome domains:

Pain (PA)
Physical functioning (PF2)
Fatigue (FA)

The simulation is evaluated under the base-case alternative treatment-effect scenario. Simulated survival curves are summarized across Monte Carlo replicates and compared with Kaplan–Meier curves estimated from the raw Study 1 data.

The survival-generating process is repeated within each PRO-domain simulation. Therefore, simulated OS and PFS results are summarized separately by PRO domain to confirm that each domain-specific simulation run reproduces the intended survival distributions.

Setup

library(dplyr)
library(tidyr)
library(haven)
library(survival)
library(survminer)
library(ggplot2)
library(ggpubr)
library(knitr)

PROJECT_ROOT <- "/Users/emmarisner/Desktop/BU/Gilead Fellowship /TTdD2026/7.13Survival_validation"

knitr::opts_knit$set(
  root.dir = PROJECT_ROOT
)

options(
  dplyr.summarise.inform = FALSE,
  scipen = 999
)

PRO_DOMAINS <- c("PA", "PF2", "FA")

DOMAIN_LABELS <- c(
  PA = "Pain",
  PF2 = "Physical Functioning",
  FA = "Fatigue"
)

VALIDATION_SCENARIO <- "base_case_alt"
VALIDATION_EFFECT <- "alternative"

KM_MAX_TIME <- 36
KM_TIME_INCREMENT <- 0.25

KM_GRID <- seq(
  from = 0,
  to = KM_MAX_TIME,
  by = KM_TIME_INCREMENT
)

OUTPUT_DIRECTORY <- file.path(
  PROJECT_ROOT,
  "simulation_outputs",
  "validation"
)

PLOT_DIRECTORY <- file.path(
  OUTPUT_DIRECTORY,
  "plots"
)

TABLE_DIRECTORY <- file.path(
  OUTPUT_DIRECTORY,
  "tables"
)

dir.create(
  PLOT_DIRECTORY,
  recursive = TRUE,
  showWarnings = FALSE
)

dir.create(
  TABLE_DIRECTORY,
  recursive = TRUE,
  showWarnings = FALSE
)

clean_id <- function(x){

  as.character(
    gsub(
      pattern = "^X",
      replacement = "",
      x = as.character(x)
    )
  )
}

clean_visit <- function(x){

  x <- toupper(
    trimws(
      as.character(x)
    )
  )

  gsub(
    pattern = "\\s+",
    replacement = " ",
    x = x
  )
}

check_required_columns <- function(
  data,
  required_columns,
  data_name
){

  missing_columns <- setdiff(
    required_columns,
    names(data)
  )

  if(length(missing_columns) > 0){

    stop(
      paste0(
        data_name,
        " is missing the following required column(s): ",
        paste(
          missing_columns,
          collapse = ", "
        )
      )
    )
  }

  invisible(TRUE)
}

Load simulation outputs

simulation_tte_path <- file.path(
  PROJECT_ROOT,
  "simulation_validation_tte.csv"
)

ice_events_all_path <- file.path(
  PROJECT_ROOT,
  "ice_events_all_all.rds"
)

ice_events_first_path <- file.path(
  PROJECT_ROOT,
  "ice_events_first_all.rds"
)

cat(
  "Simulation file exists:",
  file.exists(simulation_tte_path),
  "\n"
)

## Simulation file exists: TRUE

cat(
  "ICE all exists:",
  file.exists(ice_events_all_path),
  "\n"
)

## ICE all exists: TRUE

cat(
  "ICE first exists:",
  file.exists(ice_events_first_path),
  "\n"
)

## ICE first exists: TRUE

if(!file.exists(simulation_tte_path)){
  stop(
    paste(
      "Simulation validation file not found:",
      simulation_tte_path
    )
  )
}

if(!file.exists(ice_events_all_path)){
  stop(
    paste(
      "ICE all-events file not found:",
      ice_events_all_path
    )
  )
}

if(!file.exists(ice_events_first_path)){
  stop(
    paste(
      "ICE first-event file not found:",
      ice_events_first_path
    )
  )
}

sim_tte_all <- read.csv(
  simulation_tte_path,
  stringsAsFactors = FALSE
)

ice_events_all_all <- readRDS(
  ice_events_all_path
)

ice_events_first_all <- readRDS(
  ice_events_first_path
)

Check required simulation columns

required_simulation_columns <- c(
  "PRO_DOMAIN",
  "scenario_id",
  "effect_setting",
  "endpoint",
  "sim_id",
  "trt",
  "time",
  "event"
)

check_required_columns(
  data = sim_tte_all,
  required_columns = required_simulation_columns,
  data_name = "sim_tte_all"
)

Select validation scenario

sim_tte <- sim_tte_all %>%
  dplyr::filter(
    PRO_DOMAIN %in% PRO_DOMAINS,
    scenario_id == VALIDATION_SCENARIO,
    effect_setting == VALIDATION_EFFECT
  ) %>%
  dplyr::mutate(
    PRO_DOMAIN = as.character(PRO_DOMAIN),
    endpoint = as.character(endpoint),
    sim_id = as.integer(sim_id),
    time = as.numeric(time),
    event = as.integer(event),
    trt = factor(
      trt,
      levels = c(0, 1),
      labels = c(
        "Control",
        "Treatment"
      )
    )
  )

if(nrow(sim_tte) == 0){
  stop(
    "No simulation rows remain after applying the validation filters."
  )
}

sim_tte %>%
  dplyr::count(
    PRO_DOMAIN,
    endpoint,
    sim_id
  )

##    PRO_DOMAIN    endpoint sim_id   n
## 1          FA          OS      1 500
## 2          FA          OS      2 500
## 3          FA          OS      3 500
## 4          FA          OS      4 500
## 5          FA          OS      5 500
## 6          FA          OS      6 500
## 7          FA          OS      7 500
## 8          FA          OS      8 500
## 9          FA          OS      9 500
## 10         FA          OS     10 500
## 11         FA Progression      1 500
## 12         FA Progression      2 500
## 13         FA Progression      3 500
## 14         FA Progression      4 500
## 15         FA Progression      5 500
## 16         FA Progression      6 500
## 17         FA Progression      7 500
## 18         FA Progression      8 500
## 19         FA Progression      9 500
## 20         FA Progression     10 500
## 21         PA          OS      1 500
## 22         PA          OS      2 500
## 23         PA          OS      3 500
## 24         PA          OS      4 500
## 25         PA          OS      5 500
## 26         PA          OS      6 500
## 27         PA          OS      7 500
## 28         PA          OS      8 500
## 29         PA          OS      9 500
## 30         PA          OS     10 500
## 31         PA Progression      1 500
## 32         PA Progression      2 500
## 33         PA Progression      3 500
## 34         PA Progression      4 500
## 35         PA Progression      5 500
## 36         PA Progression      6 500
## 37         PA Progression      7 500
## 38         PA Progression      8 500
## 39         PA Progression      9 500
## 40         PA Progression     10 500
## 41        PF2          OS      1 500
## 42        PF2          OS      2 500
## 43        PF2          OS      3 500
## 44        PF2          OS      4 500
## 45        PF2          OS      5 500
## 46        PF2          OS      6 500
## 47        PF2          OS      7 500
## 48        PF2          OS      8 500
## 49        PF2          OS      9 500
## 50        PF2          OS     10 500
## 51        PF2 Progression      1 500
## 52        PF2 Progression      2 500
## 53        PF2 Progression      3 500
## 54        PF2 Progression      4 500
## 55        PF2 Progression      5 500
## 56        PF2 Progression      6 500
## 57        PF2 Progression      7 500
## 58        PF2 Progression      8 500
## 59        PF2 Progression      9 500
## 60        PF2 Progression     10 500

select PRO

sim_tte <- sim_tte_all %>%
  mutate(
    PRO_DOMAIN = as.character(PRO_DOMAIN),
    scenario_id = as.character(scenario_id),
    effect_setting = as.character(effect_setting),
    endpoint = as.character(endpoint),
    sim_id = as.integer(sim_id),
    time = as.numeric(time),
    event = as.integer(event),
    trt = factor(
      trt,
      levels = c(0, 1),
      labels = c(
        "Control",
        "Treatment"
      )
    )
  ) %>%
  filter(
    PRO_DOMAIN %in% PRO_DOMAINS,
    scenario_id == VALIDATION_SCENARIO,
    effect_setting == VALIDATION_EFFECT,
    !is.na(time),
    !is.na(event),
    !is.na(trt)
  )

if(nrow(sim_tte) == 0){
  stop(
    paste0(
      "No simulated observations were found for scenario ",
      VALIDATION_SCENARIO,
      " and effect setting ",
      VALIDATION_EFFECT,
      "."
    )
  )
}

sim output aval

simulation_availability <- sim_tte %>%
  count(
    PRO_DOMAIN,
    scenario_id,
    effect_setting,
    endpoint,
    name = "number_of_rows"
  ) %>%
  mutate(
    domain_label = unname(
      DOMAIN_LABELS[PRO_DOMAIN]
    )
  ) %>%
  dplyr::select(
    PRO_DOMAIN,
    domain_label,
    scenario_id,
    effect_setting,
    endpoint,
    number_of_rows
  ) %>%
  arrange(
    PRO_DOMAIN,
    endpoint
  )

kable(
  simulation_availability,
  caption = paste(
    "Available simulation data for",
    VALIDATION_SCENARIO,
    "under the",
    VALIDATION_EFFECT,
    "effect setting"
  )
)

Available simulation data for base_case_alt under the alternative effect setting
PRO_DOMAIN	domain_label	scenario_id	effect_setting	endpoint	number_of_rows
FA	Fatigue	base_case_alt	alternative	OS	5000
FA	Fatigue	base_case_alt	alternative	Progression	5000
PA	Pain	base_case_alt	alternative	OS	5000
PA	Pain	base_case_alt	alternative	Progression	5000
PF2	Physical Functioning	base_case_alt	alternative	OS	5000
PF2	Physical Functioning	base_case_alt	alternative	Progression	5000

missing_domains <- setdiff(
  PRO_DOMAINS,
  unique(sim_tte$PRO_DOMAIN)
)

if(length(missing_domains) > 0){
  warning(
    paste0(
      "The following requested PRO domain(s) are missing from ",
      "the selected simulation output: ",
      paste(
        missing_domains,
        collapse = ", "
      )
    )
  )
}

available_domains <- intersect(
  PRO_DOMAINS,
  unique(sim_tte$PRO_DOMAIN)
)

end and censoring checks

simulation_event_check <- sim_tte %>%
  group_by(
    PRO_DOMAIN,
    endpoint,
    sim_id,
    trt
  ) %>%
  summarise(
    n = n(),
    events = sum(
      event == 1,
      na.rm = TRUE
    ),
    censored = sum(
      event == 0,
      na.rm = TRUE
    ),
    event_percentage = 100 * events / n,
    censoring_percentage = 100 * censored / n,
    .groups = "drop"
  )

simulation_event_summary <- simulation_event_check %>%
  group_by(
    PRO_DOMAIN,
    endpoint,
    trt
  ) %>%
  summarise(
    n_simulations = n_distinct(sim_id),
    mean_sample_size = mean(n),
    minimum_sample_size = min(n),
    maximum_sample_size = max(n),
    mean_events = mean(events),
    mean_censored = mean(censored),
    mean_event_percentage = mean(event_percentage),
    mean_censoring_percentage = mean(censoring_percentage),
    minimum_events = min(events),
    maximum_events = max(events),
    .groups = "drop"
  ) %>%
  mutate(
    domain_label = unname(
      DOMAIN_LABELS[PRO_DOMAIN]
    )
  ) %>%
  dplyr::select(
    PRO_DOMAIN,
    domain_label,
    endpoint,
    trt,
    everything()
  ) %>%
  arrange(
    PRO_DOMAIN,
    endpoint,
    trt
  )

kable(
  simulation_event_summary,
  digits = 2,
  caption = "Simulated event and censoring summary across Monte Carlo replicates"
)

Simulated event and censoring summary across Monte Carlo replicates
PRO_DOMAIN	domain_label	endpoint	trt	n_simulations	mean_sample_size	minimum_sample_size	maximum_sample_size	mean_events	mean_censored	mean_event_percentage	mean_censoring_percentage	minimum_events	maximum_events
FA	Fatigue	OS	Control	10	250	250	250	149.9	100.1	59.96	40.04	139	159
FA	Fatigue	OS	Treatment	10	250	250	250	149.6	100.4	59.84	40.16	135	164
FA	Fatigue	Progression	Control	10	250	250	250	191.0	59.0	76.40	23.60	179	201
FA	Fatigue	Progression	Treatment	10	250	250	250	191.2	58.8	76.48	23.52	178	201
PA	Pain	OS	Control	10	250	250	250	151.0	99.0	60.40	39.60	134	160
PA	Pain	OS	Treatment	10	250	250	250	155.5	94.5	62.20	37.80	146	168
PA	Pain	Progression	Control	10	250	250	250	193.5	56.5	77.40	22.60	187	200
PA	Pain	Progression	Treatment	10	250	250	250	194.4	55.6	77.76	22.24	185	208
PF2	Physical Functioning	OS	Control	10	250	250	250	151.2	98.8	60.48	39.52	141	162
PF2	Physical Functioning	OS	Treatment	10	250	250	250	157.9	92.1	63.16	36.84	150	170
PF2	Physical Functioning	Progression	Control	10	250	250	250	192.5	57.5	77.00	23.00	187	198
PF2	Physical Functioning	Progression	Treatment	10	250	250	250	193.2	56.8	77.28	22.72	183	202

write.csv(
  simulation_event_check,
  file.path(
    TABLE_DIRECTORY,
    "simulation_event_check_by_replicate.csv"
  ),
  row.names = FALSE
)

write.csv(
  simulation_event_summary,
  file.path(
    TABLE_DIRECTORY,
    "simulation_event_summary_all_domains.csv"
  ),
  row.names = FALSE
)

##load study 1

PROJECT_ROOT <- "/Users/emmarisner/Desktop/BU/Gilead Fellowship /TTdD2026/7.13Survival_validation"

knitr::opts_knit$set(
  root.dir = PROJECT_ROOT
)

if(!exists("adtte1")){

  adtte1_path <- path.expand(
    "~/Desktop/adtte_study_1.sas7bdat"
  )

  if(!file.exists(adtte1_path)){
    stop(
      paste(
        "Raw Study 1 ADTTE file not found:",
        adtte1_path
      )
    )
  }

  adtte1 <- haven::read_sas(
    adtte1_path
  )
}

if(!exists("adsl1")){

  adsl1_path <- path.expand(
    "~/Desktop/adsl_study_1.csv"
  )

  if(!file.exists(adsl1_path)){
    stop(
      paste(
        "Raw Study 1 ADSL file not found:",
        adsl1_path
      )
    )
  }

  adsl1 <- read.csv(
    adsl1_path,
    stringsAsFactors = FALSE
  )
}

if(!exists("adqs1")){

  adqs1_path <- path.expand(
    "~/Desktop/adqs_study_1.csv"
  )

  if(file.exists(adqs1_path)){

    adqs1 <- read.csv(
      adqs1_path,
      stringsAsFactors = FALSE
    )

  } else {

    warning(
      paste(
        "Raw Study 1 ADQS file not found:",
        adqs1_path
      )
    )
  }
}

veify event coding

raw_event_coding_check <- adtte1 %>%
  filter(
    PARAMCD %in% c(
      "OS",
      "PFS1"
    )
  ) %>%
  count(
    PARAMCD,
    CNSR,
    EVNTDESC,
    name = "number_of_records"
  ) %>%
  arrange(
    PARAMCD,
    CNSR,
    EVNTDESC
  )

kable(
  raw_event_coding_check,
  caption = "Raw Study 1 event and censoring coding check"
)

Raw Study 1 event and censoring coding check
PARAMCD	CNSR	EVNTDESC	number_of_records
OS	0	NULL	153
OS	1	DEATH PRIOR TO DATE OF CUTOFF	390
PFS1	0	DEATH AFTER 2 OR MORE CONSECUTIVE MISSING VISITS	28
PFS1	0	DEATH AFTER STARTING NEW ANTI-CANCER THERAPY	31
PFS1	0	NO BASELINE OR NO POSTBASELINE	15
PFS1	0	NO PD AND NO DEATH	25
PFS1	0	PD AFTER STARTING NEW ANTI-CANCER THERAPY	2
PFS1	1	DEATH	36
PFS1	1	RADIOGRAPHIC DISEASE PROGRESSION	406

OS and PFS endpints

raw_trt <- adsl1 %>%
  dplyr::transmute(
    id = clean_id(USUBJID),
    trt = dplyr::case_when(
      TRT01PN == 2 ~ "Control",
      TRT01PN == 1 ~ "Treatment",
      TRUE ~ NA_character_
    )
  ) %>%
  dplyr::mutate(
    trt = factor(
      trt,
      levels = c("Control", "Treatment")
    )
  ) %>%
  dplyr::filter(
    !is.na(id),
    !is.na(trt)
  ) %>%
  dplyr::distinct(
    id,
    .keep_all = TRUE
  )

raw_os <- adtte1 %>%
  filter(
    PARAMCD == "OS"
  ) %>%
  transmute(
    id = clean_id(USUBJID),
    endpoint = "OS",
    time = as.numeric(AVAL),
    event = as.integer(CNSR == 1)
  ) %>%
  left_join(
    raw_trt,
    by = "id"
  ) %>%
  filter(
    !is.na(time),
    !is.na(event),
    !is.na(trt)
  )

raw_pfs <- adtte1 %>%
  filter(
    PARAMCD == "PFS1"
  ) %>%
  transmute(
    id = clean_id(USUBJID),
    endpoint = "PFS",
    time = as.numeric(AVAL),
  event = as.integer(
  CNSR == 1 &
    EVNTDESC == "RADIOGRAPHIC DISEASE PROGRESSION")) %>%
  left_join(
    raw_trt,
    by = "id"
  ) %>%
  filter(
    !is.na(time),
    !is.na(event),
    !is.na(trt)
  )

raw_event_summary <- bind_rows(
  raw_os,
  raw_pfs
) %>%
  group_by(
    endpoint,
    trt
  ) %>%
  summarise(
    n = n(),
    events = sum(
      event == 1,
      na.rm = TRUE
    ),
    censored = sum(
      event == 0,
      na.rm = TRUE
    ),
    event_percentage = 100 * events / n,
    censoring_percentage = 100 * censored / n,
    .groups = "drop"
  )

kable(
  raw_event_summary,
  digits = 2,
  caption = "Raw Study 1 event and censoring summary"
)

Raw Study 1 event and censoring summary
endpoint	trt	n	events	censored	event_percentage	censoring_percentage
OS	Control	271	199	72	73.43	26.57
OS	Treatment	272	191	81	70.22	29.78
PFS	Control	271	190	81	70.11	29.89
PFS	Treatment	272	216	56	79.41	20.59

write.csv(
  raw_event_summary,
  file.path(
    TABLE_DIRECTORY,
    "raw_study_1_event_summary.csv"
  ),
  row.names = FALSE
)

KM func

get_km_one <- function(
  data,
  times = KM_GRID
){

  if(nrow(data) == 0){

    return(
      data.frame(
        time = numeric(0),
        surv = numeric(0)
      )
    )
  }

  fit <- survfit(
    Surv(time, event) ~ 1,
    data = data
  )

  fit_summary <- summary(
    fit,
    times = times,
    extend = TRUE
  )

  data.frame(
    time = fit_summary$time,
    surv = fit_summary$surv
  )
}

raw KM plot

make_raw_plot <- function(
  data,
  title,
  y_label
){

  ggsurvplot(
    fit = survfit(
      Surv(time, event) ~ trt,
      data = data
    ),
    data = data,
    conf.int = TRUE,
    pval = TRUE,
    risk.table = FALSE,
    title = title,
    xlab = "Months from baseline",
    ylab = y_label,
    legend.title = "Treatment",
    legend.labs = c(
      "Control",
      "Treatment"
    ),
    ggtheme = theme_bw()
  )$plot
}

sim plot

sim_km_by_rep <- sim_tte %>%
  dplyr::group_by(
    PRO_DOMAIN,
    endpoint,
    sim_id,
    trt
  ) %>%
  dplyr::group_modify(
    ~ get_km_one(.x)
  ) %>%
  dplyr::ungroup()

sim_km_mc <- sim_km_by_rep %>%
  group_by(
    PRO_DOMAIN,
    endpoint,
    trt,
    time
  ) %>%
  summarise(
    n_sims = sum(
      !is.na(surv)
    ),
    mean_surv = mean(
      surv,
      na.rm = TRUE
    ),
    mc_sd = sd(
      surv,
      na.rm = TRUE
    ),
    mc_se = mc_sd / sqrt(n_sims),
    lower_mc = pmax(
      0,
      mean_surv - 1.96 * mc_se
    ),
    upper_mc = pmin(
      1,
      mean_surv + 1.96 * mc_se
    ),
    lower_prediction = pmax(
      0,
      mean_surv - 1.96 * mc_sd
    ),
    upper_prediction = pmin(
      1,
      mean_surv + 1.96 * mc_sd
    ),
    .groups = "drop"
  )

km_replicate_check <- sim_km_mc %>%
  group_by(
    PRO_DOMAIN,
    endpoint,
    trt
  ) %>%
  summarise(
    minimum_number_of_curves = min(
      n_sims,
      na.rm = TRUE
    ),
    maximum_number_of_curves = max(
      n_sims,
      na.rm = TRUE
    ),
    .groups = "drop"
  )

kable(
  km_replicate_check,
  caption = "Number of simulated Kaplan–Meier curves contributing to the Monte Carlo summaries"
)

Number of simulated Kaplan–Meier curves contributing to the Monte Carlo summaries
PRO_DOMAIN	endpoint	trt	minimum_number_of_curves	maximum_number_of_curves
FA	OS	Control	10	10
FA	OS	Treatment	10	10
FA	Progression	Control	10	10
FA	Progression	Treatment	10	10
PA	OS	Control	10	10
PA	OS	Treatment	10	10
PA	Progression	Control	10	10
PA	Progression	Treatment	10	10
PF2	OS	Control	10	10
PF2	OS	Treatment	10	10
PF2	Progression	Control	10	10
PF2	Progression	Treatment	10	10

write.csv(
  sim_km_by_rep,
  file.path(
    TABLE_DIRECTORY,
    "simulation_km_estimates_by_replicate.csv"
  ),
  row.names = FALSE
)

write.csv(
  sim_km_mc,
  file.path(
    TABLE_DIRECTORY,
    "simulation_km_monte_carlo_summary.csv"
  ),
  row.names = FALSE
)

make_sim_mc_plot <- function(
  domain_name,
  endpoint_name,
  title,
  y_label,
  band = c(
    "prediction",
    "mc_se"
  )
){

  band <- match.arg(band)

  plot_data <- sim_km_mc %>%
    dplyr::filter(
      PRO_DOMAIN == domain_name,
      endpoint == endpoint_name
    )

  if(nrow(plot_data) == 0){

    stop(
      paste0(
        "No simulated Kaplan–Meier data were found for domain ",
        domain_name,
        " and endpoint ",
        endpoint_name,
        "."
      )
    )
  }

  if(band == "mc_se"){

    plot_data <- plot_data %>%
      dplyr::mutate(
        ymin = lower_mc,
        ymax = upper_mc
      )

  } else {

    plot_data <- plot_data %>%
      dplyr::mutate(
        ymin = lower_prediction,
        ymax = upper_prediction
      )
  }

  ggplot(
    plot_data,
    aes(
      x = time,
      y = mean_surv,
      color = trt,
      fill = trt
    )
  ) +
    geom_ribbon(
      aes(
        ymin = ymin,
        ymax = ymax
      ),
      alpha = 0.18,
      color = NA
    ) +
    geom_step(
      linewidth = 1
    ) +
    coord_cartesian(
      xlim = c(
        0,
        KM_MAX_TIME
      ),
      ylim = c(
        0,
        1
      )
    ) +
    labs(
      title = title,
      subtitle = paste(
        "Simulation domain:",
        DOMAIN_LABELS[[domain_name]]
      ),
      x = "Months from baseline",
      y = y_label,
      color = "Treatment",
      fill = "Treatment"
    ) +
    theme_bw() +
    theme(
      plot.title = element_text(
        face = "bold"
      )
    )
}

raw_os_plot <- make_raw_plot(
  data = raw_os,
  title = "Raw Study 1 Overall Survival",
  y_label = "Overall survival probability"
)

raw_pfs_plot <- make_raw_plot(
  data = raw_pfs,
  title = "Raw Study 1 Progression-Free Survival",
  y_label = "Progression-free survival probability"
)

compare all pro domains

os_comparison_plots <- list()
pfs_comparison_plots <- list()

for(domain_name in available_domains){

  simulated_os_plot <- make_sim_mc_plot(
    domain_name = domain_name,
    endpoint_name = "OS",
    title = paste(
      "Monte Carlo Mean Simulated OS:",
      DOMAIN_LABELS[[domain_name]]
    ),
    y_label = "Overall survival probability",
    band = "prediction"
  )

  simulated_pfs_plot <- make_sim_mc_plot(
    domain_name = domain_name,
    endpoint_name = "Progression",
    title = paste(
      "Monte Carlo Mean Simulated PFS:",
      DOMAIN_LABELS[[domain_name]]
    ),
    y_label = "Progression-free survival probability",
    band = "prediction"
  )

  os_comparison_plots[[domain_name]] <- ggarrange(
  raw_os_plot,
  simulated_os_plot,
  ncol = 2,
  common.legend = FALSE
)

pfs_comparison_plots[[domain_name]] <- ggarrange(
  raw_pfs_plot,
  simulated_pfs_plot,
  ncol = 2,
  common.legend = FALSE
)
}

pain validation

if("PA" %in% names(os_comparison_plots)){
  os_comparison_plots[["PA"]]
}

if("PA" %in% names(pfs_comparison_plots)){
  pfs_comparison_plots[["PA"]]
}

survfit(Surv(time, event) ~ trt, data = raw_pfs)

## Call: survfit(formula = Surv(time, event) ~ trt, data = raw_pfs)
## 
##                 n events median 0.95LCL 0.95UCL
## trt=Control   271    190   3.71    2.83    4.27
## trt=Treatment 272    216   4.63    4.14    5.59

PF validaion

if("PF2" %in% names(os_comparison_plots)){
  os_comparison_plots[["PF2"]]
}

if("PF2" %in% names(pfs_comparison_plots)){
  pfs_comparison_plots[["PF2"]]
}

fatigue validation

if("FA" %in% names(os_comparison_plots)){
  os_comparison_plots[["FA"]]
}

if("FA" %in% names(pfs_comparison_plots)){
  pfs_comparison_plots[["FA"]]
}

save plots

for(domain_name in available_domains){

  ggsave(
    filename = file.path(
      PLOT_DIRECTORY,
      paste0(
        domain_name,
        "_raw_vs_simulated_OS.png"
      )
    ),
    plot = os_comparison_plots[[domain_name]],
    width = 12,
    height = 6,
    dpi = 300
  )

  ggsave(
    filename = file.path(
      PLOT_DIRECTORY,
      paste0(
        domain_name,
        "_raw_vs_simulated_PFS.png"
      )
    ),
    plot = pfs_comparison_plots[[domain_name]],
    width = 12,
    height = 6,
    dpi = 300
  )
}

raw summary

extract_survival_summary <- function(data){

  fit <- survival::survfit(
    survival::Surv(time, event) ~ trt,
    data = data
  )

  fit_table <- as.data.frame(
    summary(fit)$table
  )

  fit_table$trt <- rownames(fit_table)
  rownames(fit_table) <- NULL

  fit_table %>%
    dplyr::transmute(
      trt = gsub(
        pattern = "^trt=",
        replacement = "",
        x = trt
      ),
      n = records,
      events = events,
      median = median,
      lower_95 = `0.95LCL`,
      upper_95 = `0.95UCL`
    )
}

raw_os_summary <- extract_survival_summary(
  raw_os
) %>%
  mutate(
    data_source = "Raw Study 1",
    PRO_DOMAIN = NA_character_,
    endpoint = "OS"
  )

raw_pfs_summary <- extract_survival_summary(
  raw_pfs
) %>%
  mutate(
    data_source = "Raw Study 1",
    PRO_DOMAIN = NA_character_,
    endpoint = "PFS"
  )

raw_survival_summary <- bind_rows(
  raw_os_summary,
  raw_pfs_summary
) %>%
  dplyr::select(
    data_source,
    PRO_DOMAIN,
    endpoint,
    trt,
    n,
    events,
    median,
    lower_95,
    upper_95
  )

kable(
  raw_survival_summary,
  digits = 2,
  caption = "Raw Study 1 Kaplan–Meier survival summaries"
)

Raw Study 1 Kaplan–Meier survival summaries
data_source	PRO_DOMAIN	endpoint	trt	n	events	median	lower_95	upper_95
Raw Study 1	NA	OS	Control	271	199	11.20	10.35	12.78
Raw Study 1	NA	OS	Treatment	272	191	14.42	13.11	16.00
Raw Study 1	NA	PFS	Control	271	190	3.71	2.83	4.27
Raw Study 1	NA	PFS	Treatment	272	216	4.63	4.14	5.59

sim summary

 get_simulated_survival_summary <- function(data){

  fit <- survfit(
    Surv(time, event) ~ 1,
    data = data
  )

  fit_table <- summary(fit)$table

  median_value <- unname(
    fit_table["median"]
  )

  lower_value <- unname(
    fit_table["0.95LCL"]
  )

  upper_value <- unname(
    fit_table["0.95UCL"]
  )

  data.frame(
    n = nrow(data),
    events = sum(
      data$event == 1,
      na.rm = TRUE
    ),
    censored = sum(
      data$event == 0,
      na.rm = TRUE
    ),
    median = as.numeric(
      median_value
    ),
    lower_95 = as.numeric(
      lower_value
    ),
    upper_95 = as.numeric(
      upper_value
    )
  )
}

sim_survival_by_replicate <- sim_tte %>%
  group_by(
    PRO_DOMAIN,
    endpoint,
    sim_id,
    trt
  ) %>%
  group_modify(
    ~ get_simulated_survival_summary(.x)
  ) %>%
  ungroup()

sim_survival_summary <- sim_survival_by_replicate %>%
  group_by(
    PRO_DOMAIN,
    endpoint,
    trt
  ) %>%
  summarise(
    n_simulations = n_distinct(sim_id),
    mean_sample_size = mean(
      n,
      na.rm = TRUE
    ),
    mean_events = mean(
      events,
      na.rm = TRUE
    ),
    mean_censored = mean(
      censored,
      na.rm = TRUE
    ),
    mean_median_survival = mean(
      median,
      na.rm = TRUE
    ),
    median_of_median_survival = median(
      median,
      na.rm = TRUE
    ),
    sd_median_survival = sd(
      median,
      na.rm = TRUE
    ),
    empirical_lower_95 = quantile(
      median,
      probs = 0.025,
      na.rm = TRUE,
      names = FALSE
    ),
    empirical_upper_95 = quantile(
      median,
      probs = 0.975,
      na.rm = TRUE,
      names = FALSE
    ),
    .groups = "drop"
  ) %>%
  mutate(
    domain_label = unname(
      DOMAIN_LABELS[PRO_DOMAIN]
    )
  ) %>%
  dplyr::select(
    PRO_DOMAIN,
    domain_label,
    endpoint,
    trt,
    everything()
  ) %>%
  arrange(
    PRO_DOMAIN,
    endpoint,
    trt
  )

kable(
  sim_survival_summary,
  digits = 2,
  caption = "Simulated survival summaries across Monte Carlo replicates"
)

Simulated survival summaries across Monte Carlo replicates
PRO_DOMAIN	domain_label	endpoint	trt	n_simulations	mean_sample_size	mean_events	mean_censored	mean_median_survival	median_of_median_survival	sd_median_survival	empirical_lower_95	empirical_upper_95
FA	Fatigue	OS	Control	10	250	149.9	100.1	19.74	18.84	2.32	17.22	22.95
FA	Fatigue	OS	Treatment	10	250	149.6	100.4	22.11	21.39	3.17	18.08	27.02
FA	Fatigue	Progression	Control	10	250	191.0	59.0	4.51	4.45	1.48	2.49	7.00
FA	Fatigue	Progression	Treatment	10	250	191.2	58.8	5.91	5.58	1.61	3.93	9.14
PA	Pain	OS	Control	10	250	151.0	99.0	21.37	21.74	2.68	17.50	25.99
PA	Pain	OS	Treatment	10	250	155.5	94.5	22.08	22.70	1.85	18.70	24.50
PA	Pain	Progression	Control	10	250	193.5	56.5	4.42	4.56	1.04	3.10	5.97
PA	Pain	Progression	Treatment	10	250	194.4	55.6	5.96	5.87	1.02	4.58	7.33
PF2	Physical Functioning	OS	Control	10	250	151.2	98.8	21.84	22.18	2.36	17.73	24.53
PF2	Physical Functioning	OS	Treatment	10	250	157.9	92.1	20.67	20.55	1.88	17.45	23.59
PF2	Physical Functioning	Progression	Control	10	250	192.5	57.5	5.50	6.00	1.23	3.65	6.73
PF2	Physical Functioning	Progression	Treatment	10	250	193.2	56.8	6.60	6.60	1.33	4.69	8.50

write.csv(
  raw_survival_summary,
  file.path(
    TABLE_DIRECTORY,
    "raw_study_1_survival_summary.csv"
  ),
  row.names = FALSE
)

write.csv(
  sim_survival_by_replicate,
  file.path(
    TABLE_DIRECTORY,
    "simulated_survival_summary_by_replicate.csv"
  ),
  row.names = FALSE
)

write.csv(
  sim_survival_summary,
  file.path(
    TABLE_DIRECTORY,
    "simulated_survival_summary_all_domains.csv"
  ),
  row.names = FALSE
)

raw vs sim mediam surv

raw_median_for_comparison <- raw_survival_summary %>%
  transmute(
    endpoint,
    trt,
    raw_n = n,
    raw_events = events,
    raw_median_survival = median,
    raw_lower_95 = lower_95,
    raw_upper_95 = upper_95
  )

sim_median_for_comparison <- sim_survival_summary %>%
  mutate(
    endpoint = case_when(
      endpoint == "Progression" ~ "PFS",
      TRUE ~ endpoint
    )
  )

raw_vs_simulated_survival_comparison <- sim_median_for_comparison %>%
  left_join(
    raw_median_for_comparison,
    by = c(
      "endpoint",
      "trt"
    )
  ) %>%
  mutate(
    absolute_difference_in_median = 
      mean_median_survival -
      raw_median_survival
  ) %>%
  dplyr::select(
    PRO_DOMAIN,
    domain_label,
    endpoint,
    trt,
    raw_n,
    raw_events,
    raw_median_survival,
    raw_lower_95,
    raw_upper_95,
    n_simulations,
    mean_sample_size,
    mean_events,
    mean_censored,
    mean_median_survival,
    median_of_median_survival,
    sd_median_survival,
    empirical_lower_95,
    empirical_upper_95,
    absolute_difference_in_median
  ) %>%
  arrange(
    endpoint,
    PRO_DOMAIN,
    trt
  )

kable(
  raw_vs_simulated_survival_comparison,
  digits = 2,
  caption = "Raw versus simulated median survival comparison"
)

Raw versus simulated median survival comparison
PRO_DOMAIN	domain_label	endpoint	trt	raw_n	raw_events	raw_median_survival	raw_lower_95	raw_upper_95	n_simulations	mean_sample_size	mean_events	mean_censored	mean_median_survival	median_of_median_survival	sd_median_survival	empirical_lower_95	empirical_upper_95	absolute_difference_in_median
FA	Fatigue	OS	Control	271	199	11.20	10.35	12.78	10	250	149.9	100.1	19.74	18.84	2.32	17.22	22.95	8.54
FA	Fatigue	OS	Treatment	272	191	14.42	13.11	16.00	10	250	149.6	100.4	22.11	21.39	3.17	18.08	27.02	7.69
PA	Pain	OS	Control	271	199	11.20	10.35	12.78	10	250	151.0	99.0	21.37	21.74	2.68	17.50	25.99	10.17
PA	Pain	OS	Treatment	272	191	14.42	13.11	16.00	10	250	155.5	94.5	22.08	22.70	1.85	18.70	24.50	7.65
PF2	Physical Functioning	OS	Control	271	199	11.20	10.35	12.78	10	250	151.2	98.8	21.84	22.18	2.36	17.73	24.53	10.64
PF2	Physical Functioning	OS	Treatment	272	191	14.42	13.11	16.00	10	250	157.9	92.1	20.67	20.55	1.88	17.45	23.59	6.24
FA	Fatigue	PFS	Control	271	190	3.71	2.83	4.27	10	250	191.0	59.0	4.51	4.45	1.48	2.49	7.00	0.79
FA	Fatigue	PFS	Treatment	272	216	4.63	4.14	5.59	10	250	191.2	58.8	5.91	5.58	1.61	3.93	9.14	1.27
PA	Pain	PFS	Control	271	190	3.71	2.83	4.27	10	250	193.5	56.5	4.42	4.56	1.04	3.10	5.97	0.70
PA	Pain	PFS	Treatment	272	216	4.63	4.14	5.59	10	250	194.4	55.6	5.96	5.87	1.02	4.58	7.33	1.32
PF2	Physical Functioning	PFS	Control	271	190	3.71	2.83	4.27	10	250	192.5	57.5	5.50	6.00	1.23	3.65	6.73	1.79
PF2	Physical Functioning	PFS	Treatment	272	216	4.63	4.14	5.59	10	250	193.2	56.8	6.60	6.60	1.33	4.69	8.50	1.97

write.csv(
  raw_vs_simulated_survival_comparison,
  file.path(
    TABLE_DIRECTORY,
    "raw_vs_simulated_survival_comparison_all_domains.csv"
  ),
  row.names = FALSE
)

sim_tte %>%
  dplyr::count(
    PRO_DOMAIN,
    endpoint
  )

##   PRO_DOMAIN    endpoint    n
## 1         FA          OS 5000
## 2         FA Progression 5000
## 3         PA          OS 5000
## 4         PA Progression 5000
## 5        PF2          OS 5000
## 6        PF2 Progression 5000

sim_tte %>%
  dplyr::filter(
    PRO_DOMAIN == "PA",
    endpoint == "Progression"
  ) %>%
  dplyr::summarise(
    n = dplyr::n(),
    mean_time = mean(time, na.rm = TRUE),
    median_time = median(time, na.rm = TRUE),
    event_rate = mean(event, na.rm = TRUE)
  )

##      n mean_time median_time event_rate
## 1 5000  10.02974    3.130628     0.7758

nrow(raw_pfs)

## [1] 543

table(raw_pfs$event)

## 
##   0   1 
## 137 406

summary(raw_pfs$time)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##  0.03285  1.41273  2.85832  4.38722  5.76591 29.07598

TTD Simulation Validation: Survival Curve Checks for All PRO Domains

Emma Risner, M.S., Gilead Sciences

July 14, 2026

Purpose

Setup

Setup

Load simulation outputs

Check required simulation columns

Select validation scenario

select PRO

sim output aval

end and censoring checks

veify event coding

OS and PFS endpints

KM func

raw KM plot

sim plot

compare all pro domains

pain validation

PF validaion

fatigue validation

save plots

raw summary

sim summary

raw vs sim mediam surv