Suggestion: DLT profile plot and Prob(DLT) plot with shaded loss areas #783

marlenesg · 2024-01-25T14:33:41Z

Summary
The suggested additions include three functions:

get_results(): prepares the posterior samples for the following two functions
plot_DLT_profile(): returns a facet plot with one facet per cohort, showing DLT yes/no/not evaluable for each patient
plot_DLT_probs(): returns a plot for Prob(DLT) with CI for each dose of the dose grid. The different loss areas are shaded in different colors

get_results()

#' Get results
#'
#' Takes a \code{crmPack::Data} object which is obtained by function \code{Data},
#' the posterior samples produced by function \code{mcmc} and the recommended dose
#' provided by function \code{nextBest} and returns a data.frame that contains
#' various quantiles of the DLT probability for each dose, as well as the values
#' of the loss function.
#'
#' @param data A \code{Data} object, which includes all evaluable subjects.
#' @param posterior_samples A \code{crmPack} object produced by function \code{mcmc}.
#' @param recommended_dose A \code{crmPack} object produced by function \code{nextBest}.
#' @return A data.frame containing various quantiles of the DLT probability
#' for each dose, as well as the values of the loss function.
#' @examples
#' \dontrun{
#' dose_grid <- c(10, 15, 30, 45, 60, 80)
#' unit <- "mg/kg"
#' ## Specify the observed data
#' data <- Data(
#'   x = c(rep(10, 3), rep(15, 3), rep(30, 3)),
#'   y = c(rep(0, 3), rep(0, 2), 1, rep(0, 2), 1),
#'   cohort = c(rep(1, 3), rep(2, 3), rep(3, 3)),
#'   doseGrid = dose_grid,
#'   ID = 1:9
#' )
#' mcmc_options <- McmcOptions(
#'   burnin = 5000,
#'   step = 2,
#'   samples = 100000,
#'   rng_kind = "Wichmann-Hill",
#'   rng_seed = 1
#' )
#' model_bcrm <- LogisticLogNormal(
#'   mean = c(-0.8, -0.09),
#'   cov = matrix(c(1.5^2, 0, 0, 1.2^2), nrow = 2),
#'   ref_dose = 45
#' )
#' }

#' increment1 <- IncrementsRelative(
#'  intervals = c(0, 60),
#'  increments = c(1, 0.67)
#' )

#' increment2 <- IncrementsRelativeDLT(
#'  dlt_intervals = c(0, 1),
#'  increments = c(1, 0.50)
#' )
#'
#' combIncrement <- IncrementsMin(increments_list = list(increment1, increment2))

#' ncrm_loss <- NextBestNCRMLoss(
#'  target = c(0.2, 0.35),
#'  overdose = c(0.35, 0.6),
#'  unacceptable = c(0.6, 1),
#'  max_overdose_prob = 0.25,
#'  losses = c(1, 0, 2, 3)
#' )

#' postSamples <- mcmc(data, model_bcrm, mcmc_options)

#' dose_rec <- nextBest(ncrm_loss,
#'  doselimit = maxDose(combIncrement, data),
#'  postSamples, model_bcrm, data
#' )
#'
#' get_results(
#'  data = data,
#'  posterior_samples = postSamples,
#'  recommended_dose = dose_rec)
#' }
#' @source This function uses \code{ggplot} function from \code{ggplot2}
#' R package.
#' @author Burak Kuersad Guenhan, Marianna Grinberg, Marlene Schulte-Goebel
#' @seealso \code{ggplot2::ggplot}
#' @export
get_results <- function(data, posterior_samples, recommended_dose) {
  prob_samples_mat <- matrix(
    nrow = size(posterior_samples@options),
    ncol = data@nGrid
  )

  # evaluate the probs, for all samples
  for (i in seq_len(data@nGrid)) {
    prob_samples_mat[, i] <- prob(
      dose = data@doseGrid[i],
      model_bcrm,
      posterior_samples
    )
  }

  pq025 <- apply(prob_samples_mat, 2, function(x) stats::quantile(x, 0.025))
  pq25 <- apply(prob_samples_mat, 2, function(x) stats::quantile(x, 0.25))
  pq50 <- apply(prob_samples_mat, 2, function(x) stats::quantile(x, 0.50))
  pq75 <- apply(prob_samples_mat, 2, function(x) stats::quantile(x, 0.75))
  pq95 <- apply(prob_samples_mat, 2, function(x) stats::quantile(x, 0.95))
  pq975 <- apply(prob_samples_mat, 2, function(x) stats::quantile(x, 0.975))


  results <- data.frame(
    dose = data@doseGrid,
    Quant025 = pq025 * 100,
    Quant25 = pq25 * 100,
    Quant50 = pq50 * 100,
    Quant75 = pq75 * 100,
    Quant975 = pq975 * 100,
    loss.values = as.vector(recommended_dose$probs[, "posterior_loss"]),
    p.overdose = (as.vector(recommended_dose$probs[, "excessive"]) + recommended_dose$probs[, "unacceptable"]) * 100,
    dose.rec = recommended_dose$value
  )

  return(results)
}

