R/PPC_norm.R

In michelecianfriglia/SampleSizeWass: Wasserstein-based Sample Size Calculations

#' Predictive Probability Criterion: Normal Model
#'
#' @description \loadmathjax{} Implements the predictive probability criterion in the normal model using the \mjseqn{\ell_{2}} Wasserstein distance of order two.
#'
#' @usage PPC_norm(
#'   n,
#'   theta_1,
#'   n_1,
#'   theta_2,
#'   n_2,
#'   theta_D,
#'   n_D,
#'   sigma,
#'   xi,
#'   v,
#'   plot = FALSE
#' )
#'
#' @param n The sample size. Must be a vector of positive integers arranged in ascending order.
#' @param theta_1,n_1 The parameters of the first normal prior. Must be a finite value and a non-negative value, respectively.
#' @param theta_2,n_2 The parameters of the second normal prior. Must be a finite value and a non-negative value, respectively.
#' @param theta_D,n_D The parameters of the design normal prior. Must be a finite value and a positive value, respectively.
#' @param sigma A constant used to define the variance of the priors. Must be a positive value.
#' @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 normal priors, whereas the design normal prior must be proper. \cr
#'
#' If the first and second normal priors are equal, the function stops with an error. \cr
#'
#' The prior variances are given by \code{sigma} squared over the prior sample sizes, that is \code{n_1}, \code{n_2}, and \code{n_D}.
#'
#' @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.
#'
#' @examples
#' # Parameters of the first normal prior
#' prior_1 <- c(15, 25)
#'
#' # Parameters of the second normal prior
#' prior_2 <- c(10, 15)
#'
#' # Parameters of the design normal prior
#' prior_D <- c(13, 10)
#'
#' output <- PPC_norm(n = 1:800,
#'                    theta_1 = prior_1[1], n_1 = prior_1[2],
#'                    theta_2 = prior_2[1], n_2 = prior_2[2],
#'                    theta_D = prior_D[1], n_D = prior_D[2],
#'                    sigma = 2, xi = 0.05, v = 0.1)
#' @export
#'
#' @importFrom crayon italic
#' @importFrom graphics abline legend par
#' @importFrom mathjaxr preview_rd
#' @importFrom stats pchisq

PPC_norm <- function(n, theta_1, n_1, theta_2, n_2, theta_D, n_D, sigma, xi, v, plot = FALSE) {

  check_PPC_norm(n = n, theta_1 = theta_1, n_1 = n_1, theta_2 = theta_2, n_2 = n_2, theta_D = theta_D, n_D = n_D, sigma = sigma, xi = xi, v = v, plot = plot)

  if (theta_1 == theta_2 & n_1 == n_2) {
    stop(paste0("The predictive probability is zero for the selected parameters. Please choose (", italic("theta_1"), ", ", italic("n_1"), ") different from (", italic("theta_2"), ", ", italic("n_2"), ")."), call. = FALSE)
  }

  if (n_1 == n_2) {
    stop(paste0("The predictive probability is zero for the selected parameters. Please choose ", italic("n_1"), " different from ", italic("n_2"), "."), call. = FALSE)
  }

  w_1 <- n_1 / (n + n_1)
  w_2 <- n_2 / (n + n_2)
  w_n <- w_2 - w_1

  bures_dist <- (sigma / sqrt(n + n_1) - sigma / sqrt(n + n_2))^2
  denom <- (w_n * sigma)^2 * (1 / n + 1 / n_D)

  p_n <- 1 - pchisq(q = (xi - bures_dist) / denom, df = 1, ncp = (w_1 * theta_1 - w_2 * theta_2 + w_n * theta_D)^2 / denom)
  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))

}