R/ss_best_binomial.R

In ssutil: Sample Size Calculation Tools

#' Power to Correctly Select the Best Group in a Binomial Test
#'
#' Computes the exact probability of correctly identifying the best group
#' when the outcome follows a binomial distribution. It assumes that \code{p1}
#' is the probability of success in the best group, and that the success
#' probability in all other groups is lower by a fixed difference \code{dif}.
#'
#' The formula is based on the exact method described by Sobel and Huyett (1957).
#'
#' @param p1 Numeric. Probability of success in the best group (must be in \[0, 1\]).
#' @param dif Numeric. Difference in success probability between the best group and the next best (must be > 0).
#' @param ngroups Integer. Number of groups (must be greater than 1).
#' @param npergroup Integer. Number of subjects per group (must be positive).
#'
#' @return A numeric value representing the probability of correctly identifying the best group.
#'
#' @importFrom stats pbinom dbinom
#' @export
#'
#' @examples
#' power_best_binomial(p1 = 0.8, dif = 0.2, ngroups = 4, npergroup = 50)
#'
#' @references
#' Sobel, M., & Huyett, M. J. (1957). Selecting the Best One of Several Binomial Populations.
#' *Bell System Technical Journal*, 36(2), 537–576. \doi{10.1002/j.1538-7305.1957.tb02411.x}
power_best_binomial <- function(p1, dif, ngroups, npergroup) {
  stopifnot("p1 must be between 0 and 1" = p1 >= 0 & p1 <= 1)
  stopifnot("p1 must be greater than dif" = p1 >= dif)
  stopifnot("p1 - dif must be > 0" = (p1 - dif) > 0)
  stopifnot("ngroups must be greater than 1" = ngroups > 1)
  stopifnot("ngroups must be an integer" = ngroups %% 1 == 0)
  stopifnot("npergroup must be > 0" = npergroup > 0)
  stopifnot("npergroup must be an integer" = npergroup %% 1 == 0)

  d <- dif
  k <- ngroups
  n <- npergroup

  b1j <- function(j, n, p1) dbinom(j, n, p1)
  b2j <- function(j, n, p1, d) dbinom(j, n, p1 - d)
  B2j <- function(j, n, p1, d) pbinom(j, n, p1 - d)

  Pcs_sum <- 0
  for (j in 0:n) {
    b1j_val <- b1j(j, n, p1)
    inner_sum <- 0
    for (i in 0:(k - 1)) {
      coeff <- choose(k - 1, i) / (1 + i)
      b2j_val <- b2j(j, n, p1, d)^i
      # Bug correction
      #B2j_val <- B2j(j - 1 - i , n, p1, d)^(k - 1 - i)
      B2j_val <- B2j(j - 1, n, p1, d)^(k - 1 - i)
      inner_sum <- inner_sum + coeff * b2j_val * B2j_val
    }
    Pcs_sum <- Pcs_sum + b1j_val * inner_sum
  }

  return(Pcs_sum)
}

#' Sample Size to Select the Best Group in a Binomial Test
#'
#' Computes the minimum sample size per group required to achieve a target probability
#' of correctly selecting the best group in a binomial test. The best group is assumed
#' to have success probability \code{p1}, and the other groups have \code{p1 - dif}.
#'
#' The function searches for the smallest \code{npergroup} such that the power from
#' \code{\link{power_best_binomial}} is at least the target \code{power}.
#'
#' @param power Numeric. Desired probability of correctly selecting the best group (in \[0, 1\]).
#' @param p1 Numeric. Probability of success in the best group (in \[0, 1\]).
#' @param dif Numeric. Difference in success probability with the next best group (> 0).
#' @param ngroups Integer. Number of groups (must be > 1).
#' @param max_n Integer. Maximum sample size to evaluate (default is 1000).
#'
#' @return An integer representing the minimum sample size per group required to reach the specified power.
#'
#' @export
#'
#' @examples
#' ss_best_binomial(power = 0.9, p1 = 0.8, dif = 0.2, ngroups = 4)
ss_best_binomial <- function(power, p1, dif, ngroups, max_n = 1000) {
  stopifnot("power must be between 0 and 1" = power >= 0 & power <= 1)
  stopifnot("p1 must be between 0 and 1" = p1 >= 0 & p1 <= 1)
  stopifnot("p1 must be greater than dif" = p1 >= dif)
  stopifnot("p1 - dif must be > 0" = (p1 - dif) > 0)

  n <- 1
  while (n <= max_n) {
    ps <- power_best_binomial(p1, dif, ngroups, n)
    if (ps >= power) return(n)
    n <- n + 1
  }

  stop(max_n, ": max_n limit reached without achieving desired power.")
}


