16: clean postprobDist

R/dbetabinom.R

@@ -19,8 +19,6 @@
 #'  whether to return the log density value (not default).
 #' @return The density values of the beta-binomial distribution at `x`.
 #'
-#' @note `x`, `a` and `b` can be vectors.
-#'
 #' @example examples/dbetabinom.R
 #' @export
 dbetabinom <- function(x, m, a, b, log = FALSE) {
@@ -57,8 +55,6 @@ dbetabinom <- function(x, m, a, b, log = FALSE) {
 #'  whether to return the log density value (not default).
 #' @return The (log) density values of the mixture of beta-binomial distributions at `x`.
 #'
-#' @note `x` can be a vector.
-#'
 #' @example examples/dbetabinomMix.R
 #' @export
 dbetabinomMix <- function(x, m, par, weights, log = FALSE) {
@@ -89,6 +85,7 @@ dbetabinomMix <- Vectorize(dbetabinomMix, vectorize.args = "x")
 #'
 #' @example examples/getBetamixPost.R
 #' @export
+
 getBetamixPost <- function(x, n, par, weights) {
   ## check the format
   stopifnot(
@@ -125,6 +122,7 @@ getBetamixPost <- function(x, n, par, weights) {
 }
 
 
+
 #' Beta-mixture density function
 #'
 #' Note that `x` can be a vector.

R/postprob.R

@@ -35,7 +35,6 @@ postprobBeta <- function(x, n, p, a = 1, b = 1) {
   stats::pbeta(p, a + x, b + n - x, lower.tail = FALSE)
 }
 
-
 #' Posterior Probability of Efficacy Given Beta-Mixture Prior
 #'
 #' @description `r lifecycle::badge("experimental")`

R/postprobDist.R

@@ -2,135 +2,194 @@
 #' @include postprob.R
 NULL
 
-#' Compute the posterior probability with beta prior on SOC
-#'
-#' Using the approach by Thall and Simon (Biometrics, 1994), evaluate the
-#' posterior probability of having Pr(P_E > P_S + delta | data) (but see below
-#' for relative delta margin). Both for the new treatment E as well as for the
-#' SOC S data might be available. However the default is that no data is
-#' available for the SOC, corresponding to the single arm trial situation. Note
-#' that a uniform prior is the useful default for the treatment proportion,
-#' while in the single arm trial an informative prior on the SOC proportion is
-#' useful.
-#'
-#' Beta mixture prior can be specified for the treatment (\code{parE}
-#' and \code{weights} parameters) and control proportion (\code{parS} and
-#' \code{weightsS} parameters), see \code{\link{postprob}} for details. Note
-#' that being able to specify a beta mixture prior also on the control
-#' treatment is e.g. important for the futility decision making (see the
-#' \code{\link{oc2}} code).
-#'
-#' @param x number of successes (in the treatment group). Note that \code{x}
-#' can be a vector.
-#' @param n number of patients (in the treatment group)
-#' @param xS number of successes in the SOC group (default: 0)
-#' @param nS number of patients in the SOC group (default: 0)
-#' @param delta margin by which the response rate in the treatment group should
-#' be better than in the SOC group (default: 0)
-#' @param relativeDelta should the delta be relative? (not default). If this is
-#' \code{TRUE}, then a relative delta is used. This means we want to have
-#' response at least in delta proportion of the SOC non-responding patients.
-#' Non-responding patients rate is 1 - P_S, and we want to have P_S + (1 - P_S)
-#' * delta response rate (at least) in the treatment. That is, we evaluate the
-#' posterior probability Pr(P_E > P_S + (1 - P_S) * delta | data).
-#' @param parE the beta parameters matrix, with K rows and 2 columns,
-#' corresponding to the beta parameters of the K components. default is a
-#' uniform prior.
-#' @param weights the mixture weights of the beta mixture prior. Default are
-#' uniform weights across mixture components.
-#' @param parS beta parameters for the SOC group (default: uniform)
-#' @param weightsS weights for the SOC group (default: uniform)
-#' @return the posterior probability
+#' The Posterior Beta Mixture Integrand when `delta` is relative
+#'
+#' The helper function to generate Integrand function when `relative Delta = TRUE`.
+#'
+#' @typed delta : number
+#'  margin by which the response rate in the treatment group should
+#'  be better than in the SOC group. Must be >= `0`. See note.
+#' @typed p_s : number
+#'  probability of success or response rate of standard of care or `SOC` group.
+#' @typed activeBetamixPost : list
+#'  a list of posterior parameters of a beta-mixture-binomial distribution with generic names
+#'  `par` and `weights`. See `[getBetaMix()]`.
+#' @typed controlBetamixPost : list
+#'  a list of posterior parameters of a beta-mixture-binomial distribution with generic names
+#'  `par` and `weights`. See `[getBetaMix()]`.
+#'
+#' @return Function that is an argument for `[stats::integrate()]`.
+#'
+#' @keywords internal
+#'
+h_integrand_relDelta <- function(p_s, delta, activeBetamixPost, controlBetamixPost) {
+  cdf <- postprob(
+    x = 0, # Needed for Vectorize()
+    p = (1 - p_s) * delta + p_s,
+    betamixPost = activeBetamixPost
+  )
+  pdf <- with(
+    controlBetamixPost,
+    dbetaMix(x = p_s, par = par, weights = weights)
+  )
+  cdf * pdf
+}
+
+#' The Posterior Beta Mixture Integrand when Delta is absolute
+#'
+#' The helper function to generate Integrand function when `relative Delta = FALSE`,
+#' a default setting.
+#'
+#' @inheritParams h_integrand_relDelta
+#'
+#' @return Function that is an argument for `[stats::integrate()]`.
+#'
+#' @keywords internal
+#'
+h_integrand <- function(p_s, delta, activeBetamixPost, controlBetamixPost) {
+  cdf <- postprob(
+    x = 0, # Needed for Vectorize()
+    p = p_s + delta,
+    betamixPost = activeBetamixPost
+  )
+  pdf <- with(
+    controlBetamixPost,
+    dbetaMix(x = p_s, par = par, weights = weights)
+  )
+  cdf * pdf
+}
+
+#' Generating bounds for the Integration of Beta Mixture Posterior
+#'
+#' Using the quantile of the Beta Mixture Distribution from parameters given by standard of care `SOC` or
+#' experimental group `E` to determine bounds as inputs to `[stats::integrate()]`.
+#'
+#' @inheritParams h_integrand_relDelta
+#' @return Integrand function
+#'
+#' @keywords internal
+#'
+h_get_bounds <- function(controlBetamixPost) {
+  epsilon <- .Machine$double.xmin
+  with(
+    controlBetamixPost,
+    qbetaMix(
+      p = c(epsilon, 1 - epsilon),
+      par = par,
+      weights = weights
+    )
+  )
+}
+
+#' Compute the Posterior Probability with Beta Prior on `SOC`
+#'
+#' @description `r lifecycle::badge("experimental")`
+#'
+#' Using the approach by Thall and Simon (Biometrics, 1994), this evaluates the
+#' posterior probability of achieving superior response rate in the treatment group compared to standard of care (SOC).
+#' See note below for two formulations of the difference in response rates.
+#'
+#' @inheritParams h_integrand_relDelta
+#' @typed x : numeric
+#'  number of success counts in the treatment group.
+#' @typed n : number
+#'  number of patients in the treatment group.
+#' @typed xS : number
+#'  number of success counts in the SOC group.
+#' @typed nS : number
+#'  number of patients in the SOC group.
+#' @typed relativeDelta : flag
+#'  If `TRUE`, then a `relativeDelta` is used. Represents that a minimum
+#'  response rate in magnitude of `delta` of the SOC non-responding patients. See note.
+#' @typed parE : "`numeric` or `matrix`"
+#'  parameters for beta distribution. If it is a matrix, it needs to have 2 columns,
+#'  and each row corresponds to each component of a beta-mixture distribution
+#'  for the `E` group. See details.
+#' @typed weights : numeric
+#'  the non-negative mixture weights of the beta mixture prior for group `E`. See details.
+#' @typed parS : "`numeric` or `matrix`"
+#'  parameters for beta distribution. If it is a matrix, it needs to have 2 columns,
+#'  and each row corresponds to each component of a beta-mixture distribution
+#'  for the `S` group. See details.
+#' @typed weightsS : numeric
+#'  the non-negative mixture weights of the beta mixture prior for group `S`. See details.
+#' @typed epsilon : number
+#'  the smallest non-negative floating number to represent the lower bound for
+#'  the interval of integration.
+#' @return The posterior probability
+#'
+#' @note
+#'
+#' ## Delta :
+#'
+#' The desired improvement is denoted as `delta`. There are two options in using `delta`.
+#' The absolute case when `relativeDelta = FALSE` and relative as when `relativeDelta = TRUE`.
+#'
+#' 1. The absolute case is when we define an absolute delta, greater than `P_S`,
+#' the response rate of the `SOC` group such that
+#' the posterior is `Pr(P_E > P_S + delta | data)`.
+#'
+#' 2. In the relative case, we suppose that the treatment group's
+#' response rate is assumed to be greater than `P_S + (1-P_S) * delta` such that
+#' the posterior is `Pr(P_E > P_S + (1 - P_S) * delta | data)`.
+#'
+#' @details
+#'
+#' The beta mixture prior for the `E` arm requires argument `parE` and `weights`.
+#' The beta mixture prior for the `S` arm requires argument `parS` and `weightsS`.
+#' See `[postprob()]` for details.
+#'
+#' If a beta-mixture is used, by default, the weights are uniform across the components.
+#' Weights can exceed 1, to which the algorithm will normalize the weights such that all weights sum to 1.
 #'
 #' @example examples/postprobDist.R
 #' @export
-postprobDist <- function(x, n,
-                         xS = 0, nS = 0,
-                         delta = 0,
+postprobDist <- function(x,
+                         n,
+                         xS = 0,
+                         nS = 0,
+                         delta,
                          relativeDelta = FALSE,
                          parE = c(1, 1),
                          weights,
                          parS = c(1, 1),
                          weightsS) {
-  ## if parE is a vector => situation where there is only one component
   if (is.vector(parE)) {
-    ## check that it has exactly two entries
-    stopifnot(identical(length(parE), 2L))
-
-    ## and transpose to matrix with one row
+    assert_true(identical(length(parE), 2L))
     parE <- t(parE)
   }
-
-  ## if prior weights of the beta mixture are not supplied
-  if (missing(weights)) {
-    weights <- rep(1, nrow(parE))
-  }
-
-  ## if parS is a vector => situation where there is only one component
   if (is.vector(parS)) {
-    ## check that it has exactly two entries
-    stopifnot(identical(length(parS), 2L))
-
-    ## and transpose to matrix with one row
+    assert_true(identical(length(parS), 2L))
     parS <- t(parS)
   }
-
-  ## if prior weights of the beta mixture are not supplied
+  if (missing(weights)) {
+    weights <- rep(1, nrow(parE))
+  }
   if (missing(weightsS)) {
     weightsS <- rep(1, nrow(parS))
   }
-
-  ## compute updated beta parameters
+  assert_numeric(x, lower = 0, upper = n, finite = TRUE)
+  assert_number(n, lower = x, finite = TRUE)
+  assert_number(xS, lower = 0, upper = n, finite = TRUE)
+  assert_numeric(nS, lower = 0, finite = TRUE)
+  assert_number(delta, lower = 0, finite = TRUE)
+  assert_flag(relativeDelta)
+  assert_numeric(weights, lower = 0, finite = TRUE)
+  assert_numeric(weightsS, lower = 0, finite = TRUE)
+  assert_numeric(parE, lower = 0, finite = TRUE)
+  assert_numeric(parS, lower = 0, finite = TRUE)
   activeBetamixPost <- getBetamixPost(x = x, n = n, par = parE, weights = weights)
   controlBetamixPost <- getBetamixPost(x = xS, n = nS, par = parS, weights = weightsS)
-
-  ## use numerical integration to compute this probability, as given on p.338
-  ## in the article by Thall and Simon (1994):
-  integrand <-
-    if (relativeDelta) {
-      function(p) {
-        cdf <- postprob(
-          x = x,
-          p = (1 - delta) * p + delta,
-          betamixPost = activeBetamixPost
-        )
-
-        pdf <- with(
-          controlBetamixPost,
-          dbetaMix(x = p, par = par, weights = weights)
-        )
-
-        return(cdf * pdf)
-      }
-    } else {
-      function(p) {
-        cdf <- postprob(
-          x = x,
-          p = p + delta,
-          betamixPost = activeBetamixPost
-        )
-
-        pdf <- with(
-          controlBetamixPost,
-          dbetaMix(x = p, par = par, weights = weights)
-        )
-
-        return(cdf * pdf)
-      }
-    }
-
-  ## do the integration. be careful to cover the region where there can
-  ## really be any non-zero values. I.e. only integrate over the region where
-  ## the beta density of the control is non-zero.
-  epsilon <- 1e-13
-  bounds <- with(
-    controlBetamixPost,
-    qbetaMix(
-      p = c(epsilon, 1 - epsilon),
-      par = par,
-      weights = weights
-    )
-  )
+  assert_names(names(activeBetamixPost), identical.to = c("par", "weights"))
+  assert_names(names(controlBetamixPost), identical.to = c("par", "weights"))
+  if (relativeDelta) {
+    epsilon <- .Machine$double.xmin
+    integrand <- h_integrand_relDelta
+  } else {
+    epsilon <- .Machine$double.xmin
+    integrand <- h_integrand
+  }
+  bounds <- h_get_bounds(controlBetamixPost = controlBetamixPost)
   intRes <- integrate(
     f = integrand,
     lower =
@@ -142,11 +201,13 @@ postprobDist <- function(x, n,
       min(
         ifelse(relativeDelta, 1, 1 - delta),
         bounds[2]
-      )
+      ),
+    delta = delta,
+    activeBetamixPost = activeBetamixPost,
+    controlBetamixPost = controlBetamixPost
   )
-
   if (intRes$message == "OK") {
-    return(intRes$value)
+    intRes$value
   } else {
     stop(intRes$message)
   }