Updated the gs_cp() function via incoporating both original design an…

…d updated design; Also updated the test-developer-gs_cp file for the updated gs_cp() function.

R/gs_cp.R

@@ -21,12 +21,13 @@
 #' @description  \code{gs_cp()} computes conditional boundary crossing probabilities at future planned analyses
 #' for a given group sequential design assuming an interim z-statistic at a specified interim analysis.
 #' @param x An object of type \code{gsDesign2}
-#' @param i analysis at which interim z-value is given; must be from 1 to max(x$analysis$analysis)
-#' @param zi interim z-value at analysis i (scalar)
-#' @param theta_hat Scalar or vector of estimated theta, could be output of ahr_blinded() or arbitrary value(s)
+#' @param x_updated An object of type \code{gsDesign2}
+#' @param i Analysis at which interim z-value is given; must be from 1 to max(x$analysis$analysis)
+#' @param zi Interim z-value at analysis i (scalar)
+#' @param local_alternative A TRUE or FALSE value that indicates if we apply local approximation theorem
 #' @return A list with: theta -- a list containing theta0 (theta under H0, i.e., 0,) and theta1 (a vector of theta under H1)
 #'                      upper_bound -- the upper bound value return by any gs_design_ahr function
-#'                      upper_prob -- a list contains the conditional probability given zi under H0 and H1, respectively
+#'                      upper_prob -- a list contains the conditional probability given zi under H0, H1 and estimated treatment effect, respectively
 #'
 #' @examples
 #' library(gsDesign2)
@@ -83,68 +84,85 @@
 #'   ustime = ustime,
 #'   observed_data = list(observed_data_ia, NULL, NULL))
 #'
