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

[marlenesg]

Summary
The suggested additions include three functions:

1. get_results(): prepares the posterior samples for the following two functions
2. plot_DLT_profile(): returns a facet plot with one facet per cohort, showing DLT yes/no/not evaluable for each patient
3. 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
}

[marlenesg]

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]

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