Develop gs_cp_ahr following gsDesign::gsCP

NAMESPACE

@@ -30,7 +30,8 @@ export(fixed_design_rd)
 export(fixed_design_rmst)
 export(gs_b)
 export(gs_bound_summary)
-export(gs_cp_npe)
+export(gs_cp_npe1)
+export(gs_cp_npe2)
 export(gs_create_arm)
 export(gs_design_ahr)
 export(gs_design_combo)

R/gs_cp_npe.R

@@ -17,67 +17,133 @@
 #  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 #' Conditional power computation with non-constant effect size
-#'
-#' @details
-#' We assume \eqn{Z_1} and \eqn{Z_2} are the z-values at an interim analysis and later analysis, respectively.
-#' We assume further \eqn{Z_1} and \eqn{Z_2} are bivariate normal with standard group sequential assumptions
-#' on independent increments where for \eqn{i=1,2}
+#' We assume \eqn{Z_i, i = 1, ..., K} are the z-values at an interim analysis i, respectively.
+#' We assume further \eqn{Z_i, i = 1, ..., K} follows multivariate normal distribution
 #' \deqn{E(Z_i) = \theta_i\sqrt{I_i}}
-#' \deqn{Var(Z_i) = 1/I_i}
-#' \deqn{Cov(Z_1, Z_2) = t \equiv I_1/I_2}
-#' where \eqn{\theta_1, \theta_2} are real values and \eqn{0<I_1<I_2}.
+#' \deqn{Cov(Z_i, Z_j) = \sqrt{t_i/t_j}}
 #' See https://merck.github.io/gsDesign2/articles/story-npe-background.html for assumption details.
-#' Returned value is
-#' \deqn{P(Z_2 > b \mid Z_1 = a) = 1 - \Phi\left(\frac{b - \sqrt{t}a - \sqrt{I_2}(\theta_2 - \theta_1\sqrt{t})}{\sqrt{1 - t}}\right)}
+#' For simple conditional power, the returned value is
+#' \deqn{P(Z_j > b \mid Z_i = c) = 1 - \Phi\left(\frac{b - \sqrt{t}c - \sqrt{I_j}(\theta_j - \theta_i\sqrt{t})}{\sqrt{1 - t}}\right)}
+#' For non-simple conditional power, the returned value is the list of
+#' \deqn{P(\{Z_j \geq b_j\} \& \{\cap_{m=i+1}^{j-1} a_m \leq Z_m < b_m\} \mid Z_i = c)}.
+#' \deqn{P(\{Z_j \geq b_j\} \& \{\cap_{m=i+1}^{j-1} Z_m < b_m\} \mid Z_i = c)}.
+#' \deqn{P(\{Z_j \leq b_j\} \& \{\cap_{m=i+1}^{j-1} a_m \leq Z_m < b_m\} \mid Z_i = c)}.
 #'
-#' @param theta A vector of length two, which specifies the natural parameter for treatment effect.
+#' @param theta For simple conditional power, `theta` is a vector of length two, which specifies the natural parameter for treatment effect.
 #'              The first element of `theta` is the treatment effect of an interim analysis i.
 #'              The second element of `theta` is the treatment effect of a future analysis j.
-#' @param info A vector of two, which specifies the statistical information under the treatment effect `theta`.
-#' @param a Interim z-value at analysis i (scalar).
-#' @param b Future target z-value at analysis j (scalar).
-#' @return A scalar with the conditional power \eqn{P(Z_2>b\mid Z_1=a)}.
-#' @export
 #'
+#'              For non-simple conditional power, `theta` is a vector of j-i+1, which specifies the natural parameter for treatment effect.
+#'              The first element of `theta` is the treatment effect of an interim analysis i.
+#'              The second element of `theta` is the treatment effect of an interim analysis i+1.
+#'              ...
+#'              The last element of `theta` is the treatment effect of a future analysis j.
+#' @param t For non-simple conditional power, `t` is a vector of j-i+1, which specifies the information fraction under the treatment effect `theta`.
+#' @param info For simple conditional power, `info` is a vector of length two, which specifies the statistical information under the treatment effect `theta`.
+#'             For non-simple conditional power, `info` is a vector of j-i+1, which specifies the statistical information under the treatment effect `theta`.
+#' @param a For non-simple conditional power, `a` is a vector of length j-i, which specifies the futility bounds from analysis i+1 to analysis j.
+#' @param b For simple conditional power, `b` is a scalar which specifies the future target z-value at analysis j.
+#'          For non-simple conditional power, `b` is a vector of length j-i, which specifies the efficacy bounds from analysis i+1 to analysis j.
+#' @param c For both simple and non-simple conditional power, `c` is the interim z-value at analysis i (scalar).
+#' @return For simple conditional power, the function returns a scalar with the conditional power \eqn{P(Z_2 > b\mid Z_1 = c)}.
+#'         For non-simple conditional power, the function returns a list of conditional powers:
+#'         prob_alpha is a vector of c(alpha_i,i+1, ..., alpha_i,j-1, alpha_i,j), where alpha_i,j = \eqn{P(\{Z_j \geq b_j\} \& \{\cap_{m=i+1}^{j-1} a_m \leq Z_m < b_m\} \mid Z_i = c)}.
+#'                                       prob_alpha_plus is a vector of c(alpha^+_i,i+1, ..., alpha^+_i,j-1, alpha^+_i,j), where alpha^+_i,j = \eqn{P(\{Z_j \geq b_j\} \& \{\cap_{m=i+1}^{j-1} Z_m < b_m\} \mid Z_i = c)}.
+#'                                       prob_beta is a vector of c(beta_i,i+1, ..., beta_i,j-1, beta_i,j) where beta_i,j = \eqn{P(\{Z_j \leq b_j\} \& \{\cap_{m=i+1}^{j-1} a_m \leq Z_m < b_m\} \mid Z_i = c)}.
 #' @examples
 #' library(gsDesign2)
