Prototype code that could go into simulations table generation

Result can look like this: 

<img width="1371" height="577" alt="Image" src="https://github.com/user-attachments/assets/9773645d-34c6-4f0d-8a2a-778b1b1c0cb3" />

Code:

```r

library(crmPack)
library(knitr)
library(tidyverse)
library(kableExtra)
library(huxtable)
library(labelled)
library(parallel)
library(future.apply)

outdir <- "./03-reporting/CRM"
datadir <- "data/other"

saved.image <- file.path(datadir, "myfile.RData")

if (file.exists(saved.image)) {
  load(saved.image)
  
} else {
  #SPECIFY TRIAL DETAILS----------------------------------------------------------
  
  # General Parameters
  doses <- c(0.5, 1.5, 3.25, 6.5) #Doses
  
  coarseGrid <- doses
  
  #coarseGrid <- seq(min(doses),max(doses),by=1)
  
  #coarseGrid <- seq(min(doses),max(doses),by=0.5)
  
  DE.DoseGrid  <- c(0.05, doses)  # Dose grid over which Dose Escalation will proceed - added an 'escape dose' of 0.05mg to avoid stopping after a DLT at the lowest dose
  
  #DE.DoseGrid <- doses #no escape dose
  
  emptyData <- Data(doseGrid = DE.DoseGrid, placebo = F)
  
  refDose <- 1.5
  
  #mySize <- CohortSizeConst(6)
  
  mySize <- CohortSizeRange(intervals=c(0,6.5),cohort_size = c(3,6))
  
  nTrialSims <- 100
  
  burning <- 10
  steps <- 1
  samples <- 1000
  
  myOptions <- McmcOptions(burnin = burning,
                           step = steps,
                           samples = samples)
  
  # Dose escalation criteria
  DE.Incr <- IncrementsDoseLevels(levels = 1, basis_level = "last")
  
  T.low <- 0.05   # Lower target dose level
  T.upp <- 0.2   # Upper target dose level
  OD.low <- 0.2  # over-dose lower limit
  maxODprob <- 0.25
  
  myNextBest <- NextBestNCRM(
    target = c(T.low, T.upp),
    overdose = c(OD.low, 1),
    max_overdose_prob = maxODprob
  )
  
  # Stopping criteria
  myStop1 <- StoppingMinPatients(24)   # Not really needed and is arbitrary here, remove? This is a futility threshold around finding the MTD.
  myStop2 <- StoppingTargetProb(target = c(T.low, T.upp), prob = 0.99)   # prob is high so likely to drive the stopping criteria
  myStop3 <- StoppingPatientsNearDose(nPatients = 9, percentage = 0)
 # DE.Stop <- myStop1 | (myStop2 & myStop3)
  myStop4 <- StoppingMissingDose()
  
  DE.Stop <- myStop1 | myStop4
  
  # Prior
  threshmin <- 0.1 # Upper toxicity bound at lowest dose
  probmin <- 0.37 # Probability of toxicity above threshmin for lowest dose
  threshmax <- 0.1 #lowest toxicity bound at highest dose
  probmax <- 0.47 # Probability of toxicity below threshmax for highst dose
  
  minInfModel <- MinimalInformative(
    dosegrid = coarseGrid,
    refDose = refDose,
    #logNormal = TRUE,
    threshmin = threshmin,
    probmin = probmin,
    threshmax = threshmax,
    probmax = probmax,
    seed = 100,
    control = list(max.time = 20)
  )
  
  prior_samples <- mcmc(data = emptyData,
                        model = minInfModel$model,
                        options = myOptions)
  
  #PRIOR TOXICITY CURVE---------------------------------------------------------
  prior.plot <- plot(prior_samples, minInfModel$model, emptyData) +
    ggtitle("Figure 1. Prior for safety analysis") +
    labs(x = "Dose (mg)", y = paste0("DLAE Rate")) +
    geom_hline(yintercept = c(T.low, T.upp)*100, lty = 1, color = "grey")+
    theme_bw()
  
  myModelParams <- minInfModel$parameters
  myModelmeancov <- minInfModel$model
  alpha0 <- unname(myModelParams[1])
  alpha1 <- unname(myModelParams[2])
  cov1 <- myModelmeancov@params@cov[1, 1]
  cov2 <- myModelmeancov@params@cov[1, 2]
  cov3 <- myModelmeancov@params@cov[2, 1]
  cov4 <- myModelmeancov@params@cov[2, 2]
  
  # Design
  design <- Design(
    model = minInfModel$model,
    nextBest = myNextBest,
    stopping = DE.Stop,
    increments = DE.Incr,
    cohort_size = mySize,
    data = emptyData,
    startingDose = min(coarseGrid)
  )
  
  # Simulation parameters
  alpha1fun <- function(DT, DO, pT, pO) {
    return((log(pT / (1 - pT)) - log(pO / (1 - pO))) / (log(DT / refDose) - log(DO /
                                                                                  refDose)))
  }
  
  alpha0fun <- function(DT, pT, alpha1funRes) {
    return((log(pT / (1 - pT)) - alpha1funRes * log(DT / refDose)))
  }
  
  ## Truth 1: Optimistic - 1.5 mg dose achieves target toxicity and OD level (35%) not reached until the max dose (>6.5 mg)
  DT1  <- 50
  DO1  <- 100
  pT1  <- (T.upp + T.low) / 2
  pO1 <- OD.low
  
  alpha1.1 <- alpha1fun(DT1, DO1, pT1, pO1)
  alpha0.1 <- alpha0fun(DT1, pT1, alpha1.1)
  
  optTruth <- function(dose) {
    z <- exp(alpha0.1 + alpha1.1 * log(dose / refDose))
    return(z / (1 + z))
  }
  
  ## Truth 2: Realistic - 0.8 mg dose achieves target toxicity and OD level (35%) reached at 3 mg
  DT2 <- 3.25
  DO2 <- 6.5
  pT2 <- (T.upp + T.low) / 2
  pO2 <- OD.low
  
  alpha1.2 <- alpha1fun(DT2, DO2, pT2, pO2)
  alpha0.2 <- alpha0fun(DT2, pT2, alpha1.2)
  
  relTruth <- function(dose) {
    z <- exp(alpha0.2 + alpha1.2 * log(dose / refDose))
    return(z / (1 + z))
  }
  
  ## Truth 3: Pessimistic - 0.2 mg achieves target toxicity and OD (80%) reached at 5 mg - make more pessimistic
  DT3  <- 1.5
  DO3  <- 3.25
  pT3  <- (T.upp + T.low) / 2
  pO3 <- 0.6
  
  alpha1.3 <- alpha1fun(DT3, DO3, pT3, pO3)
  alpha0.3 <- alpha0fun(DT3, pT3, alpha1.3)
  
  pesTruth <- function(dose) {
    z <- exp(alpha0.3 + alpha1.3 * log(dose / refDose))
    return(z / (1 + z))
  }
  
 # pesTruth <- probFunction(minInfModel$model, alpha0 = alpha0.3, alpha1 = alpha1.3)
  curve(pesTruth(x), from = 0, to = 10, ylim = c(0, 1))
  
  ## Truth 4: Prior - uses the parameters from the prior model
  priorTruth <- function(dose) {
    #return(minInfModel$model@prob(dose, alpha0=alpha0, alpha1=exp(alpha1)))
    z <- exp(alpha0 + alpha1 * log(dose / refDose))
    return(z / (1 + z))
  }
  
  
  png(filename=file.path(outdir,"scenarios.png"))
  curve(
    relTruth(x),
    from = 0,
    to = 10,
    ylim = c(0, 1),
    main = "Figure 2. Simulation Scenarios",
    xlab = "Dose (mg)",
    ylab = "Probability of DLAE",
    col = "blue",
    lwd = 2
  )
  curve(
    optTruth(x),
    from = 0,
    to = 10,
    ylim = c(0, 1),
    add = TRUE,
    col = "green",
    lwd = 2
  )
  curve(
    pesTruth(x),
    from = 0,
    to = 10,
    ylim = c(0, 1),
    add = TRUE,
    col = "red",
    lwd = 2
  )
  curve(
    priorTruth(x),
    from = 0,
    to = 10,
    ylim = c(0, 1),
    add = TRUE,
    col = "black",
    lwd = 2
  )
  abline(
    h = c(0.2,0.05),
    lty = 2,
    lwd = 1.5,
    col = "grey"
  )
  legend(
    "topleft",
    legend = c("Pessimistic", "Realistic", "Optimistic", "Prior-based"),
    lty = 1,
    col = c("red", "blue", "green", "black"),
    lwd = 2
  )
  dev.off()
  
  # SIMULATIONS-----------------------------------------------------------------
  set.seed(41)
  out.examine <- examine(design)
  
  column_headers1 <- data.frame(
    dose = "Dose (mg)",
    DLTs = "DLAEs",
    nextDose = "Next Dose (mg)",
    stop = "Stop?",
    increment = "Increment (%)",
    stringsAsFactors = FALSE
  )
  
  
  footnote <- paste(
    "The lowest dose of ",
    min(DE.DoseGrid,na.rm=TRUE),
    "mg is an arbitrary dose below the starting dose. It has been included to allow the option for the model to continue if a toxicity is observed at the lowest tested dose."
  )
  
  ht <- as_hux(rbind(column_headers1, out.examine), add_colnames = FALSE) %>%
    set_caption("Table 1. Examination of the Chosen Safety Decision Algoritm") %>%
    set_na_string("0") %>%
    set_all_borders() %>%
    set_align(everywhere, value="center") %>%
    huxtable::add_footnote(footnote) #did not work
  
  ht
  
  # Simulate to obtain operating characteristics
  obtainOC <- function(truth) {
    #T1 <- cbind(DE.DoseGrid, truth(dose = DE.DoseGrid))
    time <- system.time(
      mySims <- crmPack::simulate(
        design,
        seed = 652,
        truth = truth,
        #args = DE.DoseGrid,
        #args=data.frame(dose=DE.DoseGrid),
        nsim = nTrialSims,
        mcmcOptions = myOptions,
        parallel = TRUE
      )
    )
    T1_sum <- list(mySims, summary(mySims, truth = truth))
  }
  
  # Testing individual truths
  # x <- obtainOC(pesTruth) #it does not run
  # x <- obtainOC(optTruth)
  # x <- obtainOC(relTruth)
  # x <- obtainOC(priorTruth)
  
  #  
  # pesTruth(dose=DE.DoseGrid)
  # relTruth(dose=DE.DoseGrid)
  # optTruth(dose=DE.DoseGrid)
  # priorTruth(dose=DE.DoseGrid)
  
  # names(results) <- c("Pessimistic", "Realistic", "Optimistic", "Prior-based")
  # 
  # simSum <- summary(results[["Pessimistic"]][[1]], truth = pesTruth)
  # 
  # print(plot(simSum))
  # 
  # simSum <- summary(results[["Realistic"]][[1]], truth = relTruth)
  # 
  # print(plot(simSum))
  # 
  # simSum <- summary(results[["Optimistic"]][[1]], truth = optTruth)
  # 
  # print(plot(simSum))
  # 
  # simSum <- summary(results[["Prior-based"]][[1]], truth = priorTruth)
  # 
  # print(plot(simSum))
  
  #results <- lapply(list(pesTruth, relTruth, optTruth, priorTruth), obtainOC)
  
  #results <- mclapply(list(pesTruth, relTruth, optTruth, priorTruth), obtainOC, mc.cores=getOption("mc.cores",2L))
  
  
  plan(multisession, workers = 32)
  
  results <- future_lapply(list(pesTruth, relTruth, optTruth, priorTruth), obtainOC,future.seed=TRUE)
  
  # Prep for the percentile row of Table 3
  results2 <- lapply(c(
    results[[1]][[1]]@data,
    results[[2]][[1]]@data,
    results[[3]][[1]]@data,
    results[[4]][[1]]@data
  ), function(z)
    tibble(
      N = z@nObs,
      Dose = z@x,
      Tox = z@y
    )) %>% bind_rows(.id = "Simulation")
  
  results2 <- results2 %>%  mutate(
    Simulation = as.numeric(Simulation),
    truthID = case_when(
      Simulation <= nTrialSims ~ "pes",
      Simulation > nTrialSims &
        Simulation <= 2 * nTrialSims ~ "rel",
      Simulation > 2 * nTrialSims &
        Simulation <= 3 * nTrialSims ~ "opt",
      Simulation > 3 * nTrialSims &
        Simulation <= 4 * nTrialSims ~ "prior",
      TRUE ~ "unknown"
    )
  )
  
  rowx <- results2 %>% group_by(truthID) %>%
    summarise(
      q50 = quantile(Dose, probs = 0.5),
      q80 = quantile(Dose, probs = 0.8),
      q90 = quantile(Dose, probs = 0.9),
      q95 = quantile(Dose, probs = 0.95)
    )
  
  
  # Produce summary table
  row1 <- data.frame(row1 = c(
    pes = paste(pT3 * 100, "% at ", DT3, "mg", sep = ""),
    rel = paste(pT2 * 100, "% at ", DT2, "mg", sep = ""),
    opt = paste(pT1 * 100, "% at ", DT1, "mg", sep = "")
  ))
  
  row2 <- data.frame(row2 = c(
    pes = paste(pO3 * 100, "% at ", DO3, "mg", sep = ""),
    rel = paste(pO2 * 100, "% at ", DO2, "mg", sep = ""),
    opt = paste(pO1 * 100, "% at ", DO1, "mg", sep = "")
  ))
  
  row3 <- data.frame(row3 = c(
    pes = paste(
      median(results[[1]][[2]]@n_obs),
      " (",
      round(quantile(results[[1]][[2]]@n_obs, 0.1), 0),
      ", ",
      round(quantile(results[[1]][[2]]@n_obs, 0.9), 0),
      ")",
      sep = ""
    ),
    rel = paste(
      median(results[[2]][[2]]@n_obs),
      " (",
      round(quantile(results[[2]][[2]]@n_obs, 0.1), 0),
      ", ",
      round(quantile(results[[2]][[2]]@n_obs, 0.9), 0),
      ")",
      sep = ""
    ),
    opt = paste(
      median(results[[3]][[2]]@n_obs),
      " (",
      round(quantile(results[[3]][[2]]@n_obs, 0.1), 0),
      ", ",
      round(quantile(results[[3]][[2]]@n_obs, 0.9), 0),
      ")",
      sep = ""
    )
  ))
  
  row4 <- data.frame(row4 = c(
    pes = paste(
      round(median(results[[1]][[2]]@prop_dlts) * 100, 0),
      "% (",
      round(quantile(results[[1]][[2]]@prop_dlts, 0.1) * 100, 0),
      "%, ",
      round(quantile(results[[1]][[2]]@prop_dlts, 0.9) * 100, 0),
      "%)",
      sep = ""
    ),
    rel = paste(
      round(median(results[[2]][[2]]@prop_dlts) * 100, 0),
      "% (",
      round(quantile(results[[2]][[2]]@prop_dlts, 0.1) * 100, 0),
      "%, ",
      round(quantile(results[[2]][[2]]@prop_dlts, 0.9) * 100, 0),
      "%)",
      sep = ""
    ),
    opt = paste(
      round(median(results[[3]][[2]]@prop_dlts) * 100, 0),
      "% (",
      round(quantile(results[[3]][[2]]@prop_dlts, 0.1) * 100, 0),
      "%, ",
      round(quantile(results[[3]][[2]]@prop_dlts, 0.9) * 100, 0),
      "%)",
      sep = ""
    )
  ))
  
  row5 <- data.frame(row5 = c(
    pes = paste(
      round(median((results[[1]][[2]]@n_above_target / results[[1]][[2]]@n_obs) *
                     100
      ), 0),
      "% (",
      round(quantile((results[[1]][[2]]@n_above_target / results[[1]][[2]]@n_obs) *
                       100, 0.1
      ), 0),
      "%, ",
      round(quantile((results[[1]][[2]]@n_above_target / results[[1]][[2]]@n_obs) *
                       100, 0.9
      ), 0),
      "%)",
      sep = ""
    ),
    rel = paste(
      round(median((results[[2]][[2]]@n_above_target / results[[2]][[2]]@n_obs) *
                     100
      ), 0),
      "% (",
      round(quantile((results[[2]][[2]]@n_above_target / results[[2]][[2]]@n_obs) *
                       100, 0.1
      ), 0),
      "%, ",
      round(quantile((results[[2]][[2]]@n_above_target / results[[2]][[2]]@n_obs) *
                       100, 0.9
      ), 0),
      "%)",
      sep = ""
    ),
    opt = paste(
      round(median((results[[3]][[2]]@n_above_target / results[[3]][[2]]@n_obs) *
                     100
      ), 0),
      "% (",
      round(quantile((results[[3]][[2]]@n_above_target / results[[3]][[2]]@n_obs) *
                       100, 0.1
      ), 0),
      "%, ",
      round(quantile((results[[3]][[2]]@n_above_target / results[[3]][[2]]@n_obs) *
                       100, 0.9
      ), 0),
      "%)",
      sep = ""
    )
  ))
  
  # for row6 it is the mode and not the median displayed
  row6 <- data.frame(row6 = c(
    pes = paste(
      results[[1]][[2]]@dose_most_selected,
      " (",
      round(quantile(results[[1]][[2]]@dose_selected, 0.1), 1),
      ", ",
      round(quantile(results[[1]][[2]]@dose_selected, 0.9), 1),
      ")",
      sep = ""
    ),
    rel = paste(
      results[[2]][[2]]@dose_most_selected,
      " (",
      round(quantile(results[[2]][[2]]@dose_selected, 0.1), 1),
      ", ",
      round(quantile(results[[2]][[2]]@dose_selected, 0.9), 1),
      ")",
      sep = ""
    ),
    opt = paste(
      results[[3]][[2]]@dose_most_selected,
      " (",
      round(quantile(results[[3]][[2]]@dose_selected, 0.1), 1),
      ", ",
      round(quantile(results[[3]][[2]]@dose_selected, 0.9), 1),
      ")",
      sep = ""
    )
  ))
  
  row7 <- data.frame(row7 = c(
    pes = paste(
      round(results[[1]][[2]]@obs_tox_rate_at_dose_most_selected * 100, 0),
      "%",
      sep = ""
    ),
    rel = paste(
      round(results[[2]][[2]]@obs_tox_rate_at_dose_most_selected * 100, 0),
      "%",
      sep = ""
    ),
    opt = paste(
      round(results[[3]][[2]]@obs_tox_rate_at_dose_most_selected * 100, 0),
      "%",
      sep = ""
    )
  )) #%>% mutate(row7 = replace_na(row7, 0))
  
  rounding <- 1
  
  row8 <- data.frame(row8 = c(
    pes = paste(
      round(rowx$q50[rowx$truthID == "pes"], rounding),
      round(rowx$q80[rowx$truthID == "pes"], rounding),
      round(rowx$q90[rowx$truthID == "pes"], rounding),
      round(rowx$q95[rowx$truthID == "pes"], rounding),
      sep = ", "
    ),
    rel = paste(
      round(rowx$q50[rowx$truthID == "rel"], rounding),
      round(rowx$q80[rowx$truthID == "rel"], rounding),
      round(rowx$q90[rowx$truthID == "rel"], rounding),
      round(rowx$q95[rowx$truthID == "rel"], rounding),
      sep = ", "
    ),
    opt = paste(
      round(rowx$q50[rowx$truthID == "opt"], rounding),
      round(rowx$q80[rowx$truthID == "opt"], rounding),
      round(rowx$q90[rowx$truthID == "opt"], rounding),
      round(rowx$q95[rowx$truthID == "opt"], rounding),
      sep = ", "
    )
  ))
  
  row_all <- cbind.data.frame(row1, row2, row3, row4, row5, row6, row7, row8)
  
  row_all <- rownames_to_column(row_all, "scenario") %>%
    rename(
      "Assumed DLAE Rates" = row1,
      "Overall Number of Participants\n enrolled in Study" = row3,
      "Proportion of Participants\n who experience a DLAE" = row4,
      "Proportion of Participants\n Treated at Doses in Overdose Range" = row5,
      "Final Dose Selected (mg)" = row6,
      "Fitted DLAE rate at Dose\n most often selected" = row7,
      "Percentiles of Doses across all\n Simulations (50th, 80th, 90th, 95th)" = row8
    )
  
  
  row_all2 <- row_all %>% pivot_longer(cols = -scenario)
  
  row_all3 <- row_all2 %>% pivot_wider(
    names_from = scenario,
    values_from = value,
    names_repair = "universal"
  )
  
  column_headers2 <- data.frame(
    name = " ",
    pes = "Pessimistic",
    rel = "Realistic",
    opt = "Optimistic",
    stringsAsFactors = FALSE
  )
  
  # Final Table
  ht2 <- as_hux(rbind(column_headers2, row_all3), add_colnames = FALSE) %>%
    set_caption("Table 2. Summary of Operating Characteristics") %>%
    merge_cells(2:3, 1) %>%
    set_align(everywhere, -1, "center") %>%
    set_bold(1, 1:ncol(row_all3), TRUE) %>%
    set_right_border(0:nrow(row_all3) + 1, 1:ncol(row_all3), 1) %>%
    set_outer_borders(0:nrow(row_all3) + 1, 1:ncol(row_all3), 1) %>%
    set_bottom_border(1, 1:ncol(row_all3), 1) %>%
    set_bottom_border(3:nrow(row_all3), 1:ncol(row_all3), 1) %>%
    huxtable::add_footnote(
      "*Median (10th and 90th percentile) for all rows except for 'Final Dose Selected', where the mode is displayed instead of the median, and for 'Percentiles of Doses across all Simulations', where the 50th, 80th, 90th and 95th percentiles are displayed.*"
    )
  
  ht2
  
  #SAVING SESSION---------------------------------------------------------------
  sesInfo <- sessionInfo()
  
  objects <- "sesInfo"
  
  objects <-
    c(
      objects,
      grep(
        "*dir|*.ls|*.df|*.plot|*.lm|*.tb|ht2|nTrialSims|*.examine|toxicities|doses|coarseGrid|refdose|myNextBest|myStopping1|myStopping2|myStopping3|myStopping|mySize|mySizePL|myIncrements|options|model|design|SimDesign|truthList|trialData|nsim|results|outdir|alpha0|alpha1|cov*|threshmin|probmin|threshmax|probmax|burning|steps|samples",
        x = ls(),
        value = TRUE
      )
    )
  
  save(list = objects, file = saved.image)
  
}

```

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Prototype code that could go into simulations table generation #878

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Prototype code that could go into simulations table generation #878

Description

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions