add exact derivations for power and type-i error in single-arm studies

Author: NightlordTW

NAMESPACE

@@ -40,7 +40,9 @@ export(power_fisher)
 export(run_bcts_app)
 export(sim_rct_normal)
 export(singlearm_beta_power)
+export(singlearm_beta_power_exact)
 export(singlearm_beta_type1)
+export(singlearm_beta_type1_exact)
 export(ssre)
 import(ggplot2)
 importFrom(Rcpp,evalCpp)

R/RcppExports.R

@@ -5,6 +5,49 @@ beta_prob_gt <- function(a, b, M, n_draws) {
     .Call(`_bcts_beta_prob_gt`, a, b, M, n_draws)
 }
 
+#' @title Analytic Power Calculation for Single-Arm Beta-Binomial Design
+#'
+#' @description Computes the exact Bayesian power for a single-arm binomial trial with a
+#' conjugate Beta prior. Power is defined as the probability that the posterior probability
+#' that the true response rate exceeds a threshold \code{M} is greater than or equal to
+#' a prespecified cutoff \code{threshold}, under a fixed true response rate \code{p_t}.
+#'
+#' This function avoids simulation and uses a deterministic summation over all possible
+#' outcomes of the binomial distribution.
+#'
+#' @param p_t Numeric in \[0, 1\]. True response probability for the treatment arm.
+#' @param n_t Integer. Sample size of the treatment arm.
+#' @param M Numeric in \[0, 1\]. Threshold on the response rate for decision-making,
+#' e.g., \code{M = 0.6}.
+#' @param threshold Numeric in \[0, 1\]. Posterior probability cutoff for declaring success,
+#' e.g., \code{0.95}.
+#' @param prior Character string. Either \code{"flat"} for a Beta(1,1) prior or
+#' \code{"beta"} to specify a custom prior using \code{a_base} and \code{b_base}.
+#' @param a_base Numeric. Alpha parameter for the Beta prior (used only if \code{prior = "beta"}).
+#' @param b_base Numeric. Beta parameter for the Beta prior (used only if \code{prior = "beta"}).
+#'
+#' @return A named list with:
+#' \describe{
+#'   \item{\code{estimate}}{The exact Bayesian power (a number between 0 and 1).}
+#'   \item{\code{mc_se}}{\code{NA_real_}, returned for compatibility (Monte Carlo SE not applicable).}
+#'   \item{\code{B}}{\code{NA_integer_}, returned for compatibility with simulation version.}
+#'   \item{\code{successes}}{\code{NA_integer_}, returned for compatibility.}
+#' }
+#'
+#' @examples
+#' singlearm_beta_power_exact(
+#'   p_t = 0.75, n_t = 35, M = 0.6,
+#'   threshold = 0.95, prior = "flat"
+#' )
+#'
+#' @seealso \code{\link{singlearm_beta_power}} for the simulation-based version.
+#'
+#' @author Thomas Debray \email{[email protected]}
+#' @export
+singlearm_beta_power_exact <- function(p_t, n_t, M, threshold, prior = "flat", a_base = 1, b_base = 1) {
+    .Call(`_bcts_singlearm_beta_power_exact`, p_t, n_t, M, threshold, prior, a_base, b_base)
+}
+
 #' @title Compute p-values for a t-distribution with Fixed Degrees of Freedom
 #'
 #' @description Simulates a single-arm binomial trial with a conjugate Beta prior,
