R/PPC_binom_bound.R

In michelecianfriglia/SampleSizeWass: Wasserstein-based Sample Size Calculations

#' Predictive Probability Criterion: Binomial Model
#'
#' @description \loadmathjax{} Implements the predictive probability criterion in the binomial model using the upper bound of the \mjseqn{\ell_{1}} Wasserstein distance of order one.
#'
#' @usage PPC_binom_bound(
#'   n,
#'   alpha_1,
#'   beta_1,
#'   alpha_2,
#'   beta_2,
#'   alpha_D,
#'   beta_D,
#'   xi,
#'   v,
#'   plot = FALSE
#' )
#'
#' @param n The sample size. Must be a vector of positive integers arranged in ascending order.
#' @param alpha_1,beta_1 The parameters of the first beta prior. Must be non-negative values.
#' @param alpha_2,beta_2 The parameters of the second beta prior. Must be non-negative values.
#' @param alpha_D,beta_D The parameters of the design beta prior. Must be positive values.
#' @param xi A constant used to compute the predictive probability. Must be a positive value.
#' @param v A constant used to determine the optimal sample size. Must be a value in \mjseqn{(0, 1)}.
#' @param plot Logical. If \code{TRUE}, a plot shows the behavior of the predictive probability as a function of the sample size.
#'
#' @details Users can use non-informative improper priors for the first and second beta priors, whereas the design beta prior must be proper. \cr
#'
#' If the first and second beta priors are equal, the function stops with an error.
#'
#' @return A list with the following components: \cr
#'
#' \item{p_n}{The predictive probability.}
#' \item{t_opt}{The optimal threshold.}
#' \item{n_opt}{The optimal sample size.}
#'
#' @author Michele Cianfriglia \email{michele.cianfriglia@@uniroma1.it}
#'
#' @references Cianfriglia, M., Padellini, T., and Brutti, P. (2023). Wasserstein consensus for Bayesian sample size determination.
#'
#' @seealso [PPC_binom()]
#'
#' @examples
#' # Parameters of the first beta prior
#' prior_1 <- c(51, 42)
#'
#' # Parameters of the second beta prior
#' prior_2 <- c(55, 29)
#'
#' # Parameters of the design beta prior
#' prior_D <- c(23, 15)
#'
#' output <- PPC_binom_bound(n = 1:200,
#'                           alpha_1 = prior_1[1], beta_1 = prior_1[2],
#'                           alpha_2 = prior_2[1], beta_2 = prior_2[2],
#'                           alpha_D = prior_D[1], beta_D = prior_D[2],
#'                           xi = 0.05, v = 0.1)
#' @export
#'
#' @importFrom crayon italic
#' @importFrom extraDistr pbbinom
#' @importFrom graphics abline legend par
#' @importFrom mathjaxr preview_rd

PPC_binom_bound <- function(n, alpha_1, beta_1, alpha_2, beta_2, alpha_D, beta_D, xi, v, plot = FALSE) {

  check_PPC_binom_bound(n = n, alpha_1 = alpha_1, beta_1 = beta_1, alpha_2 = alpha_2, beta_2 = beta_2, alpha_D = alpha_D, beta_D = beta_D, xi = xi, v = v, plot = plot)

  if (alpha_1 == alpha_2 & beta_1 == beta_2) {
    stop(paste0("The predictive probability is zero for the selected parameters. Please choose (", italic("alpha_1"), ", ", italic("beta_1"), ") different from (", italic("alpha_2"), ", ", italic("beta_2"), ")."), call. = FALSE)
  }

  if (abs(alpha_1 - alpha_2) == abs(beta_1 - beta_2)) {
    stop(paste0("The predictive probability is zero for the selected parameters. Please choose (", italic("alpha_1"), ", ", italic("beta_1"), ") different from (", italic("alpha_2"), ", ", italic("beta_2"), ")."), call. = FALSE)
  }

  if (abs(alpha_1 - alpha_2) < abs(beta_1 - beta_2)) {

    k_1 <- (xi * (n + alpha_1 + beta_1) * (n + alpha_2 + beta_2) - abs(alpha_1 - alpha_2) * (n + beta_2) - alpha_2 * abs(beta_1 - beta_2)) / (abs(beta_1 - beta_2) - abs(alpha_1 - alpha_2))
    p_n <- 1 - pbbinom(q = floor(k_1), size = n, alpha = alpha_D, beta = beta_D)
    n_max <- which.max(p_n)
    t_opt <- v * p_n[n_max]
    n_opt <- which.max(p_n[n_max:max(n)] < t_opt) + n_max - 1

    if (n_opt == 1) {
      if (all(p_n == 0) == TRUE) {
        warning(paste("The predictive probability is zero for the selected parameters. The value of", italic("xi"), "is probably too large."), call. = FALSE)
      } else {
        warning("The optimal sample size is equal to one. The selected sample size is probably too small.", call. = FALSE)
      }
    }

    if (plot) {
      par(pty = "s")
      plot(x = n, y = p_n, type = "l", lwd = 2, xlab = "n", ylab = "predictive probability")
      abline(h = t_opt, lwd = 2, col = "red")
      abline(v = n_opt, lwd = 2, col = "blue")
      legend("top", c(paste("optimal threshold =", formatC(x = t_opt, format = "e", digits = 2)), paste("optimal sample size =", n_opt)), lty = c(1, 1), lwd = c(2, 2), col = c("red", "blue"), bty = "n", inset = c(0, -0.25), xpd = TRUE)
      par(pty = "m")
    }

    return(list("p_n" = p_n, "t_opt" = t_opt, "n_opt" = n_opt))

  } else {

    k_2 <- (abs(alpha_1 - alpha_2) * (n + beta_2) + alpha_2 * abs(beta_1 - beta_2) - xi * (n + alpha_1 + beta_1) * (n + alpha_2 + beta_2)) / (abs(alpha_1 - alpha_2) - abs(beta_1 - beta_2))
    p_n <- pbbinom(q = floor(k_2) - 1, size = n, alpha = alpha_D, beta = beta_D)
    n_max <- which.max(p_n)
    t_opt <- v * p_n[n_max]
    n_opt <- which.max(p_n[n_max:max(n)] < t_opt) + n_max - 1

    if (n_opt == 1) {
      if (all(p_n == 0) == TRUE) {
        warning(paste("The predictive probability is zero for the selected parameters. The value of", italic("xi"), "is probably too large."), call. = FALSE)
      } else {
        warning("The optimal sample size is equal to one. The selected sample size is probably too small.", call. = FALSE)
      }
    }

    if (plot) {
      par(pty = "s")
      plot(x = n, y = p_n, type = "l", lwd = 2, xlab = "n", ylab = "predictive probability")
      abline(h = t_opt, lwd = 2, col = "red")
      abline(v = n_opt, lwd = 2, col = "blue")
      legend("top", c(paste("optimal threshold =", formatC(x = t_opt, format = "e", digits = 2)), paste("optimal sample size =", n_opt)), lty = c(1, 1), lwd = c(2, 2), col = c("red", "blue"), bty = "n", inset = c(0, -0.25), xpd = TRUE)
      par(pty = "m")
    }

    return(list("p_n" = p_n, "t_opt" = t_opt, "n_opt" = n_opt))

  }

}