From 1c76ef440cccbefeec6d054aedeb140be481fe9c Mon Sep 17 00:00:00 2001 From: LittleBeannie Date: Fri, 23 May 2025 17:11:28 -0400 Subject: [PATCH 01/29] update version number in DESCRIPTION --- DESCRIPTION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index c0c0294b..bad089e3 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: gsDesign2 Title: Group Sequential Design with Non-Constant Effect -Version: 1.1.3.5 +Version: 1.1.3.6 Authors@R: c( person("Keaven", "Anderson", email = "keaven_anderson@merck.com", role = c("aut")), person("Yilong", "Zhang", email = "elong0527@gmail.com", role = c("aut")), From c047205877837ce5af1b37cdd037d82f9cf67ac1 Mon Sep 17 00:00:00 2001 From: LittleBeannie Date: Fri, 23 May 2025 17:11:42 -0400 Subject: [PATCH 02/29] update pkgdown.yml --- _pkgdown.yml | 1 + 1 file changed, 1 insertion(+) diff --git a/_pkgdown.yml b/_pkgdown.yml index 111ec926..b0aa19e0 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -77,6 +77,7 @@ reference: desc: > Functions for conditional power. contents: + - gs_cp_ahr - gs_cp_npe - title: "Input definition" desc: > From a5a4d1a3f73e09f7ea25c46cb4f7e0a3abb0cd17 Mon Sep 17 00:00:00 2001 From: LittleBeannie Date: Fri, 23 May 2025 17:12:37 -0400 Subject: [PATCH 03/29] add a new function `gs_cp_ahr` --- R/gs_cp_ahr.R | 123 ++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 123 insertions(+) create mode 100644 R/gs_cp_ahr.R diff --git a/R/gs_cp_ahr.R b/R/gs_cp_ahr.R new file mode 100644 index 00000000..1fcb1d76 --- /dev/null +++ b/R/gs_cp_ahr.R @@ -0,0 +1,123 @@ +# 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 for group sequential designs using the average hazard ratio method. +#' +#' @param x A design object of type `gs_design_ahr()`, `gs_power_ahr()`, or `gs_update_ahr()`. +#' @param i Analysis at which interim z-value is given; must be from 1 to `nrow(x$analysis) - 1`. +#' @param z_i Interim z-value at analysis `i` (scalar). +#' @param hr_future A scalar or vector, which specifies the natural parameter for treatment effect at the future analyses. +#' The default gives you a range of possible HRs from 0.6 to 1. +#' +#' @return A data frame showing the conditional power based on the input design, interim test statistic, and possible future HR. +#' @export +#' +#' @examples +#' library(gsDesign) +#' library(gsDesign2) +#' +#' # Original design +#' x <- gs_design_ahr(info_frac = c(0.5, 0.8, 1)) +#' +#' # Updated design +#' xu <- gs_update_ahr( +#' x = x, +#' ustime = c(170, x$analysis$event[2:3] )/ max(x$analysis$event), +#' event_tbl = data.frame(analysis = c(1, 1), event = c(100, 70))) +#' +#' # We assume IA1 p-value of 0.04, one-sided, and future HR is possible from 0.6 to 1. +#' gs_cp_ahr(x = xu, +#' i = 1, +#' z_i = -qnorm(0.04), +#' hr_future = seq(0.6, 1, by = 0.1)) +#' +#' # We assume the test statistics z_i is derived by the function hrn2z(), which translates +#' # a HR and number of events into an approximate corresponding Z-value, +#' # using the Schoenfeld approximation. +#' # The future HR is possible from 0.6 to 1. +#' gs_cp_ahr(x = xu, +#' i = 1, +#' z_i = -hrn2z(hr = 0.7, n = xu$analysis$event[1], ratio = 1), +#' hr_future = seq(0.6, 1, by = 0.1)) +#' +gs_cp_ahr <- function(x = NULL, + i = 1, + z_i = 0, + hr_i = NULL, + hr_future = NULL) { + + # ----------------------------------------- # + # initialization # + # ----------------------------------------- # + n_analysis <- nrow(x$analysis) + ratio <- x$input$x$input$ratio + q_e <- ratio / (1 + ratio) + event_i <- x$analysis$event[i] + + # ----------------------------------------- # + # input checking # + # ----------------------------------------- # + if (is.null(x)) { + stop("Please provide the design generated by gs_design_ahr, gs_power_ahr or gs_update_ahr.") + } + + if (length(z_i) > 1) { + stop("The z_i is required to be a scalers.") + } + + theta_i <- z_i / sqrt(event_i) + + # mutate hr_future if it is null + if (is.null(hr_future)) { + hr_future <- seq(0.5, 1, 0.1) + } + theta_future <- -log(hr_future) * sqrt(ratio / (1 + ratio)^2) + + # -------------------------------------------- # + # calculate conditional power under H0 theta # + # -------------------------------------------- # + ans <- NULL + + for (k in (i+1):n_analysis) { + + ans_k <- NULL + event_k <- x$analysis$event[k] + + for (j in seq_along(theta_future)) { + cp <- gs_cp_npe(theta = c(theta_i, theta_future[j]), + info = c(event_i, event_k) * q_e * (1 - q_e), + a = z_i, + b = x$bound$z[x$bound$analysis == k & x$bound$bound == "upper"]) + + ans_k_new <- data.frame(current_analysis = i, + future_analysis = k, + future_hr = hr_future[j], + conditional_power = cp) + + ans_k <- rbind(ans_k, ans_k_new) + } + + ans <- rbind(ans, ans_k) + } + + + # -------------------------------------------- # + # output # + # -------------------------------------------- # + return(ans) +} From 9f53d7debe1a6899d61acf3bb686e06f3ce5d386 Mon Sep 17 00:00:00 2001 From: LittleBeannie Date: Fri, 23 May 2025 17:12:55 -0400 Subject: [PATCH 04/29] add `input` to the output of `gs_update_ahr` --- R/gs_update_ahr.R | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/R/gs_update_ahr.R b/R/gs_update_ahr.R index 56ac0164..2d45148a 100644 --- a/R/gs_update_ahr.R +++ b/R/gs_update_ahr.R @@ -304,6 +304,13 @@ gs_update_ahr <- function( # ----------------------------------- # ans <- list() + ans$input <- list() + ans$input$x <- x + ans$input$alpha <- alpha + ans$input$ustime <- ustime + ans$input$lstime <- lstime + ans$input$event_tbl <- event_tbl + ans$enroll_rate <- x$enroll_rate ans$fail_rate <- x$fail_rate From 922d96e3a30643142c12fa6263ddeec40310c947 Mon Sep 17 00:00:00 2001 From: LittleBeannie Date: Fri, 23 May 2025 17:22:12 -0400 Subject: [PATCH 05/29] add Rd files --- man/gs_cp_ahr.Rd | 53 ++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 53 insertions(+) create mode 100644 man/gs_cp_ahr.Rd diff --git a/man/gs_cp_ahr.Rd b/man/gs_cp_ahr.Rd new file mode 100644 index 00000000..5a21cba5 --- /dev/null +++ b/man/gs_cp_ahr.Rd @@ -0,0 +1,53 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/gs_cp_ahr.R +\name{gs_cp_ahr} +\alias{gs_cp_ahr} +\title{Conditional power computation for group sequential designs using the average hazard ratio method.} +\usage{ +gs_cp_ahr(x = NULL, i = 1, z_i = 0, hr_i = NULL, hr_future = NULL) +} +\arguments{ +\item{x}{A design object of type \code{gs_design_ahr()}, \code{gs_power_ahr()}, or \code{gs_update_ahr()}.} + +\item{i}{Analysis at which interim z-value is given; must be from 1 to \code{nrow(x$analysis) - 1}.} + +\item{z_i}{Interim z-value at analysis \code{i} (scalar).} + +\item{hr_future}{A scalar or vector, which specifies the natural parameter for treatment effect at the future analyses. +The default gives you a range of possible HRs from 0.6 to 1.} +} +\value{ +A data frame showing the conditional power based on the input design, interim test statistic, and possible future HR. +} +\description{ +Conditional power computation for group sequential designs using the average hazard ratio method. +} +\examples{ +library(gsDesign) +library(gsDesign2) + +# Original design +x <- gs_design_ahr(info_frac = c(0.5, 0.8, 1)) + +# Updated design +xu <- gs_update_ahr( + x = x, + ustime = c(170, x$analysis$event[2:3] )/ max(x$analysis$event), + event_tbl = data.frame(analysis = c(1, 1), event = c(100, 70))) + +# We assume IA1 p-value of 0.04, one-sided, and future HR is possible from 0.6 to 1. +gs_cp_ahr(x = xu, + i = 1, + z_i = -qnorm(0.04), + hr_future = seq(0.6, 1, by = 0.1)) + +# We assume the test statistics z_i is derived by the function hrn2z(), which translates +# a HR and number of events into an approximate corresponding Z-value, +# using the Schoenfeld approximation. +# The future HR is possible from 0.6 to 1. +gs_cp_ahr(x = xu, + i = 1, + z_i = -hrn2z(hr = 0.7, n = xu$analysis$event[1], ratio = 1), + hr_future = seq(0.6, 1, by = 0.1)) + +} From e79482a332a20f744faf88f31e72d6c845406200 Mon Sep 17 00:00:00 2001 From: LittleBeannie Date: Fri, 23 May 2025 17:26:36 -0400 Subject: [PATCH 06/29] fix ci/cd check --- NAMESPACE | 1 + R/gs_cp_ahr.R | 1 - man/gs_cp_ahr.Rd | 2 +- 3 files changed, 2 insertions(+), 2 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index b144d025..85858776 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -28,6 +28,7 @@ export(fixed_design_rd) export(fixed_design_rmst) export(gs_b) export(gs_bound_summary) +export(gs_cp_ahr) export(gs_cp_npe) export(gs_create_arm) export(gs_design_ahr) diff --git a/R/gs_cp_ahr.R b/R/gs_cp_ahr.R index 1fcb1d76..bf318652 100644 --- a/R/gs_cp_ahr.R +++ b/R/gs_cp_ahr.R @@ -58,7 +58,6 @@ gs_cp_ahr <- function(x = NULL, i = 1, z_i = 0, - hr_i = NULL, hr_future = NULL) { # ----------------------------------------- # diff --git a/man/gs_cp_ahr.Rd b/man/gs_cp_ahr.Rd index 5a21cba5..ccbf69ff 100644 --- a/man/gs_cp_ahr.Rd +++ b/man/gs_cp_ahr.Rd @@ -4,7 +4,7 @@ \alias{gs_cp_ahr} \title{Conditional power computation for group sequential designs using the average hazard ratio method.} \usage{ -gs_cp_ahr(x = NULL, i = 1, z_i = 0, hr_i = NULL, hr_future = NULL) +gs_cp_ahr(x = NULL, i = 1, z_i = 0, hr_future = NULL) } \arguments{ \item{x}{A design object of type \code{gs_design_ahr()}, \code{gs_power_ahr()}, or \code{gs_update_ahr()}.} From b1961e30d5a11d5ff49b6941e8355d7cdd9ef442 Mon Sep 17 00:00:00 2001 From: Shiyu Zhang Date: Mon, 4 Aug 2025 17:23:59 -0400 Subject: [PATCH 07/29] remove the gs_cp_ahr function as this is not correct --- R/gs_cp_ahr.R | 122 -------------------------------------------------- 1 file changed, 122 deletions(-) delete mode 100644 R/gs_cp_ahr.R diff --git a/R/gs_cp_ahr.R b/R/gs_cp_ahr.R deleted file mode 100644 index bf318652..00000000 --- a/R/gs_cp_ahr.R +++ /dev/null @@ -1,122 +0,0 @@ -# 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 for group sequential designs using the average hazard ratio method. -#' -#' @param x A design object of type `gs_design_ahr()`, `gs_power_ahr()`, or `gs_update_ahr()`. -#' @param i Analysis at which interim z-value is given; must be from 1 to `nrow(x$analysis) - 1`. -#' @param z_i Interim z-value at analysis `i` (scalar). -#' @param hr_future A scalar or vector, which specifies the natural parameter for treatment effect at the future analyses. -#' The default gives you a range of possible HRs from 0.6 to 1. -#' -#' @return A data frame showing the conditional power based on the input design, interim test statistic, and possible future HR. -#' @export -#' -#' @examples -#' library(gsDesign) -#' library(gsDesign2) -#' -#' # Original design -#' x <- gs_design_ahr(info_frac = c(0.5, 0.8, 1)) -#' -#' # Updated design -#' xu <- gs_update_ahr( -#' x = x, -#' ustime = c(170, x$analysis$event[2:3] )/ max(x$analysis$event), -#' event_tbl = data.frame(analysis = c(1, 1), event = c(100, 70))) -#' -#' # We assume IA1 p-value of 0.04, one-sided, and future HR is possible from 0.6 to 1. -#' gs_cp_ahr(x = xu, -#' i = 1, -#' z_i = -qnorm(0.04), -#' hr_future = seq(0.6, 1, by = 0.1)) -#' -#' # We assume the test statistics z_i is derived by the function hrn2z(), which translates -#' # a HR and number of events into an approximate corresponding Z-value, -#' # using the Schoenfeld approximation. -#' # The future HR is possible from 0.6 to 1. -#' gs_cp_ahr(x = xu, -#' i = 1, -#' z_i = -hrn2z(hr = 0.7, n = xu$analysis$event[1], ratio = 1), -#' hr_future = seq(0.6, 1, by = 0.1)) -#' -gs_cp_ahr <- function(x = NULL, - i = 1, - z_i = 0, - hr_future = NULL) { - - # ----------------------------------------- # - # initialization # - # ----------------------------------------- # - n_analysis <- nrow(x$analysis) - ratio <- x$input$x$input$ratio - q_e <- ratio / (1 + ratio) - event_i <- x$analysis$event[i] - - # ----------------------------------------- # - # input checking # - # ----------------------------------------- # - if (is.null(x)) { - stop("Please provide the design generated by gs_design_ahr, gs_power_ahr or gs_update_ahr.") - } - - if (length(z_i) > 1) { - stop("The z_i is required to be a scalers.") - } - - theta_i <- z_i / sqrt(event_i) - - # mutate hr_future if it is null - if (is.null(hr_future)) { - hr_future <- seq(0.5, 1, 0.1) - } - theta_future <- -log(hr_future) * sqrt(ratio / (1 + ratio)^2) - - # -------------------------------------------- # - # calculate conditional power under H0 theta # - # -------------------------------------------- # - ans <- NULL - - for (k in (i+1):n_analysis) { - - ans_k <- NULL - event_k <- x$analysis$event[k] - - for (j in seq_along(theta_future)) { - cp <- gs_cp_npe(theta = c(theta_i, theta_future[j]), - info = c(event_i, event_k) * q_e * (1 - q_e), - a = z_i, - b = x$bound$z[x$bound$analysis == k & x$bound$bound == "upper"]) - - ans_k_new <- data.frame(current_analysis = i, - future_analysis = k, - future_hr = hr_future[j], - conditional_power = cp) - - ans_k <- rbind(ans_k, ans_k_new) - } - - ans <- rbind(ans, ans_k) - } - - - # -------------------------------------------- # - # output # - # -------------------------------------------- # - return(ans) -} From 9c3ad1298604e42cb6e63695d6cb8be4ad230aa7 Mon Sep 17 00:00:00 2001 From: Shiyu Zhang Date: Mon, 4 Aug 2025 17:24:59 -0400 Subject: [PATCH 08/29] add gs_cp_npe2 function --- R/gs_cp_npe2.R | 191 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 191 insertions(+) create mode 100644 R/gs_cp_npe2.R diff --git a/R/gs_cp_npe2.R b/R/gs_cp_npe2.R new file mode 100644 index 00000000..1c55bb75 --- /dev/null +++ b/R/gs_cp_npe2.R @@ -0,0 +1,191 @@ +# 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-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{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 = 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 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-1, which specifies the futility bounds from analysis i+1 to analysis j-1. +#' @param b A vector of length j-i, which specifies the efficacy bounds from analysis i+1 to analysis j. +#' @param c Interim z-value at analysis i (scalar). +#' @return A list of conditional powers: prob_alpha = \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 = \eqn{P(\{Z_j \geq b_j\} \& \{\cap_{m=i+1}^{j-1} Z_m < b_m\} \mid Z_i = c)}. +#' prob_beta = eqn{P(\{Z_j \leq b_j\} \& \{\cap_{m=i+1}^{j-1} a_m \leq Z_m < b_m\} \mid Z_i = c)}. +#' @export +#' +#' @examples +#' library(gsDesign2) +#' library(dplyr) +#' library(mvtnorm) +#' # 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), +#' b = c(1.8, 2.1), +#' c = 1.5) +gs_cp_npe2 <- function(theta = NULL, + t = NULL, + info = NULL, + a = NULL, + b = NULL, + c = NULL){ + # Input check + 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(a) != (length(b) - 1)) + stop("The vector of lower bound should have the length of vector of lower bound - 1.") + + # Build mean vector of length M = j-i + i <- 1 + j <- length(b) + dim <- j - i # number of increments + + # Build mean vector + mu <- sapply(seq_len(dim), function(k){ + idx <- i + k + theta[idx] * sqrt(t[idx] * info[idx]) - theta[i] * sqrt(t[i] * info[i]) + }) + + # Build covariance matrix Sigma (dim x dim) + Sigma <- matrix(0, nrow = dim, ncol = dim) + + for(k in seq_len(dim)) { + for(l in seq_len(dim)) { + Sigma[k, l] <- t[i + min(k,l)] - t[i] + } + } + + # --------------------------------------------------------------------- + # Calculate Lower/upper bounds for alpha + # first cross efficacy bound at analysis j given no efficacy/futility + # bound crossing at analysis i+1 to j-1 + # --------------------------------------------------------------------- + # 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, j - i) + for(m in 1:(j-i-1)){ + lower_alpha[m] <- a[m] * sqrt(t[m]) - c * sqrt(t[i]) + } + lower_alpha[j-i] <- b[j-i] * sqrt(t[j]) - c * sqrt(t[i]) + + # upper bound + upper_alpha <- rep(0, j - i) + for(m in 1:(j-i-1)){ + upper_alpha[m] <- b[m] * sqrt(t[m]) - c * sqrt(t[i]) + } + upper_alpha[j-i] <- Inf + + # --------------------------------------------------------------------- + # Calculate Lower/upper bounds for alpha_plus + # first cross efficacy bound at analysis j given no efficacy bound + # crossing at analysis i+1 to j-1 + # --------------------------------------------------------------------- + # Integration limits: D_m = B_m - B_i + # for D_{i+1},...,D_{j-1} use (-Inf, b); for D_j use [b_j, +Inf) + # lower bound + lower_alpha_plus <- rep(-Inf, j - i) + lower_alpha_plus[j-i] <- b[j-i] * sqrt(t[j]) - c * sqrt(t[i]) + + # upper bound + upper_alpha_plus <- rep(0, j - i) + for(m in 1:(j-i-1)){ + upper_alpha_plus[m] <- b[m] * sqrt(t[m]) - c * sqrt(t[i]) + } + upper_alpha_plus[j-i] <- Inf + + + + # --------------------------------------------------------------------- + # Calculate Lower/upper bounds for beta + # Not cross efficacy bound at analysis j given no efficacy/futility + # bound crossing at analysis i+1 to j-1 + # --------------------------------------------------------------------- + # Integration limits: D_m = B_m - B_i + # for D_{i+1},...,D_{j-1} use [a, b); for D_j use [-Inf, b_j) + # lower bound + lower_beta <- rep(0, j - i) + for(m in 1:(j-i-1)){ + lower_beta[m] <- a[m] * sqrt(t[m]) - c * sqrt(t[i]) + } + lower_beta[j-i] <- -Inf + + # upper bound + upper_beta <- rep(0, j - i) + for(m in 1:(j-i)){ + upper_beta[m] <- b[m] * sqrt(t[m]) - c * sqrt(t[i]) + } + + + # Compute multivariate normal probability + prob_alpha <- mvtnorm::pmvnorm( + lower = lower_alpha, + upper = upper_alpha, + mean = mu, + sigma = Sigma)[1] + + prob_alpha_plus <- mvtnorm::pmvnorm( + lower=lower_alpha_plus, + upper=upper_alpha_plus, + mean=mu, + sigma=Sigma)[1] + + prob_beta <- mvtnorm::pmvnorm( + lower=lower_beta, + upper=upper_beta, + mean=mu, + sigma=Sigma)[1] + + + return(prob = list(prob_alpha = prob_alpha, + prob_alpha_plus = prob_alpha_plus, + prob_beta = prob_beta)) + + +} + + + + + + + + From 3f0518de1e668682be0f636df1f470dbcc73178f Mon Sep 17 00:00:00 2001 From: LittleBeannie Date: Wed, 10 Sep 2025 17:14:29 -0400 Subject: [PATCH 09/29] review code --- NAMESPACE | 2 +- R/gs_cp_npe2.R | 103 +++++++++++++++++++++++++++++++++------------- _pkgdown.yml | 2 +- man/gs_cp_ahr.Rd | 53 ------------------------ man/gs_cp_npe2.Rd | 79 +++++++++++++++++++++++++++++++++++ 5 files changed, 156 insertions(+), 83 deletions(-) delete mode 100644 man/gs_cp_ahr.Rd create mode 100644 man/gs_cp_npe2.Rd diff --git a/NAMESPACE b/NAMESPACE index 85858776..7f5edfdc 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -28,8 +28,8 @@ export(fixed_design_rd) export(fixed_design_rmst) export(gs_b) export(gs_bound_summary) -export(gs_cp_ahr) export(gs_cp_npe) +export(gs_cp_npe2) export(gs_create_arm) export(gs_design_ahr) export(gs_design_combo) diff --git a/R/gs_cp_npe2.R b/R/gs_cp_npe2.R index 1c55bb75..b88751d9 100644 --- a/R/gs_cp_npe2.R +++ b/R/gs_cp_npe2.R @@ -48,6 +48,7 @@ #' library(gsDesign2) #' 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. @@ -57,59 +58,105 @@ #' a = c(-0.5), #' b = c(1.8, 2.1), #' c = 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.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() +#' # update design based on IA1 observed data ---- +#' # assume at IA1, we observed 150 events (planned of 140) during the non-effect period, and 180 events (planned of 170) afterwards. +#' event_tbl <- data.frame(analysis = c(1, 1), event = c(150, 180)) +#' ustime <- c(150+180, x$analysis$event[2:3]) / max(x$analysis$event) +#' xu <- gs_update_ahr(x, alpha = NULL, ustime = ustime, lstime = NULL, event_tbl = event_tbl) +#' # calculate conditional power at IA1 based on non-constant treatment effect ---- +#' gs_cp_npe2(theta = x$analysis$theta, +#' t = x$analysis$info_frac, +#' info = x$analysis$info, +#' a = rep(-Inf, 3), +#' b = x$bound$z[x$bound$bound == "upper"], +#' c = -gsDesign::hrn2z(hr = 0.8, n = 150+180, ratio = 1)) gs_cp_npe2 <- function(theta = NULL, t = NULL, info = NULL, a = NULL, b = NULL, c = NULL){ - # Input check + # ------------------------------ # + # 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(a) != (length(b) - 1)) - stop("The vector of lower bound should have the length of vector of lower bound - 1.") - - # Build mean vector of length M = j-i + # ------------------------------ # + # Initialization + # ------------------------------ # i <- 1 + # the conditional power is calculated from analysis i to analysis j + # the analysis j is decided by the length of b (efficacy bound) j <- length(b) - dim <- j - i # number of increments + # let D_m = B_m - B_i, where m = i+1, i+2, ..., j + dim <- j - i + + # ------------------------------ # + # Build the asymptotic + # mean of B_j - B_i + # vector of length dim + # ------------------------------ # - # Build mean vector mu <- sapply(seq_len(dim), function(k){ idx <- i + k theta[idx] * sqrt(t[idx] * info[idx]) - theta[i] * sqrt(t[i] * info[i]) }) - # Build covariance matrix Sigma (dim x dim) + # ------------------------------ # + # Build the asymptotic + # covariance of B_j - B_i + # matrix of (dim x dim) + # ------------------------------ # Sigma <- matrix(0, nrow = dim, ncol = dim) for(k in seq_len(dim)) { for(l in seq_len(dim)) { - Sigma[k, l] <- t[i + min(k,l)] - t[i] + Sigma[k, l] <- t[i + min(k, l)] - t[i] } } - # --------------------------------------------------------------------- - # Calculate Lower/upper bounds for alpha - # first cross efficacy bound at analysis j given no efficacy/futility - # bound crossing at analysis i+1 to j-1 - # --------------------------------------------------------------------- - # 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 + # ------------------------------ # + # Calculate integration limit + # for alpha(i,j) + # ------------------------------ # + # alpha(i,j) is the conditional probability of first cross efficacy bound + # at analysis j given no efficacy/futility bound crossing between analysis i amd j + # D_m = B_m - B_i, where m = i+1, i+2, ..., j + # for D_{i+1},...,D_{j-1} use [a, b); + # for D_j use [b_j, +Inf) + # integration lower bound lower_alpha <- rep(0, j - i) - for(m in 1:(j-i-1)){ + for (m in 1:(j - i - 1)) { + ## ?? here looks wrong: when m = 1, it is a[1]*sqrt(t[1]) - c*sqrt(t[i]) where i is always 1 -- which is wrong + ## when m = 1, it is B2 - B1, which should be a[2]*sqrt(t[2]) - c*sqrt(t[i]) lower_alpha[m] <- a[m] * sqrt(t[m]) - c * sqrt(t[i]) } lower_alpha[j-i] <- b[j-i] * sqrt(t[j]) - c * sqrt(t[i]) - # upper bound + # integration upper bound upper_alpha <- rep(0, j - i) - for(m in 1:(j-i-1)){ + for(m in 1:(j - i - 1)){ + ## ?? same as above upper_alpha[m] <- b[m] * sqrt(t[m]) - c * sqrt(t[i]) } upper_alpha[j-i] <- Inf @@ -163,16 +210,16 @@ gs_cp_npe2 <- function(theta = NULL, sigma = Sigma)[1] prob_alpha_plus <- mvtnorm::pmvnorm( - lower=lower_alpha_plus, - upper=upper_alpha_plus, - mean=mu, - sigma=Sigma)[1] + lower = lower_alpha_plus, + upper = upper_alpha_plus, + mean = mu, + sigma = Sigma)[1] prob_beta <- mvtnorm::pmvnorm( - lower=lower_beta, - upper=upper_beta, - mean=mu, - sigma=Sigma)[1] + lower = lower_beta, + upper = upper_beta, + mean = mu, + sigma = Sigma)[1] return(prob = list(prob_alpha = prob_alpha, diff --git a/_pkgdown.yml b/_pkgdown.yml index b0aa19e0..3b0bbb37 100644 --- a/_pkgdown.yml +++ b/_pkgdown.yml @@ -77,8 +77,8 @@ reference: desc: > Functions for conditional power. contents: - - gs_cp_ahr - gs_cp_npe + - gs_cp_npe2 - title: "Input definition" desc: > Helper functions to define inputs for study design. diff --git a/man/gs_cp_ahr.Rd b/man/gs_cp_ahr.Rd deleted file mode 100644 index ccbf69ff..00000000 --- a/man/gs_cp_ahr.Rd +++ /dev/null @@ -1,53 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/gs_cp_ahr.R -\name{gs_cp_ahr} -\alias{gs_cp_ahr} -\title{Conditional power computation for group sequential designs using the average hazard ratio method.} -\usage{ -gs_cp_ahr(x = NULL, i = 1, z_i = 0, hr_future = NULL) -} -\arguments{ -\item{x}{A design object of type \code{gs_design_ahr()}, \code{gs_power_ahr()}, or \code{gs_update_ahr()}.} - -\item{i}{Analysis at which interim z-value is given; must be from 1 to \code{nrow(x$analysis) - 1}.} - -\item{z_i}{Interim z-value at analysis \code{i} (scalar).} - -\item{hr_future}{A scalar or vector, which specifies the natural parameter for treatment effect at the future analyses. -The default gives you a range of possible HRs from 0.6 to 1.} -} -\value{ -A data frame showing the conditional power based on the input design, interim test statistic, and possible future HR. -} -\description{ -Conditional power computation for group sequential designs using the average hazard ratio method. -} -\examples{ -library(gsDesign) -library(gsDesign2) - -# Original design -x <- gs_design_ahr(info_frac = c(0.5, 0.8, 1)) - -# Updated design -xu <- gs_update_ahr( - x = x, - ustime = c(170, x$analysis$event[2:3] )/ max(x$analysis$event), - event_tbl = data.frame(analysis = c(1, 1), event = c(100, 70))) - -# We assume IA1 p-value of 0.04, one-sided, and future HR is possible from 0.6 to 1. -gs_cp_ahr(x = xu, - i = 1, - z_i = -qnorm(0.04), - hr_future = seq(0.6, 1, by = 0.1)) - -# We assume the test statistics z_i is derived by the function hrn2z(), which translates -# a HR and number of events into an approximate corresponding Z-value, -# using the Schoenfeld approximation. -# The future HR is possible from 0.6 to 1. -gs_cp_ahr(x = xu, - i = 1, - z_i = -hrn2z(hr = 0.7, n = xu$analysis$event[1], ratio = 1), - hr_future = seq(0.6, 1, by = 0.1)) - -} diff --git a/man/gs_cp_npe2.Rd b/man/gs_cp_npe2.Rd new file mode 100644 index 00000000..0fdbb627 --- /dev/null +++ b/man/gs_cp_npe2.Rd @@ -0,0 +1,79 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/gs_cp_npe2.R +\name{gs_cp_npe2} +\alias{gs_cp_npe2} +\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_npe2(theta = NULL, t = NULL, info = NULL, a = NULL, b = NULL, c = NULL) +} +\arguments{ +\item{theta}{A vector of j-i+1, 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 an interim analysis i+1. +... +The last element of \code{theta} is the treatment effect of a future analysis j.} + +\item{t}{A vector of j-i+1, which specifies the information fraction under the treatment effect \code{theta}.} + +\item{info}{A vector of j-i+1, which specifies the statistical information under the treatment effect \code{theta}.} + +\item{a}{A vector of length j-i-1, which specifies the futility bounds from analysis i+1 to analysis j-1.} + +\item{b}{A vector of length j-i, which specifies the efficacy bounds from analysis i+1 to analysis j.} + +\item{c}{Interim z-value at analysis i (scalar).} +} +\value{ +A list of conditional powers: prob_alpha = \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 = \eqn{P(\{Z_j \geq b_j\} \& \{\cap_{m=i+1}^{j-1} Z_m < b_m\} \mid Z_i = c)}. +prob_beta = eqn{P(\{Z_j \leq b_j\} \& \{\cap_{m=i+1}^{j-1} a_m \leq Z_m < b_m\} \mid Z_i = c)}. +} +\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{ + +} +\examples{ +library(gsDesign2) +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), + b = c(1.8, 2.1), + c = 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.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() +# update design based on IA1 observed data ---- +# assume at IA1, we observed 150 events (planned of 140) during the non-effect period, and 180 events (planned of 170) afterwards. +event_tbl <- data.frame(analysis = c(1, 1), event = c(150, 180)) +ustime <- c(150+180, x$analysis$event[2:3]) / max(x$analysis$event) +xu <- gs_update_ahr(x, alpha = NULL, ustime = ustime, lstime = NULL, event_tbl = event_tbl) +# calculate conditional power at IA1 based on non-constant treatment effect ---- +gs_cp_npe2(theta = x$analysis$theta, + t = x$analysis$info_frac, + info = x$analysis$info, + a = rep(-Inf, 3), + b = x$bound$z[x$bound$bound == "upper"], + c = -gsDesign::hrn2z(hr = 0.8, n = 150+180, ratio = 1)) +} From 1c6a86aafb10ca1f84bf6bceb7c0cd384d75e648 Mon Sep 17 00:00:00 2001 From: Shiyu Zhang Date: Thu, 18 Sep 2025 18:17:30 -0400 Subject: [PATCH 10/29] added the incremental probability --- R/gs_cp_npe2.R | 232 ++++++++++++++++++++++++++----------------------- 1 file changed, 121 insertions(+), 111 deletions(-) diff --git a/R/gs_cp_npe2.R b/R/gs_cp_npe2.R index b88751d9..dcc5a654 100644 --- a/R/gs_cp_npe2.R +++ b/R/gs_cp_npe2.R @@ -39,9 +39,9 @@ #' @param a A vector of length j-i-1, which specifies the futility bounds from analysis i+1 to analysis j-1. #' @param b A vector of length j-i, which specifies the efficacy bounds from analysis i+1 to analysis j. #' @param c Interim z-value at analysis i (scalar). -#' @return A list of conditional powers: prob_alpha = \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 = \eqn{P(\{Z_j \geq b_j\} \& \{\cap_{m=i+1}^{j-1} Z_m < b_m\} \mid Z_i = c)}. -#' prob_beta = eqn{P(\{Z_j \leq b_j\} \& \{\cap_{m=i+1}^{j-1} a_m \leq Z_m < b_m\} \mid Z_i = c)}. +#' @return 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)}. #' @export #' #' @examples @@ -101,132 +101,142 @@ gs_cp_npe2 <- function(theta = NULL, 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) != 1) # note that the input of a could be NULL + stop("The length of b should equal to length of a + 1. ") # ------------------------------ # # Initialization # ------------------------------ # - i <- 1 + # the conditional power is calculated from analysis i to analysis j # the analysis j is decided by the length of b (efficacy bound) - j <- length(b) # let D_m = B_m - B_i, where m = i+1, i+2, ..., j - dim <- j - i - - # ------------------------------ # - # Build the asymptotic - # mean of B_j - B_i - # vector of length dim - # ------------------------------ # - - mu <- sapply(seq_len(dim), function(k){ - idx <- i + k - theta[idx] * sqrt(t[idx] * info[idx]) - theta[i] * sqrt(t[i] * info[i]) - }) - - # ------------------------------ # - # Build the asymptotic - # covariance of B_j - B_i - # matrix of (dim x dim) - # ------------------------------ # - Sigma <- matrix(0, nrow = dim, ncol = dim) - - for(k in seq_len(dim)) { - for(l in seq_len(dim)) { - Sigma[k, l] <- t[i + min(k, l)] - t[i] + dim <- length(b) # = j-i + prob_alpha <- rep(0, dim) + prob_alpha_plus <- rep(0, dim) + prob_beta <- rep(0, dim) + + for(x in 1:dim){ + # x ranges from 1 to j-i, represents cases for alpha_i,i+1, ..., alpha_i,j-1, alpha_i,j + + # ------------------------------ # + # 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 of B_x - B_i + # matrix of (x x x) + # ------------------------------ # + Sigma <- matrix(0, nrow = x, ncol = x) + + for(k in seq_len(x)) { + for(l in seq_len(x)) { + Sigma[k, l] <- t[1 + min(k, l)] - t[1] + } } - } - - # ------------------------------ # - # Calculate integration limit - # for alpha(i,j) - # ------------------------------ # - # alpha(i,j) is the conditional probability of first cross efficacy bound - # at analysis j given no efficacy/futility bound crossing between analysis i amd j - # D_m = B_m - B_i, where m = i+1, i+2, ..., j - # for D_{i+1},...,D_{j-1} use [a, b); - # for D_j use [b_j, +Inf) - # integration lower bound - lower_alpha <- rep(0, j - i) - for (m in 1:(j - i - 1)) { - ## ?? here looks wrong: when m = 1, it is a[1]*sqrt(t[1]) - c*sqrt(t[i]) where i is always 1 -- which is wrong - ## when m = 1, it is B2 - B1, which should be a[2]*sqrt(t[2]) - c*sqrt(t[i]) - lower_alpha[m] <- a[m] * sqrt(t[m]) - c * sqrt(t[i]) - } - lower_alpha[j-i] <- b[j-i] * sqrt(t[j]) - c * sqrt(t[i]) - - # integration upper bound - upper_alpha <- rep(0, j - i) - for(m in 1:(j - i - 1)){ - ## ?? same as above - upper_alpha[m] <- b[m] * sqrt(t[m]) - c * sqrt(t[i]) - } - upper_alpha[j-i] <- Inf - - # --------------------------------------------------------------------- - # Calculate Lower/upper bounds for alpha_plus - # first cross efficacy bound at analysis j given no efficacy bound - # crossing at analysis i+1 to j-1 - # --------------------------------------------------------------------- - # Integration limits: D_m = B_m - B_i - # for D_{i+1},...,D_{j-1} use (-Inf, b); for D_j use [b_j, +Inf) - # lower bound - lower_alpha_plus <- rep(-Inf, j - i) - lower_alpha_plus[j-i] <- b[j-i] * sqrt(t[j]) - c * sqrt(t[i]) - - # upper bound - upper_alpha_plus <- rep(0, j - i) - for(m in 1:(j-i-1)){ - upper_alpha_plus[m] <- b[m] * sqrt(t[m]) - c * sqrt(t[i]) - } - upper_alpha_plus[j-i] <- Inf + # -------------------------------------------- # + # Calculate integration limit + # for alpha(i, x), x = i + 1, ..., j + # -------------------------------------------- # + # alpha(i, x) is the conditional probability of first cross efficacy bound + # at analysis x given no efficacy/futility bound crossing before + # D_m = B_m - B_i, where m = i+1, i+2, ..., j + # for D_{i+1},...,D_{j-1} use [a, b); for D_j use [b_j, +Inf) + # integration lower bound + lower_alpha <- rep(0, x) + for (m in 1:(x - 1)) { + ## ?? here looks wrong: when m = 1, it is a[1]*sqrt(t[1]) - c*sqrt(t[i]) where i is always 1 -- which is wrong + ## when m = 1, it is B2 - B1, which should be a[2]*sqrt(t[2]) - c*sqrt(t[i]) + ## corrected - since 'a' is vector of length (j-i-1) that specifies futility bounds from analysis i+1 to analysis j-1. + ## for example, j = 4 and i = 2, then 'a' contains futility bound at analysis #3, and t[1] here is the IF for analysis i + ## 't' is a vector of length j-i+1, indicating the IF for analysis i to j. + lower_alpha[m] <- a[m] * sqrt(t[m + 1]) - c * sqrt(t[1]) + } + lower_alpha[x] <- b[x] * sqrt(t[dim + 1]) - c * sqrt(t[1]) + + # integration upper bound + # Note: should we consider the case where j = 4, i = 3, in this case, dim = 1, length of a is 0 --> we simply integrate from bj (the only element in b) to Inf + upper_alpha <- rep(0, x) + for(m in 1:(x - 1)){ + ## ?? same as above + upper_alpha[m] <- b[m] * sqrt(t[m + 1]) - c * 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 = Sigma)[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, b); for D_j use [b_j, +Inf) + # lower bound + lower_alpha_plus <- rep(-Inf, x - 1) + lower_alpha_plus[x] <- b[x] * sqrt(t[x + 1]) - c * sqrt(t[1]) + + # upper bound + upper_alpha_plus <- rep(0, x) + for(m in 1:(x - 1)){ + upper_alpha_plus[m] <- b[m] * sqrt(t[m + 1]) - c * 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 = Sigma)[1] + + # --------------------------------------------------------------------------------------- # + # Calculate Lower/upper bounds for beta + # Not 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 [-Inf, b_j) + # lower bound + lower_beta <- rep(0, x) + for(m in 1:(x - 1)){ + lower_beta[m] <- a[m] * sqrt(t[m + 1]) - c * sqrt(t[1]) + } + lower_beta[x] <- -Inf + # upper bound + upper_beta <- rep(0, x) + for(m in 1:x){ + upper_beta[m] <- b[m] * sqrt(t[m + 1]) - c * sqrt(t[1]) + } - # --------------------------------------------------------------------- - # Calculate Lower/upper bounds for beta - # Not cross efficacy bound at analysis j given no efficacy/futility - # bound crossing at analysis i+1 to j-1 - # --------------------------------------------------------------------- - # Integration limits: D_m = B_m - B_i - # for D_{i+1},...,D_{j-1} use [a, b); for D_j use [-Inf, b_j) - # lower bound - lower_beta <- rep(0, j - i) - for(m in 1:(j-i-1)){ - lower_beta[m] <- a[m] * sqrt(t[m]) - c * sqrt(t[i]) - } - lower_beta[j-i] <- -Inf + prob_beta[x] <- mvtnorm::pmvnorm( + lower = lower_beta, + upper = upper_beta, + mean = mu, + sigma = Sigma)[1] - # upper bound - upper_beta <- rep(0, j - i) - for(m in 1:(j-i)){ - upper_beta[m] <- b[m] * sqrt(t[m]) - c * sqrt(t[i]) } - # Compute multivariate normal probability - prob_alpha <- mvtnorm::pmvnorm( - lower = lower_alpha, - upper = upper_alpha, - mean = mu, - sigma = Sigma)[1] - - prob_alpha_plus <- mvtnorm::pmvnorm( - lower = lower_alpha_plus, - upper = upper_alpha_plus, - mean = mu, - sigma = Sigma)[1] - - prob_beta <- mvtnorm::pmvnorm( - lower = lower_beta, - upper = upper_beta, - mean = mu, - sigma = Sigma)[1] - - return(prob = list(prob_alpha = prob_alpha, prob_alpha_plus = prob_alpha_plus, prob_beta = prob_beta)) - } From ba017715741ab735272bde9ec9a51c571e2dda83 Mon Sep 17 00:00:00 2001 From: LittleBeannie Date: Mon, 29 Sep 2025 13:30:05 -0400 Subject: [PATCH 11/29] minor update of matrix index --- R/gs_cp_npe2.R | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/R/gs_cp_npe2.R b/R/gs_cp_npe2.R index dcc5a654..9c674b7a 100644 --- a/R/gs_cp_npe2.R +++ b/R/gs_cp_npe2.R @@ -119,14 +119,14 @@ gs_cp_npe2 <- function(theta = NULL, prob_beta <- rep(0, dim) for(x in 1:dim){ - # x ranges from 1 to j-i, represents cases for alpha_i,i+1, ..., alpha_i,j-1, alpha_i,j + # 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. @@ -149,9 +149,9 @@ gs_cp_npe2 <- function(theta = NULL, # Calculate integration limit # for alpha(i, x), x = i + 1, ..., j # -------------------------------------------- # - # alpha(i, x) is the conditional probability of first cross efficacy bound - # at analysis x given no efficacy/futility bound crossing before - # D_m = B_m - B_i, where m = i+1, i+2, ..., j + # alpha_{i, i+x} is the conditional probability of first cross efficacy bound + # at analysis i+x given no efficacy/futility bound crossing before. + # Let's denote D_m = B_m - B_i, where m = i+1, i+2, ..., i+x # for D_{i+1},...,D_{j-1} use [a, b); for D_j use [b_j, +Inf) # integration lower bound lower_alpha <- rep(0, x) @@ -161,16 +161,16 @@ gs_cp_npe2 <- function(theta = NULL, ## corrected - since 'a' is vector of length (j-i-1) that specifies futility bounds from analysis i+1 to analysis j-1. ## for example, j = 4 and i = 2, then 'a' contains futility bound at analysis #3, and t[1] here is the IF for analysis i ## 't' is a vector of length j-i+1, indicating the IF for analysis i to j. - lower_alpha[m] <- a[m] * sqrt(t[m + 1]) - c * sqrt(t[1]) + lower_alpha[m] <- a[m] * sqrt(t[m]) - c * sqrt(t[1]) } - lower_alpha[x] <- b[x] * sqrt(t[dim + 1]) - c * sqrt(t[1]) + lower_alpha[x] <- b[x] * sqrt(t[x]) - c * sqrt(t[1]) # integration upper bound # Note: should we consider the case where j = 4, i = 3, in this case, dim = 1, length of a is 0 --> we simply integrate from bj (the only element in b) to Inf upper_alpha <- rep(0, x) for(m in 1:(x - 1)){ ## ?? same as above - upper_alpha[m] <- b[m] * sqrt(t[m + 1]) - c * sqrt(t[1]) + upper_alpha[m] <- b[m] * sqrt(t[m]) - c * sqrt(t[1]) } upper_alpha[x] <- Inf From 1d0333845f375ffb20171c21f33d6d24f5da6ff2 Mon Sep 17 00:00:00 2001 From: LittleBeannie Date: Mon, 29 Sep 2025 14:21:45 -0400 Subject: [PATCH 12/29] fix gt testing error --- tests/testthat/_snaps/independent_as_gt.md | 84 +++++++++++----------- 1 file changed, 42 insertions(+), 42 deletions(-) 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} \\ From 9b502436ebddf67ec8558d9983d574c929a10b70 Mon Sep 17 00:00:00 2001 From: LittleBeannie Date: Mon, 29 Sep 2025 14:21:59 -0400 Subject: [PATCH 13/29] revert some index back --- R/gs_cp_npe2.R | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/R/gs_cp_npe2.R b/R/gs_cp_npe2.R index 9c674b7a..e5af9690 100644 --- a/R/gs_cp_npe2.R +++ b/R/gs_cp_npe2.R @@ -138,7 +138,6 @@ gs_cp_npe2 <- function(theta = NULL, # matrix of (x x x) # ------------------------------ # Sigma <- matrix(0, nrow = x, ncol = x) - for(k in seq_len(x)) { for(l in seq_len(x)) { Sigma[k, l] <- t[1 + min(k, l)] - t[1] @@ -161,15 +160,14 @@ gs_cp_npe2 <- function(theta = NULL, ## corrected - since 'a' is vector of length (j-i-1) that specifies futility bounds from analysis i+1 to analysis j-1. ## for example, j = 4 and i = 2, then 'a' contains futility bound at analysis #3, and t[1] here is the IF for analysis i ## 't' is a vector of length j-i+1, indicating the IF for analysis i to j. - lower_alpha[m] <- a[m] * sqrt(t[m]) - c * sqrt(t[1]) + lower_alpha[m] <- a[m] * sqrt(t[m + 1]) - c * sqrt(t[1]) } - lower_alpha[x] <- b[x] * sqrt(t[x]) - c * sqrt(t[1]) + lower_alpha[x] <- b[x] * sqrt(t[dim + 1]) - c * sqrt(t[1]) # integration upper bound # Note: should we consider the case where j = 4, i = 3, in this case, dim = 1, length of a is 0 --> we simply integrate from bj (the only element in b) to Inf upper_alpha <- rep(0, x) for(m in 1:(x - 1)){ - ## ?? same as above upper_alpha[m] <- b[m] * sqrt(t[m]) - c * sqrt(t[1]) } upper_alpha[x] <- Inf From 289a333b3562c7b6673aa11711d4a955c6ca5153 Mon Sep 17 00:00:00 2001 From: LittleBeannie Date: Mon, 29 Sep 2025 14:25:53 -0400 Subject: [PATCH 14/29] fix the example --- R/gs_cp_npe2.R | 4 ++-- man/gs_cp_npe2.Rd | 10 +++++----- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/R/gs_cp_npe2.R b/R/gs_cp_npe2.R index e5af9690..68cb46cc 100644 --- a/R/gs_cp_npe2.R +++ b/R/gs_cp_npe2.R @@ -83,8 +83,8 @@ #' gs_cp_npe2(theta = x$analysis$theta, #' t = x$analysis$info_frac, #' info = x$analysis$info, -#' a = rep(-Inf, 3), -#' b = x$bound$z[x$bound$bound == "upper"], +#' a = -Inf, +#' b = x$bound$z[x$bound$bound == "upper"][2:3], #' c = -gsDesign::hrn2z(hr = 0.8, n = 150+180, ratio = 1)) gs_cp_npe2 <- function(theta = NULL, t = NULL, diff --git a/man/gs_cp_npe2.Rd b/man/gs_cp_npe2.Rd index 0fdbb627..2fa2bafd 100644 --- a/man/gs_cp_npe2.Rd +++ b/man/gs_cp_npe2.Rd @@ -24,9 +24,9 @@ The last element of \code{theta} is the treatment effect of a future analysis j. \item{c}{Interim z-value at analysis i (scalar).} } \value{ -A list of conditional powers: prob_alpha = \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 = \eqn{P(\{Z_j \geq b_j\} \& \{\cap_{m=i+1}^{j-1} Z_m < b_m\} \mid Z_i = c)}. -prob_beta = eqn{P(\{Z_j \leq b_j\} \& \{\cap_{m=i+1}^{j-1} a_m \leq Z_m < b_m\} \mid Z_i = c)}. +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)}. } \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 @@ -73,7 +73,7 @@ xu <- gs_update_ahr(x, alpha = NULL, ustime = ustime, lstime = NULL, event_tbl = gs_cp_npe2(theta = x$analysis$theta, t = x$analysis$info_frac, info = x$analysis$info, - a = rep(-Inf, 3), - b = x$bound$z[x$bound$bound == "upper"], + a = -Inf, + b = x$bound$z[x$bound$bound == "upper"][2:3], c = -gsDesign::hrn2z(hr = 0.8, n = 150+180, ratio = 1)) } From ed16761241651ffdad5ec4bb2e12763bbba5ff46 Mon Sep 17 00:00:00 2001 From: Shiyu Zhang Date: Thu, 2 Oct 2025 15:01:43 -0400 Subject: [PATCH 15/29] fix index issues --- R/gs_cp_npe2.R | 100 +++++++++++++++++++++++++++++-------------------- 1 file changed, 59 insertions(+), 41 deletions(-) diff --git a/R/gs_cp_npe2.R b/R/gs_cp_npe2.R index 68cb46cc..a56d00b0 100644 --- a/R/gs_cp_npe2.R +++ b/R/gs_cp_npe2.R @@ -83,8 +83,8 @@ #' gs_cp_npe2(theta = x$analysis$theta, #' t = x$analysis$info_frac, #' info = x$analysis$info, -#' a = -Inf, -#' b = x$bound$z[x$bound$bound == "upper"][2:3], +#' a = rep(-Inf, 3), +#' b = x$bound$z[x$bound$bound == "upper"], #' c = -gsDesign::hrn2z(hr = 0.8, n = 150+180, ratio = 1)) gs_cp_npe2 <- function(theta = NULL, t = NULL, @@ -119,14 +119,14 @@ gs_cp_npe2 <- function(theta = NULL, prob_beta <- rep(0, dim) for(x in 1:dim){ - # 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 + # x ranges from 1 to j-i, represents cases for alpha_i,i+1, ..., alpha_i,j-1, alpha_i,j # ------------------------------ # # 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. @@ -138,39 +138,47 @@ gs_cp_npe2 <- function(theta = NULL, # matrix of (x x x) # ------------------------------ # Sigma <- matrix(0, nrow = x, ncol = x) + for(k in seq_len(x)) { for(l in seq_len(x)) { Sigma[k, l] <- t[1 + min(k, l)] - t[1] } } - # -------------------------------------------- # - # Calculate integration limit - # for alpha(i, x), x = i + 1, ..., j - # -------------------------------------------- # - # alpha_{i, i+x} is the conditional probability of first cross efficacy bound - # at analysis i+x given no efficacy/futility bound crossing before. - # Let's denote D_m = B_m - B_i, where m = i+1, i+2, ..., i+x + # ----------------------------------------------------------------------------------------- # + # Calculate Lower/upper bounds for alpha + # first cross efficacy bound at analysis x given no efficacy/futility bound crossing before + # ----------------------------------------------------------------------------------------- # + # D_m = B_m - B_i # for D_{i+1},...,D_{j-1} use [a, b); for D_j use [b_j, +Inf) # integration lower bound lower_alpha <- rep(0, x) - for (m in 1:(x - 1)) { - ## ?? here looks wrong: when m = 1, it is a[1]*sqrt(t[1]) - c*sqrt(t[i]) where i is always 1 -- which is wrong - ## when m = 1, it is B2 - B1, which should be a[2]*sqrt(t[2]) - c*sqrt(t[i]) - ## corrected - since 'a' is vector of length (j-i-1) that specifies futility bounds from analysis i+1 to analysis j-1. - ## for example, j = 4 and i = 2, then 'a' contains futility bound at analysis #3, and t[1] here is the IF for analysis i - ## 't' is a vector of length j-i+1, indicating the IF for analysis i to j. - lower_alpha[m] <- a[m] * sqrt(t[m + 1]) - c * sqrt(t[1]) + if(x == 1){ + lower_alpha[x] <- b[x] * sqrt(t[x + 1]) - c * sqrt(t[1]) + }else{ + for (m in 1:(x - 1)) { + ## ?? here looks wrong: when m = 1, it is a[1]*sqrt(t[1]) - c*sqrt(t[i]) where i is always 1 -- which is wrong + ## when m = 1, it is B2 - B1, which should be a[2]*sqrt(t[2]) - c*sqrt(t[i]) + ## corrected - since 'a' is vector of length (j-i-1) that specifies futility bounds from analysis i+1 to analysis j-1. + ## for example, j = 4 and i = 2, then 'a' contains futility bound at analysis #3, and t[1] here is the IF for analysis i + ## 't' is a vector of length j-i+1, indicating the IF for analysis i to j. + lower_alpha[m] <- a[m] * sqrt(t[m + 1]) - c * sqrt(t[1]) + } + lower_alpha[x] <- b[x] * sqrt(t[x + 1]) - c * sqrt(t[1]) } - lower_alpha[x] <- b[x] * sqrt(t[dim + 1]) - c * sqrt(t[1]) # integration upper bound # Note: should we consider the case where j = 4, i = 3, in this case, dim = 1, length of a is 0 --> we simply integrate from bj (the only element in b) to Inf upper_alpha <- rep(0, x) - for(m in 1:(x - 1)){ - upper_alpha[m] <- b[m] * sqrt(t[m]) - c * sqrt(t[1]) + if(x == 1){ + upper_alpha[x] = Inf + }else{ + for(m in 1:(x - 1)){ + ## ?? same as above + upper_alpha[m] <- b[m] * sqrt(t[m + 1]) - c * sqrt(t[1]) + } + upper_alpha[x] <- Inf } - upper_alpha[x] <- Inf # Compute multivariate normal probability prob_alpha[x] <- mvtnorm::pmvnorm( @@ -184,17 +192,26 @@ gs_cp_npe2 <- function(theta = NULL, # 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, b); for D_j use [b_j, +Inf) + # for D_{i+1},...,D_{j-1} use (-Inf, bm); for D_j use [b_j, +Inf) # lower bound - lower_alpha_plus <- rep(-Inf, x - 1) - lower_alpha_plus[x] <- b[x] * sqrt(t[x + 1]) - c * sqrt(t[1]) + lower_alpha_plus <- rep(0, x) + if(x == 1){ + lower_alpha_plus[x] = b[x] * sqrt(t[x + 1]) - c * sqrt(t[1]) + }else{ + lower_alpha_plus <- rep(-Inf, x - 1) + lower_alpha_plus[x] <- b[x] * sqrt(t[x + 1]) - c * sqrt(t[1]) + } # upper bound upper_alpha_plus <- rep(0, x) - for(m in 1:(x - 1)){ - upper_alpha_plus[m] <- b[m] * sqrt(t[m + 1]) - c * sqrt(t[1]) + 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]) - c * sqrt(t[1]) + } + upper_alpha_plus[x] <- Inf } - upper_alpha_plus[x] <- Inf prob_alpha_plus[x] <- mvtnorm::pmvnorm( @@ -211,17 +228,26 @@ gs_cp_npe2 <- function(theta = NULL, # for D_{i+1},...,D_{j-1} use [a, b); for D_j use [-Inf, b_j) # lower bound lower_beta <- rep(0, x) - for(m in 1:(x - 1)){ - lower_beta[m] <- a[m] * sqrt(t[m + 1]) - c * sqrt(t[1]) + if(x == 1){ + lower_beta[x] <- -Inf + }else{ + for(m in 1:(x - 1)){ + lower_beta[m] <- a[m] * sqrt(t[m + 1]) - c * sqrt(t[1]) + } + lower_beta[x] <- -Inf } - lower_beta[x] <- -Inf # upper bound upper_beta <- rep(0, x) - for(m in 1:x){ - upper_beta[m] <- b[m] * sqrt(t[m + 1]) - c * sqrt(t[1]) + if(x == 1){ + upper_beta[x] <- b[x] * sqrt(t[x + 1]) - c * sqrt(t[1]) + }else{ + for(m in 1:x){ + upper_beta[m] <- b[m] * sqrt(t[m + 1]) - c * sqrt(t[1]) + } } + prob_beta[x] <- mvtnorm::pmvnorm( lower = lower_beta, upper = upper_beta, @@ -236,11 +262,3 @@ gs_cp_npe2 <- function(theta = NULL, prob_beta = prob_beta)) } - - - - - - - - From a080273a5ad59803afe80e63af744104fdc2e387 Mon Sep 17 00:00:00 2001 From: Shiyu Zhang Date: Thu, 2 Oct 2025 15:53:38 -0400 Subject: [PATCH 16/29] update for beta case, at jth analysis, crossing the futility bound. --- R/gs_cp_npe2.R | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/R/gs_cp_npe2.R b/R/gs_cp_npe2.R index a56d00b0..42856da1 100644 --- a/R/gs_cp_npe2.R +++ b/R/gs_cp_npe2.R @@ -119,7 +119,8 @@ gs_cp_npe2 <- function(theta = NULL, prob_beta <- rep(0, dim) for(x in 1:dim){ - # x ranges from 1 to j-i, represents cases for alpha_i,i+1, ..., alpha_i,j-1, alpha_i,j + # 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 From ef450ec490cce99bad0d431a82a6a9af27530cea Mon Sep 17 00:00:00 2001 From: LittleBeannie Date: Thu, 2 Oct 2025 17:10:15 -0400 Subject: [PATCH 17/29] add 3 cases as an example --- R/gs_cp_npe2.R | 55 +++++++++++++++++++++++++++++++++++++++++--------- 1 file changed, 46 insertions(+), 9 deletions(-) diff --git a/R/gs_cp_npe2.R b/R/gs_cp_npe2.R index 42856da1..aa91fbb6 100644 --- a/R/gs_cp_npe2.R +++ b/R/gs_cp_npe2.R @@ -36,7 +36,7 @@ #' 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-1, which specifies the futility bounds from analysis i+1 to analysis j-1. +#' @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 c Interim z-value at analysis i (scalar). #' @return 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)}. @@ -65,7 +65,7 @@ #' 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.8, 1), analysis_time = 30, +#' 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), @@ -74,17 +74,54 @@ #' h1_spending = TRUE, #' test_lower = FALSE, #' info_scale = "h0_h1_info") |> to_integer() +#' #' # update design based on IA1 observed data ---- #' # assume at IA1, we observed 150 events (planned of 140) during the non-effect period, and 180 events (planned of 170) afterwards. #' event_tbl <- data.frame(analysis = c(1, 1), event = c(150, 180)) -#' ustime <- c(150+180, x$analysis$event[2:3]) / max(x$analysis$event) +#' ustime <- c(150+180, x$analysis$event[2:4]) / max(x$analysis$event) #' xu <- gs_update_ahr(x, alpha = NULL, ustime = ustime, lstime = NULL, event_tbl = event_tbl) -#' # calculate conditional power at IA1 based on non-constant treatment effect ---- -#' gs_cp_npe2(theta = x$analysis$theta, -#' t = x$analysis$info_frac, -#' info = x$analysis$info, -#' a = rep(-Inf, 3), -#' b = x$bound$z[x$bound$bound == "upper"], +#' +#' # 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 +#' c = -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 +#' c = -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:3], +#' # IA1, IA2, IA3 and FA's information fraction +#' t = x$analysis$info_frac[1:3], +#' # IA1, IA2, IA3 and FA's statistical information +#' info = x$analysis$info[1:3], +#' # IA2, IA3 and FA's futility bound +#' a = c(-Inf, Inf) +#' # IA2, IA3 and FA's efficacy bound +#' b = x$bound$z[x$bound$bound == "upper" & x$bound$analysis %in% c(2, 3)], +#' # IA1's Z-score #' c = -gsDesign::hrn2z(hr = 0.8, n = 150+180, ratio = 1)) gs_cp_npe2 <- function(theta = NULL, t = NULL, From c58d68fc87f9af3105634e4d3349529931de2017 Mon Sep 17 00:00:00 2001 From: Shiyu Zhang Date: Wed, 8 Oct 2025 15:35:58 -0400 Subject: [PATCH 18/29] fix the data check --- R/gs_cp_npe2.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/gs_cp_npe2.R b/R/gs_cp_npe2.R index aa91fbb6..6d6db578 100644 --- a/R/gs_cp_npe2.R +++ b/R/gs_cp_npe2.R @@ -141,8 +141,8 @@ gs_cp_npe2 <- function(theta = NULL, if(length(theta) - length(b) != 1) stop("The length of theta should equal to length of b + 1. ") - if(length(b) - length(a) != 1) # note that the input of a could be NULL - stop("The length of b should equal to length of a + 1. ") + if(length(b) != length(a)) + stop("The length of b should equal to length of a. ") # ------------------------------ # # Initialization # ------------------------------ # From ebc064ca243f66c35696ee74fbb33f572f60ebc4 Mon Sep 17 00:00:00 2001 From: LittleBeannie Date: Tue, 14 Oct 2025 11:57:01 -0400 Subject: [PATCH 19/29] fix cicd error --- R/gs_cp_npe2.R | 4 ++-- man/gs_cp_npe2.Rd | 53 ++++++++++++++++++++++++++++++++++++++++------- 2 files changed, 47 insertions(+), 10 deletions(-) diff --git a/R/gs_cp_npe2.R b/R/gs_cp_npe2.R index 6d6db578..cd86b955 100644 --- a/R/gs_cp_npe2.R +++ b/R/gs_cp_npe2.R @@ -104,7 +104,7 @@ #' # IA1, IA2, and IA3's statistical information #' info = x$analysis$info[1:3], #' # IA2 and IA3's futility bound -#' a = c(-Inf, Inf) +#' 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 @@ -118,7 +118,7 @@ #' # IA1, IA2, IA3 and FA's statistical information #' info = x$analysis$info[1:3], #' # IA2, IA3 and FA's futility bound -#' a = c(-Inf, Inf) +#' a = c(-Inf, Inf), #' # IA2, IA3 and FA's efficacy bound #' b = x$bound$z[x$bound$bound == "upper" & x$bound$analysis %in% c(2, 3)], #' # IA1's Z-score diff --git a/man/gs_cp_npe2.Rd b/man/gs_cp_npe2.Rd index 2fa2bafd..8c3b1e21 100644 --- a/man/gs_cp_npe2.Rd +++ b/man/gs_cp_npe2.Rd @@ -17,7 +17,7 @@ The last element of \code{theta} is the treatment effect of a future analysis j. \item{info}{A vector of j-i+1, which specifies the statistical information under the treatment effect \code{theta}.} -\item{a}{A vector of length j-i-1, which specifies the futility bounds from analysis i+1 to analysis j-1.} +\item{a}{A vector of length j-i, which specifies the futility bounds from analysis i+1 to analysis j.} \item{b}{A vector of length j-i, which specifies the efficacy bounds from analysis i+1 to analysis j.} @@ -55,7 +55,7 @@ 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.8, 1), analysis_time = 30, + 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), @@ -64,16 +64,53 @@ x <- gs_design_ahr(enroll_rate = enroll_rate, fail_rate = fail_rate, h1_spending = TRUE, test_lower = FALSE, info_scale = "h0_h1_info") |> to_integer() + # update design based on IA1 observed data ---- # assume at IA1, we observed 150 events (planned of 140) during the non-effect period, and 180 events (planned of 170) afterwards. event_tbl <- data.frame(analysis = c(1, 1), event = c(150, 180)) -ustime <- c(150+180, x$analysis$event[2:3]) / max(x$analysis$event) +ustime <- c(150+180, x$analysis$event[2:4]) / max(x$analysis$event) xu <- gs_update_ahr(x, alpha = NULL, ustime = ustime, lstime = NULL, event_tbl = event_tbl) -# calculate conditional power at IA1 based on non-constant treatment effect ---- -gs_cp_npe2(theta = x$analysis$theta, - t = x$analysis$info_frac, - info = x$analysis$info, + +# 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, - b = x$bound$z[x$bound$bound == "upper"][2:3], + # IA2's efficacy bound + b = x$bound$z[x$bound$bound == "upper" & x$bound$analysis == 2], + # IA1's Z-score + c = -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 + c = -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:3], + # IA1, IA2, IA3 and FA's information fraction + t = x$analysis$info_frac[1:3], + # IA1, IA2, IA3 and FA's statistical information + info = x$analysis$info[1:3], + # IA2, IA3 and FA's futility bound + a = c(-Inf, Inf), + # IA2, IA3 and FA's efficacy bound + b = x$bound$z[x$bound$bound == "upper" & x$bound$analysis \%in\% c(2, 3)], + # IA1's Z-score c = -gsDesign::hrn2z(hr = 0.8, n = 150+180, ratio = 1)) } From b3a0f43d74cfd68563f1edc0bd3b97a0a1db5ee3 Mon Sep 17 00:00:00 2001 From: LittleBeannie Date: Wed, 15 Oct 2025 14:13:16 -0400 Subject: [PATCH 20/29] add developer tests to `gs_cp_npe2` --- tests/testthat/test-developer-gs_cp_npe2.R | 98 ++++++++++++++++++++++ 1 file changed, 98 insertions(+) create mode 100644 tests/testthat/test-developer-gs_cp_npe2.R 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..b689a9fb --- /dev/null +++ b/tests/testthat/test-developer-gs_cp_npe2.R @@ -0,0 +1,98 @@ +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 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 # + # --------------------------------------------------- # + 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 + c = -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]) + + # 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 + c = -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]) +}) From 1e68dc8b21eca4b1f9b1d723f1cb528f09b41019 Mon Sep 17 00:00:00 2001 From: LittleBeannie Date: Wed, 15 Oct 2025 14:15:02 -0400 Subject: [PATCH 21/29] delete some unnecessary comments --- R/gs_cp_npe2.R | 11 ----------- man/gs_cp_npe2.Rd | 6 ------ 2 files changed, 17 deletions(-) diff --git a/R/gs_cp_npe2.R b/R/gs_cp_npe2.R index cd86b955..0eb43181 100644 --- a/R/gs_cp_npe2.R +++ b/R/gs_cp_npe2.R @@ -75,12 +75,6 @@ #' test_lower = FALSE, #' info_scale = "h0_h1_info") |> to_integer() #' -#' # update design based on IA1 observed data ---- -#' # assume at IA1, we observed 150 events (planned of 140) during the non-effect period, and 180 events (planned of 170) afterwards. -#' event_tbl <- data.frame(analysis = c(1, 1), event = c(150, 180)) -#' ustime <- c(150+180, x$analysis$event[2:4]) / max(x$analysis$event) -#' xu <- gs_update_ahr(x, alpha = NULL, ustime = ustime, lstime = NULL, event_tbl = event_tbl) -#' #' # calculate conditional power #' # case 1: currently at IA1, compute conditional power at IA2 #' gs_cp_npe2(# IA1 and IA2's theta @@ -195,11 +189,6 @@ gs_cp_npe2 <- function(theta = NULL, lower_alpha[x] <- b[x] * sqrt(t[x + 1]) - c * sqrt(t[1]) }else{ for (m in 1:(x - 1)) { - ## ?? here looks wrong: when m = 1, it is a[1]*sqrt(t[1]) - c*sqrt(t[i]) where i is always 1 -- which is wrong - ## when m = 1, it is B2 - B1, which should be a[2]*sqrt(t[2]) - c*sqrt(t[i]) - ## corrected - since 'a' is vector of length (j-i-1) that specifies futility bounds from analysis i+1 to analysis j-1. - ## for example, j = 4 and i = 2, then 'a' contains futility bound at analysis #3, and t[1] here is the IF for analysis i - ## 't' is a vector of length j-i+1, indicating the IF for analysis i to j. lower_alpha[m] <- a[m] * sqrt(t[m + 1]) - c * sqrt(t[1]) } lower_alpha[x] <- b[x] * sqrt(t[x + 1]) - c * sqrt(t[1]) diff --git a/man/gs_cp_npe2.Rd b/man/gs_cp_npe2.Rd index 8c3b1e21..1cb2b2f0 100644 --- a/man/gs_cp_npe2.Rd +++ b/man/gs_cp_npe2.Rd @@ -65,12 +65,6 @@ x <- gs_design_ahr(enroll_rate = enroll_rate, fail_rate = fail_rate, test_lower = FALSE, info_scale = "h0_h1_info") |> to_integer() -# update design based on IA1 observed data ---- -# assume at IA1, we observed 150 events (planned of 140) during the non-effect period, and 180 events (planned of 170) afterwards. -event_tbl <- data.frame(analysis = c(1, 1), event = c(150, 180)) -ustime <- c(150+180, x$analysis$event[2:4]) / max(x$analysis$event) -xu <- gs_update_ahr(x, alpha = NULL, ustime = ustime, lstime = NULL, event_tbl = event_tbl) - # calculate conditional power # case 1: currently at IA1, compute conditional power at IA2 gs_cp_npe2(# IA1 and IA2's theta From 5c86264fff04904603112f6699573b88b120e514 Mon Sep 17 00:00:00 2001 From: LittleBeannie Date: Wed, 15 Oct 2025 14:31:53 -0400 Subject: [PATCH 22/29] fix cicd error on Rd files --- R/gs_cp_npe2.R | 8 ++++++-- man/gs_cp_npe2.Rd | 8 ++++++-- 2 files changed, 12 insertions(+), 4 deletions(-) diff --git a/R/gs_cp_npe2.R b/R/gs_cp_npe2.R index 0eb43181..e3810bf2 100644 --- a/R/gs_cp_npe2.R +++ b/R/gs_cp_npe2.R @@ -61,8 +61,12 @@ #' #' # 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)) +#' 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, diff --git a/man/gs_cp_npe2.Rd b/man/gs_cp_npe2.Rd index 1cb2b2f0..a304e630 100644 --- a/man/gs_cp_npe2.Rd +++ b/man/gs_cp_npe2.Rd @@ -51,8 +51,12 @@ gs_cp_npe2(theta = c(0.1, 0.2, 0.3), # 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)) +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, From c1603be459ee8d7d773ea6c8f1c621db3405c7de Mon Sep 17 00:00:00 2001 From: LittleBeannie Date: Wed, 15 Oct 2025 14:33:46 -0400 Subject: [PATCH 23/29] fix the first example error in the help page --- R/gs_cp_npe2.R | 2 +- man/gs_cp_npe2.Rd | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/R/gs_cp_npe2.R b/R/gs_cp_npe2.R index e3810bf2..f2d89a18 100644 --- a/R/gs_cp_npe2.R +++ b/R/gs_cp_npe2.R @@ -55,7 +55,7 @@ #' 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), +#' a = c(-0.5, -0.5), #' b = c(1.8, 2.1), #' c = 1.5) #' diff --git a/man/gs_cp_npe2.Rd b/man/gs_cp_npe2.Rd index a304e630..74447aa7 100644 --- a/man/gs_cp_npe2.Rd +++ b/man/gs_cp_npe2.Rd @@ -45,7 +45,7 @@ library(mvtnorm) 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), + a = c(-0.5, -0.5), b = c(1.8, 2.1), c = 1.5) From 287388d4e08d7cddb5615494c46562c2f0d55be9 Mon Sep 17 00:00:00 2001 From: Shiyu Zhang Date: Fri, 17 Oct 2025 12:16:37 -0400 Subject: [PATCH 24/29] 1. Update dim and Sigma 2. clean up the comments 3. fixed the details braces issue 4. update the beta the last upper bound, which should be the futility bound from a, not b --- R/gs_cp_npe2.R | 44 ++++++++++++++++++++------------------------ 1 file changed, 20 insertions(+), 24 deletions(-) diff --git a/R/gs_cp_npe2.R b/R/gs_cp_npe2.R index f2d89a18..942d420f 100644 --- a/R/gs_cp_npe2.R +++ b/R/gs_cp_npe2.R @@ -22,7 +22,7 @@ #' 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{Cov(Z_i, Z_j) = \sqrt{t_i/t_j} +#' \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 = c)}. @@ -141,19 +141,19 @@ gs_cp_npe2 <- function(theta = NULL, if(length(b) != length(a)) stop("The length of b should equal to length of a. ") - # ------------------------------ # + # ------------------------------ #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 - dim <- length(b) # = j-i - prob_alpha <- rep(0, dim) - prob_alpha_plus <- rep(0, dim) - prob_beta <- rep(0, dim) + 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:dim){ + 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 @@ -173,11 +173,11 @@ gs_cp_npe2 <- function(theta = NULL, # covariance of B_x - B_i # matrix of (x x x) # ------------------------------ # - Sigma <- matrix(0, nrow = x, ncol = x) + cov_matrix <- matrix(0, nrow = x, ncol = x) for(k in seq_len(x)) { for(l in seq_len(x)) { - Sigma[k, l] <- t[1 + min(k, l)] - t[1] + cov_matrix[k, l] <- t[1 + min(k, l)] - t[1] } } @@ -185,9 +185,9 @@ gs_cp_npe2 <- function(theta = NULL, # Calculate Lower/upper bounds for alpha # first cross efficacy bound at analysis x given no efficacy/futility bound crossing before # ----------------------------------------------------------------------------------------- # - # D_m = B_m - B_i + # 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) - # integration lower bound + # lower bound lower_alpha <- rep(0, x) if(x == 1){ lower_alpha[x] <- b[x] * sqrt(t[x + 1]) - c * sqrt(t[1]) @@ -198,14 +198,12 @@ gs_cp_npe2 <- function(theta = NULL, lower_alpha[x] <- b[x] * sqrt(t[x + 1]) - c * sqrt(t[1]) } - # integration upper bound - # Note: should we consider the case where j = 4, i = 3, in this case, dim = 1, length of a is 0 --> we simply integrate from bj (the only element in b) to Inf + # upper bound upper_alpha <- rep(0, x) if(x == 1){ upper_alpha[x] = Inf }else{ for(m in 1:(x - 1)){ - ## ?? same as above upper_alpha[m] <- b[m] * sqrt(t[m + 1]) - c * sqrt(t[1]) } upper_alpha[x] <- Inf @@ -216,7 +214,7 @@ gs_cp_npe2 <- function(theta = NULL, lower = lower_alpha, upper = upper_alpha, mean = mu, - sigma = Sigma)[1] + sigma = cov_matrix)[1] # --------------------------------------------------------------------------------- # # Calculate Lower/upper bounds for alpha_plus @@ -244,19 +242,18 @@ gs_cp_npe2 <- function(theta = NULL, upper_alpha_plus[x] <- Inf } - prob_alpha_plus[x] <- mvtnorm::pmvnorm( lower = lower_alpha_plus, upper = upper_alpha_plus, mean = mu, - sigma = Sigma)[1] + sigma = cov_matrix)[1] - # --------------------------------------------------------------------------------------- # + # ---------------------------------------------------------------------------------------------- # # Calculate Lower/upper bounds for beta - # Not cross efficacy bound at analysis x given no efficacy/futility bound crossing before - # --------------------------------------------------------------------------------------- # + # 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, b_j) + # 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){ @@ -274,16 +271,15 @@ gs_cp_npe2 <- function(theta = NULL, upper_beta[x] <- b[x] * sqrt(t[x + 1]) - c * sqrt(t[1]) }else{ for(m in 1:x){ - upper_beta[m] <- b[m] * sqrt(t[m + 1]) - c * sqrt(t[1]) + upper_beta[m] <- a[m] * sqrt(t[m + 1]) - c * sqrt(t[1]) } } - prob_beta[x] <- mvtnorm::pmvnorm( lower = lower_beta, upper = upper_beta, mean = mu, - sigma = Sigma)[1] + sigma = cov_matrix)[1] } From 8abb050e0c13d9adec9a21c0c9fc668022eb44d2 Mon Sep 17 00:00:00 2001 From: Shiyu Zhang Date: Tue, 21 Oct 2025 16:22:44 -0400 Subject: [PATCH 25/29] update case 3 & some other typos --- R/gs_cp_npe2.R | 23 +++++++++++------------ 1 file changed, 11 insertions(+), 12 deletions(-) diff --git a/R/gs_cp_npe2.R b/R/gs_cp_npe2.R index 942d420f..f31f1d3c 100644 --- a/R/gs_cp_npe2.R +++ b/R/gs_cp_npe2.R @@ -71,9 +71,9 @@ #' 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", +#' 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, @@ -110,15 +110,15 @@ #' #' # 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:3], +#' theta = x$analysis$theta[1:4], #' # IA1, IA2, IA3 and FA's information fraction -#' t = x$analysis$info_frac[1:3], +#' t = x$analysis$info_frac[1:4], #' # IA1, IA2, IA3 and FA's statistical information -#' info = x$analysis$info[1:3], +#' info = x$analysis$info[1:4], #' # IA2, IA3 and FA's futility bound -#' a = c(-Inf, Inf), +#' 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)], +#' b = x$bound$z[x$bound$bound == "upper" & x$bound$analysis %in% c(2, 3, 4)], #' # IA1's Z-score #' c = -gsDesign::hrn2z(hr = 0.8, n = 150+180, ratio = 1)) gs_cp_npe2 <- function(theta = NULL, @@ -168,11 +168,10 @@ gs_cp_npe2 <- function(theta = NULL, 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 of B_x - B_i - # matrix of (x x x) - # ------------------------------ # + # covariance matrix of B_x - B_i + # ---------------------------------- # cov_matrix <- matrix(0, nrow = x, ncol = x) for(k in seq_len(x)) { From 3eacea8cf52a15f12ed354653592d5ecf43d4eba Mon Sep 17 00:00:00 2001 From: LittleBeannie Date: Tue, 21 Oct 2025 16:26:53 -0400 Subject: [PATCH 26/29] add one more test --- tests/testthat/test-developer-gs_cp_npe2.R | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/tests/testthat/test-developer-gs_cp_npe2.R b/tests/testthat/test-developer-gs_cp_npe2.R index b689a9fb..e0738837 100644 --- a/tests/testthat/test-developer-gs_cp_npe2.R +++ b/tests/testthat/test-developer-gs_cp_npe2.R @@ -64,10 +64,18 @@ test_that("Compare the gs_cp_npe2 with gsDesign::gsCP", { # IA2 and FA's efficacy bound b = x_gsd$upper$bound[2:3] ) + + gsDesign2_simple_cp <- gs_cp_npe(theta = c(-log(0.8), -log(0.8)), + info = x_gsd$n.I[1:2] / 4, + a = -qnorm(0.04), + b = 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]) From e6774809819c371b221e076e80c42778faac1a3c Mon Sep 17 00:00:00 2001 From: Shiyu Zhang Date: Tue, 11 Nov 2025 12:23:25 -0500 Subject: [PATCH 27/29] 1. update gs_cp_npe to gs_cp_npe1 2. Create gs_cp_npe to call gs_cp_npe1 and gs_cp_npe2 --- R/gs_cp_npe.R | 148 +++++++++++++++++++++++++++++++++++-------------- R/gs_cp_npe1.R | 81 +++++++++++++++++++++++++++ 2 files changed, 188 insertions(+), 41 deletions(-) create mode 100644 R/gs_cp_npe1.R diff --git a/R/gs_cp_npe.R b/R/gs_cp_npe.R index cf0f2a93..84c58b44 100644 --- a/R/gs_cp_npe.R +++ b/R/gs_cp_npe.R @@ -17,67 +17,133 @@ # along with this program. If not, see . #' 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 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) } + + + + diff --git a/R/gs_cp_npe1.R b/R/gs_cp_npe1.R new file mode 100644 index 00000000..f3316467 --- /dev/null +++ b/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 . + +#' 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 = 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) +} From bb3437338ca29b041b1b0ac7f18688047c7fbed3 Mon Sep 17 00:00:00 2001 From: Shiyu Zhang Date: Tue, 11 Nov 2025 13:27:08 -0500 Subject: [PATCH 28/29] update the test: gs_cp_npe is updated to gs_cp_npe2 --- ...gs_cp_npe.R => test-developer-gs_cp_npe1.R} | 18 +++++++++--------- tests/testthat/test-developer-gs_cp_npe2.R | 4 ++-- 2 files changed, 11 insertions(+), 11 deletions(-) rename tests/testthat/{test-developer-gs_cp_npe.R => test-developer-gs_cp_npe1.R} (89%) diff --git a/tests/testthat/test-developer-gs_cp_npe.R b/tests/testthat/test-developer-gs_cp_npe1.R similarity index 89% rename from tests/testthat/test-developer-gs_cp_npe.R rename to tests/testthat/test-developer-gs_cp_npe1.R index 959370ec..94ae168f 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 # # ------------------------------ # @@ -63,26 +63,26 @@ test_that("Compare the gs_cp_npe with gsDesign::gsCP", { # ----------------------------------------- # # conditional power by gs_cp_npe under H0 # # ----------------------------------------- # - cp12_0 <- gs_cp_npe(theta = c(0,0), + cp12_0 <- gs_cp_npe1(theta = c(0,0), info = x$analysis$info0[c(1,2)], - a = -qnorm(0.04), + c = -qnorm(0.04), b = 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), + c = -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)], + cp12_1 <- gs_cp_npe1(theta = x$analysis$theta[c(1,2)], info = x$analysis$info[c(1,2)], - a = -qnorm(0.04), + c = -qnorm(0.04), b = 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), + c = -qnorm(0.04), b = x$bound$z[x$bound$bound == "upper" & x$bound$analysis == 3]) # ------------------------------ # diff --git a/tests/testthat/test-developer-gs_cp_npe2.R b/tests/testthat/test-developer-gs_cp_npe2.R index e0738837..0b75b4f2 100644 --- a/tests/testthat/test-developer-gs_cp_npe2.R +++ b/tests/testthat/test-developer-gs_cp_npe2.R @@ -65,9 +65,9 @@ test_that("Compare the gs_cp_npe2 with gsDesign::gsCP", { b = x_gsd$upper$bound[2:3] ) - gsDesign2_simple_cp <- gs_cp_npe(theta = c(-log(0.8), -log(0.8)), + gsDesign2_simple_cp <- gs_cp_npe1(theta = c(-log(0.8), -log(0.8)), info = x_gsd$n.I[1:2] / 4, - a = -qnorm(0.04), + c = -qnorm(0.04), b = x_gsd$upper$bound[2]) # IA2's CP given IA1 expect_equal(gsDesign_cp$upper$prob[1], From f2f90bb607bd7821064200c3592da5838d2c6e66 Mon Sep 17 00:00:00 2001 From: Shiyu Zhang Date: Tue, 11 Nov 2025 13:55:06 -0500 Subject: [PATCH 29/29] update RD file --- NAMESPACE | 2 +- man/gs_cp_npe.Rd | 87 ++++++++++++++++++++++++++++++++++------------- man/gs_cp_npe1.Rd | 47 +++++++++++++++++++++++++ man/gs_cp_npe2.Rd | 26 +++++++++----- 4 files changed, 129 insertions(+), 33 deletions(-) create mode 100644 man/gs_cp_npe1.Rd diff --git a/NAMESPACE b/NAMESPACE index 83cdcfdf..95e765a8 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -30,7 +30,7 @@ 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) diff --git a/man/gs_cp_npe.Rd b/man/gs_cp_npe.Rd index b4823984..adb75aa3 100644 --- a/man/gs_cp_npe.Rd +++ b/man/gs_cp_npe.Rd @@ -2,46 +2,87 @@ % 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} +\title{Conditional power computation with non-constant effect size +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{Cov(Z_i, Z_j) = \sqrt{t_i/t_j}} +See https://merck.github.io/gsDesign2/articles/story-npe-background.html for assumption details. +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)}.} \usage{ -gs_cp_npe(theta = NULL, info = NULL, a = NULL, b = NULL) +gs_cp_npe( + theta = NULL, + t = NULL, + info = NULL, + a = NULL, + b = NULL, + c = NULL, + simple_cp = FALSE +) } \arguments{ -\item{theta}{A vector of length two, which specifies the natural parameter for treatment effect. +\item{theta}{For simple conditional power, \code{theta} is 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.} +The second element of \code{theta} is the treatment effect of a future analysis j. + +\if{html}{\out{
}}\preformatted{ 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. +}\if{html}{\out{
}}} + +\item{t}{For non-simple conditional power, \code{t} is a vector of j-i+1, which specifies the information fraction under the treatment effect \code{theta}.} -\item{info}{A vector of two, which specifies the statistical information under the treatment effect \code{theta}.} +\item{info}{For simple conditional power, \code{info} is a vector of length two, which specifies the statistical information under the treatment effect \code{theta}. +For non-simple conditional power, \code{info} is a vector of j-i+1, which specifies the statistical information under the treatment effect \code{theta}.} -\item{a}{Interim z-value at analysis i (scalar).} +\item{a}{For non-simple conditional power, \code{a} is a vector of length j-i, which specifies the futility bounds from analysis i+1 to analysis j.} -\item{b}{Future target z-value at analysis j (scalar).} +\item{b}{For simple conditional power, \code{b} is a scalar which specifies the future target z-value at analysis j. +For non-simple conditional power, \code{b} is a vector of length j-i, which specifies the efficacy bounds from analysis i+1 to analysis j.} + +\item{c}{For both simple and non-simple conditional power, \code{c} is the interim z-value at analysis i (scalar).} } \value{ -A scalar with the conditional power \eqn{P(Z_2>b\mid Z_1=a)}. +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)}. } \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} +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 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)}. } \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) } diff --git a/man/gs_cp_npe1.Rd b/man/gs_cp_npe1.Rd new file mode 100644 index 00000000..6d16c503 --- /dev/null +++ b/man/gs_cp_npe1.Rd @@ -0,0 +1,47 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/gs_cp_npe1.R +\name{gs_cp_npe1} +\alias{gs_cp_npe1} +\title{Conditional power computation with non-constant effect size} +\usage{ +gs_cp_npe1(theta = NULL, info = NULL, c = 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 length two, which specifies the statistical information under the treatment effect \code{theta}.} + +\item{c}{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 = c)}. +} +\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 = c) = 1 - \Phi\left(\frac{b - \sqrt{t}c - \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_npe1(theta = c(.1, .2), + info = c(15, 35), + c = 1.5, b = 1.96) +} diff --git a/man/gs_cp_npe2.Rd b/man/gs_cp_npe2.Rd index 74447aa7..92949eda 100644 --- a/man/gs_cp_npe2.Rd +++ b/man/gs_cp_npe2.Rd @@ -32,7 +32,15 @@ prob_beta is a vector of c(beta_i,i+1, ..., beta_i,j-1, beta_i,j) where beta_i,j 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-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{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 = 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)}. } \examples{ library(gsDesign2) @@ -61,9 +69,9 @@ 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", + 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, @@ -100,15 +108,15 @@ gs_cp_npe2(# IA1, IA2 and IA3's theta # 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:3], + theta = x$analysis$theta[1:4], # IA1, IA2, IA3 and FA's information fraction - t = x$analysis$info_frac[1:3], + t = x$analysis$info_frac[1:4], # IA1, IA2, IA3 and FA's statistical information - info = x$analysis$info[1:3], + info = x$analysis$info[1:4], # IA2, IA3 and FA's futility bound - a = c(-Inf, Inf), + 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)], + b = x$bound$z[x$bound$bound == "upper" & x$bound$analysis \%in\% c(2, 3, 4)], # IA1's Z-score c = -gsDesign::hrn2z(hr = 0.8, n = 150+180, ratio = 1)) }