add visuals to interpret power and type-i error

Author: NightlordTW

NAMESPACE

@@ -20,7 +20,6 @@ export(bcts)
 export(bcts_calibrate_betaBinom_conj)
 export(bcts_one_interim)
 export(bcts_power_betaBinom_conj)
-export(bcts_run_betaBinom_conj)
 export(bcts_two_interim)
 export(bcts_type1_betaBinom_conj)
 export(dose_selection)
@@ -37,14 +36,15 @@ export(plot_prior_weight)
 export(posterior)
 export(power)
 export(power_fisher)
+export(rct_beta_calibrate)
 export(rct_beta_power)
 export(rct_power_beta_binom_cpp_vec)
 export(run_bcts_app)
+export(sat_betabinom_power)
+export(sat_betabinom_power_exact)
+export(sat_betabinom_type1)
+export(sat_betabinom_type1_exact)
 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

@@ -1,6 +1,64 @@
 # Generated by using Rcpp::compileAttributes() -> do not edit by hand
 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393
 
+#' 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)
+}
+
 beta_prob_gt <- function(a, b, M, n_draws) {
     .Call(`_bcts_beta_prob_gt`, a, b, M, n_draws)
 }
@@ -21,8 +79,11 @@ beta_prob_gt <- function(a, b, M, n_draws) {
 #' 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 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 (used only if \code{prior = "beta"}).
 #' @param b_base Numeric. Beta parameter for the Beta prior (used only if \code{prior = "beta"}).
 #'
@@ -35,17 +96,17 @@ beta_prob_gt <- function(a, b, M, n_draws) {
 #' }
 #'
 #' @examples
-#' singlearm_beta_power_exact(
+#' sat_betabinom_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.
+#' @seealso \code{\link{sat_betabinom_power}} for the simulation-based version.
 #'
 #' @author Thomas Debray \email{tdebray@fromdatatowisdom.com}
 #' @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)
+sat_betabinom_power_exact <- function(p_t, n_t, M, threshold, prior = "flat", a_base = 1, b_base = 1) {
+    .Call(`_bcts_sat_betabinom_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
@@ -61,24 +122,27 @@ singlearm_beta_power_exact <- function(p_t, n_t, M, threshold, prior = "flat", a
 #' @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{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 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 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(
+#' sat_betabinom_power(
 #'   B = 1000, p_t = 0.75, n_t = 35, M = 0.60,
 #'   threshold = 0.95, prior = "flat"
 #' )
 #'
 #' @author Thomas Debray \email{tdebray@fromdatatowisdom.com}
 #' @export
-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)
+sat_betabinom_power <- function(B, p_t, n_t, M, threshold, prior = "flat", a_base = 1, b_base = 1, show_progress = TRUE) {
+    .Call(`_bcts_sat_betabinom_power`, B, p_t, n_t, M, threshold, prior, a_base, b_base, show_progress)
 }
 
 #' @title Estimate Type-I Error for Single-Arm Trial
@@ -91,23 +155,27 @@ singlearm_beta_power <- function(B, p_t, n_t, M, threshold, prior = "flat", a_ba
 #' @param n_t Integer. Sample size of the treatment arm.
 #' @param M Numeric. Decision threshold for θ (on probability scale, e.g., 0.6).
 #' @param threshold Posterior probability threshold γ (e.g., 0.95).
-#' @param prior "flat" or "beta".
+#' @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 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(
+#' sat_betabinom_type1(
 #'   B = 1000, n_t = 35, M = 0.6,
 #'   threshold = 0.95, prior = "flat"
 #' )
 #'
 #' @author Thomas Debray \email{tdebray@fromdatatowisdom.com}
 #' @export
-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)
+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)
 }
 
 #' @title Exact Type-I Error for Single-Arm Trial (Beta-Binomial)
@@ -124,8 +192,11 @@ singlearm_beta_type1 <- function(B, n_t, M, threshold, prior = "flat", a_base =
 #' @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 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 p_null Optional. True response probability under the null hypothesis (e.g., \code{p_null = 0.6}).
@@ -141,75 +212,17 @@ singlearm_beta_type1 <- function(B, n_t, M, threshold, prior = "flat", a_base =
 #'
 #' @examples
 #' # Type-I error under flat prior at boundary
-#' singlearm_beta_type1_exact(n_t = 40, M = 0.65, threshold = 0.9, prior = "flat")
+#' sat_betabinom_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,
+#' sat_betabinom_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.
+#' @seealso \code{\link{sat_betabinom_type1}} for the simulation-based version.
 #'
 #' @author Thomas Debray \email{tdebray@fromdatatowisdom.com}
 #' @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)
-}
-
-#' 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)
+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)
 }
 

R/rct_beta_calibrate.R

@@ -8,12 +8,15 @@
 #' @param p_c,p_t True response probabilities in control and treatment (numeric in `[0,1]`).
 #' @param n_c,n_t Sample sizes in control and treatment arms (integers).
 #' @param NI_margin Non-inferiority margin on the risk-difference scale (numeric).
-#' @param prior Character. Either `"flat"` (Beta(1,1) both arms) or `"power"`
-#'   (power prior on control); see `prior_args`.
-#' @param prior_args List of hyperparameters when `prior = "power"`:
-#'   - `a0`: discount factor in `[0,1]`
-#'   - `y_0`, `n_0`: historical control responders and sample size
-#'   - `a_base`, `b_base`: baseline Beta parameters (defaults `1`, `1`)
+#' @param prior Character. Specifies the prior distribution. Options:
+#'   - `"flat"`: Flat Beta(1,1) prior for both arms (non-informative),
+#'   - `"jeffreys"`: Jeffreys prior, i.e., Beta(0.5, 0.5) for both arms,
+#'   - `"power"`: Power prior on the control arm, requires specification in `prior_args`.
+#' @param prior_args List of additional prior hyperparameters. Required when `prior = "power"`:
+#'   - `a0`: Discount factor for historical control data (numeric in `[0,1]`),
+#'   - `y_0`, `n_0`: Historical control responders and total sample size (integers),
+#'   - `a_base`, `b_base`: Baseline Beta parameters (default is `1`, `1`).
+#'   Ignored when `prior` is `"flat"` or `"jeffreys"`.
 #' @param B_cal Number of simulated trials for calibration of `gamma`.
 #' @param B_power Number of simulated trials for power estimation.
 #' @param n_draws Posterior Monte Carlo draws per trial for evaluating `Pr(NI)`.
@@ -39,7 +42,7 @@
 #'   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,
-#'   method = "cpp
+#'   method = "cpp"
 #' )
 #' res$calibration$gamma
 #' res$power
@@ -50,7 +53,7 @@ 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,
-                              prior = c("flat","power"),
+                              prior = c("flat", "jeffreys", "power"),
                               prior_args = list(),
                               B_cal = 2000, B_power = 1000,
                               n_draws = 2000,
@@ -94,3 +97,69 @@ rct_beta_calibrate <- function(alpha = 0.10,
   class(out) <- "bayesNI"
   out
 }
+
+
+#' Plot calibration trace for calibrated Bayesian NI threshold
+#'
+#' Visualizes the estimated Type-I error at different posterior thresholds
+#' during the calibration procedure.
+#'
+#' @param x Output from [bcts_calibrate_betaBinom_conj()]
+#' @param ... Not used.
+#'
+#' @return A [ggplot2::ggplot()] object.
+#' @export
+#'
+#' @examples
+#' # res <- bcts_calibrate_betaBinom_conj(...)  # Run separately
+#' # plot(res)
+plot.bcts_calibration <- function(x, ...) {
+  stopifnot(!is.null(x$trace), !is.null(x$type1))
+
+  tr     <- x$trace
+  alpha  <- x$alpha
+  gamma  <- x$gamma
+  type1  <- x$type1
+
+  # Small vertical offset above CI
+  offset <- 0.01
+
+  ggplot2::ggplot(tr, ggplot2::aes(x = .data$gamma_try, y = .data$type1)) +
+    ggplot2::geom_errorbar(
+      ggplot2::aes(ymin = .data$ci_lower, ymax = .data$ci_upper),
+      width = 0.002, alpha = 0.4
+    ) +
+    ggplot2::geom_point(size = 2) +
+    ggplot2::geom_line() +
+    ggplot2::geom_point(
+      x = gamma,
+      y = type1["estimate"],
+      color = "blue",
+      size = 3
+    ) +
+    ggplot2::annotate(
+      "text",
+      x = gamma,
+      y = type1["ci_upper"] + offset,
+      label = sprintf("gamma = %.3f\nType-I = %.1f%%",
+                      gamma, 100 * type1["estimate"]),
+      hjust = 0.5,
+      size = 3.5
+    ) +
+    ggplot2::geom_hline(
+      yintercept = alpha,
+      linetype = "dotted",
+      color = "red"
+    ) +
+    ggplot2::labs(
+      x = expression(gamma),
+      y = "Estimated Type-I error",
+      subtitle = sprintf(
+        "Calibrated gamma = %.3f after %d steps; dotted = target alpha = %.2f",
+        gamma, x$iters, alpha
+      )
+    ) +
+    ggplot2::theme_minimal(base_size = 12)+
+    ggplot2::scale_y_continuous(labels = scales::label_percent(accuracy = 1))+
+    ggplot2::scale_x_continuous(labels = scales::label_percent(accuracy = 1))
+}