Add minimum risk weighting strategy to gs_design_rd and gs_power_rd (#613)

* Development
   - Add the MR weighting strategy to `gs_info_rd`
   - Update the `weight = ...` in `gs_design_rd` and `gs_power_rd`
   - Update Rd files
   - Add developer tests to `gs_power_rd` and `gs_design_rd` regarding the MR weights

* Reviewers, Co-authored with John
   - Add `alpha_per_k_per_s` and `beta_per_k_per_s` to the globals.R
   - Add an example of MR weight in gs_info_rd.R
   - Update news.md

Author: LittleBeannie

NEWS.md

@@ -1,3 +1,10 @@
+# gsDesign2 1.1.9
+
+## New features
+
+- The minimal risk weighting strategy has been added to `gs_design_rd()` and `gs_power_rd()` for risk difference design (#611, thanks to @LittleBeannie).
+- The `sequential_pval()` function has been added to calculate the sequential p-value for a AHR group sequential design (#605, thanks to @LittleBeannie).
+
 # gsDesign2 1.1.8
 
 ## New features

R/globals.R

@@ -52,7 +52,8 @@ utils::globalVariables(
       c(
         "prevalence", "z", "info", "theta", "rd", "info0", "n",
         "probability", "probability0", "info_frac0",
-        "~risk difference at bound", "nominal p"
+        "~risk difference at bound", "nominal p",
+        "alpha_per_k_per_s", "beta_per_k_per_s"
       ),
       # From `gs_design_wlr()`
       c(

R/gs_design_rd.R

@@ -135,7 +135,7 @@ gs_design_rd <- function(p_c = tibble::tibble(stratum = "All", rate = .2),
                          beta = .1,
                          ratio = 1,
                          stratum_prev = NULL,
-                         weight = c("unstratified", "ss", "invar"),
+                         weight = c("unstratified", "ss", "invar", "mr"),
                          upper = gs_b,
                          lower = gs_b,
                          upar = gsDesign(k = 3, test.type = 1, sfu = sfLDOF, sfupar = NULL)$upper$bound,

R/gs_info_rd.R

@@ -92,7 +92,7 @@
 #' )
 #'
 #' # Example 5 ----
-#' # stratified case under sample size weighting and H0: rd0 != 0
+#' # stratified case under minimal risk weighting and H0: rd0 = 0
 #' gs_info_rd(
 #'   p_c = tibble::tibble(
 #'     stratum = c("S1", "S2", "S3"),
@@ -107,13 +107,13 @@
 #'     analysis = rep(1:3, 3),
 #'     n = c(50, 100, 200, 40, 80, 160, 60, 120, 240)
 #'   ),
-#'   rd0 = 0.02,
+#'   rd0 = 0,
 #'   ratio = 1,
-#'   weight = "ss"
+#'   weight = "mr"
 #' )
 #'
 #' # Example 6 ----
-#' # stratified case under inverse variance weighting and H0: rd0 != 0
+#' # stratified case under sample size weighting and H0: rd0 != 0
 #' gs_info_rd(
 #'   p_c = tibble::tibble(
 #'     stratum = c("S1", "S2", "S3"),
@@ -130,12 +130,13 @@
 #'   ),
 #'   rd0 = 0.02,
 #'   ratio = 1,
-#'   weight = "invar"
+#'   # users can switch to either "invar" or "mr" weighting as well
+#'   weight = "ss"
 #' )
 #'
 #' # Example 7 ----
 #' # stratified case under inverse variance weighting and H0: rd0 != 0 and
-#' # rd0 difference for different statum
+#' # rd0 difference for different stratum
 #' gs_info_rd(
 #'   p_c = tibble::tibble(
 #'     stratum = c("S1", "S2", "S3"),
@@ -173,7 +174,7 @@ gs_info_rd <- function(
     ),
     rd0 = 0,
     ratio = 1,
-    weight = c("unstratified", "ss", "invar")) {
+    weight = c("unstratified", "ss", "invar", "mr")) {
   n_analysis <- max(n$analysis)
   weight <- match.arg(weight)
 
@@ -237,6 +238,24 @@ gs_info_rd <- function(
         mutate(weight_per_k_per_s = 1 / sigma2_H1_per_k_per_s / sum_inv_var_per_s) |>
         select(-sum_inv_var_per_s)
     )
+  } else if (weight == "mr") {
+    suppressMessages(
+      tbl <- tbl |>
+        left_join(
+          tbl |>
+            group_by(analysis) |>
+            summarize(sum_inv_var_per_s = sum(1 / sigma2_H1_per_k_per_s))
+        ) |>
+        ungroup() |>
+        group_by(analysis) |>
+        mutate(alpha_per_k_per_s = (p_c - p_e) * sum_inv_var_per_s - sum((p_c - p_e) / sigma2_H1_per_k_per_s),
+               beta_per_k_per_s = 1/sigma2_H1_per_k_per_s * (1 + alpha_per_k_per_s * sum((p_c - p_e) * n / max(n))),
+               weight_per_k_per_s = beta_per_k_per_s / sum_inv_var_per_s -
+                                       alpha_per_k_per_s / sigma2_H1_per_k_per_s / (sum_inv_var_per_s + sum(alpha_per_k_per_s * (p_c - p_e) / sigma2_H1_per_k_per_s)) *
+                                       sum((p_c - p_e) * beta_per_k_per_s) / sum_inv_var_per_s
+        ) |>
+        select(-c(sum_inv_var_per_s, alpha_per_k_per_s, beta_per_k_per_s))
+    )
   }
 
   # Pool the strata together ----

R/gs_power_rd.R

@@ -243,7 +243,7 @@ gs_power_rd <- function(
     ),
     rd0 = 0,
     ratio = 1,
-    weight = c("unstratified", "ss", "invar"),
+    weight = c("unstratified", "ss", "invar", "mr"),
     upper = gs_b,
     lower = gs_b,
     upar = gsDesign(k = 3, test.type = 1, sfu = sfLDOF, sfupar = NULL)$upper$bound,

man/gs_design_rd.Rd

@@ -28,3 +28,39 @@ test_that("fixed design", {
 
   expect_equal(x1, x2$analysis$n)
 })
+
+test_that("Stratified GSD: if RD is constant across strata, then MR weights are equal to the INVAR weights", {
+  # Reference: Section 3 of Mehrotra, Devan V., and Radha Railkar.
+  # "Minimum risk weights for comparing treatments in stratified binomial trials." Statistics in Medicine 19.6 (2000): 811-825.
+  x_invar <- gs_design_rd(
+    p_c = tibble::tibble(stratum = c("biomarker positive", "biomarker negative"),
+                         rate = c(.2, .25)),
+    p_e = tibble::tibble(stratum = c("biomarker positive", "biomarker negative"),
+                         rate = c(.15, .20)),
+    rd0 = 0, info_frac = c(0.7, 1),
+    alpha = .025, beta = .1, ratio = 1,
+    stratum_prev = tibble::tibble(stratum = c("biomarker positive", "biomarker negative"),
+                                  prevalence = c(.4, .6)),
+    weight = "invar",
+    upper = gs_spending_bound, lower = gs_b,
+    upar = list(sf = gsDesign::sfLDOF, total_spend = 0.025, param = NULL, timing = NULL),
+    lpar = rep(-Inf, 2))
+
+  x_mr <- gs_design_rd(
+    p_c = tibble::tibble(stratum = c("biomarker positive", "biomarker negative"),
+                         rate = c(.2, .25)),
+    p_e = tibble::tibble(stratum = c("biomarker positive", "biomarker negative"),
+                         rate = c(.15, .20)),
+    rd0 = 0, info_frac = c(0.7, 1),
+    alpha = .025, beta = .1, ratio = 1,
+    stratum_prev = tibble::tibble(stratum = c("biomarker positive", "biomarker negative"),
+                                  prevalence = c(.4, .6)),
+    weight = "mr",
+    upper = gs_spending_bound, lower = gs_b,
+    upar = list(sf = gsDesign::sfLDOF, total_spend = 0.025, param = NULL, timing = NULL),
+    lpar = rep(-Inf, 2))
+
+  expect_equal(x_invar$analysis, x_mr$analysis)
+  expect_equal(x_invar$bound, x_mr$bound)
+})
+

man/gs_info_rd.Rd

@@ -29,3 +29,39 @@ test_that("fixed design", {
 
   expect_equal(x1, x2$bound$probability)
 })
+
+test_that("Stratified GSD: if RD is constant across strata, then MR weights are equal to the INVAR weights", {
+  # Reference: Section 3 of Mehrotra, Devan V., and Radha Railkar.
+  # "Minimum risk weights for comparing treatments in stratified binomial trials." Statistics in Medicine 19.6 (2000): 811-825.
+  x_invar <- gs_power_rd(
+    p_c = tibble::tibble(stratum = c("biomarker positive", "biomarker negative"),
+                         rate = c(.2, .25)),
+    p_e = tibble::tibble(stratum = c("biomarker positive", "biomarker negative"),
+                         rate = c(.15, .20)),
+    n = tibble::tibble(stratum = rep(c("biomarker positive", "biomarker negative"), each = 2),
+                       n = c(1000, 1500, 1000, 1500),
+                       analysis = c(1, 2, 1, 2)),
+    rd0 = 0, ratio = 1,
+    weight = "invar",
+    upper = gs_spending_bound, lower = gs_b,
+    upar = list(sf = gsDesign::sfLDOF, total_spend = 0.025, param = NULL, timing = NULL),
+    lpar = rep(-Inf, 2))
+
+  x_mr <- gs_power_rd(
+    p_c = tibble::tibble(stratum = c("biomarker positive", "biomarker negative"),
+                         rate = c(.2, .25)),
+    p_e = tibble::tibble(stratum = c("biomarker positive", "biomarker negative"),
+                         rate = c(.15, .20)),
+    n = tibble::tibble(stratum = rep(c("biomarker positive", "biomarker negative"), each = 2),
+                       n = c(1000, 1500, 1000, 1500),
+                       analysis = c(1, 2, 1, 2)),
+    rd0 = 0, ratio = 1,
+    weight = "mr",
+    upper = gs_spending_bound, lower = gs_b,
+    upar = list(sf = gsDesign::sfLDOF, total_spend = 0.025, param = NULL, timing = NULL),
+    lpar = rep(-Inf, 2))
+
+  expect_equal(x_invar$analysis, x_mr$analysis)
+  expect_equal(x_invar$bound, x_mr$bound)
+})
+