plot_DLT_profile()

#' Plot DLT profiles
#'
#' Takes a \code{crmPack::Data} object which is obtained by function \code{Data} and plots
#' DLT profiles, showing dose levels for each subjects, and whether DLT is observed
#' or not.
#'
#' @param data A \code{Data} object, which includes all evaluable subjects.
#' @param data_NE A data.frame, which includes all non-evaluable subjects.
#' @param unit A string specifying the dose unit used in the trials, for example "mg/kg".
#' @return The return value is invisible \code{NULL}.
#' @examples
#' \dontrun{
#' dose_grid <- c(10, 15, 30, 45, 60, 80)
#' unit <- "mg/kg"
#' ## Specify the observed data
#' data <- Data(
#'   x = c(rep(10, 3), rep(15, 3), rep(30, 3)),
#'   y = c(rep(0, 3), rep(0, 2), 1, rep(0, 2), 1),
#'   cohort = c(rep(1, 3), rep(2, 3), rep(3, 3)),
#'   doseGrid = dose_grid,
#'   ID = 1:9
#' )
#' ## Specify the NON-evaluable data
#' ## (if none, set data_NE <- NULL)
#' data_NE <- data.frame(
#'   IDs_NE = 10,
#'   doses_NE = 30,
#'   dlts_NE = 2,
#'   cohorts_NE = 3
#' )
#'
#' plot_DLT_profile(
#'   data = data,
#'   data_ne = data_NE,
#'   unit = unit
#' )
#' }
#' @source This function uses \code{ggplot} function from \code{ggplot2}
#' R package.
#' @author Burak Kuersad Guenhan, Marianna Grinberg, Marlene Schulte-Goebel
#' @seealso \code{ggplot2::ggplot}
#' @export
plot_DLT_profile <- function(data = data,
                             data_NE = data_NE,
                             unit) {
  # combine evaluable patients in 'data' with not-evaluable patients in 'data_NE'
  ID <- c(data@ID, data_NE$IDs_NE)
  dose <- c(data@x, data_NE$doses_NE)
  cohort <- c(data@cohort, data_NE$cohorts_NE)
  DLT <- c(data@y, data_NE$dlts_NE)

  ## build data frame
  dat <- data.frame(
    "ID" = ID,
    "dose" = dose,
    "DLT" = factor(DLT,
      levels = c("0", "1", "2"),
      labels = c("DLT No", "DLT Yes", "Not evaluable")
    ),
    "cohort" = paste("Cohort", cohort)
  )
  dat <- dat[order(dat$ID), ]

  ## DLT profile plot
  p <- ggplot(data = dat) +
    geom_point(aes(x = factor(ID, levels = unique(ID[order(cohort, ID)])), y = dose, shape = DLT, color = DLT), size = 2) +
    scale_shape_manual(values = c(
      "DLT No" = 19, "DLT Yes" = 17,
      "Not evaluable" = 0
    ), drop = FALSE) +
    scale_color_manual(
      values = c(
        "DLT No" = "black",
        "DLT Yes" = "red",
        "Not evaluable" = "blue"
      ),
      drop = FALSE
    ) +
    scale_x_discrete(breaks = dat$ID, labels = sort(dat$ID)) +
    scale_y_continuous(
      limits = c(0, max(data@doseGrid)),
      breaks = data@doseGrid,
      labels = data@doseGrid
    ) +
    facet_wrap(. ~ cohort, strip.position = "bottom", scales = "free_x") +
    ggtitle("DLT Profile Plot") +
    xlab("Subject IDs") +
    ylab(paste0("Dose (", unit, ")")) +
    theme_bw() +
    theme(plot.title = element_text(hjust = 0.5)) +
    theme(legend.title = element_blank()) +
    theme(
      legend.background = element_blank(),
      legend.box.background = element_rect(colour = "black"),
      axis.text.x = element_text(angle = 45, hjust = 1)
    )

  p
}

plot_DLT_probs()

#' Plot DLT probabilities
#'
#' Takes a \code{crmPack::Data} object which is obtained by function \code{Data} and plots
#' DLT profiles, showing dose levels for each subjects, and whether DLT is observed
#' or not.
#'
#' @param results A data.frame obtained by \code{get_results}.
#' @param loss_intervals A vector containing the loss interval bounds,
#' for example c(0.2, 0.35, 0.6, 1) which represents the four intervals (0, 0.2),
#' (0.2, 0.35), (0.35, 0.6), (0.6, 1).
#' @param loss_values A vector containing the loss values (as strings) for each
#' interval, for example c("1", "0", "1", "2").
#' @param unit A string specifying the dose unit used in the trials, for example "mg/kg".
#' @return The return value is invisible \code{NULL}.
#' @examples
#' \dontrun{
#' dose_grid <- c(10, 15, 30, 45, 60, 80)
#' unit <- "mg/kg"
#' ## Specify the observed data
#' data <- Data(
#'   x = c(rep(10, 3), rep(15, 3), rep(30, 3)),
#'   y = c(rep(0, 3), rep(0, 2), 1, rep(0, 2), 1),
#'   cohort = c(rep(1, 3), rep(2, 3), rep(3, 3)),
#'   doseGrid = dose_grid,
#'   ID = 1:9
#' )
#' mcmc_options <- McmcOptions(
#'   burnin = 5000,
#'   step = 2,
#'   samples = 100000,
#'   rng_kind = "Wichmann-Hill",
#'   rng_seed = 1
#' )
#' model_bcrm <- LogisticLogNormal(
#'   mean = c(-0.8, -0.09),
#'   cov = matrix(c(1.5^2, 0, 0, 1.2^2), nrow = 2),
#'   ref_dose = 45
#' )
#' }

#' increment1 <- IncrementsRelative(
#'  intervals = c(0, 60),
#'  increments = c(1, 0.67)
#' )

#' increment2 <- IncrementsRelativeDLT(
#'  dlt_intervals = c(0, 1),
#'  increments = c(1, 0.50)
#' )
#'
#' combIncrement <- IncrementsMin(increments_list = list(increment1, increment2))

#' ncrm_loss <- NextBestNCRMLoss(
#'  target = c(0.2, 0.35),
#'  overdose = c(0.35, 0.6),
#'  unacceptable = c(0.6, 1),
#'  max_overdose_prob = 0.25,
#'  losses = c(1, 0, 2, 3)
#' )

#' postSamples <- mcmc(data, model_bcrm, mcmc_options)

