R/PPC_binom.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 \mjseqn{\ell_{1}} Wasserstein distance of order \mjseqn{q}.
#'
#' @usage PPC_binom(
#'   n,
#'   alpha_1,
#'   beta_1,
#'   alpha_2,
#'   beta_2,
#'   alpha_D,
#'   beta_D,
#'   xi,
#'   v,
#'   q = 1,
#'   M = 3000,
#'   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 q The order of the Wasserstein distance. Must be a value in \mjseqn{[1, \infty)}.
#' @param M The number of Monte Carlo replications. Must be a positive integer.
#' @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. \cr
#'
#' The [PPC_binom_bound()] function should be used when \mjseqn{q = 1}.
#'
#' @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_bound()]
#'
#' @examples
#' \dontrun{
#' # 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)
#'
#' # Seed
#' set.seed(1234)
#'
#' output <- PPC_binom(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 graphics abline legend par
#' @importFrom mathjaxr preview_rd
#' @importFrom pbapply pbmapply

PPC_binom <- function(n, alpha_1, beta_1, alpha_2, beta_2, alpha_D, beta_D, xi, v, q = 1, M = 3000, plot = FALSE) {

  check_PPC_binom(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, q = q, M = M, 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)
  }

  p_n <- pbmapply(FUN = pred_prob_binom, n = n, MoreArgs = list("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, "q" = q, "M" = M))
  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))

}