updated spending functions to include WT

Author: xidongdxi

NAMESPACE

@@ -56,5 +56,6 @@ export(spending_linear)
 export(spending_of)
 export(spending_pocock)
 export(spending_with_time)
+export(spending_wt)
 export(three_doses_two_primary_two_secondary)
 export(two_doses_two_primary_two_secondary)

R/graph_test_shortcut_gsd.R

@@ -44,8 +44,11 @@
 #'     `NA` padding, `info_frac` must be a matrix with `NA` in the same
 #'     positions as `p`.
 #'
-#'   Non-`NA` values must be in (0, 1] and monotonically non-decreasing per
-#'   hypothesis. The last non-`NA` value does not need to be 1, allowing
+#'   Non-`NA` values must be positive and monotonically non-decreasing per
+#'   hypothesis. Values greater than 1 are allowed (e.g., when more
+#'   information is collected than planned). The spending functions cap
+#'   the cumulative spending at `alpha` for information fractions at or
+#'   above 1. The last non-`NA` value does not need to be 1, allowing
 #'   the procedure to be applied up to an interim analysis.
 #' @param spending_fn Spending function(s) for computing group sequential
 #'   boundaries. Can be:
@@ -562,7 +565,7 @@ gsd_test <- function(graph, p, alpha, info_frac, spending_fn, look_back,
     if (test_values) {
       tv_details[[k]] <- gsd_test_values_details(
         step_graph, p, k, alpha, info_frac, spending_fn,
-        newly_in_order, hyp_names, rejected, has_data_k
+        newly_in_order, hyp_names, rejected, active_at_k
       )
 
       # Add Look_back column (FALSE for all standard rows)
@@ -587,9 +590,14 @@ gsd_test <- function(graph, p, alpha, info_frac, spending_fn, look_back,
                            tv_details[[k]]$Analysis == k)
 
           if (length(hyp_row) > 0) {
-            # Has data at analysis k: set Reject to FALSE and insert
-            # look_back rows after it
-            tv_details[[k]]$Reject[hyp_row] <- FALSE
+            # Has data at analysis k: check if the nominal p-value at
+            # analysis k also crosses the boundary. If not, set Reject
+            # to FALSE (the rejection is only via look_back).
+            p_at_k <- tv_details[[k]]$p[hyp_row]
+            b_at_k <- tv_details[[k]]$Boundary[hyp_row]
+            if (is.na(p_at_k) || p_at_k > b_at_k) {
+              tv_details[[k]]$Reject[hyp_row] <- FALSE
+            }
             before <- tv_details[[k]][seq_len(hyp_row), , drop = FALSE]
             after <- if (hyp_row < nrow(tv_details[[k]])) {
               tv_details[[k]][(hyp_row + 1):nrow(tv_details[[k]]), ,
@@ -822,8 +830,8 @@ gsd_input_val <- function(graph, p, alpha, info_frac, spending_fn, look_back,
     "Information fractions must have the same number of columns as p" =
       ncol(info_frac) == num_analyses,
     "Information fractions must be numeric" = is.numeric(info_frac),
-    "Non-NA information fractions must be in (0, 1]" =
-      length(if_non_na) == 0 || all(if_non_na > 0 & if_non_na <= 1),
+    "Non-NA information fractions must be positive" =
+      length(if_non_na) == 0 || all(if_non_na > 0),
     "Spending functions must be a list of functions" =
       is.list(spending_fn) &&
       all(vapply(spending_fn, is.function, logical(1))),

R/print.gsd_graph_report.R

@@ -145,10 +145,15 @@ print.gsd_graph_report <- function(x, ..., precision = 6, indent = 2) {
     check.names = FALSE
   )
   names(df_summary)[[1]] <- formatC("Hypothesis", width = hyp_width)
-  names(df_summary)[[2]] <- "Adj.P-value"
+  names(df_summary)[[2]] <- "Adj.P-value*"
 
   print(df_summary, row.names = FALSE)
 
+  cat(pad, "(*) Adjusted p-values account for both the group sequential",
+      " design and the\n", pad, "    graphical multiple comparison procedure.",
+      " Based on repeated p-values when\n", pad, "    look_back = FALSE,",
+      " and sequential p-values when look_back = TRUE.\n", sep = "")
+
   # Rejection sequence
   rej_seq <- x$outputs$rejection_sequence
   if (length(rej_seq) > 0) {
@@ -201,9 +206,10 @@ print.gsd_graph_report <- function(x, ..., precision = 6, indent = 2) {
 
       # Print footnote for look_back hypotheses
       if (has_look_back) {
-        cat(pad, "(*) Rejected via look_back: nominal p-value did not cross",
-            " the boundary at the\n", pad, "    current analysis, but",
-            " crossed the boundary at an earlier analysis.\n", sep = "")
+        cat(pad, "(*) Rejected via look_back: the nominal p-value crossed",
+            " the boundary at an\n", pad, "    earlier analysis with the",
+            " hypothesis weight updated via graph propagation.\n",
+            sep = "")
       }
       cat("\n")
     }

R/spending_functions.R

@@ -16,8 +16,8 @@
 #' @param alpha A numeric scalar of the total significance level to be spent.
 #'   Must be between 0 and 1.
 #' @param info_frac A numeric scalar or vector of information fractions. Values
-#'   must be in \[0, 1\]. When `info_frac = 0`, the spending is 0. When
-#'   `info_frac = 1`, the spending equals `alpha`.
+#'   must be non-negative. When `info_frac = 0`, the spending is 0. When
+#'   `info_frac >= 1`, the spending is capped at `alpha`.
 #' @param gamma A numeric scalar for the gamma parameter of the
 #'   Hwang-Shih-DeCani spending function. Common choices are `gamma = -4`
 #'   (approximates O'Brien-Fleming), `gamma = 1` (approximates Pocock), and
@@ -89,10 +89,12 @@
 #' spending_linear(0.025, c(1/3, 2/3, 1))
 spending_of <- function(alpha, info_frac) {
   stopifnot(
-    "info_frac must be in [0, 1]" = all(info_frac >= 0 & info_frac <= 1)
+    "info_frac must be non-negative" = all(info_frac >= 0),
+    "At most one info_frac value can be >= 1" = sum(info_frac >= 1) <= 1
   )
   result <- 2 * (1 - stats::pnorm(stats::qnorm(1 - alpha / 2) / sqrt(info_frac)))
   result[info_frac == 0] <- 0
+  result <- pmin(result, alpha)
   result
 }
 
@@ -103,10 +105,12 @@ spending_of <- function(alpha, info_frac) {
 #' spending_pocock(0.025, 0.5)
 spending_pocock <- function(alpha, info_frac) {
   stopifnot(
-    "info_frac must be in [0, 1]" = all(info_frac >= 0 & info_frac <= 1)
+    "info_frac must be non-negative" = all(info_frac >= 0),
+    "At most one info_frac value can be >= 1" = sum(info_frac >= 1) <= 1
   )
   result <- alpha * log(1 + (exp(1) - 1) * info_frac)
   result[info_frac == 0] <- 0
+  result <- pmin(result, alpha)
   result
 }
 
@@ -119,14 +123,16 @@ spending_pocock <- function(alpha, info_frac) {
 #' spending_hsd(0.025, 0.5, gamma = 0)
 spending_hsd <- function(alpha, info_frac, gamma = -4) {
   stopifnot(
-    "info_frac must be in [0, 1]" = all(info_frac >= 0 & info_frac <= 1)
+    "info_frac must be non-negative" = all(info_frac >= 0),
+    "At most one info_frac value can be >= 1" = sum(info_frac >= 1) <= 1
   )
   if (gamma == 0) {
     result <- alpha * info_frac
   } else {
     result <- alpha * (1 - exp(-gamma * info_frac)) / (1 - exp(-gamma))
   }
   result[info_frac == 0] <- 0
+  result <- pmin(result, alpha)
   result
 }
 
@@ -137,9 +143,10 @@ spending_hsd <- function(alpha, info_frac, gamma = -4) {
 #' spending_linear(0.025, 0.5)
 spending_linear <- function(alpha, info_frac) {
   stopifnot(
-    "info_frac must be in [0, 1]" = all(info_frac >= 0 & info_frac <= 1)
+    "info_frac must be non-negative" = all(info_frac >= 0),
+    "At most one info_frac value can be >= 1" = sum(info_frac >= 1) <= 1
   )
-  alpha * info_frac
+  pmin(alpha * info_frac, alpha)
 }
 
 
@@ -207,12 +214,163 @@ spending_with_time <- function(spending_fn, spending_time) {
   stopifnot(
     "spending_fn must be a function" = is.function(spending_fn),
     "spending_time must be a numeric vector" = is.numeric(spending_time),
-    "spending_time must be in [0, 1]" =
-      all(spending_time >= 0 & spending_time <= 1)
+    "spending_time must be non-negative" =
+      all(spending_time >= 0),
+    "At most one spending_time value can be >= 1" =
+      sum(spending_time >= 1) <= 1
   )
 
   function(alpha, info_frac) {
     st <- spending_time[seq_along(info_frac)]
     spending_fn(alpha, st)
   }
 }
+
+
+#' Wang-Tsiatis spending function
+#'
+#' @description
+#' Computes the implied cumulative alpha spending from the Wang-Tsiatis family
+#' of group sequential boundaries. The Wang-Tsiatis boundaries at analysis
+#' \eqn{k} with information fraction \eqn{t_k} are defined as:
+#' \deqn{c_k = C \cdot t_k^{\Delta - 0.5},}
+#' where \eqn{\Delta} is the shape parameter and \eqn{C} is a constant
+#' calibrated so that the overall Type I error equals \eqn{\alpha}.
+#'
+#' Special cases:
+#' * \eqn{\Delta = 0.5}: Pocock boundaries (equal Z-scale boundaries across
+#'   analyses).
+#' * \eqn{\Delta = 0}: O'Brien-Fleming boundaries (very conservative at
+#'   early analyses).
+#' * \eqn{0 < \Delta < 0.5}: intermediate between O'Brien-Fleming and Pocock.
+#'
+#' Unlike the Lan-DeMets approximations ([spending_of()], [spending_pocock()]),
+#' this function computes the **exact** boundaries from the Wang-Tsiatis
+#' family and derives the implied spending. It is computationally more
+#' expensive because it requires root-finding and multivariate normal
+#' integration at each call.
+#'
+#' @param alpha A numeric scalar of the total significance level.
+#' @param info_frac A numeric vector of information fractions at each analysis.
+#'   Must be non-negative, with at most one value \eqn{\geq 1}.
+#' @param delta A numeric scalar for the shape parameter \eqn{\Delta}.
+#'   The default is `0.5` (Pocock). Use `0` for O'Brien-Fleming.
+#' @param maxpts An integer scalar for the maximum number of function values
+#'   for [mvtnorm::GenzBretz()]. The default is 25000.
+#' @param abseps A numeric scalar for the absolute error tolerance for
+#'   [mvtnorm::GenzBretz()]. The default is 1e-6.
+#'
+#' @return A numeric vector the same length as `info_frac` of cumulative alpha
+#'   spent at each information fraction.
+#'
+#' @seealso [spending_of()] and [spending_pocock()] for the Lan-DeMets
+#'   approximations, [gs_boundaries()] for computing boundaries from spending
+#'   functions, [graph_test_shortcut_gsd()] for the graphical procedure.
+#'
+#' @references
+#'   Wang, S. K., and Tsiatis, A. A. (1987). Approximately optimal one-parameter
+#'   boundaries for group sequential trials. \emph{Biometrics}, 43(1), 193-199.
+#'
+#' @export
+#'
+#' @examples
+#' # Exact O'Brien-Fleming (delta = 0)
+#' spending_wt(0.025, c(0.5, 1), delta = 0)
+#'
+#' # Exact Pocock (delta = 0.5)
+#' spending_wt(0.025, c(0.5, 1), delta = 0.5)
+#'
+#' # Intermediate (delta = 0.25)
+#' spending_wt(0.025, c(1/3, 2/3, 1), delta = 0.25)
+#'
+#' # Compare with Lan-DeMets approximations
+#' spending_of(0.025, c(1/3, 2/3, 1))     # Lan-DeMets OBF approximation
+#' spending_wt(0.025, c(1/3, 2/3, 1), 0)  # Exact OBF
+#'
+#' # Use in graph_test_shortcut_gsd (wrap to fix delta)
+#' \donttest{
+#' g <- graph_create(c(0.5, 0.5), rbind(c(0, 1), c(1, 0)))
+#' p <- rbind(H1 = c(0.024, 0.01), H2 = c(0.015, 0.005))
+#' graph_test_shortcut_gsd(
+#'   graph = g, p = p, alpha = 0.025,
+#'   info_frac = c(0.5, 1),
+#'   spending_fn = function(a, t) spending_wt(a, t, delta = 0.25)
+#' )
+#' }
+spending_wt <- function(alpha, info_frac, delta = 0.5,
+                        maxpts = 25000, abseps = 1e-6) {
+  stopifnot(
+    "info_frac must be non-negative" = all(info_frac >= 0),
+    "At most one info_frac value can be >= 1" = sum(info_frac >= 1) <= 1,
+    "delta must be a numeric scalar" = is.numeric(delta) && length(delta) == 1
+  )
+
+  K <- length(info_frac)
+
+  # Handle edge cases
+  if (alpha <= 0) return(rep(0, K))
+  if (K == 1) return(pmin(alpha, alpha))
+
+  # Correlation matrix
+  corr <- gs_corr(info_frac)
+
+  # Wang-Tsiatis boundary shape: c_k = C * t_k^(delta - 0.5)
+  # For info_frac = 0, the shape is Inf (or 0 depending on delta),
+  # handle by setting those boundaries to Inf (never cross)
+  shape <- ifelse(info_frac == 0, 0, info_frac^(delta - 0.5))
+
+  algo <- mvtnorm::GenzBretz(maxpts = maxpts, abseps = abseps)
+
+  # Find C such that P(cross at some k | H0) = alpha
+  # P(cross) = 1 - P(Z_1 < c_1, ..., Z_K < c_K)
+  find_C <- function(C_val) {
+    bounds_z <- C_val * shape
+    # Replace any Inf or very large bounds with 20 for numerical stability
+    bounds_z <- pmin(bounds_z, 20)
+
+    prob_no_cross <- mvtnorm::pmvnorm(
+      upper = bounds_z,
+      corr = corr,
+      algorithm = algo
+    )[[1]]
+
+    (1 - prob_no_cross) - alpha
+  }
+
+  # Search for C. Boundaries are on the Z-scale, so C is typically 1-5
+  C_root <- tryCatch(
+    stats::uniroot(find_C, interval = c(0.1, 20), tol = abseps),
+    error = function(e) {
+      # Widen search if needed
+      stats::uniroot(find_C, interval = c(0.01, 50), tol = abseps)
+    }
+  )
+  C_val <- C_root$root
+  bounds_z <- C_val * shape
+  bounds_z <- pmin(bounds_z, 20)
+
+  # Compute implied cumulative spending at each analysis k:
+  # alpha_k = P(Z_1 >= c_1 or ... or Z_k >= c_k)
+  #         = 1 - P(Z_1 < c_1, ..., Z_k < c_k)
+  cum_spending <- numeric(K)
+  for (k in seq_len(K)) {
+    if (info_frac[k] == 0) {
+      cum_spending[k] <- 0
+      next
+    }
+    if (k == 1) {
+      # Univariate case: P(Z >= c_1) = 1 - Phi(c_1)
+      cum_spending[k] <- stats::pnorm(bounds_z[1], lower.tail = FALSE)
+    } else {
+      cum_spending[k] <- 1 - mvtnorm::pmvnorm(
+        upper = bounds_z[seq_len(k)],
+        corr = corr[seq_len(k), seq_len(k)],
+        algorithm = algo
+      )[[1]]
+    }
+  }
+
+  # Cap at alpha for numerical stability
+  cum_spending <- pmin(cum_spending, alpha)
+  cum_spending
+}