R/decision1S_boundary.R

In RBesT: R Bayesian Evidence Synthesis Tools

#' Decision Boundary for 1 Sample Designs
#'
#' Calculates the decision boundary for a 1 sample design. This is the
#' critical value at which the decision function will change from 0
#' (failure) to 1 (success).
#'
#' @template args-boundary1S
#'
#' @details The specification of the 1 sample design (prior, sample
#' size and decision function, \eqn{D(y)}), uniquely defines the
#' decision boundary
#'
#' \deqn{y_c = \max_y\{D(y) = 1\},}{y_c = max_{y}{D(y) = 1},}
#'
#' which is the maximal value of \eqn{y} whenever the decision \eqn{D(y)}
#' function changes its value from 1 to 0 for a decision function
#' with `lower.tail=TRUE` (otherwise the definition is \eqn{y_c =
#' \max_{y}\{D(y) = 0\}}{y_c = max_{y}{D(y) = 0}}). The decision
#' function may change at most at a single critical value as only
#' one-sided decision functions are supported. Here,
#' \eqn{y} is defined for binary and Poisson endpoints as the sufficient
#' statistic \eqn{y = \sum_{i=1}^{n} y_i} and for the normal
#' case as the mean \eqn{\bar{y} = 1/n \sum_{i=1}^n y_i}.
#'
#' The convention for the critical value \eqn{y_c} depends on whether
#' a left (`lower.tail=TRUE`) or right-sided decision function
#' (`lower.tail=FALSE`) is used. For `lower.tail=TRUE` the
#' critical value \eqn{y_c} is the largest value for which the
#' decision is 1, \eqn{D(y \leq y_c) = 1}, while for
#' `lower.tail=FALSE` then \eqn{D(y > y_c) = 1} holds. This is
#' aligned with the cumulative density function definition within R
#' (see for example [pbinom()]).
#'
#' @return Returns the critical value \eqn{y_c}.
#'
#' @family design1S
#'
#' @examples
#'
#' # non-inferiority example using normal approximation of log-hazard
#' # ratio, see ?decision1S for all details
#' s <- 2
#' flat_prior <- mixnorm(c(1, 0, 100), sigma = s)
#' nL <- 233
#' theta_ni <- 0.4
#' theta_a <- 0
#' alpha <- 0.05
#' beta <- 0.2
#' za <- qnorm(1 - alpha)
#' zb <- qnorm(1 - beta)
#' n1 <- round((s * (za + zb) / (theta_ni - theta_a))^2)
#' theta_c <- theta_ni - za * s / sqrt(n1)
#'
#' # double criterion design
#' # statistical significance (like NI design)
#' dec1 <- decision1S(1 - alpha, theta_ni, lower.tail = TRUE)
#' # require mean to be at least as good as theta_c
#' dec2 <- decision1S(0.5, theta_c, lower.tail = TRUE)
#' # combination
#' decComb <- decision1S(c(1 - alpha, 0.5), c(theta_ni, theta_c), lower.tail = TRUE)
#'
#' # critical value of double criterion design
#' decision1S_boundary(flat_prior, nL, decComb)
#'
#' # ... is limited by the statistical significance ...
#' decision1S_boundary(flat_prior, nL, dec1)
#'
#' # ... or the indecision point (whatever is smaller)
#' decision1S_boundary(flat_prior, nL, dec2)
#'
#' @export
decision1S_boundary <- function(prior, n, decision, ...)
  UseMethod("decision1S_boundary")
#' @export
decision1S_boundary.default <- function(prior, n, decision, ...)
  "Unknown density"

#' @templateVar fun decision1S_boundary
#' @template design1S-binomial
#' @export
decision1S_boundary.betaMix <- function(prior, n, decision, ...) {
  VdecisionLazy <- Vectorize(function(r) {
    decision(postmix(prior, r = r, n = n)) - 0.25
  })

  ## find decision boundary
  bounds <- VdecisionLazy(c(0, n))
  lower.tail <- attr(decision, "lower.tail")
  if (prod(bounds) > 0) {
    ## decision is always the same
    if (lower.tail) {
      crit <- ifelse(bounds[1] < 0, 0, n + 1)
    } else {
      crit <- ifelse(bounds[1] < 0, n, -1)
    }
  } else {
    crit <- uniroot_int(
      VdecisionLazy,
      c(0, n),
      f.lower = bounds[1],
      f.upper = bounds[2]
    )
  }

  ## crit is always pointing to the 0 just before the decision which
  ## is why we need a discrimination here
  if (lower.tail) {
    crit <- crit - 1
  }

  crit
}


