Extend normal outcome functions to two-sided decisions

.Rbuildignore

@@ -37,3 +37,6 @@ src/stan_files/.*\.hpp$
 ^Makefile$
 ^LICENSE$
 ^vignettes/articles$
+^\.vscode$
+^[.]?air[.]toml$
+^\.lintr$

.gitignore

@@ -3,10 +3,9 @@
 .RData
 .Ruserdata
 docs
-src/*.so
-src/*.o
-src/stan_files/*.hpp
-src/stan_files/*.o
+src/
 *.*~
 inst/sbc/sbc*html
 .vscode/
+*.dll
+.Renviron

.lintr

@@ -0,0 +1,5 @@
+linters: linters_with_defaults(
+    line_length_linter(120),
+    object_name_linter = NULL,
+    commented_code_linter = NULL
+  )

NAMESPACE

@@ -75,7 +75,9 @@ S3method(print,EMnmm)
 S3method(print,betaBinomialMix)
 S3method(print,betaMix)
 S3method(print,decision1S)
+S3method(print,decision1S_2sided)
 S3method(print,decision2S)
+S3method(print,decision2S_2sided)
 S3method(print,dlink)
 S3method(print,gMAP)
 S3method(print,gMAPpred)
@@ -135,6 +137,7 @@ export(gMAP)
 export(inv_logit)
 export(likelihood)
 export(logit)
+export(lower)
 export(mixbeta)
 export(mixcombine)
 export(mixfit)
@@ -165,6 +168,7 @@ export(rmix)
 export(rmixdiff)
 export(robustify)
 export(sigma)
+export(upper)
 export(write_mix_json)
 import(Formula)
 import(Rcpp)

NEWS.md

@@ -3,6 +3,7 @@
 ## Enhancements
 
 * Reformat R sources using `Air`.
+* Extend normal outcome functions to two-sided decisions: `decision1S` and `decision2S` now allow `lower.tail` to have as many elements as `pc` to allow the specification of two-sided decision boundaries (mixed `lower.tail` elements, i.e. some specifying "lower" and some "upper") to capture intermediate result scenarios. This is currently only respected by the normal outcome functions (methods for `normMix`), because it is work in progress.
 
 # RBesT 1.8-2 - April 25th, 2025
 

R/decision1S.R

@@ -6,9 +6,12 @@
 #' @param pc Vector of critical cumulative probabilities.
 #' @param qc Vector of respective critical values. Must match the length of `pc`.
 #' @param lower.tail Logical; if `TRUE` (default), probabilities
-#' are \eqn{P(X \leq x)}, otherwise, \eqn{P(X > x)}.
+#' are \eqn{P(X \leq x)}, otherwise, \eqn{P(X > x)}. Either length 1 or same
+#' length as `pc`.
+#' @param x Two-sided decision function.
 #'
-#' @details The function creates a one-sided decision function which
+#' @details For `lower.tail` being either `TRUE` or `FALSE`,
+#' the function creates a one-sided decision function which
 #' takes two arguments. The first argument is expected to be a mixture
 #' (posterior) distribution. This distribution is tested whether it
 #' fulfills all the required threshold conditions specified with the
@@ -25,6 +28,10 @@
 #'
 #' \deqn{\Pi_i H_i(P(\theta \leq q_{c,i}) - p_{c,i} ).}
 #'
+#' For the case of a boolen vector given to `lower.tail` the
+#' direction of each decision aligns respectively, and a two-sided
+#' decision function is created.
+#'
 #' When the second argument is set to `TRUE` a distance metric is
 #' returned component-wise per defined condition as
 #'
@@ -36,7 +43,9 @@
 #'
 #' @family design1S
 #'
-#' @return The function returns a decision function which takes two
+#' @return The function returns a decision function (of class
+#' `decision1S` for one-sided, and of class `decision1S_2sided`
+#' for two-sided decisions) which takes two
 #' arguments. The first argument is expected to be a mixture
 #' (posterior) distribution which is tested if the specified
 #' conditions are met. The logical second argument determines if the
@@ -45,6 +54,11 @@
 #' log-space, i.e. the distance is 0 at the decision boundary,
 #' negative for a 0 decision and positive for a 1 decision.
 #'
+#' For two-sided decision functions, the two components can be
+#' extracted with functions [lower()] and [upper()]. The distance
+#' as calculated by the decision function is returned as a list with
+#' components `lower` and `upper`.
+#'
 #' @references Neuenschwander B, Rouyrre N, Hollaender H, Zuber E,
 #' Branson M. A proof of concept phase II non-inferiority
 #' criterion. *Stat. in Med.*. 2011, 30:1618-1627
@@ -53,7 +67,7 @@
 #'
 #' # see Neuenschwander et al., 2011
 #'
-#' # example is for a time-to-event trial evaluating non-inferiority
+#' # example is for a time-to-event trial evaluating non-inferiority (NI)
 #' # using a normal approximation for the log-hazard ratio
 #'
 #' # reference scale
@@ -99,9 +113,52 @@
 #' # here with HR of 0.8 for 40 events
 #' decComb(postmix(flat_prior, m = log(0.8), n = 40))
 #'
+#' # A two-sided decision function can be useful to determine if
+#' # certain intermediate (i.e. neither "go" nor "stop") decisions
+#' # are to be made based on the posterior distribution.
+#' # For example, in the above situation we might have an intermediate
+#' # scenario where the trial is significant for non-inferiority but
+#' # the mean estimate is in an intermediate range, say between theta_c
+#' # theta_f:
+#' theta_f <- 0.3
+#' decCombIntermediate <- decision1S(
+#'   c(1 - alpha, 0.5, 0.8),
+#'   c(theta_ni, theta_c, theta_f),
+#'   lower.tail = c(TRUE, FALSE, TRUE)
+#' )
+#' # Not fulfilled for the prior:
+#' decCombIntermediate(flat_prior)
+#' # But for a hypothetical trial outcome with HR 1.2 and 300 events:
+#' decCombIntermediate(postmix(flat_prior, m = log(1.2), n = 300))
+#'
 #' @export
 decision1S <- function(pc = 0.975, qc = 0, lower.tail = TRUE) {
-  assert_that(length(pc) == length(qc))
+  assert_numeric(pc)
+  assert_numeric(qc, len = length(pc))
+  assert_logical(lower.tail)
+  assert_true(length(lower.tail) == 1L || length(lower.tail) == length(pc))
+  lower.tail <- scalar_if_same(lower.tail)
+
+  is_two_sided <- length(lower.tail) > 1
+
+  if (is_two_sided) {
+    return(create_decision1S_2sided(pc, qc, lower.tail))
+  } else {
+    return(create_decision1S_1sided(pc, qc, lower.tail))
+  }
+}
+
+#' @keywords internal
+scalar_if_same <- function(x) {
+  if (length(x) > 1 && all(x == x[1])) {
+    return(x[1])
+  }
+  x
+}
+
+#' Internal Constructor for 1 Sample One-sided Decision Function
+#' @keywords internal
+create_decision1S_1sided <- function(pc, qc, lower.tail) {
   lpc <- log(pc)
   fun <- function(mix, dist = FALSE) {
     test <- pmix(mix, qc, lower.tail = lower.tail, log.p = TRUE) - lpc
@@ -113,21 +170,75 @@ decision1S <- function(pc = 0.975, qc = 0, lower.tail = TRUE) {
   attr(fun, "pc") <- pc
   attr(fun, "qc") <- qc
   attr(fun, "lower.tail") <- lower.tail
+
   class(fun) <- c("decision1S", "function")
   fun
 }
 
+#' Internal Constructor for 1 Sample Two-sided Decision Function
+#' @keywords internal
+create_decision1S_2sided <- function(pc, qc, lower.tail) {
+  use_lower <- which(lower.tail)
+  use_upper <- which(!lower.tail)
+  assert_true(length(use_lower) > 0 && length(use_upper) > 0)
+
+  lower_part <- create_decision1S_1sided(pc[use_lower], qc[use_lower], TRUE)
+  upper_part <- create_decision1S_1sided(pc[use_upper], qc[use_upper], FALSE)
+
+  fun <- function(mix, dist = FALSE) {
+    dl <- lower_part(mix, dist)
+    du <- upper_part(mix, dist)
+    if (dist) {
+      return(list(lower = dl, upper = du))
+    }
+    as.numeric(all(dl > 0) && all(du > 0))
+  }
+  attr(fun, "lower") <- lower_part
+  attr(fun, "upper") <- upper_part
+
+  class(fun) <- c("decision1S_2sided", "function")
+  fun
+}
+
+#' @rdname decision1S
 #' @export
-print.decision1S <- function(x, ...) {
-  cat("1 sample decision function\n")
-  cat("Conditions for acceptance:\n")
+lower <- function(x) {
+  assert_multi_class(x, c("decision1S_2sided", "decision2S_2sided"))
+  attr(x, "lower")
+}
+
+#' @rdname decision1S
+#' @export
+upper <- function(x) {
+  assert_multi_class(x, c("decision1S_2sided", "decision2S_2sided"))
+  attr(x, "upper")
+}
+
+#' @keywords internal
+print_decision1S_1sided <- function(x) {
   qc <- attr(x, "qc")
   pc <- attr(x, "pc")
   low <- attr(x, "lower.tail")
   cmp <- ifelse(low, "<=", ">")
-  for (i in seq_along(qc)) {
-    cat(paste0("P(theta ", cmp, " ", qc[i], ") > ", pc[i], "\n"))
-  }
+  cat(paste0("P(theta ", cmp, " ", qc, ") > ", pc, "\n"), sep = "")
+}
+
+#' @export
+print.decision1S <- function(x, ...) {
+  cat("1 sample decision function\n")
+  cat("Conditions for acceptance:\n")
+  print_decision1S_1sided(x)
+  invisible(x)
+}
+
+#' @export
+print.decision1S_2sided <- function(x, ...) {
+  cat("1 sample decision function (two-sided)\n")
+  cat("Conditions for acceptance:\n")
+  cat("Lower tail conditions:\n")
+  print_decision1S_1sided(lower(x))
+  cat("Upper tail conditions:\n")
+  print_decision1S_1sided(upper(x))
   invisible(x)
 }
 

R/decision1S_boundary.R

@@ -31,7 +31,10 @@
 #' aligned with the cumulative density function definition within R
 #' (see for example [pbinom()]).
 #'
-#' @return Returns the critical value \eqn{y_c}.
+#' @return Returns the critical value \eqn{y_c}. For two-sided
+#'   decision functions a named vector with components
+#'   `lower_or_equal_than` and `higher_than` is returned, containing
+#'   the critical values for the lower and upper decision boundaries.
 #'
 #' @family design1S
 #'
@@ -157,13 +160,67 @@ decision1S_boundary.normMix <- function(
   eps = 1e-6,
   ...
 ) {
-  ## distributions of the means of the data generating distributions
-  ## for now we assume that the underlying standard deviation
-  ## matches the respective reference scales
+  # Get the default sigma if not provided already here to only message once.
   if (missing(sigma)) {
     sigma <- RBesT::sigma(prior)
     message("Using default prior reference scale ", sigma)
   }
+
+  if (is(decision, "decision1S_2sided")) {
+    decision1S_boundary_normMix_2sided(
+      prior,
+      n,
+      decision,
+      sigma,
+      eps
+    )
+  } else {
+    decision1S_boundary_normMix_1sided(
+      prior,
+      n,
+      decision,
+      sigma,
+      eps
+    )
+  }
+}
+
+#' @keywords internal
+decision1S_boundary_normMix_2sided <- function(
+  prior,
+  n,
+  decision,
+  sigma,
+  eps
+) {
+  crit_lower <- decision1S_boundary_normMix_1sided(
+    prior,
+    n,
+    lower(decision),
+    sigma,
+    eps
+  )
+  crit_upper <- decision1S_boundary_normMix_1sided(
+    prior,
+    n,
+    upper(decision),
+    sigma,
+    eps
+  )
+  c(lower_or_equal_than = crit_lower, higher_than = crit_upper)
+}
+
+#' @keywords internal
+decision1S_boundary_normMix_1sided <- function(
+  prior,
+  n,
+  decision,
+  sigma,
+  eps
+) {
+  ## distributions of the means of the data generating distributions
+  ## for now we assume that the underlying standard deviation
+  ## matches the respective reference scales
   assert_number(sigma, lower = 0)
 
   sd_samp <- sigma / sqrt(n)
@@ -172,9 +229,6 @@ decision1S_boundary.normMix <- function(
 
   ## change the reference scale of the prior such that the prior
   ## represents the distribution of the respective means
-  ## mean_prior <- prior
-  ## sigma(mean_prior) <- sd_samp
-
   m <- summary(prior, probs = c())["mean"]
 
   lim <- qnorm(p = c(eps / 2, 1 - eps / 2), mean = m, sd = sd_samp)
@@ -185,7 +239,6 @@ decision1S_boundary.normMix <- function(
   crit
 }
 
-
 #' @templateVar fun decision1S_boundary
 #' @template design1S-poisson
 #' @export