## ----setup, include=FALSE----------------------------------------------------- knitr::opts_chunk$set(collapse = TRUE, comment = "#>") library(AHTauDesign) filename <- "getting-started-results.rds" paths <- c(file.path("..", "inst", "extdata", filename), file.path("inst", "extdata", filename), system.file("extdata", filename, package = "AHTauDesign")) paths <- paths[nzchar(paths) & file.exists(paths)] if (!length(paths)) { stop("Missing precomputed vignette data. Run the one-time preparation chunk.") } example <- readRDS(paths[1]) if (!identical(example$example_version, "two-branch-202610")) { stop("The vignette RDS is from an older example. Regenerate it once.") } reference <- example$reference alternative <- example$alternative stopifnot(inherits(reference, "ahtau_planning_report"), inherits(alternative, "ahtau_planning_report")) ## ----trial-generator, eval=FALSE---------------------------------------------- # trial_generator <- function(N, A, F, rate_C, rate_T) { # arm <- rep(c(0, 1), each = N / 2) # entry <- runif(N, 0, A) # event_time <- rexp(N, ifelse(arm == 0, rate_C, rate_T)) # censor_time <- A + F - entry # data.frame(arm = arm, X = pmin(event_time, censor_time), # status = as.numeric(event_time <= censor_time)) # } # # null_generator <- function(N, A, F, rate_C) { # trial_generator(N, A, F, rate_C, rate_C) # } # # rate_C <- -log(0.60) / 36 # rate_T <- 0.55 * rate_C ## ----calibration-code, eval=FALSE--------------------------------------------- # calibration_designs <- data.frame( # scenario_id = c("small_early", "medium_early", "large_mature", "large_later"), # N = c(60, 120, 240, 360), # A = 48, F = c(0, 6, 24, 36), rate_C = -log(0.60) / 36 # ) # # rules <- data.frame( # rule_id = c("relative_only", "relative_plus_count"), # rule_family = c("A", "B"), # p_risk_threshold = c(0.10, 0.10), # n_risk_threshold = c(NA, 30), # p_support_threshold = c(0.95, 0.95), # use_absolute_n = c(FALSE, TRUE) # ) # # cal <- simulate_calibration_null( # design_space = calibration_designs, # trial_generator = null_generator, # tau_grid = c(24, 36), n_sim = 10000, seed = 202610, # rules = rules, # reliability_criteria = list(type1_tolerance = 0.06, # coverage_tolerance = 0.94, # support_tolerance = 0.95) # ) # # # Review the candidate comparison before explicitly freezing a rule. # cal$rule_evaluation$performance # frozen <- freeze_support_rule( # diagnostic = cal$diagnostic, rules = rules, # rule_id = "relative_plus_count", # rationale = "Selected from qualifying rules after UPR/APR review.", # calibration_scope = list(designs = calibration_designs, # tau_grid = c(24, 36), # null_model = "equal exponential event hazards"), # provenance = list(source = "synthetic vignette simulation", # inference_conditioning = "support_valid"), # simulation_settings = cal$settings # ) ## ----calibration-summary, echo=FALSE------------------------------------------ counts <- reference$frozen_rule$reliability_counts knitr::kable( counts[, c("n_total", "n_reliable", "n_unreliable", "n_indeterminate")], row.names = FALSE, caption = "Null-calibration reliability classifications" ) perf <- reference$frozen_rule$performance knitr::kable( perf[, c("rule_id", "n_pass", "pass_reliable", "pass_unreliable", "pass_indeterminate", "UPR_conservative", "APR_confirmed")], digits = 3, row.names = FALSE, caption = "Candidate support-rule comparison (illustrative calibration)" ) th <- reference$frozen_rule$thresholds knitr::kable(data.frame( selected_rule = reference$frozen_rule$rule_id, p_risk_min = th$p_risk_min, n_risk_min = th$n_risk_min, p_support_min = th$p_support_min ), digits = 3, row.names = FALSE, caption = "Frozen support thresholds for this synthetic example") ## ----stage1-code, eval=FALSE-------------------------------------------------- # raw_reference <- simulate_design( # trial_generator = trial_generator, # n_sim = 1000, # tau_grid = c(24, 36), # trial_args = list(N = 160, A = 48, F = 36, # rate_C = rate_C, rate_T = rate_T), # scenario_id = "reference", # seed = 202611 # ) # audit_reference <- stage1_audit( # raw_reference, thresholds = frozen$thresholds # ) ## ----stage1-summary, echo=FALSE----------------------------------------------- clin <- reference$stage1_clinical_summary knitr::kable( clin[, c("tau", "n_sim", "min_mean_p_risk", "min_mean_n_risk", "p_support", "Stage1_pass", "power_overall")], digits = 3, row.names = FALSE, caption = "Reference design: support and overall power at 36 months" ) ## ----stage2-code, eval=FALSE-------------------------------------------------- # raw_short <- simulate_design( # trial_generator = trial_generator, # n_sim = 1000, # tau_grid = c(24, 36), # trial_args = list(N = 160, A = 48, F = 6, # rate_C = rate_C, rate_T = rate_T), # scenario_id = "short_followup", # seed = 202612 # ) # audit_short <- stage1_audit(raw_short, thresholds = frozen$thresholds) # selection <- stage2_select(audit_short$scenario_summary, # tau_clin = 36, clinical_window = c(24, 36)) # # # Generate fresh trials; the selected horizon is held fixed. # raw_eval <- simulate_design( # trial_generator = trial_generator, # n_sim = 1000, # tau_grid = c(36, selection$final_decision_summary$tau_star), # trial_args = list(N = 160, A = 48, F = 6, # rate_C = rate_C, rate_T = rate_T), # scenario_id = "independent_evaluation", # seed = 202613 # ) # audit_eval <- stage1_audit(raw_eval, thresholds = frozen$thresholds) # independent_eval <- stage2_evaluate(selection, audit_eval$scenario_summary) # independent_eval$power_difference ## ----stage2-summary, echo=FALSE----------------------------------------------- cand <- alternative$stage2_result$candidate_summary knitr::kable( cand[, c("tau", "min_mean_n_risk", "p_support", "Stage1_pass", "within_clinical_window", "utility_U", "selected")], digits = 3, row.names = FALSE, caption = "Supported candidates, selection utility, and chosen horizon" ) decision <- alternative$stage2_result$final_decision_summary knitr::kable( data.frame(clinical_horizon = decision$tau_clin, chosen_horizon = decision$tau_star, decision = decision$decision_type), row.names = FALSE, caption = "Conditional Stage 2 decision" ) knitr::kable( alternative$stage2_evaluation$eval_summary[ , c("role", "tau", "Stage1_pass", "p_support", "power_overall")], digits = 3, row.names = FALSE, caption = "Independent evaluation of the fixed horizons" ) knitr::kable( data.frame(power_difference = alternative$stage2_evaluation$power_difference), digits = 3, row.names = FALSE, caption = "Power difference: selected alternative minus clinical horizon" ) ## ----part2-code, eval=FALSE--------------------------------------------------- # part2_12 <- find_n_req( # F_followup = 12, tau = 36, thresholds = frozen$thresholds, # trial_generator = trial_generator, # trial_args = list(A = 48, rate_C = rate_C, rate_T = rate_T), # N_lower = 160, N_upper = 600, coarse_step = 20, fine_step = 10, # n_coarse = 200, n_fine = 400, n_validate = 2000, # power_target = 0.80, # search_seed = 202614, validation_seed = 202615 # ) # part2_36 <- find_n_req( # F_followup = 36, tau = 36, thresholds = frozen$thresholds, # trial_generator = trial_generator, # trial_args = list(A = 48, rate_C = rate_C, rate_T = rate_T), # N_lower = 160, N_upper = 600, coarse_step = 20, fine_step = 10, # n_coarse = 200, n_fine = 400, n_validate = 2000, # power_target = 0.80, # search_seed = 202616, validation_seed = 202617 # ) ## ----part2-summary, echo=FALSE------------------------------------------------ req <- reference$part2_requirements knitr::kable( req[, c("F_followup", "tau", "N_req", "status", "power_target", "raw_power", "isotonic_power")], digits = 3, row.names = FALSE, caption = "Validated sample-size requirements for the fixed 36-month horizon" ) ## ----planning-report-code, eval=FALSE----------------------------------------- # report_reference <- planning_report( # trial_spec = list(scenario_id = "reference", N = 160, A = 48, F = 36, # tau_clin = 36, clinical_window = c(24, 36), # power_target = 0.80, survival_model = "exponential"), # frozen_rule = frozen, # stage1_result = audit_reference, # part2_results = list(F12 = part2_12, F36 = part2_36), # simulation_diagnostics = list(base_seed = 202610, # stage1_seed = 202611) # ) # report_alternative <- planning_report( # trial_spec = list(scenario_id = "short_followup", N = 160, A = 48, # F = 6, tau_clin = 36, clinical_window = c(24, 36), # survival_model = "exponential"), # frozen_rule = frozen, # stage1_result = audit_short, # stage2_result = selection, # stage2_evaluation = independent_eval, # simulation_diagnostics = list(selection_seed = 202612, # independent_evaluation_seed = 202613) # ) ## ----planning-report-summary, echo=FALSE-------------------------------------- a <- reference$decision_summary b <- alternative$decision_summary knitr::kable(data.frame( design = c("Reference", "Short follow-up"), Stage1_pass = c(a$Stage1_pass, b$Stage1_pass), power_overall = c(a$power_overall, b$power_overall), Stage2_decision = c(a$stage2_decision, b$stage2_decision), final_horizon = c(a$tau_star, b$tau_star) ), digits = 3, row.names = FALSE, caption = "Major decisions from the two synthetic planning reports") ## ----prepare-vignette-data, eval=FALSE, include=FALSE------------------------- # package_root <- if (file.exists("DESCRIPTION")) "." else ".." # stopifnot(file.exists(file.path(package_root, "DESCRIPTION"))) # base_seed <- 202610 # # trial_generator <- function(N, A, F, rate_C, rate_T) { # stopifnot(N >= 2, N %% 2 == 0) # arm <- rep(c(0, 1), each = N / 2) # entry <- runif(N, 0, A) # event_time <- rexp(N, ifelse(arm == 0, rate_C, rate_T)) # censor_time <- A + F - entry # data.frame(arm = arm, X = pmin(event_time, censor_time), # status = as.numeric(event_time <= censor_time)) # } # null_generator <- function(N, A, F, rate_C) { # trial_generator(N, A, F, rate_C, rate_C) # } # rate_C <- -log(0.60) / 36 # rate_T <- 0.55 * rate_C # # # Calibration includes small, immature settings and larger, mature settings. # calibration_designs <- data.frame( # scenario_id = c("small_early", "medium_early", "large_mature", "large_later"), # N = c(60, 120, 240, 360), A = 48, # F = c(0, 6, 24, 36), rate_C = rate_C # ) # rules <- data.frame( # rule_id = c("relative_only", "relative_plus_count"), # rule_family = c("A", "B"), # p_risk_threshold = c(0.10, 0.10), # n_risk_threshold = c(NA, 30), # p_support_threshold = c(0.95, 0.95), # use_absolute_n = c(FALSE, TRUE) # ) # cal <- simulate_calibration_null( # design_space = calibration_designs, # trial_generator = null_generator, tau_grid = c(24, 36), # n_sim = 10000, seed = base_seed, # rules = rules, # reliability_criteria = list(type1_tolerance = 0.06, # coverage_tolerance = 0.94, # support_tolerance = 0.95) # ) # print(cal$reliability_counts) # print(cal$rule_evaluation$performance[, c("rule_id", "n_pass", # "UPR_conservative", "APR_confirmed")]) # if (!any(cal$classified_diagnostic$reliability_status != "Reliable")) { # stop("Calibration did not show a non-Reliable setting; review the example.") # } # if (identical(as.vector(cal$rule_evaluation$pass_matrix[, 1]), # as.vector(cal$rule_evaluation$pass_matrix[, 2]))) { # stop("Candidate rules gave identical decisions; revise calibration space.") # } # # # Rule must qualify under the actual Monte Carlo classifications. # frozen <- freeze_support_rule( # diagnostic = cal$diagnostic, rules = rules, # rule_id = "relative_plus_count", # rationale = "Qualifying rule selected after comparing conservative UPR and confirmed APR.", # calibration_scope = list(designs = calibration_designs, # tau_grid = c(24, 36), # null_model = "equal exponential hazards"), # provenance = list(source = "synthetic vignette simulation", # inference_conditioning = "support_valid"), # simulation_settings = cal$settings # ) # # # Stage 1: supported but below target power. # raw_reference <- simulate_design( # trial_generator = trial_generator, # n_sim = 1000, # tau_grid = c(24, 36), # trial_args = list(N = 160, A = 48, F = 36, # rate_C = rate_C, rate_T = rate_T), # scenario_id = "reference", # seed = base_seed + 1 # ) # audit_reference <- stage1_audit(raw_reference, thresholds = frozen$thresholds) # clinical_reference <- subset(audit_reference$scenario_summary, tau == 36) # if (!isTRUE(clinical_reference$Stage1_pass) || # !(clinical_reference$power_overall < 0.80)) { # stop("Reference example did not exhibit Stage 1 Pass with power < 0.80.") # } # # # Independent Stage 2 selection scenario: short follow-up causes support failure. # raw_short <- simulate_design( # trial_generator = trial_generator, # n_sim = 1000, # tau_grid = c(24, 36), # trial_args = list(N = 160, A = 48, F = 6, # rate_C = rate_C, rate_T = rate_T), # scenario_id = "short_followup", # seed = base_seed + 2 # ) # audit_short <- stage1_audit(raw_short, thresholds = frozen$thresholds) # selection <- stage2_select(audit_short$scenario_summary, # tau_clin = 36, clinical_window = c(24, 36)) # if (!identical(selection$final_decision_summary$decision_type, # "Utility-driven alternative")) { # stop("Stage 2 example did not select an eligible planning alternative.") # } # # # Fresh simulations only for the two FIXED horizons. # chosen_tau <- selection$final_decision_summary$tau_star # raw_eval <- simulate_design( # trial_generator = trial_generator, # n_sim = 1000, # tau_grid = c(36, chosen_tau), # trial_args = list(N = 160, A = 48, F = 6, # rate_C = rate_C, rate_T = rate_T), # scenario_id = "independent_evaluation", # seed = base_seed + 3 # ) # audit_eval <- stage1_audit(raw_eval, thresholds = frozen$thresholds) # independent_eval <- stage2_evaluate(selection, audit_eval$scenario_summary) # # # Part II keeps the 36-month AH estimand fixed. Simulation counts are illustrative. # part2_12 <- find_n_req( # F_followup = 12, tau = 36, thresholds = frozen$thresholds, # trial_generator = trial_generator, # trial_args = list(A = 48, rate_C = rate_C, rate_T = rate_T), # N_lower = 160, N_upper = 600, coarse_step = 20, fine_step = 10, # n_coarse = 200, n_fine = 400, n_validate = 2000, # power_target = 0.80, search_seed = base_seed + 4, # validation_seed = base_seed + 5 # ) # part2_36 <- find_n_req( # F_followup = 36, tau = 36, thresholds = frozen$thresholds, # trial_generator = trial_generator, # trial_args = list(A = 48, rate_C = rate_C, rate_T = rate_T), # N_lower = 160, N_upper = 600, coarse_step = 20, fine_step = 10, # n_coarse = 200, n_fine = 400, n_validate = 2000, # power_target = 0.80, search_seed = base_seed + 6, # validation_seed = base_seed + 7 # ) # # report_reference <- planning_report( # trial_spec = list(scenario_id = "reference", N = 160, A = 48, F = 36, # allocation = "1:1", tau_clin = 36, # clinical_window = c(24, 36), power_target = 0.80, # survival_model = "exponential"), # frozen_rule = frozen, stage1_result = audit_reference, # part2_results = list(F12 = part2_12, F36 = part2_36), # simulation_diagnostics = list( # base_seed = base_seed, stage1_seed = base_seed + 1, # part2_F12_search_seed = base_seed + 4, # part2_F12_validation_seed = base_seed + 5, # part2_F36_search_seed = base_seed + 6, # part2_F36_validation_seed = base_seed + 7, # stage1_n_sim = 1000, note = "Illustrative Monte Carlo planning example") # ) # report_alternative <- planning_report( # trial_spec = list(scenario_id = "short_followup", N = 160, A = 48, F = 6, # allocation = "1:1", tau_clin = 36, # clinical_window = c(24, 36), survival_model = "exponential"), # frozen_rule = frozen, stage1_result = audit_short, # stage2_result = selection, stage2_evaluation = independent_eval, # simulation_diagnostics = list( # base_seed = base_seed, selection_seed = base_seed + 2, # evaluation_seed = base_seed + 3, # selection_n_sim = 1000, evaluation_n_sim = 1000) # ) # # # Keep only the completed, auditable summary objects. Raw null-calibration # # replicate results are deliberately not bundled with the vignette. # example <- list( # example_version = "two-branch-202610", # reference = report_reference, # alternative = report_alternative # ) # output_dir <- file.path(package_root, "inst", "extdata") # dir.create(output_dir, recursive = TRUE, showWarnings = FALSE) # output_file <- file.path(output_dir, "getting-started-results.rds") # saveRDS(example, output_file, compress = "xz") # stopifnot(file.exists(output_file)) # cat("SUCCESS: updated precomputed vignette data:", # normalizePath(output_file), "\n")