## returns a function object which is the decision boundary. That is
## the function finds at a regular grid between llim1 and ulim1 the
## roots of the decision function and returns an interpolation
## function object
solve_boundary1S_normMix <- function(decision, mix, n, lim) {
  sigma <- sigma(mix)

  cond_decisionStep <- function() {
    fn <- function(m) {
      decision(postmix(mix, m = m, se = sigma / sqrt(n))) - 0.75
    }
    Vectorize(fn)
  }

  ## ensure that at the limiting boundaries the decision function
  ## has a different sign (which must be true)
  ind_fun <- cond_decisionStep()
  dec_bounds <- ind_fun(lim)
  while (prod(dec_bounds) > 0) {
    w <- diff(lim)
    lim <- c(lim[1] - w / 2, lim[2] + w / 2)
    dec_bounds <- ind_fun(lim)
  }

  uniroot(
    ind_fun,
    interval = lim,
    f.lower = dec_bounds[1],
    f.upper = dec_bounds[2]
  )$root
}

#' @templateVar fun decision1S_boundary
#' @template design1S-normal
#' @export
decision1S_boundary.normMix <- function(
  prior,
  n,
  decision,
  sigma,
  eps = 1e-6,
  ...
) {
  ## distributions of the means of the data generating distributions
  ## for now we assume that the underlying standard deviation
  ## matches the respective reference scales
  if (missing(sigma)) {
    sigma <- RBesT::sigma(prior)
    message("Using default prior reference scale ", sigma)
  }
  assert_number(sigma, lower = 0)

  sd_samp <- sigma / sqrt(n)

  sigma(prior) <- sigma

  ## change the reference scale of the prior such that the prior
  ## represents the distribution of the respective means
  ## mean_prior <- prior
  ## sigma(mean_prior) <- sd_samp

  m <- summary(prior, probs = c())["mean"]

  lim <- qnorm(p = c(eps / 2, 1 - eps / 2), mean = m, sd = sd_samp)

  ## find the boundary of the decision function within the domain we integrate
  crit <- solve_boundary1S_normMix(decision, prior, n, lim)

  crit
}


#' @templateVar fun decision1S_boundary
#' @template design1S-poisson
#' @export
decision1S_boundary.gammaMix <- function(prior, n, decision, eps = 1e-6, ...) {
  assert_that(likelihood(prior) == "poisson")

  cond_decisionStep <- function() {
    fn <- function(m) {
      decision(postmix(prior, n = n, m = m / n)) - 0.25
    }
    Vectorize(fn)
  }

  Vdecision <- cond_decisionStep()

  m <- summary(prior, probs = c())["mean"]

  lambda_prior <- m * n

  lim <- qpois(p = c(eps / 2, 1 - eps / 2), lambda = lambda_prior)
  lim[1] <- 0

  bounds <- Vdecision(lim)
  ## make sure there is a decision somewhere
  while (prod(bounds) > 0) {
    lim[2] <- round(lim[2] * 2)
    bounds[2] <- Vdecision(lim[2])
    if (diff(lim) > 1e9) {
      break
    }
  }

  lower.tail <- attr(decision, "lower.tail")

  ## check if the decision is constantly 1 or 0
  if (prod(bounds) > 0) {
    ## decision is always the same
    if (lower.tail) {
      crit <- ifelse(bounds[1] < 0, 0, Inf)
    } else {
      crit <- ifelse(bounds[1] < 0, Inf, -1)
    }
  } else {
    ## find decision boundary
    crit <- uniroot_int(
      Vdecision,
      c(lim[1], lim[2]),
      f.lower = bounds[1],
      f.upper = bounds[2]
    )
  }

  ## crit is always pointing to the 0 just before the decision which
  ## is why we need a discrimination here
  if (lower.tail) {
    crit <- crit - 1
  }

  crit
}