-#' theta_hat <- 0.29
-#' gs_cp(x_updated, i = 1, zi = 2, j = 3, theta_hat = theta_hat)
+#' gs_cp(x = x, x_updated = x_updated, i = 1, zi = 2, j = 3)
 
 
 
 
-gs_cp <- function(x, i, zi, j, theta_hat){
+gs_cp <- function(x, x_updated, i, zi, j, local_alternative){
 
   # input check
   # Check if 'x' has an 'analysis' element and it's a matrix or data frame
   if (!is.list(x) || !"analysis" %in% names(x) || (!is.matrix(x$analysis) && !is.data.frame(x$analysis))) {
     stop("'x' should be a list containing an 'analysis' matrix or data frame.")
   }
 
+  # Check if 'x_updated' has an 'analysis' element and it's a matrix or data frame
+  if (!is.list(x_updated) || !"analysis" %in% names(x_updated) || (!is.matrix(x_updated$analysis) && !is.data.frame(x_updated$analysis))) {
+    stop("'x_updated' should be a list containing an 'analysis' matrix or data frame.")
+  }
+
   # Check if 'i' and 'j' are numeric and within the appropriate range
   if (!is.numeric(i) || !is.numeric(j)) {
     stop("'i' and 'j' should be numeric.")
   }
 
-  if (!(1 <= i && i < j && j <= dim(x$analysis)[1])) {
-    stop("Please ensure that 1 <= i < j and j <= dim(x$analysis)[1].")
+  if (!(1 <= i && i < j && j <= dim(x$analysis)[1]) |
+      !(j <= dim(x_updated$analysis)[1])) {
+    stop("Please ensure that 1 <= i < j and j <= dim(x$analysis)[1] and j <= dim(x_updated$analysis)[1].")
   }
 
   # Check the class of `x`, which should be an output from `gs_update_ahr`
-  if(!("updated_design" %in% class(x))){
-    stop("'x' should be an output from gs_update_ahr() ")
+  if(!("updated_design" %in% class(x_updated))){
+    stop("'x_updated' should be an output from gs_update_ahr() ")
   }
 
   # Check the # of upper bound is equal to # of analysis
   # ! what is efficacy is only at IA2 and FA? >= 3 .. Check if j has upper bound
-  if(length(which(x$bound$bound == "upper")) != dim(x$analysis)[1]){
-    stop("'x' should contains the same number of upper bounds as the number of analysis")
+  if(length(which(x$bound$bound == "upper")) != dim(x$analysis)[1] |
+     length(which(x_updated$bound$bound == "upper")) != dim(x_updated$analysis)[1]){
+    stop("'x' and 'x_updated' should contains the same number of upper bounds as the number of analysis")
   }
 
   # obtain necessary information from x
   info_frac <- x$analysis$info_frac
   info0 <- x$analysis$info0 # under local asymptotic, assume H0 ~= H1
-
+  info <- x$analysis$info
+  info0_hat <- x_updated$analysis$info0
   # default theta: under H0 and under H1
-  theta0 <- c(0, 0, 0)
+  theta0 <- rep(0, dim(x$analysis)[1])
   theta1 <- x$analysis$theta
-  theta_est <- c(theta_hat, theta1[-c(1:length(theta_hat))])
+  theta_est <- x_updated$analysis$theta
 
   # compute the conditional probability under H0
   eff_bound <- as.data.frame(x$bound) %>%
     filter(bound == "upper") %>%
     select(z)
 
   # compute conditional power under H0
-  # prob0 <- 1 - pnorm((eff_bound[i+1, ] - sqrt(info_frac[i]/info_frac[i+1]) * zi)/sqrt((info_frac[i+1]-info_frac[i])/info_frac[i+1]))
   prob0 <- 1 - pnorm((sqrt(info_frac[j])*eff_bound[j, ] - sqrt(info_frac[i])*zi)/sqrt(info_frac[j] - info_frac[i]))
 
 
   # compute conditional power under H1
   mu_star <- sqrt(info_frac[j])*sqrt(info0[j])*theta1[j] - sqrt(info_frac[i])*sqrt(info0[i])*theta1[i]
-  sigma2_star <- info_frac[j] - info_frac[i]
+
+  if(local_alternative){
+    sigma2_star <- info_frac[j] - info_frac[i]
+  }else{
+    sigma2_star <- info_frac[j]*info0[j]/info[j] - info_frac[i]*info0[i]/info[i]
+  }
+
   prob1 <- 1 - pnorm((sqrt(info_frac[j])*eff_bound[j, ] - sqrt(info_frac[i])*zi - mu_star)/sqrt(sigma2_star))
 
 
   # compute conditional power under estimated theta
-  mu_star_est <- sqrt(info_frac[j])*sqrt(info0[j])*theta_est[j] - sqrt(info_frac[i])*sqrt(info0[i])*theta_est[i]
-  sigma2_star_est <- info_frac[j] - info_frac[i]
+  mu_star_est <- sqrt(info_frac[j])*sqrt(info0[j])*theta1[j] - sqrt(info_frac[i])*sqrt(info0_hat[i])*theta_est[i]
+  if(local_alternative){
+    sigma2_star_est <- info_frac[j] - info_frac[i]*info0_hat[i]/info[i]
+  }else{
+    sigma2_star_est <- info_frac[j]*info0[j]/info[j] - info_frac[i]*info0_hat[i]/info[i]
+  }
+
   prob_est <- 1 - pnorm((sqrt(info_frac[j])*eff_bound[j, ] - sqrt(info_frac[i])*zi - mu_star_est)/sqrt(sigma2_star_est))
 
   # return list of results

tests/testthat/test-developer-gs_cp.R

@@ -92,9 +92,8 @@ test_that("Compare the conditional power of gsDesign and gsDesign2 under PH with
   # theta_hat <- blinded$theta
 
   theta_hat <- x2$theta[1] # observed HR
-  y2 <- gs_cp(x = y1, i = 1, zi = -qnorm(p), j = 2, theta_hat = theta_hat)
-  y3 <- gs_cp(x = y1, i = 1, zi = -qnorm(p), j = 3, theta_hat = theta_hat)
-
+  y2 <- gs_cp(x = y0, x_updated = y1, i = 1, zi = -qnorm(p), j = 2, local_alternative = TRUE)
+  y3 <- gs_cp(x = y0, x_updated = y1, i = 1, zi = -qnorm(p), j = 3, local_alternative = TRUE)
 
   # ------------------------------ #
   #       comparison               #
@@ -104,37 +103,33 @@ test_that("Compare the conditional power of gsDesign and gsDesign2 under PH with
 
   # comparison of events of the original design and updated design
   expect_equal(x0$n.I, y0$analysis$event, tolerance = 1e-5)
-  expect_equal(x1$n.I, y1_updated$analysis$event, tolerance = 1e-5)
+  expect_equal(x1$n.I, y1$analysis$event, tolerance = 1e-5)
 
   # comparison of timing of the original design and updated design
   expect_equal((x0$T) |> as.vector(), y0$analysis$time, tolerance = 1e-5)
   expect_equal((x1$timing) |> as.vector(), y1$analysis$info_frac0, tolerance = 1e-5)
 
-  # comparison of crossing probability under H1
-  expect_equal(x0$upper$prob[, 2] |> cumsum(), y0$bound$probability, tolerance = 1e-2)
-  expect_equal(x1$upper$prob[, 2] |> cumsum(), y1$bound$probability, tolerance = 1e-2)
-
   # comparison of crossing probability under H0
   expect_equal(x0$upper$prob[, 1] |> cumsum(), y0$bound$probability0, tolerance = 1e-2)
   expect_equal(x1$upper$prob[, 1] |> cumsum(), y1$bound$probability0, tolerance = 1e-2)
 
+  # comparison of crossing probability under H1
+  expect_equal(x0$upper$prob[, 2] |> cumsum(), y0$bound$probability, tolerance = 1e-2)
+  expect_equal(x1$upper$prob[, 2] |> cumsum(), y1$bound$probability, tolerance = 1e-2)
+
   # comparison of conditional power
   # theta = H0
-  expect_equal(x2$upper$prob[1, 2], y2$upper_prob$prob0, tolerance = 1e-3)
+  expect_equal(x2$upper$prob[1, 2], y2$upper_prob$prob0, tolerance = 1e-2)
   expect_equal(sum(x2$upper$prob[, 2]), y3$upper_prob$prob0, tolerance = 1e-1)
+
   # theta = H1
-  expect_equal(x2$upper$prob[1, 3], y2$upper_prob$prob1, tolerance = 1e-2)
+  expect_equal(x2$upper$prob[1, 3], y2$upper_prob$prob1, tolerance = 1e-1)
   expect_equal(sum(x2$upper$prob[, 3]), y3$upper_prob$prob1, tolerance = 1e-2)
+
   # theta = IA estimated theta
-  ## !!! Not pass the tests: 1. consider the theta_hat  2. consider info0_hat
-  #expect_equal(sum(x2$upper$prob[1, 1]), y2$upper_prob$prob_est, tolerance = 1e-2)
-  #expect_equal(sum(x2$upper$prob[, 1]), y3$upper_prob$prob_est, tolerance = 1e-2)
+  expect_equal(sum(x2$upper$prob[1, 1]), y2$upper_prob$prob_est, tolerance = 2e-1)
+  expect_equal(sum(x2$upper$prob[, 1]), y3$upper_prob$prob_est, tolerance = 2e-1)
 })
 
 
 
-
-
-
-
-