diff --git a/DESCRIPTION b/DESCRIPTION
index 327cc9ea..9775fec8 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -68,4 +68,4 @@ VignetteBuilder:
LinkingTo:
Rcpp
Roxygen: list(markdown = TRUE)
-RoxygenNote: 7.3.2
+RoxygenNote: 7.3.3
diff --git a/NAMESPACE b/NAMESPACE
index 702ca757..c56190c0 100644
--- a/NAMESPACE
+++ b/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)
+export(gs_cp_simple)
export(gs_create_arm)
export(gs_design_ahr)
export(gs_design_combo)
diff --git a/R/gs_cp.R b/R/gs_cp.R
new file mode 100644
index 00000000..1f01791e
--- /dev/null
+++ b/R/gs_cp.R
@@ -0,0 +1,262 @@
+# Copyright (c) 2026 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 .
+
+#' Conditional power computation with non-constant effect size for non-/crossing an upper boundary at analysis j given observed Z value at analysis i
+#'
+#' @details
+#' We assume \eqn{Z_i, i = 1, ..., K} are the z-statistics 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{Cov(Z_i, Z_j) = I_i/I_j}.
+#' See https://merck.github.io/gsDesign2/articles/story-npe-background.html for assumption details.
+#' Returned value is 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 = z_i)}.
+#' \deqn{P(\{Z_j \geq b_j\} \& \{\cap_{m=i+1}^{j-1} Z_m < b_m\} \mid Z_i = z_i)}.
+#' \deqn{P(\{Z_j \leq b_j\} \& \{\cap_{m=i+1}^{j-1} a_m \leq Z_m < b_m\} \mid Z_i = z_i)}.
+#'
+#' @param x An object of type gsDesign2.
+#' @param theta Optional numeric vector with length \eqn{j-i+1}, which specifies the natural parameter for treatment effect of interim analysis \eqn{i} through analysis \eqn{j}.
+#' @param i Index of current analysis, with default of 1.
+#' @param zi Numeric scalar z-value observed at analysis \eqn{i}.
+#' @return A list of conditional powers: prob_alpha is a numeric vector of (\eqn{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 = z_i)}.
+#' prob_alpha_plus is a numeric vector of (\eqn{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 = z_i)}.
+#' prob_beta is a numeric vector of (\eqn{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 = z_i)}.
+#' @export
+#'
+#' @examples
+#' library(gsDesign2)
+#' library(gsDesign)
+#' library(dplyr)
+#' library(mvtnorm)
+#' # Example 1
+#' # original design ----
+#' enroll_rate <- define_enroll_rate(duration = c(2, 2, 2, 18),
+#' rate = c(1, 2, 3, 4))
+#' fail_rate <- define_fail_rate(duration = c(3, Inf),
+#' fail_rate = log(2) / 10,
+#' dropout_rate = 0.001,
+#' hr = c(1, 0.7))
+#' x <- gs_design_ahr(enroll_rate = enroll_rate, fail_rate = fail_rate,
+#' alpha = 0.025, beta = 0.1, ratio = 1,
+#' info_frac = c(0.4, 0.6, 0.8, 1), analysis_time = 30,
+#' binding = FALSE,
+#' upper = gs_spending_bound,
+#' upar = list(sf = sfLDOF, total_spend = 0.025, param = NULL),
+#' lower = gs_spending_bound,
+#' lpar = list(sf = sfLDOF, total_spend = 0.1),
+#' h1_spending = TRUE,
+#' test_lower = TRUE,
+#' info_scale = "h0_h1_info") |> to_integer()
+#'
+#' # calculate conditional power
+#' # case 1: currently at IA1, compute conditional power at IA2, IA3 and FA
+#' gs_cp(x = x, i = 1, zi = -gsDesign::hrn2z(hr = 0.8, n = 150+180, ratio = 1))
+#'
+#' # case 2: currently at IA1, compute conditional power at IA2, IA3 and FA, with user-input theta
+#' gs_cp(x = x, i = 1, zi = -gsDesign::hrn2z(hr = 0.8, n = 150+180, ratio = 1))
+#'
+gs_cp <- function(x = NULL, theta = NULL, i = 1, zi = NULL){
+ # ------------------------------ #
+ # Input checking
+ # ------------------------------ #
+ # check x
+ if(is.null(x)){
+ stop("Please provide a design object of type gsDesign2 to calculate conditional power.")
+ }
+
+ # check theta
+ if(!is.null(theta) && length(theta) != dim(x$analysis)[1] - i + 1){
+ stop("Please provide the vector of theta from current analysis i to final analysis.")
+ }
+ # check i
+ if(i > dim(x$analysis)[1]){
+ stop("i cannot be larger than the number of analysis of design x.")
+ }
+
+ # check zi
+ if (is.null(zi) || length(zi) != 1 || !is.finite(zi)) {
+ stop("Please provide a finite scalar Z-score `zi` for analysis `i`.")
+ }
+ # ------------------------------ #a
+ # Initialization
+ # ------------------------------ #
+
+ # the conditional power is calculated from analysis i to analysis j
+ # the analysis j is decided by the length of b (efficacy bound)
+ # let D_m = B_m - B_i, where m = i+1, i+2, ..., j
+
+ k_max <- dim(x$analysis)[1]
+ n_future_analysis <- k_max - i
+ prob_alpha <- rep(0, n_future_analysis)
+ prob_alpha_plus <- rep(0, n_future_analysis)
+ prob_beta <- rep(0, n_future_analysis)
+
+ # futility bound from analysis i+1 to analysis j
+ if("lower" %in% unique(x$bound$bound)){
+ a <- x$bound$z[x$bound$bound == "lower"][(i+1):k_max]
+ }else{
+ a <- rep(-Inf, k_max)
+ }
+
+ # efficacy bound from analysis i+1 to analysis j
+ if("upper" %in% unique(x$bound$bound)){
+ b <- x$bound$z[x$bound$bound == "upper"][(i+1):k_max]
+ }else{
+ b <- rep(Inf, k_max)
+ }
+
+ # statistical information from analysis i to analysis j
+ info <- x$analysis$info[i:k_max]
+ # information fraction from analysis i to analysis j
+ t <- x$analysis$info_frac[i:k_max]
+
+ if(is.null(theta)){theta <- x$analysis$theta[i:k_max]}
+
+ for(y in seq_len(n_future_analysis)){
+ # y ranges from 1 to j-i, represents cases for alpha_{i,i+1}, ..., {alpha_i,j-1}, alpha_{i,j}
+ # y is the increment from i
+
+ # ------------------------------ #
+ # Build the asymptotic
+ # mean of B_k - B_i
+ # vector of length y
+ # ------------------------------ #
+ mu <- sapply(seq_len(y), function(k){
+ idx <- k + 1 # i.e., start from i+1
+ theta[idx] * sqrt(t[idx] * info[idx]) - theta[1] * sqrt(t[1] * info[1]) #first element of `theta` is the treatment effect of IA i.
+ })
+
+ # ---------------------------------- #
+ # Build the asymptotic
+ # covariance matrix of B_y - B_i
+ # ---------------------------------- #
+ cov_matrix <- matrix(0, nrow = y, ncol = y)
+
+ for(k in seq_len(y)) {
+ for(l in seq_len(y)) {
+ cov_matrix[k, l] <- t[1 + min(k, l)] - t[1]
+ }
+ }
+
+ # ----------------------------------------------------------------------------------------- #
+ # Calculate Lower/upper bounds for alpha
+ # first cross efficacy bound at analysis y given no efficacy/futility bound crossing before
+ # ----------------------------------------------------------------------------------------- #
+ # Integration limits: D_m = B_m - B_i
+ # for D_{i+1},...,D_{j-1} use [a, b); for D_j use [b_j, +Inf)
+ # lower bound
+ lower_alpha <- rep(0, y)
+ if(y == 1){
+ lower_alpha[y] <- b[y] * sqrt(t[y + 1]) - zi * sqrt(t[1])
+ }else{
+ for (m in 1:(y - 1)) {
+ lower_alpha[m] <- a[m] * sqrt(t[m + 1]) - zi * sqrt(t[1])
+ }
+ lower_alpha[y] <- b[y] * sqrt(t[y + 1]) - zi * sqrt(t[1])
+ }
+
+ # upper bound
+ upper_alpha <- rep(0, y)
+ if(y == 1){
+ upper_alpha[y] = Inf
+ }else{
+ for(m in 1:(y - 1)){
+ upper_alpha[m] <- b[m] * sqrt(t[m + 1]) - zi * sqrt(t[1])
+ }
+ upper_alpha[y] <- Inf
+ }
+
+ # Compute multivariate normal probability
+ prob_alpha[y] <- mvtnorm::pmvnorm(
+ lower = lower_alpha,
+ upper = upper_alpha,
+ mean = mu,
+ sigma = cov_matrix)[1]
+
+ # --------------------------------------------------------------------------------- #
+ # Calculate Lower/upper bounds for alpha_plus
+ # first cross efficacy bound at analysis y given no efficacy bound crossing before
+ # --------------------------------------------------------------------------------- #
+ # Integration limits: D_m = B_m - B_i
+ # for D_{i+1},...,D_{j-1} use (-Inf, bm); for D_j use [b_j, +Inf)
+ # lower bound
+ lower_alpha_plus <- rep(0, y)
+ if(y == 1){
+ lower_alpha_plus[y] = b[y] * sqrt(t[y + 1]) - zi * sqrt(t[1])
+ }else{
+ lower_alpha_plus <- rep(-Inf, y - 1)
+ lower_alpha_plus[y] <- b[y] * sqrt(t[y + 1]) - zi * sqrt(t[1])
+ }
+
+ # upper bound
+ upper_alpha_plus <- rep(0, y)
+ if(y == 1){
+ upper_alpha_plus[y] <- Inf
+ }else{
+ for(m in 1:(y - 1)){
+ upper_alpha_plus[m] <- b[m] * sqrt(t[m + 1]) - zi * sqrt(t[1])
+ }
+ upper_alpha_plus[y] <- Inf
+ }
+
+ prob_alpha_plus[y] <- mvtnorm::pmvnorm(
+ lower = lower_alpha_plus,
+ upper = upper_alpha_plus,
+ mean = mu,
+ sigma = cov_matrix)[1]
+
+ # ---------------------------------------------------------------------------------------------- #
+ # Calculate Lower/upper bounds for beta
+ # First crossing the futility bound at analysis y given no efficacy/futility bound crossing before
+ # ---------------------------------------------------------------------------------------------- #
+ # Integration limits: D_m = B_m - B_i
+ # for D_{i+1},...,D_{j-1} use [a, b); for D_j use [-Inf, a_j)
+ # lower bound
+ lower_beta <- rep(0, y)
+ if(y == 1){
+ lower_beta[y] <- -Inf
+ }else{
+ for(m in 1:(y - 1)){
+ lower_beta[m] <- a[m] * sqrt(t[m + 1]) - zi * sqrt(t[1])
+ }
+ lower_beta[y] <- -Inf
+ }
+
+ # upper bound
+ upper_beta <- rep(0, y)
+ if(y == 1){
+ upper_beta[y] <- b[y] * sqrt(t[y + 1]) - zi * sqrt(t[1])
+ }else{
+ for(m in 1:y){
+ upper_beta[m] <- a[m] * sqrt(t[m + 1]) - zi * sqrt(t[1])
+ }
+ }
+
+ prob_beta[y] <- mvtnorm::pmvnorm(
+ lower = lower_beta,
+ upper = upper_beta,
+ mean = mu,
+ sigma = cov_matrix)[1]
+
+ }
+
+ return(list(prob_alpha = prob_alpha,
+ prob_alpha_plus = prob_alpha_plus,
+ prob_beta = prob_beta))
+
+}
diff --git a/R/gs_cp_npe.R b/R/gs_cp_npe1.R
similarity index 66%
rename from R/gs_cp_npe.R
rename to R/gs_cp_npe1.R
index cf0f2a93..712541f1 100644
--- a/R/gs_cp_npe.R
+++ b/R/gs_cp_npe1.R
@@ -19,25 +19,24 @@
#' 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}
+#' Suppose there are \eqn{K} analyses. Let \eqn{Z_i} and \eqn{Z_j} be the Z-statistics at a current analysis \eqn{i} and a future analysis \eqn{j}, respectively.
+#' Let's denote the statistical information at the \eqn{i}-th analysis as \eqn{I_i}.
+#' We further assume \eqn{Z_i} and \eqn{Z_j} are bivariate normal with standard group sequential assumptions on independent increments, then
#' \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 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)}
+#' \deqn{P(Z_j > z_j \mid Z_i = z_i) = 1 - \Phi\left(\frac{z_j - \sqrt{t}z_i - \sqrt{I_j}(\theta_j - \theta_i\sqrt{t})}{\sqrt{1 - t}}\right)}
#'
-#' @param theta A vector of length two, which specifies the natural parameter for treatment effect.
+#' @param theta A numeric 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
+#' @param info A vector of length two, which specifies the statistical information under the treatment effect `theta`.
+#' @param zi Numeric scalar z-value observed at analysis \eqn{i}.
+#' @param zj Numeric scalar at the future analysis \eqn{j}.
+#' @return A scalar with the conditional power \eqn{P(Z_j > z_i \mid Z_i = z_i)}.
+#' @noRd
#'
#' @examples
#' library(gsDesign2)
@@ -45,12 +44,9 @@
#' # 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_npe(theta = c(.1, .2),
-#' info = c(15, 35),
-#' a = 1.5, b = 1.96)
-gs_cp_npe <- function(theta = NULL,
- info = NULL,
- a = NULL, b = NULL
+#' gs_cp_npe1(theta = c(.1, .2), info = c(15, 35), zi = 1.5, zj = 1.96)
+#'
+gs_cp_npe1 <- function(theta = NULL, info = NULL, zi = NULL, zj = NULL
) {
# ----------------------------------------- #
# input checking #
@@ -67,14 +63,12 @@ gs_cp_npe <- function(theta = NULL,
stop("Please provide info (statistical information given the treatment effect theta) to calculate conditional power.")
}
- check_info(info)
-
# ----------------------------------------- #
# calculate conditional power under theta #
# ----------------------------------------- #
t <- info[1] / info[2]
- numerator1 <- b - a * sqrt(t)
+ numerator1 <- zj - zi * 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)
diff --git a/R/gs_cp_npe2.R b/R/gs_cp_npe2.R
new file mode 100644
index 00000000..e8e54ac5
--- /dev/null
+++ b/R/gs_cp_npe2.R
@@ -0,0 +1,296 @@
+# Copyright (c) 2024 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 .
+
+#' Conditional power computation with non-constant effect size for non-/crossing an upper boundary at analysis j given observed Z value at analysis i
+#'
+#' @details
+#' We assume \eqn{Z_i, i = 1, ..., K} are the z-statistics 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{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 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 = z_i)}.
+#' \deqn{P(\{Z_j \geq b_j\} \& \{\cap_{m=i+1}^{j-1} Z_m < b_m\} \mid Z_i = z_i)}.
+#' \deqn{P(\{Z_j \leq b_j\} \& \{\cap_{m=i+1}^{j-1} a_m \leq Z_m < b_m\} \mid Z_i = z_i)}.
+#'
+#' @param theta 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 A vector of j-i+1, which specifies the information fraction under the treatment effect `theta`.
+#' @param info A vector of j-i+1, which specifies the statistical information under the treatment effect `theta`.
+#' @param a A vector of length j-i, which specifies the futility bounds from analysis i+1 to analysis j.
+#' @param b A vector of length j-i, which specifies the efficacy bounds from analysis i+1 to analysis j.
+#' @param zi Numeric scalar z-value observed at analysis \eqn{i}.
+#' @return A list of conditional powers: prob_alpha is a numeric vector of (\eqn{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 = z_i)}.
+#' prob_alpha_plus is a numeric vector of (\eqn{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 = z_i)}.
+#' prob_beta is a numeric vector of (\eqn{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 = z_i)}.
+#' @noRd
+#'
+#' @examples
+#' library(gsDesign2)
+#' library(gsDesign)
+#' 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_npe2(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),
+#' zi = 1.5)
+#'
+#' # Example 2
+#' # original design ----
+#' enroll_rate <- define_enroll_rate(duration = c(2, 2, 2, 18),
+#' rate = c(1, 2, 3, 4))
+#' fail_rate <- define_fail_rate(duration = c(3, Inf),
+#' fail_rate = log(2) / 10,
+#' dropout_rate = 0.001,
+#' hr = c(1, 0.7))
+#' x <- gs_design_ahr(enroll_rate = enroll_rate, fail_rate = fail_rate,
+#' alpha = 0.025, beta = 0.1, ratio = 1,
+#' info_frac = c(0.4, 0.6, 0.8, 1), analysis_time = 30,
+#' binding = FALSE,
+#' upper = gs_spending_bound,
+#' upar = list(sf = sfLDOF, total_spend = 0.025, param = NULL),
+#' lower = gs_b,
+#' lpar = c(-Inf, -Inf, -Inf),
+#' h1_spending = TRUE,
+#' test_lower = FALSE,
+#' info_scale = "h0_h1_info") |> to_integer()
+#'
+#' # calculate conditional power
+#' # case 1: currently at IA1, compute conditional power at IA2
+#' gs_cp_npe2(# IA1 and IA2's theta
+#' theta = x$analysis$theta[1:2],
+#' # IA1 and IA2's information fraction
+#' t = x$analysis$info_frac[1:2],
+#' # IA1 and IA2's statistical information
+#' info = x$analysis$info[1:2],
+#' # IA2's futility bound
+#' a = -Inf,
+#' # IA2's efficacy bound
+#' b = x$bound$z[x$bound$bound == "upper" & x$bound$analysis == 2],
+#' # IA1's Z-score
+#' zi = -gsDesign::hrn2z(hr = 0.8, n = 150+180, ratio = 1))
+#'
+#' # case 2: currently at IA1, compute conditional power at IA2 and IA3
+#' gs_cp_npe2(# IA1, IA2 and IA3's theta
+#' theta = x$analysis$theta[1:3],
+#' # IA1, IA2 and IA3's information fraction
+#' t = x$analysis$info_frac[1:3],
+#' # IA1, IA2, and IA3's statistical information
+#' info = x$analysis$info[1:3],
+#' # IA2 and IA3's futility bound
+#' a = c(-Inf, Inf),
+#' # IA2 and IA3's efficacy bound
+#' b = x$bound$z[x$bound$bound == "upper" & x$bound$analysis %in% c(2, 3)],
+#' # IA1's Z-score
+#' zi = -gsDesign::hrn2z(hr = 0.8, n = 150+180, ratio = 1))
+#'
+#' # case 3: currently at IA1, compute conditional power at IA2, IA3 and FA
+#' gs_cp_npe2(# IA1, IA2, IA3 and FA's theta
+#' theta = x$analysis$theta[1:4],
+#' # IA1, IA2, IA3 and FA's information fraction
+#' t = x$analysis$info_frac[1:4],
+#' # IA1, IA2, IA3 and FA's statistical information
+#' info = x$analysis$info[1:4],
+#' # IA2, IA3 and FA's futility bound
+#' a = c(-Inf, -Inf, -Inf),
+#' # IA2, IA3 and FA's efficacy bound
+#' b = x$bound$z[x$bound$bound == "upper" & x$bound$analysis %in% c(2, 3, 4)],
+#' # IA1's Z-score
+#' zi = -gsDesign::hrn2z(hr = 0.8, n = 150+180, ratio = 1))
+gs_cp_npe2 <- function(theta = NULL,
+ t = NULL,
+ info = NULL,
+ a = NULL,
+ b = NULL,
+ zi = NULL){
+ # ------------------------------ #
+ # Input checking
+ # ------------------------------ #
+ if(length(theta) != length(t))
+ stop("The input of theta should have the same length of the input of t.")
+
+ if(length(theta) != length(info))
+ stop("The input of theta should have the same length of the input of info.")
+
+ if(length(theta) - length(b) != 1)
+ stop("The length of theta should equal to length of b + 1. ")
+
+ if(length(b) != length(a))
+ stop("The length of b should equal to length of a. ")
+
+ # check zi
+ if (is.null(zi) || length(zi) != 1 || !is.finite(zi)) {
+ stop("Please provide a finite scalar Z-score `zi` for analysis `i`.")
+ }
+ # ------------------------------ #
+ # Initialization
+ # ------------------------------ #
+
+ # the conditional power is calculated from analysis i to analysis j
+ # the analysis j is decided by the length of b (efficacy bound)
+ # let D_m = B_m - B_i, where m = i+1, i+2, ..., j
+ n_future_analysis <- length(b) # = j-i
+ prob_alpha <- rep(0, n_future_analysis)
+ prob_alpha_plus <- rep(0, n_future_analysis)
+ prob_beta <- rep(0, n_future_analysis)
+
+ for(x in 1:n_future_analysis){
+ # x ranges from 1 to j-i, represents cases for alpha_{i,i+1}, ..., {alpha_i,j-1}, alpha_{i,j}
+ # x is the increment from i
+
+ # ------------------------------ #
+ # Build the asymptotic
+ # mean of B_x - B_i
+ # vector of length x
+ # ------------------------------ #
+
+ mu <- sapply(seq_len(x), function(k){
+ idx <- k + 1 # i.e., start from i+1
+ theta[idx] * sqrt(t[idx] * info[idx]) - theta[1] * sqrt(t[1] * info[1]) #first element of `theta` is the treatment effect of IA i.
+ })
+
+ # ---------------------------------- #
+ # Build the asymptotic
+ # covariance matrix of B_x - B_i
+ # ---------------------------------- #
+ cov_matrix <- matrix(0, nrow = x, ncol = x)
+
+ for(k in seq_len(x)) {
+ for(l in seq_len(x)) {
+ cov_matrix[k, l] <- t[1 + min(k, l)] - t[1]
+ }
+ }
+
+ # ----------------------------------------------------------------------------------------- #
+ # Calculate Lower/upper bounds for alpha
+ # first cross efficacy bound at analysis x given no efficacy/futility bound crossing before
+ # ----------------------------------------------------------------------------------------- #
+ # Integration limits: D_m = B_m - B_i
+ # for D_{i+1},...,D_{j-1} use [a, b); for D_j use [b_j, +Inf)
+ # lower bound
+ lower_alpha <- rep(0, x)
+ if(x == 1){
+ lower_alpha[x] <- b[x] * sqrt(t[x + 1]) - zi * sqrt(t[1])
+ }else{
+ for (m in 1:(x - 1)) {
+ lower_alpha[m] <- a[m] * sqrt(t[m + 1]) - zi * sqrt(t[1])
+ }
+ lower_alpha[x] <- b[x] * sqrt(t[x + 1]) - zi * sqrt(t[1])
+ }
+
+ # upper bound
+ upper_alpha <- rep(0, x)
+ if(x == 1){
+ upper_alpha[x] = Inf
+ }else{
+ for(m in 1:(x - 1)){
+ upper_alpha[m] <- b[m] * sqrt(t[m + 1]) - zi * sqrt(t[1])
+ }
+ upper_alpha[x] <- Inf
+ }
+
+ # Compute multivariate normal probability
+ prob_alpha[x] <- mvtnorm::pmvnorm(
+ lower = lower_alpha,
+ upper = upper_alpha,
+ mean = mu,
+ sigma = cov_matrix)[1]
+
+ # --------------------------------------------------------------------------------- #
+ # Calculate Lower/upper bounds for alpha_plus
+ # first cross efficacy bound at analysis x given no efficacy bound crossing before
+ # --------------------------------------------------------------------------------- #
+ # Integration limits: D_m = B_m - B_i
+ # for D_{i+1},...,D_{j-1} use (-Inf, bm); for D_j use [b_j, +Inf)
+ # lower bound
+ lower_alpha_plus <- rep(0, x)
+ if(x == 1){
+ lower_alpha_plus[x] = b[x] * sqrt(t[x + 1]) - zi * sqrt(t[1])
+ }else{
+ lower_alpha_plus <- rep(-Inf, x - 1)
+ lower_alpha_plus[x] <- b[x] * sqrt(t[x + 1]) - zi * sqrt(t[1])
+ }
+
+ # upper bound
+ upper_alpha_plus <- rep(0, x)
+ if(x == 1){
+ upper_alpha_plus[x] <- Inf
+ }else{
+ for(m in 1:(x - 1)){
+ upper_alpha_plus[m] <- b[m] * sqrt(t[m + 1]) - zi * sqrt(t[1])
+ }
+ upper_alpha_plus[x] <- Inf
+ }
+
+ prob_alpha_plus[x] <- mvtnorm::pmvnorm(
+ lower = lower_alpha_plus,
+ upper = upper_alpha_plus,
+ mean = mu,
+ sigma = cov_matrix)[1]
+
+ # ---------------------------------------------------------------------------------------------- #
+ # Calculate Lower/upper bounds for beta
+ # First crossing the futility bound at analysis x given no efficacy/futility bound crossing before
+ # ---------------------------------------------------------------------------------------------- #
+ # Integration limits: D_m = B_m - B_i
+ # for D_{i+1},...,D_{j-1} use [a, b); for D_j use [-Inf, a_j)
+ # lower bound
+ lower_beta <- rep(0, x)
+ if(x == 1){
+ lower_beta[x] <- -Inf
+ }else{
+ for(m in 1:(x - 1)){
+ lower_beta[m] <- a[m] * sqrt(t[m + 1]) - zi * sqrt(t[1])
+ }
+ lower_beta[x] <- -Inf
+ }
+
+ # upper bound
+ upper_beta <- rep(0, x)
+ if(x == 1){
+ upper_beta[x] <- b[x] * sqrt(t[x + 1]) - zi * sqrt(t[1])
+ }else{
+ for(m in 1:x){
+ upper_beta[m] <- a[m] * sqrt(t[m + 1]) - zi * sqrt(t[1])
+ }
+ }
+
+ prob_beta[x] <- mvtnorm::pmvnorm(
+ lower = lower_beta,
+ upper = upper_beta,
+ mean = mu,
+ sigma = cov_matrix)[1]
+
+ }
+
+
+ return(list(prob_alpha = prob_alpha,
+ prob_alpha_plus = prob_alpha_plus,
+ prob_beta = prob_beta))
+
+}
diff --git a/R/gs_cp_simple.R b/R/gs_cp_simple.R
new file mode 100644
index 00000000..e734fde6
--- /dev/null
+++ b/R/gs_cp_simple.R
@@ -0,0 +1,191 @@
+# Copyright (c) 2026 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 .
+
+#' Conditional power computation with non-constant effect size
+#'
+#' @details
+#' Suppose there are \eqn{K} analyses. Let \eqn{Z_i} and \eqn{Z_j} be the Z-statistics at a current analysis \eqn{i} and a future analysis \eqn{j}, respectively.
+#' Let's denote the statistical information at the \eqn{i}-th analysis as \eqn{I_i}.
+#' We further assume \eqn{Z_i} and \eqn{Z_j} are bivariate normal with standard group sequential assumptions on independent increments, then
+#' \deqn{E(Z_i) = \theta_i\sqrt{I_i}}
+#' \deqn{Var(Z_i) = 1/I_i}
+#' \deqn{Cov(Z_i, Z_j) = t \equiv I_i/I_j}
+#' See https://merck.github.io/gsDesign2/articles/story-npe-background.html for assumption details.
+#' Returned value is
+#' \deqn{P(Z_j > b_j \mid Z_i = z_i) = 1 - \Phi\left(\frac{b_j - \sqrt{t}z_i - \sqrt{I_j}(\theta_j - \theta_i\sqrt{t})}{\sqrt{1 - t}}\right)}
+#' where \eqn{b_j} is the efficacy bound at analysis \eqn{j}, \eqn{z_i} is the observed Z-value at analysis \eqn{i},
+#' and \eqn{\theta_i} natural parameter for treatment effect at analysis \eqn{i}.
+#'
+#' @param x An object of type gsDesign2.
+#' @param theta Optional numeric vector of natural parameter for treatment effects from analysis \eqn{i} through final analysis \eqn{K}.
+#' If `NULL`, values are taken from `x$analysis$theta`.
+#' @param i Integer index for current analysis. Default is 1.
+#' @param zi Numeric scalar z-value observed at analysis \eqn{i}.
+#' @return A numeric vector with the conditional power
+#' \eqn{P(Z_j > b_j \mid Z_i = z_i)} for \eqn{j = i+1, ..., K}.
+#' @export
+#'
+#' @examples
+#' library(gsDesign2)
+#' library(gsDesign)
+#' # Calculate conditional power with optional user-input theta (if NULL, will come from the input design)
+#' alpha <- 0.025
+#' beta <- 0.1
+#' ratio <- 1
+#' enroll_rate <- define_enroll_rate(
+#' duration = c(2, 2, 10),
+#' rate = (1:3) / 3)
+#' fail_rate <- define_fail_rate(
+#' duration = Inf, fail_rate = log(2) / 9,
+#' hr = 0.6, dropout_rate = .0001)
+#' analysis_time <- c(12, 24, 36)
+#' ratio <- 1
+#' upper <- gs_spending_bound
+#' lower <- gs_b
+#' upar <- list(sf = sfLDOF, total_spend = alpha)
+#' lpar <- rep(-Inf, 3)
+#' x <- gs_design_ahr(
+#' enroll_rate = enroll_rate, fail_rate = fail_rate,
+#' alpha = alpha, beta = beta, ratio = ratio,
+#' info_scale = "h0_h1_info",
+#' info_frac = NULL,
+#' analysis_time = c(12, 24, 36),
+#' upper = upper, upar = upar, test_upper = TRUE,
+#' lower = lower, lpar = lpar, test_lower = FALSE,
+#' ) |>
+#' to_integer()
+#'
+#' # theta is NULL case
+#' gs_cp_simple(x = x, theta = NULL, i = 1, zi = 1.5)
+#'
+#' # User-input theta case
+#' gs_cp_simple(x = x, theta = c(0.1, 0.2, 0.3), i = 1, zi = 1.5)
+#'
+gs_cp_simple <- function(x = NULL, theta = NULL, i = 1, zi = NULL) {
+ # ----------------------------------------- #
+ # input checking #
+ # ----------------------------------------- #
+
+ # check x
+ if(is.null(x)){
+ stop("Please provide a design object of type gsDesign2 to calculate conditional power.")
+ }
+
+ if (!is.list(x) || is.null(x$analysis) || is.null(x$bound)) {
+ stop("`x` must contain `analysis` and `bound` components.")
+ }
+
+ analysis <- x$analysis
+ bound <- x$bound
+ k_max <- nrow(analysis)
+
+ if (is.null(k_max) || k_max < 1) {
+ stop("`x$analysis` must have at least one analysis row.")
+ }
+
+ # check theta
+ if(!is.null(theta) && length(theta) != k_max - i + 1){
+ stop("Please provide the vector of theta from current analysis i to final analysis.")
+ }
+
+ # check i
+ if (length(i) != 1 || !is.finite(i) || i != as.integer(i)) {
+ stop("`i` must be a single finite integer.")
+ }
+
+ if (i < 1 || i > k_max) {
+ stop("`i` must be between 1 and the max number of analyses in `x$analysis`.")
+ }
+
+ # check zi
+ if (is.null(zi) || length(zi) != 1 || !is.finite(zi)) {
+ stop("Please provide a finite scalar Z-score `zi` for analysis `i`.")
+ }
+
+ # ----------------------------------------- #
+ # Simple conditional power calculation #
+ # ----------------------------------------- #
+ n_analysis <- k_max - i + 1 # number of future analysis
+ conditional_power <- rep(0, max(n_analysis - 1, 1)) # if i = FA
+ info <- analysis$info
+ info0 <- analysis$info0
+ b <- bound$z[bound$bound == "upper"]
+
+ if(is.null(theta)){theta <- analysis$theta}
+
+ # cp function under H0 (use info0)
+ cp_func_h0 <- function(k){
+
+ t <- info0[i]/info0[i + k]
+ numerator1 <- b[i + k] - zi * sqrt(t)
+
+ if(is.null(theta)){
+ theta <- analysis$theta
+ numerator2 <- theta[i + k] * sqrt(info0[i + k]) - theta[i] * sqrt(t * info0[i])
+ }else{
+ numerator2 <- theta[1 + k] * sqrt(info0[i + k]) - theta[1] * sqrt(t * info0[i]) # starting index of theta is 1, which is different from the above
+ }
+
+ denominator <- sqrt(1 - t)
+ cp <- pnorm((numerator1 - numerator2) / denominator, lower.tail = FALSE)
+
+ return(cp)
+ }
+
+ # cp function under H1 (use info)
+ cp_func_h1 <- function(k){
+
+ t <- info[i]/info[i + k]
+ numerator1 <- b[i + k] - zi * sqrt(t)
+
+
+ if(is.null(theta)){
+ theta <- analysis$theta
+ numerator2 <- theta[i + k] * sqrt(info[i + k]) - theta[i] * sqrt(t * info[i])
+ }else{
+ numerator2 <- theta[1 + k] * sqrt(info[i + k]) - theta[1] * sqrt(t * info[i]) # starting index of theta is different from the above
+ }
+
+ denominator <- sqrt(1 - t)
+ cp <- pnorm((numerator1 - numerator2) / denominator, lower.tail = FALSE)
+
+ return(cp)
+ }
+
+ # ----------------------------------------- #
+ # Output #
+ # ----------------------------------------- #
+
+ if(n_analysis == 1){ # i.e., i = last analysis of x
+
+ conditional_power <- cp_func_h1(0) # no difference of using cp_func_h1 or cp_func_h0
+
+ }else{
+ if(sum(theta) == 0){
+ conditional_power <- vapply(1:(n_analysis - 1), cp_func_h0, numeric(1)) # cp under H0
+ }else{
+ conditional_power <- vapply(1:(n_analysis - 1), cp_func_h1, numeric(1)) # cp under H1
+ }
+ }
+
+ return(conditional_power)
+
+ }
+
+
+
diff --git a/R/gs_update_ahr.R b/R/gs_update_ahr.R
index 023bc07f..c4d8d6d7 100644
--- a/R/gs_update_ahr.R
+++ b/R/gs_update_ahr.R
@@ -302,7 +302,7 @@ gs_update_ahr <- function(
# Tidy outputs #
# ----------------------------------- #
ans <- list()
-
+
ans$design <- x$design
ans$enroll_rate <- x$enroll_rate
diff --git a/_pkgdown.yml b/_pkgdown.yml
index 59f7f563..e6ba5ae1 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -77,7 +77,8 @@ reference:
desc: >
Functions for conditional power.
contents:
- - gs_cp_npe
+ - gs_cp
+ - gs_cp_simple
- title: "Input definition"
desc: >
Helper functions to define inputs for study design.
diff --git a/man/gsDesign2-package.Rd b/man/gsDesign2-package.Rd
index 22ed6d5a..c31fdc74 100644
--- a/man/gsDesign2-package.Rd
+++ b/man/gsDesign2-package.Rd
@@ -45,7 +45,7 @@ Other contributors:
\item Yalin Zhu \email{yalin.zhu@outlook.com} [contributor]
\item Shiyu Zhang \email{shiyu.zhang@merck.com} [contributor]
\item Dickson Wanjau \email{dickson.wanjau@merck.com} [contributor]
- \item Merck & Co., Inc., Rahway, NJ, USA and its affiliates (02891sr49) [copyright holder]
+ \item Merck & Co., Inc., Rahway, NJ, USA and its affiliates (\href{https://ror.org/02891sr49}{ROR}) [copyright holder]
}
}
diff --git a/man/gs_cp.Rd b/man/gs_cp.Rd
new file mode 100644
index 00000000..e9739b24
--- /dev/null
+++ b/man/gs_cp.Rd
@@ -0,0 +1,69 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/gs_cp.R
+\name{gs_cp}
+\alias{gs_cp}
+\title{Conditional power computation with non-constant effect size for non-/crossing an upper boundary at analysis j given observed Z value at analysis i}
+\usage{
+gs_cp(x = NULL, theta = NULL, i = 1, zi = NULL)
+}
+\arguments{
+\item{x}{An object of type gsDesign2.}
+
+\item{theta}{Optional numeric vector with length \eqn{j-i+1}, which specifies the natural parameter for treatment effect of interim analysis \eqn{i} through analysis \eqn{j}.}
+
+\item{i}{Index of current analysis, with default of 1.}
+
+\item{zi}{Numeric scalar z-value observed at analysis \eqn{i}.}
+}
+\value{
+A list of conditional powers: prob_alpha is a numeric vector of (\eqn{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 = z_i)}.
+prob_alpha_plus is a numeric vector of (\eqn{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 = z_i)}.
+prob_beta is a numeric vector of (\eqn{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 = z_i)}.
+}
+\description{
+Conditional power computation with non-constant effect size for non-/crossing an upper boundary at analysis j given observed Z value at analysis i
+}
+\details{
+We assume \eqn{Z_i, i = 1, ..., K} are the z-statistics 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{Cov(Z_i, Z_j) = I_i/I_j}.
+See https://merck.github.io/gsDesign2/articles/story-npe-background.html for assumption details.
+Returned value is 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 = z_i)}.
+\deqn{P(\{Z_j \geq b_j\} \& \{\cap_{m=i+1}^{j-1} Z_m < b_m\} \mid Z_i = z_i)}.
+\deqn{P(\{Z_j \leq b_j\} \& \{\cap_{m=i+1}^{j-1} a_m \leq Z_m < b_m\} \mid Z_i = z_i)}.
+}
+\examples{
+library(gsDesign2)
+library(gsDesign)
+library(dplyr)
+library(mvtnorm)
+# Example 1
+# original design ----
+enroll_rate <- define_enroll_rate(duration = c(2, 2, 2, 18),
+ rate = c(1, 2, 3, 4))
+fail_rate <- define_fail_rate(duration = c(3, Inf),
+ fail_rate = log(2) / 10,
+ dropout_rate = 0.001,
+ hr = c(1, 0.7))
+x <- gs_design_ahr(enroll_rate = enroll_rate, fail_rate = fail_rate,
+ alpha = 0.025, beta = 0.1, ratio = 1,
+ info_frac = c(0.4, 0.6, 0.8, 1), analysis_time = 30,
+ binding = FALSE,
+ upper = gs_spending_bound,
+ upar = list(sf = sfLDOF, total_spend = 0.025, param = NULL),
+ lower = gs_spending_bound,
+ lpar = list(sf = sfLDOF, total_spend = 0.1),
+ h1_spending = TRUE,
+ test_lower = TRUE,
+ info_scale = "h0_h1_info") |> to_integer()
+
+# calculate conditional power
+# case 1: currently at IA1, compute conditional power at IA2, IA3 and FA
+gs_cp(x = x, i = 1, zi = -gsDesign::hrn2z(hr = 0.8, n = 150+180, ratio = 1))
+
+# case 2: currently at IA1, compute conditional power at IA2, IA3 and FA, with user-input theta
+gs_cp(x = x, i = 1, zi = -gsDesign::hrn2z(hr = 0.8, n = 150+180, ratio = 1))
+
+}
diff --git a/man/gs_cp_npe.Rd b/man/gs_cp_npe.Rd
deleted file mode 100644
index b4823984..00000000
--- a/man/gs_cp_npe.Rd
+++ /dev/null
@@ -1,47 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/gs_cp_npe.R
-\name{gs_cp_npe}
-\alias{gs_cp_npe}
-\title{Conditional power computation with non-constant effect size}
-\usage{
-gs_cp_npe(theta = NULL, info = NULL, a = NULL, b = NULL)
-}
-\arguments{
-\item{theta}{A vector of length two, which specifies the natural parameter for treatment effect.
-The first element of \code{theta} is the treatment effect of an interim analysis i.
-The second element of \code{theta} is the treatment effect of a future analysis j.}
-
-\item{info}{A vector of two, which specifies the statistical information under the treatment effect \code{theta}.}
-
-\item{a}{Interim z-value at analysis i (scalar).}
-
-\item{b}{Future target z-value at analysis j (scalar).}
-}
-\value{
-A scalar with the conditional power \eqn{P(Z_2>b\mid Z_1=a)}.
-}
-\description{
-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 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)}
-}
-\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_npe(theta = c(.1, .2),
- info = c(15, 35),
- a = 1.5, b = 1.96)
-}
diff --git a/man/gs_cp_simple.Rd b/man/gs_cp_simple.Rd
new file mode 100644
index 00000000..8e8f3dec
--- /dev/null
+++ b/man/gs_cp_simple.Rd
@@ -0,0 +1,75 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/gs_cp_simple.R
+\name{gs_cp_simple}
+\alias{gs_cp_simple}
+\title{Conditional power computation with non-constant effect size}
+\usage{
+gs_cp_simple(x = NULL, theta = NULL, i = 1, zi = NULL)
+}
+\arguments{
+\item{x}{An object of type gsDesign2.}
+
+\item{theta}{Optional numeric vector of natural parameter for treatment effects from analysis \eqn{i} through final analysis \eqn{K}.
+If \code{NULL}, values are taken from \code{x$analysis$theta}.}
+
+\item{i}{Integer index for current analysis. Default is 1.}
+
+\item{zi}{Numeric scalar z-value observed at analysis \eqn{i}.}
+}
+\value{
+A numeric vector with the conditional power
+\eqn{P(Z_j > b_j \mid Z_i = z_i)} for \eqn{j = i+1, ..., K}.
+}
+\description{
+Conditional power computation with non-constant effect size
+}
+\details{
+Suppose there are \eqn{K} analyses. Let \eqn{Z_i} and \eqn{Z_j} be the Z-statistics at a current analysis \eqn{i} and a future analysis \eqn{j}, respectively.
+Let's denote the statistical information at the \eqn{i}-th analysis as \eqn{I_i}.
+We further assume \eqn{Z_i} and \eqn{Z_j} are bivariate normal with standard group sequential assumptions on independent increments, then
+\deqn{E(Z_i) = \theta_i\sqrt{I_i}}
+\deqn{Var(Z_i) = 1/I_i}
+\deqn{Cov(Z_i, Z_j) = t \equiv I_i/I_j}
+See https://merck.github.io/gsDesign2/articles/story-npe-background.html for assumption details.
+Returned value is
+\deqn{P(Z_j > b_j \mid Z_i = z_i) = 1 - \Phi\left(\frac{b_j - \sqrt{t}z_i - \sqrt{I_j}(\theta_j - \theta_i\sqrt{t})}{\sqrt{1 - t}}\right)}
+where \eqn{b_j} is the efficacy bound at analysis \eqn{j}, \eqn{z_i} is the observed Z-value at analysis \eqn{i},
+and \eqn{\theta_i} natural parameter for treatment effect at analysis \eqn{i}.
+}
+\examples{
+library(gsDesign2)
+library(gsDesign)
+# Calculate conditional power with optional user-input theta (if NULL, will come from the input design)
+alpha <- 0.025
+beta <- 0.1
+ratio <- 1
+enroll_rate <- define_enroll_rate(
+ duration = c(2, 2, 10),
+ rate = (1:3) / 3)
+fail_rate <- define_fail_rate(
+ duration = Inf, fail_rate = log(2) / 9,
+ hr = 0.6, dropout_rate = .0001)
+analysis_time <- c(12, 24, 36)
+ratio <- 1
+upper <- gs_spending_bound
+lower <- gs_b
+upar <- list(sf = sfLDOF, total_spend = alpha)
+lpar <- rep(-Inf, 3)
+x <- gs_design_ahr(
+ enroll_rate = enroll_rate, fail_rate = fail_rate,
+ alpha = alpha, beta = beta, ratio = ratio,
+ info_scale = "h0_h1_info",
+ info_frac = NULL,
+ analysis_time = c(12, 24, 36),
+ upper = upper, upar = upar, test_upper = TRUE,
+ lower = lower, lpar = lpar, test_lower = FALSE,
+) |>
+ to_integer()
+
+# theta is NULL case
+gs_cp_simple(x = x, theta = NULL, i = 1, zi = 1.5)
+
+# User-input theta case
+gs_cp_simple(x = x, theta = c(0.1, 0.2, 0.3), i = 1, zi = 1.5)
+
+}
diff --git a/tests/testthat/_snaps/independent_as_gt.md b/tests/testthat/_snaps/independent_as_gt.md
index 5980018c..1e83b9d5 100644
--- a/tests/testthat/_snaps/independent_as_gt.md
+++ b/tests/testthat/_snaps/independent_as_gt.md
@@ -2,9 +2,9 @@
\begin{table}[t]
\caption*{
- {\large Fixed Design under AHR Method\textsuperscript{\textit{1}}}
+ {\fontsize{20}{25}\selectfont Fixed Design under AHR Method\textsuperscript{\textit{1}}\fontsize{12}{15}\selectfont }
}
- \fontsize{12.0pt}{14.4pt}\selectfont
+ \fontsize{12.0pt}{14.0pt}\selectfont
\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}}lrrrrrrr}
\toprule
Design & N & Events & Time & AHR & Bound & alpha & Power \\
@@ -21,9 +21,9 @@
\begin{table}[t]
\caption*{
- {\large Custom Title\textsuperscript{\textit{1}}}
+ {\fontsize{20}{25}\selectfont Custom Title\textsuperscript{\textit{1}}\fontsize{12}{15}\selectfont }
}
- \fontsize{12.0pt}{14.4pt}\selectfont
+ \fontsize{12.0pt}{14.0pt}\selectfont
\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}}lrrrrrrr}
\toprule
Design & N & Events & Time & AHR & Bound & alpha & Power \\
@@ -40,9 +40,9 @@
\begin{table}[t]
\caption*{
- {\large Fixed Design under Fleming-Harrington Method\textsuperscript{\textit{1}}}
+ {\fontsize{20}{25}\selectfont Fixed Design under Fleming-Harrington Method\textsuperscript{\textit{1}}\fontsize{12}{15}\selectfont }
}
- \fontsize{12.0pt}{14.4pt}\selectfont
+ \fontsize{12.0pt}{14.0pt}\selectfont
\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}}lrrrrrrr}
\toprule
Design & N & Events & Time & AHR & Bound & alpha & Power \\
@@ -59,10 +59,10 @@
\begin{table}[t]
\caption*{
- {\large Bound summary for AHR design} \\
- {\small AHR approximations of \textasciitilde{}HR at bound}
+ {\fontsize{20}{25}\selectfont Bound summary for AHR design\fontsize{12}{15}\selectfont } \\
+ {\fontsize{14}{17}\selectfont AHR approximations of \textasciitilde{}HR at bound\fontsize{12}{15}\selectfont }
}
- \fontsize{12.0pt}{14.4pt}\selectfont
+ \fontsize{12.0pt}{14.0pt}\selectfont
\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}}lrrrrr}
\toprule
& & & & \multicolumn{2}{c}{Cumulative boundary crossing probability} \\
@@ -84,10 +84,10 @@
\begin{table}[t]
\caption*{
- {\large Bound summary for AHR design} \\
- {\small AHR approximations of \textasciitilde{}HR at bound}
+ {\fontsize{20}{25}\selectfont Bound summary for AHR design\fontsize{12}{15}\selectfont } \\
+ {\fontsize{14}{17}\selectfont AHR approximations of \textasciitilde{}HR at bound\fontsize{12}{15}\selectfont }
}
- \fontsize{12.0pt}{14.4pt}\selectfont
+ \fontsize{12.0pt}{14.0pt}\selectfont
\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}}lrrrrr}
\toprule
& & & & \multicolumn{2}{c}{Cumulative boundary crossing probability} \\
@@ -120,10 +120,10 @@
\begin{table}[t]
\caption*{
- {\large Bound summary for WLR design} \\
- {\small WLR approximation of \textasciitilde{}wHR at bound}
+ {\fontsize{20}{25}\selectfont Bound summary for WLR design\fontsize{12}{15}\selectfont } \\
+ {\fontsize{14}{17}\selectfont WLR approximation of \textasciitilde{}wHR at bound\fontsize{12}{15}\selectfont }
}
- \fontsize{12.0pt}{14.4pt}\selectfont
+ \fontsize{12.0pt}{14.0pt}\selectfont
\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}}lrrrrr}
\toprule
& & & & \multicolumn{2}{c}{Cumulative boundary crossing probability} \\
@@ -146,10 +146,10 @@
\begin{table}[t]
\caption*{
- {\large Bound summary for WLR design\textsuperscript{\textit{1}}} \\
- {\small WLR approximation of \textasciitilde{}wHR at bound}
+ {\fontsize{20}{25}\selectfont Bound summary for WLR design\textsuperscript{\textit{1}}\fontsize{12}{15}\selectfont } \\
+ {\fontsize{14}{17}\selectfont WLR approximation of \textasciitilde{}wHR at bound\fontsize{12}{15}\selectfont }
}
- \fontsize{12.0pt}{14.4pt}\selectfont
+ \fontsize{12.0pt}{14.0pt}\selectfont
\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}}lrrrrr}
\toprule
& & & & \multicolumn{2}{c}{Cumulative boundary crossing probability\textsuperscript{\textit{2}}} \\
@@ -184,10 +184,10 @@
\begin{table}[t]
\caption*{
- {\large Bound summary for MaxCombo design} \\
- {\small MaxCombo approximation}
+ {\fontsize{20}{25}\selectfont Bound summary for MaxCombo design\fontsize{12}{15}\selectfont } \\
+ {\fontsize{14}{17}\selectfont MaxCombo approximation\fontsize{12}{15}\selectfont }
}
- \fontsize{12.0pt}{14.4pt}\selectfont
+ \fontsize{12.0pt}{14.0pt}\selectfont
\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}}lrrrr}
\toprule
& & & \multicolumn{2}{c}{Cumulative boundary crossing probability} \\
@@ -220,10 +220,10 @@
\begin{table}[t]
\caption*{
- {\large Bound summary of Binary Endpoint} \\
- {\small measured by risk difference}
+ {\fontsize{20}{25}\selectfont Bound summary of Binary Endpoint\fontsize{12}{15}\selectfont } \\
+ {\fontsize{14}{17}\selectfont measured by risk difference\fontsize{12}{15}\selectfont }
}
- \fontsize{12.0pt}{14.4pt}\selectfont
+ \fontsize{12.0pt}{14.0pt}\selectfont
\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}}lrrrrr}
\toprule
& & & & \multicolumn{2}{c}{Cumulative boundary crossing probability} \\
@@ -244,10 +244,10 @@
\begin{table}[t]
\caption*{
- {\large Bound summary of Binary Endpoint} \\
- {\small measured by risk difference}
+ {\fontsize{20}{25}\selectfont Bound summary of Binary Endpoint\fontsize{12}{15}\selectfont } \\
+ {\fontsize{14}{17}\selectfont measured by risk difference\fontsize{12}{15}\selectfont }
}
- \fontsize{12.0pt}{14.4pt}\selectfont
+ \fontsize{12.0pt}{14.0pt}\selectfont
\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}}lrrrrr}
\toprule
& & & & \multicolumn{2}{c}{Cumulative boundary crossing probability} \\
@@ -278,10 +278,10 @@
\begin{table}[t]
\caption*{
- {\large Bound Summary} \\
- {\small from gs\_power\_wlr}
+ {\fontsize{20}{25}\selectfont Bound Summary\fontsize{12}{15}\selectfont } \\
+ {\fontsize{14}{17}\selectfont from gs\_power\_wlr\fontsize{12}{15}\selectfont }
}
- \fontsize{12.0pt}{14.4pt}\selectfont
+ \fontsize{12.0pt}{14.0pt}\selectfont
\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}}lrrrrr}
\toprule
& & & & \multicolumn{2}{c}{Cumulative boundary crossing probability} \\
@@ -315,10 +315,10 @@
\begin{table}[t]
\caption*{
- {\large Bound summary for WLR design} \\
- {\small WLR approximation of \textasciitilde{}wHR at bound}
+ {\fontsize{20}{25}\selectfont Bound summary for WLR design\fontsize{12}{15}\selectfont } \\
+ {\fontsize{14}{17}\selectfont WLR approximation of \textasciitilde{}wHR at bound\fontsize{12}{15}\selectfont }
}
- \fontsize{12.0pt}{14.4pt}\selectfont
+ \fontsize{12.0pt}{14.0pt}\selectfont
\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}}lrrrrr}
\toprule
& & & & \multicolumn{2}{c}{Cumulative probability to cross boundaries} \\
@@ -352,10 +352,10 @@
\begin{table}[t]
\caption*{
- {\large Bound summary for WLR design\textsuperscript{\textit{1}}} \\
- {\small WLR approximation of \textasciitilde{}wHR at bound}
+ {\fontsize{20}{25}\selectfont Bound summary for WLR design\textsuperscript{\textit{1}}\fontsize{12}{15}\selectfont } \\
+ {\fontsize{14}{17}\selectfont WLR approximation of \textasciitilde{}wHR at bound\fontsize{12}{15}\selectfont }
}
- \fontsize{12.0pt}{14.4pt}\selectfont
+ \fontsize{12.0pt}{14.0pt}\selectfont
\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}}lrrrrr}
\toprule
& & & & \multicolumn{2}{c}{Cumulative boundary crossing probability\textsuperscript{\textit{2}}} \\
@@ -390,10 +390,10 @@
\begin{table}[t]
\caption*{
- {\large Bound summary for WLR design} \\
- {\small WLR approximation of \textasciitilde{}wHR at bound}
+ {\fontsize{20}{25}\selectfont Bound summary for WLR design\fontsize{12}{15}\selectfont } \\
+ {\fontsize{14}{17}\selectfont WLR approximation of \textasciitilde{}wHR at bound\fontsize{12}{15}\selectfont }
}
- \fontsize{12.0pt}{14.4pt}\selectfont
+ \fontsize{12.0pt}{14.0pt}\selectfont
\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}}lrrrrr}
\toprule
& & & & \multicolumn{2}{c}{Cumulative boundary crossing probability} \\
@@ -424,10 +424,10 @@
\begin{table}[t]
\caption*{
- {\large Bound summary for WLR design} \\
- {\small WLR approximation of \textasciitilde{}wHR at bound}
+ {\fontsize{20}{25}\selectfont Bound summary for WLR design\fontsize{12}{15}\selectfont } \\
+ {\fontsize{14}{17}\selectfont WLR approximation of \textasciitilde{}wHR at bound\fontsize{12}{15}\selectfont }
}
- \fontsize{12.0pt}{14.4pt}\selectfont
+ \fontsize{12.0pt}{14.0pt}\selectfont
\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}}lrrrr}
\toprule
& & & \multicolumn{2}{c}{Cumulative boundary crossing probability} \\
diff --git a/tests/testthat/test-developer-gs_cp.R b/tests/testthat/test-developer-gs_cp.R
new file mode 100644
index 00000000..e6b76c38
--- /dev/null
+++ b/tests/testthat/test-developer-gs_cp.R
@@ -0,0 +1,132 @@
+library(gsDesign)
+
+test_that("Compare the gs_cp with gsDesign::gsCP", {
+ # ------------------------------ #
+ # parameters #
+ # ------------------------------ #
+ alpha <- 0.025
+ beta <- 0.1
+ ratio <- 1
+
+ # Enrollment
+ enroll_rate <- define_enroll_rate(
+ duration = c(2, 2, 10),
+ rate = (1:3) / 3)
+
+ # Failure and dropout
+ fail_rate <- define_fail_rate(
+ duration = Inf, fail_rate = log(2) / 9,
+ hr = 0.6, dropout_rate = .0001)
+
+ # IA and FA analysis time
+ analysis_time <- c(12, 24, 36)
+
+ # Randomization ratio
+ ratio <- 1
+
+ # Spending
+ upper <- gs_spending_bound
+ lower <- gs_b
+ upar <- list(sf = sfLDOF, total_spend = alpha)
+ lpar <- rep(-Inf, 3)
+ #lpar <- list(sf = sfLDOF, total_spend = beta)
+ # ------------------------------ #
+ # original design of gsDesign2 #
+ # ------------------------------ #
+ x <- gs_design_ahr(
+ enroll_rate = enroll_rate, fail_rate = fail_rate,
+ alpha = alpha, beta = beta, ratio = ratio,
+ info_scale = "h0_h1_info",
+ info_frac = 1:3/3,
+ analysis_time = 36,
+ upper = upper, upar = upar, test_upper = TRUE,
+ lower = lower, lpar = lpar, test_lower = FALSE,
+ )
+ # ------------------------------ #
+ # original design of gsDesign #
+ # ------------------------------ #
+ x_gsd <- gsSurv(k = 3, test.type = 1, alpha = alpha, beta = beta,
+ astar = 0, timing = 1:3/3,
+ sfu = sfLDOF, sfupar = 0,
+ sfl = sfLDOF, sflpar = 0,
+ lambdaC = log(2) / 9, hr = 0.6, hr0 = 1,
+ eta = fail_rate$dropout_rate |> unique(),
+ gamma = enroll_rate$rate,
+ R = enroll_rate$duration,
+ S = NULL, T = analysis_time[3],
+ minfup = analysis_time[3] - sum(enroll_rate$duration),
+ ratio = ratio)
+
+ # ---------------------------------------------------------- #
+ # case 1 #
+ # currently at IA1, compute conditional power at IA2 and FA #
+ # ---------------------------------------------------------- #
+ gsDesign_cp <- gsCP(x_gsd, i = 1, zi = -qnorm(0.04), theta = -log(0.8)/2)
+
+ gsDesign2_cp <- gs_cp(x = x, i = 1, zi = -qnorm(0.04), theta = c(-log(0.8), -log(0.8), -log(0.8)))
+
+ gsDesign2_npe2 <- gs_cp_npe2(# IA1's Z-score
+ zi = -qnorm(0.04),
+ # IA1, IA2 and FA's theta
+ theta = c(-log(0.8), -log(0.8), -log(0.8)),
+ # IA1, IA2 and FA's information fraction
+ t = x_gsd$timing[1:3],
+ # IA1, IA2 and FA's statistical information
+ info = x_gsd$n.I[1:3] / 4,
+ # IA2 and FA's futility bound
+ a = c(-Inf, -Inf),
+ # IA2 and FA's efficacy bound
+ b = x_gsd$upper$bound[2:3]
+ )
+
+ # IA2's CP given IA1
+ expect_equal(gsDesign_cp$upper$prob[1],
+ gsDesign2_cp$prob_alpha[1],
+ tolerance = 1e-2)
+
+ expect_equal(gsDesign2_cp$prob_alpha[1],
+ gsDesign2_npe2$prob_alpha[1],
+ tolerance = 1e-2)
+
+ # FA's CP given IA1
+ expect_equal(gsDesign_cp$upper$prob[2],
+ gsDesign2_cp$prob_alpha[2],
+ tolerance = 1e-2)
+
+ expect_equal(gsDesign2_cp$prob_alpha[2],
+ gsDesign2_npe2$prob_alpha[2],
+ tolerance = 1e-2)
+
+ # --------------------------------------------------- #
+ # case 2 #
+ # currently at IA2, compute conditional power at FA #
+ # --------------------------------------------------- #
+ gsDesign_cp <- gsCP(x_gsd, i = 2, zi = -qnorm(0.04), theta = -log(0.8)/2)
+
+ gsDesign2_cp <- gs_cp(x = x, i = 2, zi = -qnorm(0.04), theta = c(-log(0.8), -log(0.8)))
+
+ gsDesign2_npe2 <- gs_cp_npe2(# IA2's Z-score
+ zi = -qnorm(0.04),
+ # IA2 and FA's theta
+ theta = c(-log(0.8), -log(0.8)),
+ # IA2 and FA's information fraction
+ t = x_gsd$timing[2:3],
+ # IA2 and FA's statistical information
+ info = x_gsd$n.I[2:3] / 4,
+ # FA's futility bound
+ a = -Inf,
+ # FA's efficacy bound
+ b = x_gsd$upper$bound[3]
+ )
+
+ # FA's CP given IA2
+ expect_equal(gsDesign_cp$upper$prob[1],
+ gsDesign2_cp$prob_alpha[1],
+ tolerance = 1e-2)
+
+ expect_equal(gsDesign2_cp$prob_alpha[1],
+ gsDesign2_npe2$prob_alpha[1],
+ tolerance = 1e-2)
+
+})
+
diff --git a/tests/testthat/test-developer-gs_cp_npe.R b/tests/testthat/test-developer-gs_cp_npe1.R
similarity index 73%
rename from tests/testthat/test-developer-gs_cp_npe.R
rename to tests/testthat/test-developer-gs_cp_npe1.R
index 959370ec..e74674f4 100644
--- a/tests/testthat/test-developer-gs_cp_npe.R
+++ b/tests/testthat/test-developer-gs_cp_npe1.R
@@ -1,6 +1,6 @@
library(gsDesign)
-test_that("Compare the gs_cp_npe with gsDesign::gsCP", {
+test_that("Compare the gs_cp_npe1 with gsDesign::gsCP", {
# ------------------------------ #
# parameters #
# ------------------------------ #
@@ -41,9 +41,7 @@ test_that("Compare the gs_cp_npe with gsDesign::gsCP", {
analysis_time = c(12, 24, 36),
upper = upper, upar = upar, test_upper = TRUE,
lower = lower, lpar = lpar, test_lower = FALSE,
- ) |>
- to_integer()
-
+ )
# ------------------------------ #
# original design of gsDesign #
@@ -60,30 +58,30 @@ test_that("Compare the gs_cp_npe with gsDesign::gsCP", {
minfup = analysis_time[3] - sum(enroll_rate$duration),
ratio = ratio)
- # ----------------------------------------- #
- # conditional power by gs_cp_npe under H0 #
- # ----------------------------------------- #
- cp12_0 <- gs_cp_npe(theta = c(0,0),
+ # ------------------------------------------ #
+ # conditional power by gs_cp_npe1 under H0 #
+ # ------------------------------------------ #
+ cp12_0 <- gs_cp_npe1(theta = c(0,0),
info = x$analysis$info0[c(1,2)],
- a = -qnorm(0.04),
- b = x$bound$z[x$bound$bound == "upper" & x$bound$analysis == 2])
+ zi = -qnorm(0.04),
+ zj = x$bound$z[x$bound$bound == "upper" & x$bound$analysis == 2])
- cp13_0 <- gs_cp_npe(theta = c(0,0),
+ cp13_0 <- gs_cp_npe1(theta = c(0,0),
info = x$analysis$info0[c(1,3)],
- a = -qnorm(0.04),
- b = x$bound$z[x$bound$bound == "upper" & x$bound$analysis == 3])
- # ----------------------------------------- #
- # conditional power by gs_cp_npe under H1 #
- # ----------------------------------------- #
- cp12_1 <- gs_cp_npe(theta = x$analysis$theta[c(1,2)],
+ zi = -qnorm(0.04),
+ zj = x$bound$z[x$bound$bound == "upper" & x$bound$analysis == 3])
+ # ------------------------------------------ #
+ # conditional power by gs_cp_npe1 under H1 #
+ # ------------------------------------------ #
+ cp12_1 <- gs_cp_npe1(theta = x$analysis$theta[c(1,2)],
info = x$analysis$info[c(1,2)],
- a = -qnorm(0.04),
- b = x$bound$z[x$bound$bound == "upper" & x$bound$analysis == 2])
+ zi = -qnorm(0.04),
+ zj = x$bound$z[x$bound$bound == "upper" & x$bound$analysis == 2])
- cp13_1 <- gs_cp_npe(theta = x$analysis$theta[c(1,3)],
+ cp13_1 <- gs_cp_npe1(theta = x$analysis$theta[c(1,3)],
info = x$analysis$info[c(1,3)],
- a = -qnorm(0.04),
- b = x$bound$z[x$bound$bound == "upper" & x$bound$analysis == 3])
+ zi = -qnorm(0.04),
+ zj = x$bound$z[x$bound$bound == "upper" & x$bound$analysis == 3])
# ------------------------------ #
# conditional power by gsDesign #
diff --git a/tests/testthat/test-developer-gs_cp_npe2.R b/tests/testthat/test-developer-gs_cp_npe2.R
new file mode 100644
index 00000000..50025dfa
--- /dev/null
+++ b/tests/testthat/test-developer-gs_cp_npe2.R
@@ -0,0 +1,130 @@
+library(gsDesign)
+
+test_that("Compare the gs_cp_npe2 with gsDesign::gsCP", {
+ # ------------------------------ #
+ # parameters #
+ # ------------------------------ #
+ alpha <- 0.025
+ beta <- 0.1
+ ratio <- 1
+
+ # Enrollment
+ enroll_rate <- define_enroll_rate(
+ duration = c(2, 2, 10),
+ rate = (1:3) / 3)
+
+ # Failure and dropout
+ fail_rate <- define_fail_rate(
+ duration = Inf, fail_rate = log(2) / 9,
+ hr = 0.6, dropout_rate = .0001)
+
+ # IA and FA analysis time
+ analysis_time <- c(12, 24, 36)
+
+ # Randomization ratio
+ ratio <- 1
+
+ # Spending
+ upper <- gs_spending_bound
+ lower <- gs_b
+ upar <- list(sf = sfLDOF, total_spend = alpha)
+ lpar <- rep(-Inf, 3)
+
+ # ------------------------------ #
+ # original design of gsDesign2 #
+ # ------------------------------ #
+ x <- gs_design_ahr(
+ enroll_rate = enroll_rate, fail_rate = fail_rate,
+ alpha = alpha, beta = beta, ratio = ratio,
+ info_scale = "h0_h1_info",
+ info_frac = NULL,
+ analysis_time = c(12, 24, 36),
+ upper = upper, upar = upar, test_upper = TRUE,
+ lower = lower, lpar = lpar, test_lower = FALSE,
+ )
+
+ # ------------------------------ #
+ # original design of gsDesign #
+ # ------------------------------ #
+ x_gsd <- gsSurv(k = 3, test.type = 1, alpha = alpha, beta = beta,
+ astar = 0, timing = 1:3/3,
+ sfu = sfLDOF, sfupar = 0,
+ sfl = sfLDOF, sflpar = 0,
+ lambdaC = log(2) / 9, hr = 0.6, hr0 = 1,
+ eta = fail_rate$dropout_rate |> unique(),
+ gamma = enroll_rate$rate,
+ R = enroll_rate$duration,
+ S = NULL, T = analysis_time[3],
+ minfup = analysis_time[3] - sum(enroll_rate$duration),
+ ratio = ratio)
+
+ # ------------------------------------------------------ #
+ # case 1 #
+ # currently at IA1, compute conditional power at IA2/FA #
+ # ------------------------------------------------------ #
+ gsDesign_cp <- gsCP(x_gsd, i = 1, zi = -qnorm(0.04), theta = -log(0.8)/2)
+
+ gsDesign2_cp <- gs_cp_npe2(# IA1's Z-score
+ zi = -qnorm(0.04),
+ # IA1, IA2 and FA's theta
+ theta = c(-log(0.8), -log(0.8), -log(0.8)),
+ # IA1, IA2 and FA's information fraction
+ t = x_gsd$timing[1:3],
+ # IA1, IA2 and FA's statistical information
+ info = x_gsd$n.I[1:3] / 4,
+ # IA2 and FA's futility bound
+ a = c(-Inf, -Inf),
+ # IA2 and FA's efficacy bound
+ b = x_gsd$upper$bound[2:3]
+ )
+
+ gsDesign2_simple_cp <- gs_cp_npe1(theta = c(-log(0.8), -log(0.8)),
+ info = x_gsd$n.I[1:2] / 4,
+ zi = -qnorm(0.04),
+ zj = x_gsd$upper$bound[2])
+
+ # IA2's CP given IA1
+ expect_equal(gsDesign_cp$upper$prob[1],
+ gsDesign2_cp$prob_alpha[1])
+
+ expect_equal(gsDesign_cp$upper$prob[1],
+ gsDesign2_simple_cp)
+
+ # FA's CP given IA1
+ expect_equal(gsDesign_cp$upper$prob[2],
+ gsDesign2_cp$prob_alpha[2])
+
+ # --------------------------------------------------- #
+ # case 2 #
+ # currently at IA2, compute conditional power at FA #
+ # --------------------------------------------------- #
+ gsDesign_cp <- gsCP(x_gsd, i = 2, zi = -qnorm(0.04), theta = -log(0.8)/2)
+
+ gsDesign2_cp <- gs_cp_npe2(# IA2's Z-score
+ zi = -qnorm(0.04),
+ # IA2 and FA's theta
+ theta = c(-log(0.8), -log(0.8)),
+ # IA2 and FA's information fraction
+ t = x_gsd$timing[2:3],
+ # IA2 and FA's statistical information
+ info = x_gsd$n.I[2:3] / 4,
+ # FA's futility bound
+ a = -Inf,
+ # FA's efficacy bound
+ b = x_gsd$upper$bound[3]
+ )
+
+ gsDesign2_simple_cp <- gs_cp_npe1(theta = c(-log(0.8)),
+ info = x_gsd$n.I[2:3] / 4,
+ zi = -qnorm(0.04),
+ zj = x_gsd$upper$bound[3])
+ # FA's CP given IA2
+ expect_equal(gsDesign_cp$upper$prob[1],
+ gsDesign2_cp$prob_alpha[1])
+
+ expect_equal(gsDesign_cp$upper$prob[1],
+ gsDesign2_simple_cp)
+
+})
+
+
diff --git a/tests/testthat/test-developer-gs_cp_simple.R b/tests/testthat/test-developer-gs_cp_simple.R
new file mode 100644
index 00000000..d9d67bb3
--- /dev/null
+++ b/tests/testthat/test-developer-gs_cp_simple.R
@@ -0,0 +1,146 @@
+library(gsDesign)
+
+test_that("Compare the gs_cp_simple with gsDesign::gsCP", {
+ # ------------------------------ #
+ # parameters #
+ # ------------------------------ #
+ alpha <- 0.025
+ beta <- 0.1
+ ratio <- 1
+
+ # Enrollment
+ enroll_rate <- define_enroll_rate(
+ duration = c(2, 2, 10),
+ rate = (1:3) / 3)
+
+ # Failure and dropout
+ fail_rate <- define_fail_rate(
+ duration = Inf, fail_rate = log(2) / 9,
+ hr = 0.6, dropout_rate = .0001)
+
+ # IA and FA analysis time
+ analysis_time <- c(12, 24, 36)
+
+ # Randomization ratio
+ ratio <- 1
+
+ # Spending
+ upper <- gs_spending_bound
+ lower <- gs_b
+ upar <- list(sf = sfLDOF, total_spend = alpha)
+ lpar <- rep(-Inf, 3)
+
+ # ------------------------------ #
+ # original design of gsDesign2 #
+ # ------------------------------ #
+ x <- gs_design_ahr(
+ enroll_rate = enroll_rate, fail_rate = fail_rate,
+ alpha = alpha, beta = beta, ratio = ratio,
+ info_scale = "h0_h1_info",
+ info_frac = NULL,
+ analysis_time = c(12, 24, 36),
+ upper = upper, upar = upar, test_upper = TRUE,
+ lower = lower, lpar = lpar, test_lower = FALSE,
+ ) |>
+ to_integer()
+
+
+ # ------------------------------ #
+ # original design of gsDesign #
+ # ------------------------------ #
+ x_gsd <- gsSurv(k = 3, test.type = 1, alpha = alpha, beta = beta,
+ astar = 0, timing = x$analysis$info_frac,
+ sfu = sfLDOF, sfupar = 0,
+ sfl = sfLDOF, sflpar = 0,
+ lambdaC = log(2) / 9, hr = 0.6, hr0 = 1,
+ eta = fail_rate$dropout_rate |> unique(),
+ gamma = enroll_rate$rate,
+ R = enroll_rate$duration,
+ S = NULL, T = analysis_time[3],
+ minfup = analysis_time[3] - sum(enroll_rate$duration),
+ ratio = ratio)
+
+ # -------------------------------------------- #
+ # conditional power by gs_cp_simple under H0 #
+ # -------------------------------------------- #
+ cp_0 <- gs_cp_simple(x = x, theta = c(0,0,0), i = 1, zi = -qnorm(0.04))
+
+ # -------------------------------------------- #
+ # conditional power by gs_cp_simple under H1 #
+ # -------------------------------------------- #
+ cp_1 <- gs_cp_simple(x = x, i = 1, zi = -qnorm(0.04))
+
+ # ------------------------------------------ #
+ # conditional power by gs_cp_npe1 under H0 #
+ # ------------------------------------------ #
+ cp12_0 <- gs_cp_npe1(theta = c(0,0),
+ info = x$analysis$info0[c(1,2)],
+ zi = -qnorm(0.04),
+ zj = x$bound$z[x$bound$bound == "upper" & x$bound$analysis == 2])
+
+ cp13_0 <- gs_cp_npe1(theta = c(0,0),
+ info = x$analysis$info0[c(1,3)],
+ zi = -qnorm(0.04),
+ zj = x$bound$z[x$bound$bound == "upper" & x$bound$analysis == 3])
+ # ------------------------------------------ #
+ # conditional power by gs_cp_npe1 under H1 #
+ # ------------------------------------------ #
+ cp12_1 <- gs_cp_npe1(theta = x$analysis$theta[c(1,2)],
+ info = x$analysis$info[c(1,2)],
+ zi = -qnorm(0.04),
+ zj = x$bound$z[x$bound$bound == "upper" & x$bound$analysis == 2])
+
+ cp13_1 <- gs_cp_npe1(theta = x$analysis$theta[c(1,3)],
+ info = x$analysis$info[c(1,3)],
+ zi = -qnorm(0.04),
+ zj = x$bound$z[x$bound$bound == "upper" & x$bound$analysis == 3])
+
+ # ------------------------------ #
+ # conditional power by gsDesign #
+ # ------------------------------ #
+ xcp_gsd <- gsCP(x_gsd, i = 1, zi = -qnorm(0.04))
+
+ # ------------------------------ #
+ # comparison #
+ # ------------------------------ #
+ # under H0
+ # given IA1 assumed blinded data and compute IA2 conditional power
+ expect_equal(xcp_gsd$upper$prob[1, 2],
+ cp_0[1],
+ tolerance = 5e-2)
+
+ expect_equal(cp_0[1],
+ cp12_0,
+ tolerance = 5e-2)
+
+ # given IA1 assumed blinded data and compute FA conditional power
+ expect_equal(sum(xcp_gsd$upper$prob[, 2]),
+ cp_0[2],
+ tolerance = 5e-2)
+
+ expect_equal(cp_0[2],
+ cp13_0,
+ tolerance = 5e-2)
+ # under H1
+ # given IA1 assumed blinded data and compute IA2 conditional power
+ expect_equal(xcp_gsd$upper$prob[1, 3],
+ cp_1[1],
+ tolerance = 5e-2)
+
+ expect_equal(cp_1[1],
+ cp12_1,
+ tolerance = 5e-2)
+ # given IA1 assumed blinded data and compute FA conditional power
+ expect_equal(sum(xcp_gsd$upper$prob[, 3]),
+ cp_1[2],
+ tolerance = 5e-2)
+
+ expect_equal(cp_1[2],
+ cp13_1,
+ tolerance = 5e-2)
+
+})
+
+
+
+