Skip to content

Commit dc7c71c

Browse files
committed
Improve alpha control
1 parent aa2539f commit dc7c71c

File tree

8 files changed

+199
-30
lines changed

8 files changed

+199
-30
lines changed

R/bayesNI-conj-calibrate.R

Lines changed: 103 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -288,9 +288,9 @@ bcts_calibrate_betaBinom_conj <- function(alpha = 0.10, p_c, M, n_c, n_t,
288288
prior_args = list(),
289289
B_cal = 2000, lower = 0.80, upper = 0.999,
290290
n_draws = 2000,
291-
tol = 0.001, maxit = 30,
291+
tol = NULL, maxit = 30,
292292
seed = 11L, show_progress = TRUE,
293-
verbose = TRUE, digits = 4,
293+
verbose = FALSE, digits = 4,
294294
progress_fun = NULL,
295295
calibrate_on = c("point","upper","lower")) {
296296

@@ -304,12 +304,23 @@ bcts_calibrate_betaBinom_conj <- function(alpha = 0.10, p_c, M, n_c, n_t,
304304
if (p_c < 0 || p_c > 1) stop("`p_c` must be in [0,1].")
305305
if (n_c < 1 || n_t < 1) stop("`n_c` and `n_t` must be positive integers.")
306306

307+
# --- Check B_cal ---
308+
check_B_cal(B_cal, context = "bcts_calibrate_betaBinom_conj")
309+
307310
set.seed(seed)
308311
seeds <- sample.int(1e9, B_cal) # common random numbers
309312

313+
# Monte Carlo standard error of a binomial proportion
314+
mcse <- sqrt(alpha * (1 - alpha) / B_cal)
315+
310316
if (is.null(tol)) {
311-
tol <- 1.25 * sqrt(alpha * (1 - alpha) / B_cal)
317+
tol <- 1.25 * mcse
312318
if (verbose) message(sprintf("Auto tol set to %.4f based on B_cal=%d", tol, B_cal))
319+
} else if (tol < 0.5 * mcse) {
320+
warning(sprintf(
321+
"Specified `tol` = %.4f is too narrow relative to Monte Carlo SE = %.4f (based on B_cal = %d). Consider increasing `tol` or B_cal.",
322+
tol, mcse, B_cal
323+
))
313324
}
314325

315326
# cache full results per gamma (so we can reuse MC error, etc.)
@@ -339,19 +350,20 @@ bcts_calibrate_betaBinom_conj <- function(alpha = 0.10, p_c, M, n_c, n_t,
339350
)
340351
}
341352

353+
342354
# initial bracket diagnostics (must bracket alpha on the chosen metric)
343355
r_lo <- type1_at(lower); m_lo <- get_metric(r_lo)
344356
r_hi <- type1_at(upper); m_hi <- get_metric(r_hi)
357+
358+
# One-sided check: ensure that at least one value is below alpha
359+
check_alpha_bracketing(calibrate_on, alpha, lower, upper, r_lo, r_hi, verbose = verbose)
360+
361+
# Optional: verbose printout
345362
if (verbose) {
346363
message(sprintf(
347364
"Init: lower=%.3f -> %s=%.4f | upper=%.3f -> %s=%.4f | target alpha=%.4f",
348365
lower, calibrate_on, m_lo, upper, calibrate_on, m_hi, alpha))
349366
}
350-
if (!(m_lo >= alpha && m_hi <= alpha)) {
351-
warning(sprintf(
352-
"Alpha not bracketed for '%s': metric(lower)=%.4f, metric(upper)=%.4f",
353-
calibrate_on, m_lo, m_hi))
354-
}
355367

356368

357369
lo <- lower; hi <- upper
@@ -389,7 +401,14 @@ bcts_calibrate_betaBinom_conj <- function(alpha = 0.10, p_c, M, n_c, n_t,
389401
"Iter %02d: gamma=%.4f -> %s=%.4f (diff=%.4f)",
390402
i, mid, calibrate_on, m_mid, diff))
391403