#' dose_rec <- nextBest(ncrm_loss,
#'  doselimit = maxDose(combIncrement, data),
#'  postSamples, model_bcrm, data
#' )
#'
#' get_results(
#'  data = data,
#'  posterior_samples = postSamples,
#'  recommended_dose = dose_rec)
#' }
#'
#' loss_intervals <- c(0.2, 0.35, 0.6, 1)
#' loss_values <- c("1", "0", "1", "2")
#'
#' plot_DLT_probs(
#'  results = results,
#'  loss_intervals = loss_intervals * 100,
#'  loss_values = loss_values,
#'  unit = "mg/kg"
#'  )
#'
#' }
#' @source This function uses \code{ggplot} function from \code{ggplot2}
#' R package.
#' @author Burak Kuersad Guenhan, Marianna Grinberg, Marlene Schulte-Goebel
#' @seealso \code{ggplot2::ggplot}, \code{crmPack::get_results}
#' @export
plot_DLT_probs <- function(results,
                           loss_intervals,
                           loss_values,
                           unit) {
  dat <- data.frame(
    x = results$dose,
    y = results$Quant50,
    ylo = results$Quant025,
    yup = results$Quant975,
    yinlo = results$Quant25,
    yinup = results$Quant75
  )

  args_inner <- list(
    mapping = aes_(
      ymax = ~yinup,
      ymin = ~yinlo,
      x = ~x
    ),
    size = 1.5,
    show.legend = FALSE,
    position = position_dodge(width = 0.50),
    data = dat
  )

  layer_inner <- do.call(geom_linerange, args_inner)


  p <- ggplot() +
    geom_pointrange(
      position = position_dodge(width = 0.50),
      data = dat,
      mapping = aes(
        x = x,
        y = y,
        ymin = ylo,
        ymax = yup
      )
    ) +
    layer_inner +
    geom_hline(aes(yintercept = loss_intervals[2]),
      lty = 2
    ) +
    geom_hline(aes(yintercept = loss_intervals[3]),
      lty = 2
    ) +
    theme(
      legend.justification = c(0, 1),
      legend.position = c(0, 1)
    ) +
    theme(legend.title.align = 0.5) +
    theme_bw() +
    theme(plot.title = element_text(hjust = 0.5)) +
    scale_y_continuous(breaks = seq(0, 100, 20)) +
    scale_x_continuous(breaks = results$dose) +
    ylab("Probability of DLT (in %)") +
    xlab(paste0("Dose (", unit, ")")) +
    geom_rect(
      aes(
        ymin = loss_intervals[1], ymax = loss_intervals[2],
        xmin = -Inf, xmax = Inf, fill = loss_values[2]
      ),
      alpha = 0.01
    ) +
    geom_rect(
      aes(
        ymin = 0, ymax = loss_intervals[1],
        xmin = -Inf, xmax = Inf, fill = loss_values[1]
      ),
      alpha = 0.2
    ) +
    geom_rect(
      aes(
        ymin = loss_intervals[2], ymax = loss_intervals[3],
        xmin = -Inf, fill = loss_values[3], xmax = Inf
      ),
      alpha = 0.2
    ) +
    geom_rect(
      aes(
        ymin = loss_intervals[3], ymax = loss_intervals[4],
        xmin = -Inf, fill = loss_values[4], xmax = Inf
      ),
      alpha = 0.5
    ) +
    scale_fill_manual("Loss",
      values = rep("red", 3),
      guide = guide_legend(override.aes = list(alpha = c(0.001, 0.05, 0.2)))
    ) +
    theme(legend.box.background = ggplot2::element_rect(colour = "black")) +
    ggtitle("Results from the Bayesian model") +
    theme(text = element_text(size = 10))

  p
}

The text was updated successfully, but these errors were encountered:

marlenesg · 2024-01-25T14:42:58Z

Hi @danielinteractive, this is the example code I mentioned for the DLT-profile-plot and the Prob(DLT)-plot. Could you please let me know where in the repo this code would fit best? Naming of the functions and variables etc. isn't final, I am open for suggestions.

danielinteractive · 2024-01-26T08:46:34Z

Thanks a lot @marlenesg !
I would say on first sight:

plot_DLT_profile could become a specific method for Data objects, and therefore live in Data-methods.R
plot_DLT_probs could become a specific method for Samples objects, and then live in Samples-methods.R
get_results can become a helper function for plot_DLT_probs and be in helpers_samples.R
- could also be replaced by one/more tidy methods that we have now

marlenesg added the enhancement New feature or request label Jan 25, 2024

marlenesg linked a pull request Feb 9, 2024 that will close this issue

783 profile plots #792

Open

danielinteractive added this to To do in crmPack Board via automation Feb 26, 2024

danielinteractive assigned marlenesg Feb 26, 2024

danielinteractive moved this from To do to In progress in crmPack Board Feb 26, 2024

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Suggestion: DLT profile plot and Prob(DLT) plot with shaded loss areas #783

Suggestion: DLT profile plot and Prob(DLT) plot with shaded loss areas #783

marlenesg commented Jan 25, 2024

marlenesg commented Jan 25, 2024

danielinteractive commented Jan 26, 2024 •

edited

Loading

Suggestion: DLT profile plot and Prob(DLT) plot with shaded loss areas #783

Suggestion: DLT profile plot and Prob(DLT) plot with shaded loss areas #783

Comments

marlenesg commented Jan 25, 2024

marlenesg commented Jan 25, 2024

danielinteractive commented Jan 26, 2024 • edited Loading

danielinteractive commented Jan 26, 2024 •

edited

Loading