Add harm bound to AHR designs

.claude/skills/write-tests.md

@@ -16,11 +16,11 @@ This project uses **testit** for testing. testit assertions are plain R expressi
 tests/
 ├── test-all.R              # Runner: library(testit); test_pkg("gsDesign2")
 └── testit/
-    ├── helper.R            # Shared setup (sourced before test files)
-    ├── helper-*.R          # Additional helpers
+    ├── helper.R            # Shared setup (auto-sourced before test files)
+    ├── helper-*.R          # Additional helpers (also auto-sourced)
     ├── fixtures/           # Test data (.Rdata, .rds, etc.)
     ├── test-*.R            # Test files
-    └── test-*.md           # Snapshot files (paired with .R files)
+    └── test-*.md           # Snapshot files (standalone, no .R file needed)
 ```
 
 ## Core Pattern
@@ -131,7 +131,7 @@ assert("output structure is correct", {
 
 ## Snapshot Tests
 
-Create a `.md` file alongside the `.R` file (same base name). Format:
+A `.md` snapshot file is a standalone test — it does NOT require a paired `.R` file. The `.md` file contains both the code and the expected output. Format:
 
 ````markdown
 ## `function_name()` description
@@ -154,7 +154,7 @@ testit runs the R code block and compares output to the text block. To initializ
 3. **Use `all.equal(..., tolerance = t)` with the tightest tolerance that passes** — don't use overly loose tolerances.
 4. **Group related assertions in one `assert()` block** — each block should test one logical concept.
 5. **Use descriptive assert messages** — they appear in failure output.
-6. **Shared setup goes in `helper.R`** — it's sourced before all test files.
+6. **Shared setup goes in `helper*.R` files** — testit auto-sources all `helper*.R` files before test files. Never `source()` them manually.
 7. **Load fixture data with `load("fixtures/file.Rdata")`** — paths are relative to `tests/testit/`.
 8. **Use `all.equal()` only when exact comparison fails in CI** — typically macOS produces slightly different floating-point results while `identical()` works fine on Windows/Linux.
 

DESCRIPTION

@@ -68,4 +68,4 @@ VignetteBuilder:
 LinkingTo:
     Rcpp
 Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.3.3
+Config/roxygen2/version: 8.0.0

R/as_gt.R

@@ -111,7 +111,7 @@ as_gt.fixed_design_summary <- function(x, title = NULL, footnote = NULL, ...) {
 #'   "spanner")`; users can use the functions in the `gt` package to customize
 #'   the table. To disable footnotes, use `footnote = FALSE`.
 #' @param display_bound A vector of strings specifying the label of the bounds.
-#'   The default is `c("Efficacy", "Futility")`.
+#'   The default is `c("Efficacy", "Futility", "Harm")`.
 #' @param display_columns A vector of strings specifying the variables to be
 #'   displayed in the summary table.
 #' @param display_inf_bound Logical, whether to display the +/-inf bound.