392-
if (abs(diff) < tol) {
404+
stop_now <- switch(
405+
calibrate_on,
406+
point = abs(m_mid - alpha) < tol,
407+
upper = (r_mid$ci_upper <= alpha && r_mid$ci_upper > (alpha - tol)),
408+
lower = (r_mid$ci_lower >= alpha && r_mid$ci_lower < (alpha + tol))
409+
)
410+
411+
if (stop_now) {
393412
if (!is.null(progress_fun)) progress_fun(maxit, maxit) # <- force 100%
394413
return(list(
395414
gamma = mid,
@@ -445,3 +464,78 @@ bcts_calibrate_betaBinom_conj <- function(alpha = 0.10, p_c, M, n_c, n_t,
445464
)
446465
)
447466
}
467+
468+
#' Validate B_cal parameter
469+
#'
470+
#' @param B_cal Number of simulations
471+
#' @param min_reliable Minimum value considered reliable (default = 100)
472+
#' @param context Optional name of the calling function for better messages
473+
#' @return An integer-validated B_cal value (invisibly)
474+
check_B_cal <- function(B_cal, min_reliable = 100, context = NULL) {
475+
name <- if (!is.null(context)) paste0(" in ", context) else ""
476+
477+
if (!is.numeric(B_cal) || length(B_cal) != 1 || B_cal <= 0 || B_cal != as.integer(B_cal)) {
478+
stop(sprintf("`B_cal`%s must be a positive integer.", name), call. = FALSE)
479+
}
480+
481+
if (B_cal < min_reliable) {
482+
warning(sprintf("`B_cal`%s is very small; Monte Carlo estimates may be unreliable.", name), call. = FALSE)
483+
}
484+
485+
invisible(as.integer(B_cal))
486+
}
487+
488+
#' Check alpha bracketing for calibration of Type-I error
489+
#'
490+
#' @param calibrate_on One of "point", "upper", or "lower"
491+
#' @param alpha Target Type-I error
492+
#' @param lower,upper Gamma bounds
493+
#' @param r_lo Result of type1_at(lower)
494+
#' @param r_hi Result of type1_at(upper)
495+
#' @param verbose Print message
496+
#'
497+
#' @return Invisibly TRUE
498+
check_alpha_bracketing <- function(calibrate_on, alpha, lower, upper, r_lo, r_hi, verbose = TRUE) {
499+
if (calibrate_on == "point") {
500+
est_lo <- r_lo$estimate
501+
est_hi <- r_hi$estimate
502+
if (!(alpha >= est_hi && alpha <= est_lo)) {
503+
warning(sprintf(
504+
paste0(
505+
"Cannot bracket alpha = %.3f using point estimate.\n",
506+
"At gamma = %.3f: estimate = %.4f | at gamma = %.3f: estimate = %.4f.\n",
507+
"Try adjusting `lower`, `upper`, or `B_cal`."
508+
), alpha, lower, est_lo, upper, est_hi
509+
))
510+
}
511+
} else if (calibrate_on == "upper") {
512+
if (r_hi$ci_upper > alpha) {
513+
warning(sprintf(
514+
paste0(
515+
"No gamma value found such that the *upper bound* of Type-I error is <= alpha = %.3f.\n",
516+
"At gamma = %.3f: upper CI = %.4f.\n",
517+
"Try increasing `upper`, decreasing `lower`, or increasing `B_cal`."
518+
), alpha, upper, r_hi$ci_upper
519+
))
520+
}
521+
} else if (calibrate_on == "lower") {
522+
if (r_hi$ci_lower < alpha) {
523+
warning(sprintf(
524+
paste0(
525+
"No gamma value found such that the *lower bound* of Type-I error is >= alpha = %.3f.\n",
526+
"At gamma = %.3f: lower CI = %.4f.\n",
527+
"Try increasing `upper`, decreasing `lower`, or increasing `B_cal`."
528+
), alpha, upper, r_hi$ci_lower
529+
))
530+
}
531+
}
532+
533+
if (verbose) {
534+
message(sprintf(
535+
"Bracket check: lower=%.3f | upper=%.3f | alpha=%.3f (%s)",
536+
lower, upper, alpha, calibrate_on
537+
))
538+
}
539+
540+
invisible(TRUE)
541+
}

inst/app/app.R

Lines changed: 30 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -4,23 +4,31 @@
44
# ---- Load packages that the APP needs ----
55
library(shiny)
66
library(ggplot2)
7-
library(bcts) # <- your package is built & installed during deploy
8-
9-
# ---- Locate the installed app directory (inside the bcts package) ----
10-
app_dir <- system.file("app", package = "bcts")
11-
stopifnot(nzchar(app_dir), dir.exists(app_dir))
12-
13-
# Optional: serve static assets from inst/app/www as /www
14-
www_dir <- file.path(app_dir, "www")
15-
if (dir.exists(www_dir)) addResourcePath("www", www_dir)
7+
library(bcts)
168

179
# Helper to source R files from inst/app/R and inst/app/modules
1810
source_dir <- function(path) {
1911
if (dir.exists(path)) {
2012
fs <- list.files(path, pattern = "\\.R$", full.names = TRUE)
21-
for (f in fs) source(f, local = TRUE)
13+
for (f in fs) {
14+
message("Sourcing: ", basename(f))
15+
tryCatch({
16+
source(f, local = globalenv()) # Use globalenv for visibility
17+
}, error = function(e) {
18+
message("Error sourcing ", f, ": ", e$message)
19+
})
20+
}
2221
}
2322
}
23+
24+
USE_LOCAL <- interactive() && dir.exists(file.path("inst", "app")) # Only TRUE when running locally
25+
26+
app_dir <- if (USE_LOCAL) {
27+
file.path("inst", "app")
28+
} else {
29+
"."
30+
}
31+
2432
source_dir(file.path(app_dir, "R"))
2533
source_dir(file.path(app_dir, "modules"))
2634

@@ -387,10 +395,17 @@ server <- function(input, output, session) {
387395
pc_current = reactive(input$pc) # for pre-filling the range nicely
388396
)
389397

390-
391-
392-
393-
394398
}
395399

396-
shinyApp(ui, server)
400+
tryCatch(
401+
{
402+
shiny::shinyApp(
403+
ui = ui,
404+
server = server
405+
)
406+
},
407+
error = function(e) {
408+
message("🚨 App failed to launch: ", conditionMessage(e))
409+
stop(e)
410+
}
411+
)

inst/app/modules/mod_armpriors.R

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,9 @@
11
# modules/mod_armpriors.R
22
# Prior & posterior Beta densities per arm (faceted)
33

4+
message("Loaded mod_armpriors.R")
5+
6+
47
mod_armpriors_ui <- function(id, height = 320) {
58
ns <- NS(id)
69
tagList(

inst/app/modules/mod_designsummary.R

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,5 @@
11
# modules/mod_designsummary.R
2+
message("Loaded mod_designsummary.R")
23

34
mod_designsummary_ui <- function(id) {
45
ns <- NS(id)

inst/app/modules/mod_narrative.R

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,5 @@
1+
# modules/mod_narrative_ui
2+
message("Loaded mod_narrative_ui.R")
13

24
# ---- UI ----
35
mod_narrative_ui <- function(id) {
@@ -75,13 +77,11 @@ mod_narrative_server <- function(
7577
"<p>%s</p>
7678
<p>%s</p>
7779
<p>The analysis uses %s.</p>
78-
<p>Posterior probabilities are computed using %s draws per simulated trial, across %s simulated trials.</p>
79-
<p>Seed = %s (for reproducibility).</p>",
80+
<p>Posterior probabilities are computed using %s draws per simulated trial, across %s simulated trials.</p>",
8081
design_txt,
8182
arm_txt,
8283
prior_type,
83-
fmt_int(ndraws()), fmt_int(B()),
84-
fmt_int(seed())
84+
fmt_int(ndraws()), fmt_int(B())
8585
))
8686
})
8787
})

man/bcts_calibrate_betaBinom_conj.Rd

Lines changed: 2 additions & 2 deletions
Some generated files are not rendered by default. Learn more about customizing how changed files appear on GitHub.

man/check_B_cal.Rd

Lines changed: 21 additions & 0 deletions
Some generated files are not rendered by default. Learn more about customizing how changed files appear on GitHub.

man/check_alpha_bracketing.Rd

Lines changed: 35 additions & 0 deletions
Some generated files are not rendered by default. Learn more about customizing how changed files appear on GitHub.

0 commit comments

Comments
 (0)