@@ -22,21 +65,20 @@ beta_prob_gt <- function(a, b, M, n_draws) {
 #' \code{"beta"} to specify a custom prior using \code{a_base} and \code{b_base}.
 #' @param a_base Numeric. Alpha parameter for the Beta prior (only used if \code{prior = "beta"}).
 #' @param b_base Numeric. Beta parameter for the Beta prior (only used if \code{prior = "beta"}).
-#' @param n_draws Integer. Number of posterior draws per trial.
 #' @param show_progress Logical. If \code{TRUE}, prints a simple progress bar to console.
 #'
 #' @return A list with: estimate (power), mc_se, successes, B.
 #'
 #' @examples
 #' singlearm_beta_power(
 #'   B = 1000, p_t = 0.75, n_t = 35, M = 0.60,
-#'   threshold = 0.95, prior = "flat", n_draws = 2000
+#'   threshold = 0.95, prior = "flat"
 #' )
 #'
 #' @author Thomas Debray \email{[email protected]}
 #' @export
-singlearm_beta_power <- function(B, p_t, n_t, M, threshold, prior = "flat", a_base = 1, b_base = 1, n_draws = 2000L, show_progress = TRUE) {
-    .Call(`_bcts_singlearm_beta_power`, B, p_t, n_t, M, threshold, prior, a_base, b_base, n_draws, show_progress)
+singlearm_beta_power <- function(B, p_t, n_t, M, threshold, prior = "flat", a_base = 1, b_base = 1, show_progress = TRUE) {
+    .Call(`_bcts_singlearm_beta_power`, B, p_t, n_t, M, threshold, prior, a_base, b_base, show_progress)
 }
 
 #' @title Estimate Type-I Error for Single-Arm Trial
@@ -52,20 +94,64 @@ singlearm_beta_power <- function(B, p_t, n_t, M, threshold, prior = "flat", a_ba
 #' @param prior "flat" or "beta".
 #' @param a_base Alpha parameter for Beta prior (if prior = "beta").
 #' @param b_base Beta parameter for Beta prior (if prior = "beta").
-#' @param n_draws Number of posterior draws per trial.
 #' @param show_progress Logical. Show progress in console?
 #'
 #' @return A list with \code{estimate} (type-I error), \code{mc_se}, \code{B}, and \code{rejections}.
 #'
 #' @examples
 #' singlearm_beta_type1(
 #'   B = 1000, n_t = 35, M = 0.6,
-#'   threshold = 0.95, prior = "flat", n_draws = 2000
+#'   threshold = 0.95, prior = "flat"
 #' )
 #'
 #' @author Thomas Debray \email{[email protected]}
 #' @export
-singlearm_beta_type1 <- function(B, n_t, M, threshold, prior = "flat", a_base = 1, b_base = 1, n_draws = 2000L, show_progress = TRUE) {
-    .Call(`_bcts_singlearm_beta_type1`, B, n_t, M, threshold, prior, a_base, b_base, n_draws, show_progress)
+singlearm_beta_type1 <- function(B, n_t, M, threshold, prior = "flat", a_base = 1, b_base = 1, show_progress = TRUE) {
+    .Call(`_bcts_singlearm_beta_type1`, B, n_t, M, threshold, prior, a_base, b_base, show_progress)
+}
+
+#' @title Exact Type-I Error for Single-Arm Trial (Beta-Binomial)
+#'
+#' @description Computes the exact Type-I error for a single-arm binomial trial
+#' using a conjugate Beta prior. This is done by summing the probability of all
+#' outcomes where the posterior probability that \eqn{\theta > M} exceeds a
+#' decision threshold \code{threshold}, under the null hypothesis.
+#'
+#' By default, the null hypothesis is \eqn{\theta = M}, but a more conservative
+#' frequentist setting can be evaluated by setting \code{p_null < M}.
+#'
+#' @param n_t Integer. Sample size of the treatment arm.
+#' @param M Numeric in \[0, 1\]. Decision threshold on the response rate, e.g., \code{M = 0.6}.
+#' @param threshold Numeric in \[0, 1\]. Posterior probability cutoff for declaring success,
+#' e.g., \code{threshold = 0.95}.
+#' @param prior Character string. Either \code{"flat"} for a Beta(1,1) prior or
+#' \code{"beta"} to specify a custom prior using \code{a_base} and \code{b_base}.
+#' @param a_base Numeric. Alpha parameter for the Beta prior (only used if \code{prior = "beta"}).
+#' @param b_base Numeric. Beta parameter for the Beta prior (only used if \code{prior = "beta"}).
+#' @param p_null Optional. True response probability under the null hypothesis (e.g., \code{p_null = 0.6}).
+#' If not specified, defaults to \code{p_null = M} (boundary case).
+#'
+#' @return A named list with:
+#' \describe{
+#'   \item{\code{estimate}}{Exact Type-I error (a number between 0 and 1).}
+#'   \item{\code{mc_se}}{\code{NA_real_}, included for compatibility.}
+#'   \item{\code{B}}{\code{NA_integer_}, included for compatibility.}
+#'   \item{\code{rejections}}{\code{NA_integer_}, included for compatibility.}
+#' }
+#'
+#' @examples
+#' # Type-I error under flat prior at boundary
+#' singlearm_beta_type1_exact(n_t = 40, M = 0.65, threshold = 0.9, prior = "flat")
+#'
+#' # Type-I error under true p < M (frequentist view)
+#' singlearm_beta_type1_exact(n_t = 40, M = 0.65, threshold = 0.9,
+#'                            prior = "flat", p_null = 0.60)
+#'
+#' @seealso \code{\link{singlearm_beta_type1}} for the simulation-based version.
+#'
+#' @author Thomas Debray \email{[email protected]}
+#' @export
+singlearm_beta_type1_exact <- function(n_t, M, threshold, prior = "flat", a_base = 1, b_base = 1, p_null = -1.0) {
+    .Call(`_bcts_singlearm_beta_type1_exact`, n_t, M, threshold, prior, a_base, b_base, p_null)
 }
 

R/bcts_power_sat_betaBinom_conj.R

@@ -1,22 +1,46 @@
-#' Simulate Bayesian power for a single-arm binomial trial (via Rcpp)
+#' Estimate Bayesian power for a single-arm binomial trial
 #'
-#' @param B Number of trial simulations.
+#' @description Computes the Bayesian power for a single-arm binomial trial
+#' using either simulation or exact analytic summation over all possible outcomes.
+#'
+#' @param B Number of trial simulations (used only if \code{method = "simulate"}).
 #' @param p_t True response probability.
 #' @param n_t Sample size.
 #' @param M Threshold for success (e.g. 0.6).
 #' @param threshold Posterior probability threshold (e.g. 0.95).
 #' @param prior Prior type: "flat" or "beta".
-#' @param a_base Alpha of Beta prior (only used if prior = "beta").
-#' @param b_base Beta of Beta prior (only used if prior = "beta").
-#' @param n_draws Number of posterior samples per trial.
-#' @param show_progress Show progress bar?
+#' @param a_base Alpha of Beta prior (only used if \code{prior = "beta"}).
+#' @param b_base Beta of Beta prior (only used if \code{prior = "beta"}).
+#' @param show_progress Show progress bar? (only relevant for simulation).
+#' @param method Either \code{"simulate"} (default) or \code{"exact"}.
+#'
+#' @return A list with \code{estimate} (power), \code{mc_se}, \code{successes}, and \code{B}.
+#'
+#' @examples
+#' singlearm_beta_power(
+#'   B = 1000, p_t = 0.75, n_t = 35, M = 0.60,
+#'   threshold = 0.95, prior = "flat", method = "simulate"
+#' )
 #'
-#' @return A list with `estimate` (power), `mc_se`, `successes`, and `B`.
+#' singlearm_beta_power(
+#'   p_t = 0.75, n_t = 35, M = 0.60,
+#'   threshold = 0.95, prior = "flat", method = "exact"
+#' )
 #'
 #' @export
-singlearm_beta_power <- function(B, p_t, n_t, M, threshold,
+singlearm_beta_power <- function(B = 10000, p_t, n_t, M, threshold,
                                  prior = "flat", a_base = 1, b_base = 1,
-                                 n_draws = 2000, show_progress = TRUE) {
-  .Call(`_bcts_singlearm_beta_power`, B, p_t, n_t, M, threshold,
-        prior, a_base, b_base, n_draws, show_progress)
+                                 show_progress = TRUE,
+                                 method = c("exact", "simulate")) {
+  method <- match.arg(method)
+
+  if (method == "simulate") {
+    .Call(`_bcts_singlearm_beta_power`,
+          B, p_t, n_t, M, threshold,
+          prior, a_base, b_base, show_progress)
+  } else {
+    .Call(`_bcts_singlearm_beta_power_exact`,
+          p_t, n_t, M, threshold,
+          prior, a_base, b_base)
+  }
 }

R/bcts_type1_sat_betaBinom_conj.R

@@ -1,21 +1,64 @@
-#' Simulate Bayesian Type-I error for a single-arm binomial trial (via Rcpp)
+#' Estimate Bayesian Type-I Error for a Single-Arm Binomial Trial
 #'
-#' @param B Number of trial simulations.
-#' @param n_t Sample size.
-#' @param M Decision threshold (also used as null hypothesis value).
-#' @param threshold Posterior probability threshold (e.g. 0.95).
-#' @param prior Prior type: "flat" or "beta".
-#' @param a_base Alpha of Beta prior (only used if prior = "beta").
-#' @param b_base Beta of Beta prior (only used if prior = "beta").
-#' @param n_draws Number of posterior samples per trial.
-#' @param show_progress Show progress bar?
+#' @description Computes the Bayesian Type-I error for a single-arm binomial trial
+#' using a conjugate Beta prior. You can choose between a simulation-based approach
+#' or an exact summation over all possible outcomes under the null hypothesis.
 #'
-#' @return A list with `estimate` (type-I error), `mc_se`, `rejections`, and `B`.
+#' By default, the null hypothesis assumes \code{p_null = M} (i.e., testing at the boundary),
+#' but for a more conservative frequentist setting, you may specify \code{p_null < M}.
+#'
+#' @param B Integer. Number of simulations (only used if \code{method = "simulate"}).
+#' @param n_t Integer. Sample size of the treatment arm.
+#' @param M Numeric. Decision threshold (also used as null value if \code{p_null} not specified).
+#' @param threshold Numeric. Posterior probability threshold for declaring success, e.g., \code{0.95}.
+#' @param prior Character string. Either \code{"flat"} (Beta(1,1)) or \code{"beta"}.
+#' @param a_base Numeric. Alpha parameter of the Beta prior (used only if \code{prior = "beta"}).
+#' @param b_base Numeric. Beta parameter of the Beta prior (used only if \code{prior = "beta"}).
+#' @param p_null Optional numeric. True response rate under the null hypothesis (default is \code{M}).
+#' Only used when \code{method = "exact"}.
+#' @param n_draws Ignored. Kept for backward compatibility.
+#' @param show_progress Logical. Show progress bar during simulation (only relevant for \code{method = "simulate"}).
+#' @param method Character. Either \code{"exact"} (default) or \code{"simulate"}.
+#'
+#' @return A named list with:
+#' \describe{
+#'   \item{\code{estimate}}{Estimated Type-I error rate.}
+#'   \item{\code{mc_se}}{Monte Carlo standard error (NA for exact).}
+#'   \item{\code{rejections}}{Number of rejections (NA for exact).}
+#'   \item{\code{B}}{Number of simulations (NA for exact).}
+#' }
+#'
+#' @examples
+#' # Simulated Type-I error at boundary (p_null = M)
+#' singlearm_beta_type1(B = 10000, n_t = 40, M = 0.65, threshold = 0.9, method = "simulate")
+#'
+#' # Exact Type-I error at boundary (p_null = M)
+#' singlearm_beta_type1(n_t = 40, M = 0.65, threshold = 0.9, method = "exact")
+#'
+#' # Exact Type-I error with p_null < M (frequentist-style)
+#' singlearm_beta_type1(n_t = 40, M = 0.65, threshold = 0.9, method = "exact", p_null = 0.60)
+#'
+#' @seealso \code{\link{singlearm_beta_power}}, \code{\link{singlearm_beta_type1_exact}}
 #'
 #' @export
-singlearm_beta_type1 <- function(B, n_t, M, threshold,
+singlearm_beta_type1 <- function(B = 10000, n_t, M, threshold,
                                  prior = "flat", a_base = 1, b_base = 1,
-                                 n_draws = 2000, show_progress = TRUE) {
-  .Call(`_bcts_singlearm_beta_type1`, B, n_t, M, threshold,
-        prior, a_base, b_base, n_draws, show_progress)
+                                 p_null = NULL,
+                                 n_draws = 2000,  # ignored
+                                 show_progress = TRUE,
+                                 method = c("exact", "simulate")) {
+  method <- match.arg(method)
+
+  if (method == "simulate") {
+    .Call(`_bcts_singlearm_beta_type1`,
+          B, n_t, M, threshold,
+          prior, a_base, b_base, show_progress)
+  } else {
+    # Default p_null = M if not specified
+    if (is.null(p_null)) p_null <- M
+
+    .Call(`_bcts_singlearm_beta_type1_exact`,
+          n_t, M, threshold,
+          prior, a_base, b_base, p_null)
+  }
 }

inst/app/app.R

@@ -410,7 +410,8 @@ server <- function(input, output, session) {
         prior = prior_type,
         a_base = a_base,
         b_base = b_base,
-        n_draws = n_draws,
+        #n_draws = n_draws,
+        method = "exact",
         show_progress = FALSE
       )