fix code

Author: NightlordTW

NAMESPACE

@@ -40,6 +40,8 @@ export(rct_beta_calibrate)
 export(rct_beta_power)
 export(rct_power_beta_binom_cpp_vec)
 export(run_bcts_app)
+export(sat_betabinom_calibrate)
+export(sat_betabinom_calibrate_exact)
 export(sat_betabinom_power)
 export(sat_betabinom_power_exact)
 export(sat_betabinom_type1)

R/RcppExports.R

@@ -153,15 +153,11 @@ sat_betabinom_power <- function(B, p_t, n_t, M, threshold, prior = "flat", a_bas
 #'
 #' @param B Integer. Number of simulations.
 #' @param n_t Integer. Sample size of the treatment arm.
+#' @param p_null True response probability under H₀. Used to simulate Binomial outcomes.
 #' @param M Numeric. Decision threshold for θ (on probability scale, e.g., 0.6).
 #' @param threshold Posterior probability threshold γ (e.g., 0.95).
-#' @param prior Character string specifying the prior distribution.
-#' Options are:
-#' "flat" for a non-informative Beta(1,1) prior;
-#' "jeffreys" for the Jeffreys prior (Beta(0.5, 0.5));
-#' "beta" for a custom Beta(\code{a_base}, \code{b_base}) prior (user must provide \code{a_base} and \code{b_base}).
-#' @param a_base Alpha parameter for Beta prior (if prior = "beta").
-#' @param b_base Beta parameter for Beta prior (if prior = "beta").
+#' @param a_base Alpha parameter for Beta prior.
+#' @param b_base Beta parameter for Beta prior.
 #' @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}.
@@ -174,8 +170,8 @@ sat_betabinom_power <- function(B, p_t, n_t, M, threshold, prior = "flat", a_bas
 #'
 #' @author Thomas Debray \email{[email protected]}
 #' @export
-sat_betabinom_type1 <- function(B, n_t, M, threshold, prior = "flat", a_base = 1, b_base = 1, show_progress = TRUE) {
-    .Call(`_bcts_sat_betabinom_type1`, B, n_t, M, threshold, prior, a_base, b_base, show_progress)
+sat_betabinom_type1 <- function(B, n_t, p_null, M, threshold, a_base = 1, b_base = 1, show_progress = TRUE) {
+    .Call(`_bcts_sat_betabinom_type1`, B, n_t, p_null, M, threshold, a_base, b_base, show_progress)
 }
 
 #' @title Exact Type-I Error for Single-Arm Trial (Beta-Binomial)
@@ -192,37 +188,27 @@ sat_betabinom_type1 <- function(B, n_t, M, threshold, prior = "flat", a_base = 1
 #' @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 specifying the prior distribution.
-#' Options are:
-#' "flat" for a non-informative Beta(1,1) prior;
-#' "jeffreys" for the Jeffreys prior (Beta(0.5, 0.5));
-#' "beta" for a custom Beta(\code{a_base}, \code{b_base}) prior (user must provide \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 a_base Numeric. Alpha parameter for the Beta prior.
+#' @param b_base Numeric. Beta parameter for the Beta prior.
 #' @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{B}}{\code{NA_real_}, included for compatibility.}
 #'   \item{\code{rejections}}{\code{NA_integer_}, included for compatibility.}
 #' }
 #'
 #' @examples
-#' # Type-I error under flat prior at boundary
-#' sat_betabinom_type1_exact(n_t = 40, M = 0.65, threshold = 0.9, prior = "flat")
-#'
-#' # Type-I error under true p < M (frequentist view)
-#' sat_betabinom_type1_exact(n_t = 40, M = 0.65, threshold = 0.9,
-#'                            prior = "flat", p_null = 0.60)
+#' sat_betabinom_type1_exact(n_t = 40, M = 0.65, threshold = 0.9)
 #'
 #' @seealso \code{\link{sat_betabinom_type1}} for the simulation-based version.
 #'
 #' @author Thomas Debray \email{[email protected]}
 #' @export
-sat_betabinom_type1_exact <- function(n_t, M, threshold, prior = "flat", a_base = 1, b_base = 1, p_null = -1.0) {
-    .Call(`_bcts_sat_betabinom_type1_exact`, n_t, M, threshold, prior, a_base, b_base, p_null)
+sat_betabinom_type1_exact <- function(n_t, M, threshold, a_base = 1, b_base = 1, p_null = -1.0) {
+    .Call(`_bcts_sat_betabinom_type1_exact`, n_t, M, threshold, a_base, b_base, p_null)
 }
 

R/rct_beta_calibrate.R R/rct_betabinom_calibrate.RR/rct_beta_calibrate.R renamed to R/rct_betabinom_calibrate.R

@@ -38,6 +38,7 @@
 #'   method = "cpp", show_progress = FALSE
 #' )
 #'
+#' @author Thomas Debray \email{[email protected]}
 #' @export
 rct_beta_power <- function(B = 10000,
                            p_c, p_t, n_c, n_t,

R/rct_beta_power.R R/rct_betabinom_power.RR/rct_beta_power.R renamed to R/rct_betabinom_power.R

@@ -38,6 +38,7 @@
 #' )
 #' res$estimate
 #' res$ci
+#' @author Thomas Debray \email{[email protected]}
 #' @export
 bcts_power_betaBinom_conj <- function(B = 1000, p_c, p_t, n_c, n_t, M,
                                          threshold,

R/rct_betabinom_sim.R

@@ -0,0 +1,140 @@
+#' Calibrate posterior cutoff gamma and estimate power (Single-Arm, Beta–Binomial)
+#'
+#' Calibrates the posterior threshold `gamma` to control Type-I error at level `alpha`,
+#' and estimates Bayesian power using a conjugate Beta–Binomial model.
+#'
+#' @param alpha Target Type-I error rate (numeric in `(0,1)`).
+#' @param p_t True response probability in treatment group (numeric in `[0,1]`).
+#' @param n_t Sample size in treatment group (integer).
+#' @param M Decision threshold for theta (numeric in `[0,1]`).
+#' @param a_base,b_base Numeric. Hyperparameters of the Beta prior (used when `prior = "beta"`).
+#' @param tol Absolute tolerance for calibration convergence.
+#' @param seed RNG seed.
+#' @param show_progress Logical. Whether to display progress bars.
+#'
+#' @return A list with components:
+#'   - `gamma`: Calibrated posterior probability threshold.
+#'   - `type1`: Estimated Type-I error at `gamma`.
+#'   - `power`: Estimated Bayesian power at `gamma`.
+#'   - `settings`: Echo of input parameters.
+#'
+#' @examples
+#' \donttest{
+#' res <- sat_beta_calibrate(
+#'   alpha = 0.10, p_t = 0.80, n_t = 40, M = 0.65,
+#'   prior = "flat"
+#' )
+#' res$gamma
+#' res$power
+#' }
+#' @author Thomas Debray \email{[email protected]}
+#' @export
+sat_betabinom_calibrate <- function(alpha = 0.10,
+                              p_t = 0.80,
+                              n_t = 40,
+                              M = 0.65,
+                              prior = c("flat", "jeffreys", "beta"),
+                              a_base = 1, b_base = 1,
+                              tol = 0.002,
+                              seed = 123,
+                              show_progress = FALSE) {
+  prior <- match.arg(prior)
+
+  cal <- sat_betabinom_find_gamma(
+    alpha = alpha, p_c = M, M = M,
+    n_t = n_t,
+    prior = prior,
+    prior_args = list(a_base = a_base, b_base = b_base),
+    seed = seed,
+    show_progress = show_progress,
+    progress_fun = pf,
+    calibrate_on = calibrate_on
+  )
+
+
+  power <- sat_betabinom_power(
+    p_t = p_t, n_t = n_t, M = M,
+    threshold = cal$gamma,
+    prior = prior, a_base = a_base, b_base = b_base,
+    show_progress = show_progress,
+    method = "exact"
+  )
+
+  type1 <- sat_betabinom_type1(
+    n_t = n_t, M = M,
+    threshold = cal$gamma,
+    prior = prior, a_base = a_base, b_base = b_base,
+    show_progress = show_progress,
+    method = "exact"
+  )
+
+  out <- list(
+    gamma = cal$gamma,
+    type1 = type1,
+    power = power,
+    settings = list(
+      alpha = alpha, p_t = p_t, n_t = n_t, M = M,
+      prior = prior, a_base = a_base, b_base = b_base,
+      B_cal = B_cal, B_power = B_power,
+      tol = tol, seed = seed
+    )
+  )
+  class(out) <- "bcts_sat"
+  return(out)
+}
+
+#' Calibrate posterior threshold gamma to control Type-I error (exact, SAT)
+#'
+#' For single-arm trials with Beta–Binomial conjugate updates, computes the Type-I error
+#' exactly for all possible posterior thresholds, and selects the largest gamma such that
+#' Type-I error is below or equal to a target alpha.
+#'
+#' @param alpha Target Type-I error rate (e.g. 0.10)
+#' @param n_t Sample size in treatment arm
+#' @param M Decision margin (null threshold)
+#' @param a_base,b_base Prior Beta(a,b) parameters (default: 1,1)
+#'
+#' @return An object of class `"bcts_calibration"` with components: `gamma`, `type1`, `trace`, and `settings`.
+#' @author Thomas Debray \email{[email protected]}
+#' @export
+sat_betabinom_calibrate_exact <- function(alpha = 0.10,
+                                          n_t,
+                                          M,
+                                          a_base = 1,
+                                          b_base = 1,
+                                          grid_n = 1000) {
+
+  # Step 1: Grid of gamma values (search over [0,1])
+  gamma_vals <- seq(1, 0, length.out = grid_n)  # from most to least conservative
+
+  # Step 2: Evaluate type-I error at each gamma
+  type1_vals <- vapply(gamma_vals, function(g) {
+    .Call(`_bcts_sat_betabinom_type1_exact`,
+          n_t, M, g, a_base, b_base, -1.0)
+  }, numeric(1))
+
+  # Step 3: Build trace table
+  trace <- data.frame(
+    gamma = gamma_vals,
+    type1 = type1_vals
+  )
+
+  # Step 4: Select gamma where type-I is just below alpha
+  eligible <- trace$type1 <= alpha
+  if (!any(eligible)) {
+    stop("Calibration failed: Type-I error is above target alpha at all gamma values.")
+  }
+
+  best_idx <- which.max(trace$gamma[eligible])  # most liberal gamma ≤ alpha
+  best_gamma <- trace$gamma[eligible][best_idx]
+  best_type1 <- trace$type1[eligible][best_idx]
+
+  out <- list(
+    gamma = best_gamma,
+    type1 = best_type1,
+    trace = trace
+  )
+  class(out) <- "bcts_calibration"
+
+  return(out)
+}

R/sat_betabinom_calibrate.R

@@ -29,6 +29,7 @@
 #'   p_t = 0.75, n_t = 35, M = 0.60,
 #'   threshold = 0.95, prior = "flat", method = "exact"
 #' )
+#' @author Thomas Debray \email{[email protected]}
 #'
 #' @export
 sat_betabinom_power <- function(B = 10000, p_t, n_t, M, threshold,

R/sat_betabinom_power.R

@@ -12,14 +12,13 @@
 #' @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 Prior type: one of:
+#'   - \code{"beta"} (default) Use custom Beta(\code{a_base}, \code{b_base}) prior..
 #'   - \code{"flat"} for a non-informative Beta(1,1) prior;
 #'   - \code{"jeffreys"} for the Jeffreys prior Beta(0.5, 0.5);
-#'   - \code{"beta"} for a custom Beta(\code{a_base}, \code{b_base}) prior.
 #' @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 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"}.
 #'
@@ -43,27 +42,40 @@
 #'
 #' @seealso \code{\link{sat_betabinom_power}}, \code{\link{sat_betabinom_type1_exact}}
 #'
+#' @author Thomas Debray \email{[email protected]}
 #' @export
 sat_betabinom_type1 <- function(B = 10000, n_t, M, threshold,
-                                prior = c("flat", "jeffreys", "beta"),
+                                prior = c("beta", "flat", "jeffreys"),
                                 a_base = 1, b_base = 1,
                                  p_null = NULL,
-                                 n_draws = 2000,  # ignored
                                  show_progress = TRUE,
                                  method = c("exact", "simulate")) {
   prior <- match.arg(prior)
   method <- match.arg(method)
 
+  # Convert prior to a_base and b_base
+  if (prior == "flat") {
+    a_base <- 1
+    b_base <- 1
+  } else if (prior == "jeffreys") {
+    a_base <- 0.5
+    b_base <- 0.5
+  } else if (prior == "beta") {
+    # leave a_base and b_base as provided
+  } else {
+    stop("Invalid prior. Must be one of 'flat', 'jeffreys', or 'beta'.")
+  }
+
+  # Default p_null = M if not specified
+  if (is.null(p_null)) p_null <- M
+
   if (method == "simulate") {
     .Call(`_bcts_sat_betabinom_type1`,
-          B, n_t, M, threshold,
-          prior, a_base, b_base, show_progress)
+          B, n_t, p_null, M, threshold,
+          a_base, b_base, show_progress)
   } else {
-    # Default p_null = M if not specified
-    if (is.null(p_null)) p_null <- M
-
     .Call(`_bcts_sat_betabinom_type1_exact`,
           n_t, M, threshold,
-          prior, a_base, b_base, p_null)
+          a_base, b_base, p_null)
   }
 }

R/sat_betabinom_type1.R

@@ -127,13 +127,18 @@ mod_singlearm_server <- function(id) {
         show_progress = FALSE
       )
 
-      # Frequentist power
-      crit_val <- qbinom(gamma, size = nt, prob = M) + 1
-      freq_power <- 1 - pbinom(crit_val - 1, size = nt, prob = pt)
-
-
-      # Frequentist type-I error
-      freq_type1 <- 1 - gamma
+      # Optional frequentist comparison
+      freq_power <- freq_type1 <- NA
+      decision_mode <- input[["crit_sa-decision_mode_sa"]]
+
+      if (!is.null(decision_mode) && decision_mode == "alpha") {
+        alpha <- input[["crit_sa-alpha_sa"]] / 100  # slider is in percent
+        if (!is.null(alpha)) {
+          crit_val <- qbinom(1 - alpha, size = nt, prob = M)
+          freq_power <- 1 - pbinom(crit_val - 1, size = nt, prob = pt)
+          freq_type1 <- 1 - pbinom(crit_val - 1, size = nt, prob = M)
+        }
+      }
 
       list(
         power = power_res,
@@ -207,11 +212,11 @@ mod_singlearm_server <- function(id) {
       cat(sprintf("Estimated power: %.2f%%\n", 100 * power_res$estimate))
       cat(sprintf("Estimated Type-I error: %.2f%%\n", 100 * type1_res$estimate))
 
-      cat("\n")
-
-      cat("FREQUENTIST COMPARISON\n")
-      cat(sprintf("Frequentist power: %.2f%%\n", 100 * freq_power))
-      cat(sprintf("Frequentist Type-I error: %.2f%%\n", 100 * freq_type1))
+      if (!is.na(res$freq_power) && !is.na(res$freq_type1)) {
+        cat("\nFREQUENTIST COMPARISON\n")
+        cat(sprintf("Frequentist power: %.2f%%\n", 100 * res$freq_power))
+        cat(sprintf("Frequentist Type-I error: %.2f%%\n", 100 * res$freq_type1))
+      }
     })
 
     output$freq_binom_plot <- renderPlot({