update app

Author: NightlordTW

DESCRIPTION

@@ -12,7 +12,7 @@ Description: Sample size calculation for adaptive trials.
 License: Apache License 2.0
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.3.2
+RoxygenNote: 7.3.3
 Depends: 
     R (>= 4.1.0)
 Imports: 
@@ -36,6 +36,6 @@ Suggests:
     tinytest,
     kableExtra,
     tibble
-LinkingTo: Rcpp
+LinkingTo: Rcpp, RcppProgress
 SystemRequirements: JAGS (>= 4.0.0)
 Config/testthat/edition: 3

NAMESPACE

@@ -37,6 +37,8 @@ export(plot_prior_weight)
 export(posterior)
 export(power)
 export(power_fisher)
+export(rct_beta_power)
+export(rct_power_beta_binom_cpp_vec)
 export(run_bcts_app)
 export(sim_rct_normal)
 export(singlearm_beta_power)

R/RcppExports.R

@@ -155,3 +155,61 @@ singlearm_beta_type1_exact <- function(n_t, M, threshold, prior = "flat", a_base
     .Call(`_bcts_singlearm_beta_type1_exact`, n_t, M, threshold, prior, a_base, b_base, p_null)
 }
 
+#' Estimate Power for Bayesian RCT Using Beta-Binomial Conjugate Model (C++)
+#'
+#' Performs Monte Carlo simulation to estimate the statistical power of a
+#' Bayesian randomized controlled trial (RCT) using conjugate Beta-Binomial
+#' models for the control and treatment arms. The function supports flat priors
+#' and power priors on the control arm, and allows for fully customizable
+#' baseline Beta prior parameters for both arms.
+#'
+#' A trial is declared successful if the posterior probability
+#' \eqn{Pr(\theta_t - \theta_c > M \mid \text{data}) \ge \gamma}.
+#'
+#' @param B Integer. Number of simulated trials.
+#' @param p_c,p_t Numeric in \[0, 1\]. True response probabilities in the control and treatment arms.
+#' @param n_c,n_t Integers. Sample sizes in the control and treatment arms.
+#' @param M Numeric. Margin on the risk-difference scale:
+#'   - Negative for non-inferiority,
+#'   - Zero for equivalence,
+#'   - Positive for superiority.
+#' @param threshold Numeric in (0, 1). Posterior probability threshold \eqn{\gamma}.
+#' @param prior Character. Type of prior to use:
+#'   - `"flat"`: independent Beta priors for both arms,
+#'   - `"power"`: power prior on the control arm, flat prior on treatment.
+#' @param prior_args List of prior hyperparameters. The following elements are supported:
+#'   - `a0`: Discount factor for historical control data (only used if `prior = "power"`).
+#'   - `y_0`, `n_0`: Number of responses and total patients in the historical control data.
+#'   - `a_base_c`, `b_base_c`: Shape parameters of the Beta prior for the control arm.
+#'   - `a_base_t`, `b_base_t`: Shape parameters of the Beta prior for the treatment arm.
+#'     - Defaults for all `a_base_*` and `b_base_*` values are 1 (i.e., flat prior).
+#' @param n_draws Integer. Number of posterior draws per trial for estimating the probability.
+#' @param show_progress Logical. Show progress bar in the console.
+#'
+#' @return A logical vector of length `B`, indicating for each trial whether the decision
+#' criterion was met (i.e., trial declared successful).
+#'
+#' @details This function uses Rcpp and vectorized binomial simulation to increase speed.
+#' Posterior samples are drawn from Beta distributions parameterized using either flat priors
+#' or power priors (for control) combined with observed trial data.
+#'
+#' @examples
+#' prior_args <- list(
+#'   a0 = 0.5, y_0 = 20, n_0 = 30,
+#'   a_base_c = 1, b_base_c = 1,
+#'   a_base_t = 2, b_base_t = 2
+#' )
+#' decisions <- rct_power_beta_binom_cpp_vec(
+#'   B = 1000, p_c = 0.8, p_t = 0.8,
+#'   n_c = 25, n_t = 25,
+#'   M = -0.1, threshold = 0.9,
+#'   prior = "power", prior_args = prior_args,
+#'   n_draws = 2000, show_progress = FALSE
+#' )
+#' mean(decisions)  # Estimated power
+#'
+#' @export
+rct_power_beta_binom_cpp_vec <- function(B, p_c, p_t, n_c, n_t, M, threshold, prior, prior_args, n_draws, show_progress = TRUE) {
+    .Call(`_bcts_rct_power_beta_binom_cpp_vec`, B, p_c, p_t, n_c, n_t, M, threshold, prior, prior_args, n_draws, show_progress)
+}
+

R/bayesNI-conj-calibrate.R

@@ -38,9 +38,6 @@
 #' )
 #' res$estimate
 #' res$ci