-#'
-#' # Calculate conditional power under arbitrary theta and info
+#' library(dplyr)
+#' library(mvtnorm)
+#' # Example 1 ----
+#' # Calculate conditional power under arbitrary theta, info and lower/upper bound
 #' # In practice, the value of theta and info commonly comes from a design.
 #' # More examples are available at the pkgdown vignettes.
-#' gs_cp_npe(theta = c(.1, .2),
-#'           info = c(15, 35),
-#'           a = 1.5, b = 1.96)
+#' gs_cp_npe(theta = c(0.1, 0.2, 0.3),
+#'            t = c(0.15, 0.35, 0.6),
+#'            info = c(15, 35, 60),
+#'            a = c(-0.5, -0.5),
+#'            b = c(1.8, 2.1),
+#'            c = 1.5,
+#'            simple_cp = TRUE)
 gs_cp_npe <- function(theta = NULL,
+                      t = NULL,
                       info = NULL,
-                      a = NULL, b = NULL
-                      ) {
+                      a = NULL,
+                      b = NULL,
+                      c = NULL,
+                      simple_cp = FALSE){
+
   # ----------------------------------------- #
   #       input checking                      #
   # ----------------------------------------- #
-  # check theta
-  if (is.null(theta)) {
-    stop("Please provide theta (arbitrary treatment effect) to calculate conditional power.")
-  } else if (length(theta) == 1) {
-    theta <- rep(theta, 2)
+  if(simple_cp){
+
+    # check theta
+    if(is.null(theta) || any(is.na(theta))){
+      stop("Please provide theta (arbitrary treatment effect) to calculate conditional power.")
+    }else if(length(theta) == 1){
+      theta <- rep(theta, 2)
+    }else if(length(theta > 2)){
+      theta <- theta[1:2]
+      message("The first two elements of theta are used.")
+    }
+    # check info
+    if(is.null(info) || any(is.na(info))){
+      stop("Please provide info (statistical information given the treatment effect theta) to calculate conditional power.")
+    }
+    # check b
+    if (is.null(b) || !is.numeric(b) || length(b) != 1 || any(is.na(b))){
+      stop("Argument 'b' must be a numeric scalar and not NA.")
+    }
+    # Check c
+    if (is.null(c) || !is.numeric(c) || length(c) != 1 || any(is.na(c))){
+      stop("Argument 'c' must be a numeric scalar and not NA.")
+    }
   }
 
-  # check info
-  if (is.null(info)) {
-    stop("Please provide info (statistical information given the treatment effect theta) to calculate conditional power.")
+  if(!simple_cp){
+
+    # check theta
+    if(is.null(theta) || any(is.na(theta))){
+      stop("Please provide theta (arbitrary treatment effect) to calculate conditional power.")
+    }
+    # check t
+    if(is.null(t) || any(is.na(t))){
+      stop("Please provide t (information fraction) to calculate conditional power.")
+    }
+    # check info
+    if(is.null(info) || any(is.na(info))){
+      stop("Please provide info (statistical information given the treatment effect theta) to calculate conditional power.")
+    }
+    # check a
+    if (is.null(a) || !is.numeric(a) || any(is.na(c))){
+      stop("Argument 'a' must be numeric and not NA.")
+    }
+    # Check b
+    if (is.null(b) || !is.numeric(b) || any(is.na(b))){
+      stop("Argument 'b' must be numeric and not NA.")
+    }
+    # Check c
+    if (is.null(c) || !is.numeric(c) || length(c) != 1 || any(is.na(c))){
+      stop("Argument 'c' must be a numeric scalar and not NA.")
+    }
+
   }
 
-  check_info(info)
 
-  # ----------------------------------------- #
-  #   calculate conditional power under theta #
-  # ----------------------------------------- #
+  # --------------------- #
+  #  Call the cp function #
+  # --------------------- #
+  if(simple_cp){
+    cp = gs_cp_npe1(theta = theta, info = info, b = b, c = c)
+  }else{
+    cp = gs_cp_npe2(theta = theta, t = t, info = info, a = a, b = b, c = c)
+  }
 
-  t <- info[1] / info[2]
-  numerator1 <- b -  a * sqrt(t)
-  numerator2 <- theta[2] * sqrt(info[2]) - theta[1] * sqrt(t * info[1])
-  denominator <- sqrt(1 - t)
-  conditional_power <- pnorm((numerator1 - numerator2)  / denominator, lower.tail = FALSE)
+  return(cp)
 
-  return(conditional_power)
 }
+
+
+
+

R/gs_cp_npe1.R

@@ -0,0 +1,81 @@
+#  Copyright (c) 2025 Merck & Co., Inc., Rahway, NJ, USA and its affiliates.
+#  All rights reserved.
+#
+#  This file is part of the gsDesign2 program.
+#
+#  gsDesign2 is free software: you can redistribute it and/or modify
+#  it under the terms of the GNU General Public License as published by
+#  the Free Software Foundation, either version 3 of the License, or
+#  (at your option) any later version.
+#
+#  This program is distributed in the hope that it will be useful,
+#  but WITHOUT ANY WARRANTY; without even the implied warranty of
+#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#  GNU General Public License for more details.
+#
+#  You should have received a copy of the GNU General Public License
+#  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+#' Conditional power computation with non-constant effect size
+#'
+#' @details
+#' We assume \eqn{Z_1} and \eqn{Z_2} are the z-values at an interim analysis and later analysis, respectively.
+#' We assume further \eqn{Z_1} and \eqn{Z_2} are bivariate normal with standard group sequential assumptions
+#' on independent increments where for \eqn{i=1,2}
+#' \deqn{E(Z_i) = \theta_i\sqrt{I_i}}
+#' \deqn{Var(Z_i) = 1/I_i}
+#' \deqn{Cov(Z_1, Z_2) = t \equiv I_1/I_2}
+#' where \eqn{\theta_1, \theta_2} are real values and \eqn{0<I_1<I_2}.
+#' See https://merck.github.io/gsDesign2/articles/story-npe-background.html for assumption details.
+#' Returned value is
+#' \deqn{P(Z_2 > b \mid Z_1 = c) = 1 - \Phi\left(\frac{b - \sqrt{t}c - \sqrt{I_2}(\theta_2 - \theta_1\sqrt{t})}{\sqrt{1 - t}}\right)}
+#'
+#' @param theta A vector of length two, which specifies the natural parameter for treatment effect.
+#'              The first element of `theta` is the treatment effect of an interim analysis i.
+#'              The second element of `theta` is the treatment effect of a future analysis j.
+#' @param info A vector of length two, which specifies the statistical information under the treatment effect `theta`.
+#' @param c Interim z-value at analysis i (scalar).
+#' @param b Future target z-value at analysis j (scalar).
+#' @return A scalar with the conditional power \eqn{P(Z_2 > b \mid Z_1 = c)}.
+#' @export
+#'
+#' @examples
+#' library(gsDesign2)
+#'
+#' # Calculate conditional power under arbitrary theta and info
+#' # In practice, the value of theta and info commonly comes from a design.
+#' # More examples are available at the pkgdown vignettes.
+#' gs_cp_npe1(theta = c(.1, .2),
+#'           info = c(15, 35),
+#'           c = 1.5, b = 1.96)
+gs_cp_npe1 <- function(theta = NULL,
+                      info = NULL,
+                      c = NULL, b = NULL
+                      ) {
+  # ----------------------------------------- #
+  #       input checking                      #
+  # ----------------------------------------- #
+  # check theta
+  if (is.null(theta)) {
+    stop("Please provide theta (arbitrary treatment effect) to calculate conditional power.")
+  } else if (length(theta) == 1) {
+    theta <- rep(theta, 2)
+  }
+
+  # check info
+  if (is.null(info)) {
+    stop("Please provide info (statistical information given the treatment effect theta) to calculate conditional power.")
+  }
+
+  # ----------------------------------------- #
+  #   calculate conditional power under theta #
+  # ----------------------------------------- #
+
+  t <- info[1] / info[2]
+  numerator1 <- b -  c * sqrt(t)
+  numerator2 <- theta[2] * sqrt(info[2]) - theta[1] * sqrt(t * info[1])
+  denominator <- sqrt(1 - t)
+  conditional_power <- pnorm((numerator1 - numerator2)  / denominator, lower.tail = FALSE)
+
+  return(conditional_power)
+}