@@ -128,7 +128,22 @@ as_gt.fixed_design_summary <- function(x, title = NULL, footnote = NULL, ...) {
 #' gs_design_ahr() |>
 #'   summary() |>
 #'   as_gt()
-#'
+#' 
+#' gs_design_ahr(
+#'   analysis_time = c(12, 24, 36),
+#'   upper = gs_spending_bound,
+#'   upar = list(sf = gsDesign::sfLDOF, total_spend = 0.025),
+#'   test_upper = c(FALSE, TRUE, TRUE),
+#'   lower = gs_spending_bound,
+#'   lpar = list(sf = gsDesign::sfHSD, total_spend = 0.1, param = -2),
+#'   test_lower = c(TRUE, TRUE, FALSE),
+#'   harm = gs_spending_bound,
+#'   hpar = list(sf = gsDesign::sfHSD, total_spend = 0.2, param = -4),
+#'   test_harm = c(TRUE, TRUE, FALSE)
+#'   ) |>
+#'  summary() |>
+#'  as_gt() 
+#' 
 #' gs_power_ahr(lpar = list(sf = gsDesign::sfLDOF, total_spend = 0.1)) |>
 #'   summary() |>
 #'   as_gt()
@@ -193,7 +208,7 @@ as_gt.fixed_design_summary <- function(x, title = NULL, footnote = NULL, ...) {
 #'
 #' # Example 5 ----
 #' # Usage of display_bound = ...
-#' # to either show efficacy bound or futility bound, or both(default)
+#' # to show selected bounds
 #' gs_power_wlr(lpar = list(sf = gsDesign::sfLDOF, total_spend = 0.1)) |>
 #'   summary() |>
 #'   as_gt(display_bound = "Efficacy")
@@ -212,7 +227,7 @@ as_gt.gs_design_summary <- function(
     colname_spanner = "Cumulative boundary crossing probability",
     colname_spannersub = c("Alternate hypothesis", "Null hypothesis"),
     footnote = NULL,
-    display_bound = c("Efficacy", "Futility"),
+    display_bound = c("Efficacy", "Futility", "Harm"),
     display_columns = NULL,
     display_inf_bound = FALSE,
     ...) {
@@ -364,6 +379,7 @@ gsd_parts <- function(
   x2 <- x2[, columns]
   x2 <- subset(x2, !is.na(`Alternate hypothesis`) & !is.na(`Null hypothesis`))
   x2 <- subset(x2, Bound %in% bound)
+  x2$Bound <- factor(x2$Bound, levels = bound)
 
   i <- match(c("Alternate hypothesis", "Null hypothesis"), names(x2))
   names(x2)[i] <- spannersub
@@ -382,10 +398,10 @@ gsd_parts <- function(
   )
 
   list(
-    x = arrange(x2, Analysis),
+    x = arrange(x2, Analysis, Bound),
     title = title, subtitle = subtitle,
     footnote = if (!isFALSE(footnote)) footnote %||% gsd_footnote(method, columns),
-    alpha = max(filter(x, Bound == bound[1])[["Null hypothesis"]])
+    alpha = max(filter(x, Bound == "Efficacy")[["Null hypothesis"]])
   )
 }
 

R/gs_bound_summary.R

@@ -1,6 +1,6 @@
 #' Bound summary table
 #'
-#' Summarizes the efficacy and futility bounds for each analysis.
+#' Summarizes the efficacy, futility, and harm bounds for each analysis.
 #'
 #' @param x Design object.
 #' @param alpha Vector of alpha values to compute additional efficacy columns.
@@ -52,14 +52,25 @@ gs_bound_summary <- function(x, digits = 4, ddigits = 2, tdigits = 0, timename =
     }
   }
   out <- Reduce(cbind, outlist)
-  # Use of union() allows placement of column "Futility" at the far right, but
-  # only if it is returned by gs_bound_summary_single(). This is because
-  # one-sided designs do not produce a Futility column.
+  # Use of union() allows placement of columns "Futility" and "Harm" at the far
+  # right, but only if they are returned by gs_bound_summary_single(). This is
+  # because one-sided designs do not produce a Futility column, and designs
+  # without harm bounds do not produce a Harm column.
   column_order <- union(c("Analysis", "Value", col_efficacy_name), colnames(out))
   out <- out[, column_order]
   return(out)
 }
 
+gs_bound_summary_values <- function(bound, analysis, bound_name, columns) {
+  row_bound <- bound[
+    bound$analysis == analysis & bound$bound == bound_name,
+    columns,
+    drop = FALSE
+  ]
+  if (nrow(row_bound) == 0) return(rep(NA_real_, length(columns)))
+  as.numeric(unlist(row_bound[1, columns], use.names = FALSE))
+}
+
 gs_bound_summary_single <- function(x, col_efficacy_name = "Efficacy", digits,
                                     ddigits, tdigits, timename) {
   # Input
@@ -72,6 +83,8 @@ gs_bound_summary_single <- function(x, col_efficacy_name = "Efficacy", digits,
   col_value <- character()
   col_efficacy <- numeric()
   col_futility <- numeric()
+  col_harm <- numeric()
+  bound_columns <- c("z", "nominal p", "~hr at bound", "probability0", "probability")
 
   for (i in seq_len(nrow(analysis))) {
 
@@ -113,33 +126,31 @@ gs_bound_summary_single <- function(x, col_efficacy_name = "Efficacy", digits,
     )
 
     # Efficacy column
-    row_efficacy <- bound[
-      bound$analysis == i & bound$bound == "upper",
-      c("z", "nominal p", "~hr at bound", "probability0", "probability")
-    ]
-    col_efficacy <- c(col_efficacy, as.numeric(row_efficacy))
+    col_efficacy <- c(col_efficacy, gs_bound_summary_values(bound, i, "upper", bound_columns))
 
     # Futility column
-    row_futility <- bound[
-      bound$analysis == i & bound$bound == "lower",
-      c("z", "nominal p", "~hr at bound", "probability0", "probability")
-    ]
-    col_futility <- c(col_futility, as.numeric(row_futility))
+    col_futility <- c(col_futility, gs_bound_summary_values(bound, i, "lower", bound_columns))
+
+    # Harm column
+    col_harm <- c(col_harm, gs_bound_summary_values(bound, i, "harm", bound_columns))
   }
 
   col_efficacy <- round(col_efficacy, digits)
   col_futility <- round(col_futility, digits)
+  col_harm <- round(col_harm, digits)
 
   out <- data.frame(
     Analysis = col_analysis,
     Value = col_value,
     Efficacy = col_efficacy,
-    Futility = col_futility
+    Futility = col_futility,
+    Harm = col_harm
   )
   colnames(out)[3] <- col_efficacy_name
 
   # One-sided design should not include Futility column
   if (all(is.na(out[["Futility"]]))) out[["Futility"]] <- NULL
+  if (all(is.na(out[["Harm"]]))) out[["Harm"]] <- NULL
 
   return(out)
 }

R/gs_design_ahr.R

@@ -98,7 +98,7 @@
 #' # Example 2 ----
 #' # Single analysis
 #' gs_design_ahr(analysis_time = 40)
-#'
+#' 
 #' # Example 3 ----
 #' # Multiple analysis_time
 #' gs_design_ahr(analysis_time = c(12, 24, 36))
@@ -168,6 +168,22 @@
 #'   lpar = rep(-Inf, 3)
 #' )
 #' }
+#'
+#' # Example 8 ----
+#' # Design with an additional harm bound
+#' \donttest{
+#' gs_design_ahr(
+#'   analysis_time = c(12, 24, 36),
+#'   upper = gs_spending_bound,
+#'   upar = list(sf = gsDesign::sfLDOF, total_spend = 0.025, param = NULL, timing = NULL),
+#'   lower = gs_spending_bound,
+#'   lpar = list(sf = gsDesign::sfHSD, total_spend = 0.1, param = -2, timing = NULL),
+#'   test_lower = c(TRUE, TRUE, FALSE),
+#'   harm = gs_spending_bound,
+#'   hpar = list(sf = gsDesign::sfHSD, total_spend = 0.2, param = -4, timing = NULL),
+#'   test_harm = c(TRUE, TRUE, FALSE)
+#' )
+#' }
 gs_design_ahr <- function(
     enroll_rate = define_enroll_rate(
       duration = c(2, 2, 10),
@@ -186,9 +202,12 @@ gs_design_ahr <- function(
     upar = list(sf = gsDesign::sfLDOF, total_spend = alpha),
     lower = gs_spending_bound,
     lpar = list(sf = gsDesign::sfLDOF, total_spend = beta),
+    harm = gs_b,
+    hpar = -Inf,
     h1_spending = TRUE,
     test_upper = TRUE,
     test_lower = TRUE,
+    test_harm = FALSE,
     info_scale = c("h0_h1_info", "h0_info", "h1_info"),
     r = 18,
     tol = 1e-6,
@@ -200,6 +219,10 @@ gs_design_ahr <- function(
   info_scale <- match.arg(info_scale)
   upper <- match.fun(upper)
   lower <- match.fun(lower)
+  harm <- match.fun(harm)
+
+  # Number of analyses (including final analysis)
+  n_analysis <- max(length(analysis_time), length(info_frac))
 
   # Check inputs ----
   check_analysis_time(analysis_time)
@@ -235,9 +258,6 @@ gs_design_ahr <- function(
   final_event <- y$event[nrow(y)]
   if_alt <- y$event / final_event
 
-  # Number of analyses (including final analysis)
-  n_analysis <- max(length(analysis_time), length(info_frac))
-
   # Initialize the next_time as the study duration
   next_time <- max(analysis_time)
 
@@ -308,6 +328,7 @@ gs_design_ahr <- function(
       alpha = alpha, beta = beta, binding = binding,
       upper = upper, upar = upar, test_upper = test_upper,
       lower = lower, lpar = lpar, test_lower = test_lower,
+      harm = harm, hpar = hpar, test_harm = test_harm,
       r = r, tol = tol
     )
   )
@@ -359,7 +380,6 @@ gs_design_ahr <- function(
       spending_time_upper <- info0 / info0_final
     }
 
-
     bound$spending_time[which(bound$bound == "upper")] <- spending_time_upper
   }
 
@@ -375,6 +395,16 @@ gs_design_ahr <- function(
     bound$spending_time[which(bound$bound == "lower")] <- spending_time_lower
   }
 
+  if (identical(harm, gs_spending_bound)) {
+    if (!is.null(hpar$timing)) {
+      spending_time_harm <- hpar$timing
+    } else {
+      spending_time_harm <- info0 / info0_final
+    }
+
+    bound$spending_time[which(bound$bound == "harm")] <- spending_time_harm
+  }
+
   if (all(is.na(bound$spending_time))){
     bound$spending_time <- NULL
   }
@@ -393,7 +423,8 @@ gs_design_ahr <- function(
     info_scale = info_scale,
     upper = upper, upar = upar,
     lower = lower, lpar = lpar,
-    test_upper = test_upper, test_lower = test_lower,
+    harm = harm, hpar = hpar,
+    test_upper = test_upper, test_lower = test_lower, test_harm = test_harm,
     h1_spending = h1_spending, binding = binding,
     info_scale = info_scale, r = r, tol = tol
   )