-#'
-#' @seealso [bcts_type1_betaBinom_conj()], [bcts_calibrate_betaBinom_conj()]
-#' @export
 bcts_power_betaBinom_conj <- function(B = 1000, p_c, p_t, n_c, n_t, M,
                                          threshold,
                                          prior = c("flat","power"),
@@ -130,10 +127,7 @@ bcts_power_betaBinom_conj <- function(B = 1000, p_c, p_t, n_c, n_t, M,
 #'   prior     = "flat",
 #'   n_draws   = 200
 #' )
-#'
-#' @seealso [bcts_calibrate_betaBinom_conj()], [bcts_power_betaBinom_conj()]
 #' @family conjugate-BetaBinom
-#' @export
 bcts_type1_betaBinom_conj <- function(B = 2000, p_c, M, n_c, n_t,
                                       threshold, prior = c("flat","power"),
                                       prior_args = list(),

R/bayesNI-conj-run.R R/rct_beta_calibrate.RR/bayesNI-conj-run.R renamed to R/rct_beta_calibrate.R

@@ -21,30 +21,32 @@
 #' @param seed RNG seed for reproducibility.
 #' @param show_progress Logical. Show text progress bars during simulation.
 #' @param verbose Logical. Print iteration messages during calibration.
+#' @param method Character. Method for estimation. Either `"simulate"` or `"cpp"`.
 #'
 #' @return An object of class `"bayesNI"` with components:
 #' - `calibration`: list from [bcts_calibrate_betaBinom_conj()]
 #' - `power`: estimated Bayesian power at calibrated `gamma`
 #' - `settings`: list of all input parameters (echoed)
 #'
 #' @seealso [bcts_calibrate_betaBinom_conj()],
-#'   [bcts_power_betaBinom_conj()],
+#'   [rct_beta_power()],
 #'   [bcts_type1_betaBinom_conj()]
 #' @family conjugate-NI
 #'
 #' @examples
 #' \donttest{
-#' res <- bcts_run_betaBinom_conj(
+#' res <- rct_beta_calibrate(
 #'   alpha = 0.10, p_c = .85, p_t = .85,
 #'   n_c = 29, n_t = 29, NI_margin = -0.20,
-#'   prior = "flat", B_cal = 1000, B_power = 500
+#'   prior = "flat", B_cal = 1000, B_power = 500,
+#'   method = "cpp
 #' )
 #' res$calibration$gamma
 #' res$power
 #' }
 #'
 #' @export
-bcts_run_betaBinom_conj <- function(alpha = 0.10,
+rct_beta_calibrate <- function(alpha = 0.10,
                               p_c = 0.85, p_t = 0.85,
                               n_c = 29, n_t = 29,
                               NI_margin = -0.20,
@@ -54,9 +56,12 @@ bcts_run_betaBinom_conj <- function(alpha = 0.10,
                               n_draws = 2000,
                               tol = 0.002,
                               seed = 123,
-                              show_progress = FALSE, verbose = TRUE) {
+                              show_progress = FALSE, verbose = TRUE,
+                              method = c("cpp", "simulate")) {
   prior <- match.arg(prior)
 
+  method <- match.arg(method)
+
   # 1) Calibrate gamma using conjugate machinery
   cal <- bcts_calibrate_betaBinom_conj(alpha = alpha, p_c = p_c, M = NI_margin,
                               n_c = n_c, n_t = n_t,
@@ -66,12 +71,14 @@ bcts_run_betaBinom_conj <- function(alpha = 0.10,
                               show_progress = show_progress, verbose = verbose)
 
   # 2) Simulate Bayesian "power" at calibrated gamma
-  power <- bcts_power_betaBinom_conj(B = B_power, p_c = p_c, p_t = p_t,
-                                     n_c = n_c, n_t = n_t, M = NI_margin,
-                                     threshold = cal$gamma,
-                                     prior = prior, prior_args = prior_args,
-                                     n_draws = n_draws,
-                                     seed = seed, show_progress = show_progress)
+  power <- rct_beta_power(B = B_power, p_c = p_c, p_t = p_t,
+                          n_c = n_c, n_t = n_t, M = NI_margin,
+                          threshold = cal$gamma,
+                          prior = prior, prior_args = prior_args,
+                          n_draws = n_draws,
+                          show_progress = show_progress,
+                          method = method, seed = seed)
+
 
   out <- list(
     calibration = cal,

R/rct_beta_power.R

@@ -0,0 +1,103 @@
+#' Estimate Bayesian Power for a 2-arm Binomial RCT (Beta-Binomial Conjugate)
+#'
+#' @description Computes the power of a Bayesian randomized controlled trial
+#' using conjugate Beta-Binomial models. Can use exact simulation in R or fast
+#' simulation in C++.
+#'
+#' @param B Integer. Number of trial simulations (ignored if \code{method = "exact"}).
+#' @param p_c,p_t Numeric in \[0,1\]. True response probabilities in control and treatment arms.
+#' @param n_c,n_t Integers. Sample sizes in control and treatment arms.
+#' @param M Numeric. Margin on the risk-difference scale. Default is 0 (superiority/equivalence).
+#' @param threshold Posterior probability cutoff (e.g. 0.9).
+#' @param prior Character. Prior type: \code{"flat"} or \code{"power"}.
+#' @param prior_args List with prior hyperparameters. Optional elements:
+#'   - \code{a0}, \code{y_0}, \code{n_0} for power prior
+#'   - \code{a_base_c}, \code{b_base_c} for control arm prior
+#'   - \code{a_base_t}, \code{b_base_t} for treatment arm prior
+#' @param n_draws Integer. Number of posterior samples per trial.
+#' @param show_progress Logical. Show text progress bar?
+#' @param method Character. Either \code{"simulate"} (pure R) or \code{"cpp"} (fast C++).
+#' @param seed Optional integer. If provided, sets the RNG seed for reproducibility.
+#'
+#' @return A list with:
+#' \item{estimate}{Estimated power (probability of success)}
+#' \item{mc_se}{Monte Carlo standard error}
+#' \item{successes}{Number of trials with successful outcome}
+#' \item{B}{Number of simulations}
+#'
+#' @examples
+#' # Using C++ backend
+#' rct_beta_power(
+#'   B = 1000, p_c = 0.85, p_t = 0.85,
+#'   n_c = 30, n_t = 30, M = -0.1,
+#'   threshold = 0.9, prior = "flat",
+#'   method = "cpp", show_progress = FALSE
+#' )
+#'
+#' @export
+rct_beta_power <- function(B = 10000,
+                           p_c, p_t, n_c, n_t,
+                           M = 0,
+                           threshold,
+                           prior = c("flat", "power"),
+                           prior_args = list(),
+                           n_draws = 2000,
+                           show_progress = TRUE,
+                           method = c("cpp", "simulate"),
+                           seed = NULL) {
+  prior <- match.arg(prior)
+  method <- match.arg(method)
+
+  if (!is.null(seed)) set.seed(seed)
+
+  if (method == "simulate") {
+    res <- bcts_power_betaBinom_conj(
+      B = B,
+      p_c = p_c,
+      p_t = p_t,
+      n_c = n_c,
+      n_t = n_t,
+      M = M,
+      threshold = threshold,
+      prior = prior,
+      prior_args = prior_args,
+      n_draws = n_draws,
+      show_progress = show_progress,
+      seed = seed
+    )
+    return(res)
+  }
+
+  if (method == "cpp") {
+    # Ensure required C++ function is available
+    dec_vec <- rct_power_beta_binom_cpp_vec(
+      B = B,
+      p_c = p_c,
+      p_t = p_t,
+      n_c = n_c,
+      n_t = n_t,
+      M = M,
+      threshold = threshold,
+      prior = prior,
+      prior_args = prior_args,
+      n_draws = n_draws,
+      show_progress = show_progress
+    )
+    k <- sum(dec_vec)
+    p <- mean(dec_vec)
+    mc_se <- sqrt(p * (1 - p) / B)
+
+    ci <- binom::binom.confint(x=k, n=B, conf.level = 0.95, methods = "wilson")
+
+    return(list(
+      estimate = p,
+      ci_lower = ci["lower"],
+      ci_upper = ci["upper"],
+      mc_se = mc_se,
+      successes = k,
+      B = B
+    ))
+  }
+
+  stop("Unsupported method. Choose either 'simulate' or 'cpp'.")
+}

inst/app/R/sim-core.R

@@ -29,10 +29,11 @@ run_all_oc <- function(pt, nt, pc, nc, M, mode, gamma, alpha, calibrate_on,
     n_draws = ndraws, show_progress = FALSE
   )
 
-  pw <- bcts_power_betaBinom_conj(
+  pw <- rct_beta_power(
     B = B, p_c = pc, p_t = pt, n_c = nc, n_t = nt, M = M,
     threshold = gamma_used, prior = prior, prior_args = prior_args,
-    n_draws = ndraws, seed = seed, show_progress = FALSE
+    n_draws = ndraws, seed = seed, show_progress = FALSE,
+    method = "cpp"
   )
 
   list(

inst/app/app.R

@@ -3,6 +3,7 @@
 
 # ---- Load packages that the APP needs ----
 library(shiny)
+library(shinyjs)
 library(ggplot2)
 library(bcts)
 
@@ -322,6 +323,7 @@ server <- function(input, output, session) {
     "sens",
     sim        = sim,                           # your eventReactive from primary analysis
     pt         = reactive(input$pt / 100),
+    pc         = reactive(input$pc / 100),
     nt         = reactive(input$nt),
     nc         = reactive(input$nc),
     M          = reactive(input$M / 100),