#' Worst‐Case Scenario Power for the Best Binomial Group
#'
#' Searches for the probability in the best‐performing group that yields the lowest statistical power,
#' given an indifference zone specification, a number of groups, and a number of subjects per group.
#'
#' @param dif Numeric. Indifference zone specification (difference threshold).
#' @param ngroups Integer. Number of groups to compare.
#' @param npergroup Integer. Number of subjects per group.
#' @return A named list with components:
#' \describe{
#'   \item{p1}{Numeric. Probability in the best group that yields the minimum power.}
#'   \item{minimum_power}{Numeric. The minimum power achieved at \code{p1}.}
#' }
#' @details
#' Defines an internal function \code{fx} that wraps \code{\link{power_best_binomial}}
#' with the supplied parameters, then uses \code{\link[stats]{optimize}} over the interval \[0,1\]
#' to find the probability \code{p1} that minimizes the resulting power.
#'
#' @seealso
#' \code{\link{power_best_binomial}}, \code{\link[stats]{optimize}}
#'
#' @examples
#' wcs_power_best_binomial(dif = 0.1, ngroups = 3, npergroup = 50)
#'
#' @export
wcs_power_best_binomial <- function(dif, ngroups, npergroup) {
  fx <- function(x) {
    power_best_binomial(x, dif, ngroups, npergroup)
  }
  # correct bug as searching from 0 to 1 may probabilities lower than dif
  res <- stats::optimize(fx, interval = c(dif + 1e-6, 1))
  names(res) <- c("p1", "minimum_power")
  res
}



#' Monte Carlo power for indifference-zone "best" binomial selection
#'
#' Group 1 is the true best (p1 + d); groups 2:ngroups have probability p1.
#' Power = P(group 1 has the highest observed count), ties broken at random.
#'
#' @param p1 probability of the best group
#' @param d difference with the rest of the groups (probability = p1-d)
#' @param ngroups number of groups (k)
#' @param npergroup trials per group (n)
#' @param nsim number of Monte Carlo simulations
#' @param seed optional seed for reproducibility
#' @return list(power, se, ci_95, nsim)
#' @export
# sim_power_best_binomial <- function(p1, d, ngroups, npergroup, nsim, seed = NULL) {
#   
#   if (p1 < 0 || p1 > 1) stop("p1 must be in [0, 1]")
#   if ((p1 - d) < 0 || (p1 - d) > 1) stop("p1 - d must be in [0, 1]")
#   if (ngroups < 2) stop("ngroups must be >= 2")
#   
#   if (!is.null(seed)) set.seed(seed)
#   
#   # Group 1 is the TRUE BEST (probability p1)
#   # Groups 2:ngroups are inferior (probability p1 - d)
#   probs <- c(p1, rep(p1 - d, ngroups - 1))
#   
#   counts <- matrix(
#     rbinom(nsim * ngroups, size = npergroup, prob = rep(probs, each = nsim)),
#     nrow = nsim, ncol = ngroups
#   )
#   
#   row_max <- apply(counts, 1, max)
#   is_max  <- counts == row_max
#   
#   selected <- apply(is_max, 1, function(row) {
#     idx <- which(row)
#     if (length(idx) == 1L) idx else sample(idx, 1L)
#   })
#   
#   successes <- selected == 1L
#   power_hat <- mean(successes)
#   
#   se <- sqrt(power_hat * (1 - power_hat) / nsim)
#   ci <- c(
#     lower = max(0, power_hat - 1.96 * se),
#     upper = min(1, power_hat + 1.96 * se)
#   )
#   
#   list(power = power_hat, se = se, ci_95 = ci, nsim = nsim)
# }