1- # ' Bayesian trial power estimation with Beta– Binomial conjugate model
1+ # ' Bayesian trial power estimation with Beta- Binomial conjugate model
22# '
3- # ' Uses Monte Carlo simulation with conjugate Beta– Binomial updates
3+ # ' Uses Monte Carlo simulation with conjugate Beta- Binomial updates
44# ' (no MCMC) to estimate the probability of declaring trial success
55# ' under user-specified assumptions. A trial is declared successful
66# ' if \eqn{Pr(\theta_t - \theta_c > M | data) \ge \gamma}.
@@ -80,10 +80,10 @@ bcts_power_betaBinom_conj <- function(B = 1000, p_c, p_t, n_c, n_t, M,
8080 B = B , successes = k )
8181}
8282
83- # ' Empirical Type-I Error for Bayesian Trial (Beta– Binomial, conjugate)
83+ # ' Empirical Type-I Error for Bayesian Trial (Beta- Binomial, conjugate)
8484# '
8585# ' Estimates the Type-I error rate for a Bayesian non-inferiority (or
86- # ' superiority) test using conjugate Beta– Binomial updates. The test declares
86+ # ' superiority) test using conjugate Beta- Binomial updates. The test declares
8787# ' success if the posterior probability \eqn{Pr(\theta_t - \theta_c > M | data)}
8888# ' exceeds a fixed threshold `threshold`. Simulation is conducted under the
8989# ' least-favourable null (LFN), where \eqn{p_t = p_c + M}.
@@ -197,12 +197,12 @@ bcts_type1_betaBinom_conj <- function(B = 2000, p_c, M, n_c, n_t,
197197 )
198198}
199199
200- # ' Calibrate the posterior threshold gamma (Beta– Binomial, conjugate)
200+ # ' Calibrate the posterior threshold gamma (Beta- Binomial, conjugate)
201201# '
202202# ' Finds the posterior probability cutoff \eqn{\gamma} such that the Bayesian
203203# ' decision rule \code{Pr(theta_t - theta_c > M | data) >= gamma} achieves a
204204# ' target Type-I error \code{alpha} under the least-favourable null, using
205- # ' conjugate Beta– Binomial updates (no MCMC). A bisection search is performed
205+ # ' conjugate Beta- Binomial updates (no MCMC). A bisection search is performed
206206# ' over \code{[lower, upper]} with caching and common random numbers.
207207# '
208208# ' @param alpha Target Type-I error in `(0,1)`, e.g., `0.10`.
@@ -216,7 +216,7 @@ bcts_type1_betaBinom_conj <- function(B = 2000, p_c, M, n_c, n_t,
216216# ' @param lower,upper Initial bracketing interval for `gamma`, values in `(0,1)`.
217217# ' @param n_draws Posterior Monte Carlo draws per trial.
218218# ' @param tol Absolute tolerance for `|TypeI - alpha|`. If `NULL`, a Monte
219- # ' Carlo– driven default is used.
219+ # ' Carlo- driven default is used.
220220# ' @param maxit Maximum number of bisection iterations.
221221# ' @param seed RNG seed for generating common random numbers across thresholds.
222222# ' @param show_progress Show a progress bar inside Type-I evaluations.
@@ -227,9 +227,9 @@ bcts_type1_betaBinom_conj <- function(B = 2000, p_c, M, n_c, n_t,
227227# ' `incProgress(1/maxit, detail = sprintf("Iteration %d of %d", iter, maxit))`.
228228# ' @param calibrate_on Character string indicating which Type-I error quantity
229229# ' to match during calibration. One of:
230- # ' * `"point"` — use the Monte Carlo **point estimate** of the Type-I error
231- # ' * `"upper"` — use the **upper bound** of the 95% Monte Carlo CI
232- # ' * `"lower"` — use the **lower bound** of the 95% Monte Carlo CI
230+ # ' * `"point"` - use the Monte Carlo **point estimate** of the Type-I error
231+ # ' * `"upper"` - use the **upper bound** of the 95% Monte Carlo CI
232+ # ' * `"lower"` - use the **lower bound** of the 95% Monte Carlo CI
233233# ' Default is `"point"`.
234234# '
235235# ' @return A list with:
@@ -410,11 +410,14 @@ bcts_calibrate_betaBinom_conj <- function(alpha = 0.10, p_c, M, n_c, n_t,
410410
411411 if (stop_now ) {
412412 if (! is.null(progress_fun )) progress_fun(maxit , maxit ) # <- force 100%
413- return (list (
413+ return (structure( list (
414414 gamma = mid ,
415- type1 = r_mid $ estimate ,
416- mc_se = r_mid $ mc_se ,
417- mc_ci = r_mid $ mc_ci ,
415+ type1 = c(
416+ estimate = r_mid $ estimate ,
417+ mc_se = r_mid $ mc_se ,
418+ ci_lower = r_mid $ ci_lower ,
419+ ci_upper = r_mid $ ci_upper
420+ ),
418421 bracket = c(lower = lo , upper = hi ),
419422 iters = i ,
420423 B_cal = B_cal ,
@@ -427,7 +430,7 @@ bcts_calibrate_betaBinom_conj <- function(alpha = 0.10, p_c, M, n_c, n_t,
427430 n_draws = n_draws , tol = tol , maxit = maxit ,
428431 seed = seed , calibrate_on = calibrate_on
429432 )
430- ))
433+ ), class = " bcts_calibration " ) )
431434 }
432435 if (diff > 0 ) lo <- mid else hi <- mid
433436 }
@@ -445,11 +448,14 @@ bcts_calibrate_betaBinom_conj <- function(alpha = 0.10, p_c, M, n_c, n_t,
445448 type1 = r_mid $ estimate , ci_lower = r_mid $ ci_lower , ci_upper = r_mid $ ci_upper ,
446449 metric = m_mid , lo = lo , hi = hi , diff = m_mid - alpha ))
447450
448- list (
451+ structure( list (
449452 gamma = mid ,
450- type1 = r_mid $ estimate ,
451- mc_se = r_mid $ mc_se ,
452- mc_ci = r_mid $ mc_ci ,
453+ type1 = c(
454+ estimate = r_mid $ estimate ,
455+ mc_se = r_mid $ mc_se ,
456+ ci_lower = r_mid $ ci_lower ,
457+ ci_upper = r_mid $ ci_upper
458+ ),
453459 bracket = c(lower = lo , upper = hi ),
454460 iters = maxit ,
455461 B_cal = B_cal ,
@@ -462,7 +468,7 @@ bcts_calibrate_betaBinom_conj <- function(alpha = 0.10, p_c, M, n_c, n_t,
462468 n_draws = n_draws , tol = tol , maxit = maxit ,
463469 seed = seed , calibrate_on = calibrate_on
464470 )
465- )
471+ ), class = " bcts_calibration " )
466472}
467473
468474# ' Validate B_cal parameter
@@ -539,3 +545,68 @@ check_alpha_bracketing <- function(calibrate_on, alpha, lower, upper, r_lo, r_hi
539545
540546 invisible (TRUE )
541547}
548+
549+ # ' Plot calibration trace for calibrated Bayesian NI threshold
550+ # '
551+ # ' Visualizes the estimated Type-I error at different posterior thresholds
552+ # ' during the calibration procedure.
553+ # '
554+ # ' @param x Output from [bcts_calibrate_betaBinom_conj()]
555+ # ' @param ... Not used.
556+ # '
557+ # ' @return A [ggplot2::ggplot()] object.
558+ # ' @export
559+ # '
560+ # ' @examples
561+ # ' # res <- bcts_calibrate_betaBinom_conj(...) # Run separately
562+ # ' # plot(res)
563+ plot.bcts_calibration <- function (x , ... ) {
564+ stopifnot(! is.null(x $ trace ), ! is.null(x $ type1 ))
565+
566+ tr <- x $ trace
567+ alpha <- x $ alpha
568+ gamma <- x $ gamma
569+ type1 <- x $ type1
570+
571+ # Small vertical offset above CI
572+ offset <- 0.01
573+
574+ ggplot2 :: ggplot(tr , ggplot2 :: aes(x = .data $ gamma_try , y = .data $ type1 )) +
575+ ggplot2 :: geom_errorbar(
576+ ggplot2 :: aes(ymin = .data $ ci_lower , ymax = .data $ ci_upper ),
577+ width = 0.002 , alpha = 0.4
578+ ) +
579+ ggplot2 :: geom_point(size = 2 ) +
580+ ggplot2 :: geom_line() +
581+ ggplot2 :: geom_point(
582+ x = gamma ,
583+ y = type1 [" estimate" ],
584+ color = " blue" ,
585+ size = 3
586+ ) +
587+ ggplot2 :: annotate(
588+ " text" ,
589+ x = gamma ,
590+ y = type1 [" ci_upper" ] + offset ,
591+ label = sprintf(" gamma = %.3f\n Type-I = %.1f%%" ,
592+ gamma , 100 * type1 [" estimate" ]),
593+ hjust = 0.5 ,
594+ size = 3.5
595+ ) +
596+ ggplot2 :: geom_hline(
597+ yintercept = alpha ,
598+ linetype = " dotted" ,
599+ color = " red"
600+ ) +
601+ ggplot2 :: labs(
602+ x = expression(gamma ),
603+ y = " Estimated Type-I error" ,
604+ subtitle = sprintf(
605+ " Calibrated gamma = %.3f after %d steps; dotted = target alpha = %.2f" ,
606+ gamma , x $ iters , alpha
607+ )
608+ ) +
609+ ggplot2 :: theme_minimal(base_size = 12 )+
610+ ggplot2 :: scale_y_continuous(labels = scales :: label_percent(accuracy = 1 ))+
611+ ggplot2 :: scale_x_continuous(labels = scales :: label_percent(accuracy = 1 ))
612+ }
0 commit comments