14 oc post prob (#16)

* cleaning up dbetabinom * removed rounding, not needed * rmarkdown language and assert for dbetabinom * changed x and m to number instead of numeric * test * pbetaMix function and test * added new tests and changed documentation * upper case for test_that comments and more tests for qbeta * test_thats for pbetaMix and qbetaMix * tidy syntax and last checks * tidy syntax and last checks * Update R/dbetabinom.R Co-authored-by: Daniel Sabanes Bove <[email protected]> * [skip actions] Roxygen Man Pages Auto Update * Update R/dbetabinom.R Co-authored-by: Daniel Sabanes Bove <[email protected]> * Update R/dbetabinom.R Co-authored-by: Daniel Sabanes Bove <[email protected]> * Update R/dbetabinom.R Co-authored-by: Daniel Sabanes Bove <[email protected]> * Update R/dbetabinom.R Co-authored-by: Daniel Sabanes Bove <[email protected]> * Update R/dbetabinom.R Co-authored-by: Daniel Sabanes Bove <[email protected]> * Update R/dbetabinom.R Co-authored-by: Daniel Sabanes Bove <[email protected]> * Update R/dbetabinom.R Co-authored-by: Daniel Sabanes Bove <[email protected]> * Update examples/pbetaMix.R Co-authored-by: Daniel Sabanes Bove <[email protected]> * p and q inputs * q,p changes * Update R/dbetabinom.R Co-authored-by: Daniel Sabanes Bove <[email protected]> * [skip actions] Roxygen Man Pages Auto Update * documentation to reflect p and q changes * amending qt * intro x-> q * amented postprob file * amended postprob file * fixed mistake on x * qbetaMix was in postprobDist * trying to find R CMD check error * [skip actions] Roxygen Man Pages Auto Update * manual of ocPostprobDist * fixing roxygen type to p * see if this fixes the CMD check * see if fixes CMD check error * xi removed * example post prob works * clean * assert_number * in postprobDist, I changed p-> q * change p to q * postprobold * postprobOld * clean post prob * nice titles * documentations * clean * clean * clean * clean * clean * clean * some examples and tests * clean * fix R CMD check errors * fix R CMD check errors for PR * Update R/ocPostprob.R Co-authored-by: Daniel Sabanes Bove <[email protected]> * Update R/ocPostprob.R Co-authored-by: Daniel Sabanes Bove <[email protected]> * [skip actions] Roxygen Man Pages Auto Update * Update R/ocPostprob.R Co-authored-by: Daniel Sabanes Bove <[email protected]> * Update R/ocPostprob.R Co-authored-by: Daniel Sabanes Bove <[email protected]> * [skip actions] Roxygen Man Pages Auto Update * Update R/ocPostprob.R Co-authored-by: Daniel Sabanes Bove <[email protected]> * PR feedback * [skip actions] Roxygen Man Pages Auto Update * clean * [skip actions] Roxygen Man Pages Auto Update * clean * Revert "[skip actions] Roxygen Man Pages Auto Update" This reverts commit 35126ab. * Revert "clean" This reverts commit 55939e9. * clean * clean * Update R/ocPostprob.R Co-authored-by: Daniel Sabanes Bove <[email protected]> * Update R/ocPostprob.R Co-authored-by: Daniel Sabanes Bove <[email protected]> * Update R/ocPostprob.R Co-authored-by: Daniel Sabanes Bove <[email protected]> * Update R/ocPostprob.R Co-authored-by: Daniel Sabanes Bove <[email protected]> * Update R/postprob.R Co-authored-by: Daniel Sabanes Bove <[email protected]> * Update R/postprob.R Co-authored-by: Daniel Sabanes Bove <[email protected]> * Update tests/testthat/test-ocPostprob.R Co-authored-by: Daniel Sabanes Bove <[email protected]> * clean * [skip actions] Roxygen Man Pages Auto Update * clean * clean * clean * Update tests/testthat/test-ocPostprob.R Co-authored-by: Daniel Sabanes Bove <[email protected]> * clean * [skip actions] Roxygen Man Pages Auto Update * clean --------- Co-authored-by: Daniel Sabanes Bove <[email protected]> Co-authored-by: 27856297+dependabot-preview[bot]@users.noreply.github.com <27856297+dependabot-preview[bot]@users.noreply.github.com>
Genentech · Oct 13, 2023 · b1ad258 · b1ad258
1 parent 6ee771d
commit b1ad258
Show file tree

Hide file tree

Showing 17 changed files with 697 additions and 4,310 deletions.
diff --git a/R/ocPostprob.R b/R/ocPostprob.R
@@ -1,114 +1,263 @@
 #' @include postprob.R
 NULL
 
-#' Calculate operating characteristics for posterior probability method
+#' Generating random distance in given looks for sample sizes for Efficacy and Futility.
 #'
-#' The trial is stopped for efficacy if the posterior probability to be
-#' above p1 is larger than tU, and stopped for futility if the posterior
-#' probability to be below p0 is larger than tL.
+#' A helper function for `ocPostprob` to generate random distance's wiggle room around looks `nn`.
+#' Numeric looks `nn` must be of minimum two elements and will generate `length(nn)-1` distances.
 #'
-#' Returned operating characteristics in a matrix include:
-#' ExpectedN: expected number of patients in the trials
-#' PrStopEarly: probability to stop the trial early (before reaching the
+#' @typed nn : number or numeric
+#'  the union of `nnE` and `nnF` (if futility analysis or looks exists) supplied.
+#'
+#' @return A numeric with `length(nn)-1` elements.
+#'
+#' @keywords internal
+#'
+h_get_distance <- function(nn) {
+  assert_numeric(nn, unique = TRUE, sorted = TRUE, min.len = 1)
+  dist0 <- floor(min(nn - c(0, nn[-length(nn)])) / 2)
+  assert_numeric(dist0, sorted = TRUE)
+  sample(-dist0:dist0,
+    size = length(nn) - 1,
+    replace = TRUE,
+    prob = 2^(-abs(-dist0:dist0) / 2)
+  )
+}
+
+#' Generating looks from random distance
+#'
+#' A helper function for `ocPostprob` that applies the numeric element of `dist` to looks `nn`.
+#'
+#' @typed dist : numeric or logical
+#'  distance for random looks around the look locations in `nn`,
+#'  where `dist`is generated from `h_get_distance` in a numeric of at least one element.
+#'  If `NULL`, only one location look will be set at `nnE` or `nnF`.
+#' @typed nnE : numeric
+#'  sample size or sizes where study can be stopped for Efficacy decision. If different for Futility decision,
+#'  specify in `nnF`.
+#' @typed nnF : numeric
+#'  sample size or sizes where study can be stopped for Futility decision if different from Efficacy decision.
+#'
+#' @return Uses distance from `h_get_distance` to add to looks, creating wiggled looks:
+#'  - `nnrE` is the result for Efficacy looks with random distance added.
+#'  - `nnrF` is the result for Futility looks with random distance added.
+#'
+#' @keywords internal
+#'
+h_get_looks <- function(dist, nnE, nnF) {
+  assert_numeric(nnE)
+  assert_numeric(nnF)
+  nn <- unique(c(nnE, nnF))
+  assert_numeric(nn)
+  assert_numeric(dist)
+  nnr <- nn + c(dist, 0)
+  list(
+    nnrE = nnr[nn %in% nnE],
+    nnrF = nnr[nn %in% nnF]
+  )
+}
+
+#' Generating random decision and sample size looks.
+#'
+#' A helper function for `ocPostprob` to generate numeric of decisions `decisions` and random looks `all_sizes`.
+#'
+#' @inheritParams h_get_looks
+#' @typed nnr : numeric
+#'  union of `nnE`and `nnF`.
+#' @typed response : numeric
+#'  a numeric of Bernoulli successes based on `size_look`.
+#' @typed truep : number
+#'  assumed true response rate or true rate (scenario).
+#' @typed p0 : number
+#'  lower Futility threshold of response rate.
+#' @typed p1 : number
+#'  upper Efficacy threshold of response rate.
+#' @typed tL : number
+#'  posterior probability threshold for being below `p0`.
+#' @typed tU : number
+#'  posterior probability threshold for being above `p1`.
+#' @typed parE : numeric
+#'  alpha and beta parameters for the prior on the treatment proportion.
+#'  Default set at alpha = 1, beta = 1, or uniform prior.
+#'
+#' @return A list of the following objects :
+#'  - `decision` : resulting numeric of decision, one of `TRUE` for Go, `FALSE` for Stop, `NA` for Gray zone.
+#'  - `all_sizes` : resulting numeric of look size, anything below maximum
+#'   look size is an indicated interim, Futility or Efficacy or both.
+#'
+#' @keywords internal
+#'
+h_get_decision <- function(nnr, response, truep, p0, p1, parE = c(1, 1), nnE, nnF, tL, tU) {
+  index_look <- 1
+  assert_numeric(nnr)
+  size_look <- nnr[index_look]
+  all_sizes <- decision <- NA
+  response <- stats::rbinom(max(nnr), size = 1, truep)
+  assert_numeric(response, lower = 0, upper = 1)
+  while (is.na(decision) && index_look <= length(nnr)) {
+    if (size_look %in% nnF) {
+      qL <- 1 - postprob(x = sum(response[1:size_look]), n = size_look, p = p0, parE = parE)
+      assert_number(qL, lower = 0, upper = 1)
+      decision <- ifelse(qL >= tL, FALSE, NA)
+    }
+    if (size_look %in% nnE) {
+      qU <- postprob(x = sum(response[1:size_look]), n = size_look, p = p1, parE = parE)
+      assert_number(qU, lower = 0, upper = 1)
+      decision <- ifelse(qU < tU, decision, TRUE)
+    }
+    all_sizes <- size_look
+    index_look <- index_look + 1
+    size_look <- nnr[index_look]
+  }
+  list(
+    decision = decision,
+    all_sizes = all_sizes
+  )
+}
+
+#' Creating list for operating characteristics.
+#'
+#' Generates operating characteristics.
+#'
+#' @inheritParams h_get_looks
+#' @inheritParams h_get_decision
+#' @typed nnrE : numeric
+#'  looks with random distance, if applied on `nnE`.
+#' @typed nnrF : numeric
+#'  looks with random distance, if applied on `nnF`.
+#' @typed all_sizes : numeric
+#'  sample sizes of all looks simulated `length(sim)` times if `dist` applied.
+#' @typed decision : numeric
+#'  Go, Stop or Gray Zone decisions of all looks simulated `length(sim)` times.
+#'
+#' @return A list of results containing :
+#'
+#' - `ExpectedN`: expected number of patients in the trials
+#' - `PrStopEarly`: probability to stop the trial early (before reaching the
 #' maximum sample size)
-#' PrEarlyEff: probability to decide for efficacy early
-#' PrEarlyFut: probability to decide for futility early
-#' PrEfficacy: probability to decide for efficacy
-#' PrFutility: probability to decide for futility
-#' PrGrayZone: probability of no decision at the end ("gray zone")
-#'
-#' @param nn vector of look locations for efficacy
-#' (if futility looks should be different, please specify also \code{nnF})
-#' @param p true rate (scenario)
-#' @param p0 lower efficacy threshold
-#' @param p1 upper efficacy threshold
-#' @param tL probability threshold for being below p0
-#' @param tU probability threshold for being above p1
-#' @param parE beta parameters for the prior on the treatment proportion
-#' @param ns number of simulations
-#' @param nr generate random look locations? (not default)
-#' @param d distance for random looks around the look locations in \code{nn}
-#' @param nnF vector of look locations for futility
-#' (default: same as efficacy)
+#' - `PrEarlyEff`: probability of Early Go decision
+#' - `PrEarlyFut`: probability of for Early Stop decision
+#' - `PrEfficacy`: probability of Go decision
+#' - `PrFutility`: probability of Stop decision
+#' - `PrGrayZone`: probability between Go and Stop ,"Evaluate" or Gray decision zone
+#'
+#' @keywords internal
+#'
+h_get_oc <- function(all_sizes, nnr, decision, nnrE, nnrF) {
+  sim <- length(all_sizes)
+  assert_logical(decision, len = sim)
+  assert_numeric(all_sizes)
+  assert_numeric(nnrE, lower = 0, upper = max(nnrE))
+  assert_numeric(nnrF, lower = 0, upper = max(nnrF))
+  data.frame(
+    ExpectedN = mean(all_sizes, na.rm = TRUE),
+    PrStopEarly = mean(all_sizes < max(nnrF), na.rm = TRUE),
+    PrEarlyEff = sum(decision * (all_sizes < max(nnrE)), na.rm = TRUE) / sim,
+    PrEarlyFut = sum((1 - decision) * (all_sizes < max(nnrF)), na.rm = TRUE) / sim,
+    PrEfficacy = sum(decision, na.rm = TRUE) / sim,
+    PrFutility = sum(1 - decision, na.rm = TRUE) / sim,
+    PrGrayZone = sum(is.na(decision)) / sim
+  )
+}
+
+#' Operating Characteristics for Posterior Probability method
+#'
+#' @description `r lifecycle::badge("experimental")`
+#'
+#' Calculate operating characteristics for posterior probability method.
+#'
+#' It is assumed that the true response rate is `truep`.
+#' The trial is stopped for Efficacy if the posterior probability to be
+#' above `p1` is larger than `tU`, and stopped for Futility if the posterior
+#' probability to be below `p0` is larger than `tL`:
+#'
+#' Stop criteria for Efficacy :
+#'
+#' `P_E(p > p1) > tU`
+#'
+#' Stop criteria for Futility :
+#'
+#' `P_E(p < p0) > tL`
+#'
+#' Resulting operating characteristics include the following:
+#'
+#' - `ExpectedN`: expected number of patients in the trials
+#' - `PrStopEarly`: probability to stop the trial early (before reaching the
+#' maximum sample size)
+#' - `PrEarlyEff`: probability of Early Go decision
+#' - `PrEarlyFut`: probability of for Early Stop decision
+#' - `PrEfficacy`: probability of Go decision
+#' - `PrFutility`: probability of Stop decision
+#' - `PrGrayZone`: probability between Go and Stop ,"Evaluate" or Gray decision zone
+#'
+#' @inheritParams h_get_looks
+#' @inheritParams h_get_decision
+#' @typed sim : number
+#'  number of simulations.
+#' @typed wiggle : logical
+#'  generate random look locations (not default).
+#'  if `TRUE`, optional to specify `dist` (see @details).
+#' @typed randomdist : logical
+#'  Random distance added to looks. if `NULL`, and `wiggle = TRUE`, function will
+#'  generate and add a random distance within range of the closest looks.
+#'
 #' @return A list with the following elements:
-#' oc: matrix with operating characteristics (see Details section)
-#' Decision: vector of the decisions made in the simulated trials
-#' (\code{TRUE} for success, \code{FALSE} for failure, \code{NA} for no
-#' decision)
-#' SampleSize: vector of the sample sizes in the simulated trials
-#' nn: vector of look locations that was supplied
-#' nnE: vector of efficacy look locations
-#' nnF: vector of futility look locations
-#' params: multiple parameters
+#'
+#' - `oc`: matrix with operating characteristics (see @details section)
+#' - `nn`: vector of look locations that was supplied
+#' - `nnE`: vector of Efficacy look locations
+#' - `nnF`: vector of Futility look locations
+#' - `params`: multiple parameters
+#'
+#' @details
+#' `ExpectedN` is an average of the simulated sample sizes.
+#'  If `wiggle = TRUE`, one can specify `dist`, though the algorithm will generate it if `dist = NULL`.
+#'  If `nnF = NULL`, no Futility or decision to Stop will be analysed. Note that `nnF = c(0)` is equivalent.
+#'  As default, `nnF` is set to the identical looks of `nnE`, and if `wiggle = TRUE`, all looks are the same, e.g.
+#'  `nnE = nnF` when wiggle and distance is applied.
 #'
 #' @example examples/ocPostprob.R
 #' @export
-ocPostprob <- function(nn, p, p0, p1, tL, tU, parE = c(1, 1),
-                       ns = 10000, nr = FALSE, d = NULL, nnF = nn) {
-  # Calculate operating characteristics via simulation
-  # nn: vector of look locations
-  # s: decision reject H0 (TRUE) or fail to reject (FALSE)
-  #    during trial if continuing (NA)
-
-  ## copy nn to nnE:
-  nnE <- sort(nn)
-  nnF <- sort(nnF)
-  s <- rep(NA, ns)
-  n <- s
+ocPostprob <- function(nnE, truep, p0, p1, tL, tU, parE = c(1, 1),
+                       sim = 50000, wiggle = FALSE, randomdist = NULL, nnF = nnE) {
   nn <- sort(unique(c(nnF, nnE)))
-  nL <- length(nn)
-  Nstart <- nn[1]
-  Nmax <- nn[nL]
-  if (nr && is.null(d)) {
-    # set parameter d for randomly generating look locations
-    d <- floor(min(nn - c(0, nn[-nL])) / 2)
+  assert_number(sim, lower = 1, finite = TRUE)
+  if (sim < 50000) {
+    warning("Advise to use sim >= 50000 to achieve convergence")
   }
-  nnr <- nn
-  nnrE <- nnE
-  nnrF <- nnF
-  for (k in 1:ns) {
-    # simulate a clinical trial ns times
-    if (nr && (d > 0)) {
-      # randomly generate look locations
-      dd <- sample(-d:d,
-        size = nL - 1, replace = TRUE,
-        prob = 2^(c(-d:0, rev(-d:(-1))) / 2)
-      )
-      nnr <- nn + c(dd, 0)
-
-      nnrE <- nnr[nn %in% nnE]
-      nnrF <- nnr[nn %in% nnF]
-    }
-    x <- stats::rbinom(Nmax, 1, p)
-    j <- 1
-    i <- nnr[j]
-    while (is.na(s[k]) && (j <= length(nnr))) {
-      if (i %in% nnrF) {
-        qL <- 1 - postprob(x = sum(x[1:i]), n = i, p = p0, parE = parE)
-        s[k] <- ifelse(qL >= tL, FALSE, NA)
-      }
-
-      if (i %in% nnrE) {
-        qU <- postprob(x = sum(x[1:i]), n = i, p = p1, parE = parE)
-        s[k] <- ifelse(qU < tU, s[k], TRUE)
-      }
-
-      n[k] <- i
-      j <- j + 1
-      i <- nnr[j]
+  decision <- vector(length = sim)
+  all_sizes <- vector(length = sim)
+  for (k in seq_len(sim)) {
+    if (length(nn) != 1 && wiggle && is.null(randomdist)) {
+      dist <- h_get_distance(nn = nn)
+      nnr <- h_get_looks(dist = dist, nnE = nnE, nnF = nnF)
+      nnrE <- nnr$nnrE
+      nnrF <- nnr$nnrF
+    } else {
+      nnrE <- nnE
+      nnrF <- nnF
     }
+    nnr <- unique(c(nnrE, nnrF))
+    tmp <- h_get_decision(
+      nnr = nnr, response = response,
+      truep = truep, p0 = p0, p1 = p1,
+      parE = c(1, 1), nnE = nnrE,
+      nnF = nnrF, tL = tL, tU = tU
+    )
+    decision[k] <- tmp$decision
+    all_sizes[k] <- tmp$all_sizes
   }
-  oc <- cbind(
-    ExpectedN = mean(n), PrStopEarly = mean(n < Nmax),
-    PrEarlyEff = sum(s * (n < Nmax), na.rm = TRUE) / ns,
-    PrEarlyFut = sum((1 - s) * (n < Nmax), na.rm = TRUE) / ns,
-    PrEfficacy = sum(s, na.rm = TRUE) / ns,
-    PrFutility = sum(1 - s, na.rm = TRUE) / ns,
-    PrGrayZone = sum(is.na(s) / ns)
-  )
-  return(list(
-    oc = oc, Decision = s, SampleSize = n,
-    nn = nn, nnE = nnE, nnF = nnF,
+  oc <- h_get_oc(all_sizes = all_sizes, nnr = nnr, decision = decision, nnrE = nnrE, nnrF = nnrF)
+  list(
+    oc = oc,
+    Decision = decision,
+    SampleSize = all_sizes,
+    union_nn = nnr,
+    input_nnE = nnE,
+    input_nnF = nnF,
+    wiggled_Eff_n = nnrE,
+    wiggled_Fut_n = nnrF,
+    wiggle_dist = dist,
     params = as.list(match.call(expand.dots = FALSE))
-  ))
+  )
 }
diff --git a/R/plotOc.R b/R/plotOc.R
@@ -13,7 +13,7 @@
 #' @keywords graphics
 plotOc <- function(z) {
   ## plot function for oc.predprob or oc.postprob, or the dist versions of them
-  graphics::barplot(table(z$Decision, z$SampleSize) / z$params$ns, beside = TRUE)
+  graphics::barplot(table(z$Decision, z$SampleSize) / z$params$sim, beside = TRUE)
 
   ## get the parameter
   parDat <- lapply(z$params, deparse)