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)
+}
diff --git a/R/gs_cp_npe2.R b/R/gs_cp_npe2.R
new file mode 100644
index 00000000..f31f1d3c
--- /dev/null
+++ b/R/gs_cp_npe2.R
@@ -0,0 +1,290 @@
+# 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, 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)}.
+#' 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
+#' 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, -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.6, 0.8, 1), analysis_time = 30,
+#' binding = FALSE,
+#' upper = gs_spending_bound,
+#' upar = list(sf = sfLDOF, total_spend = 0.025, param = NULL),
+#' lower = gs_b,
+#' lpar = c(-Inf, -Inf, -Inf),
+#' h1_spending = TRUE,
+#' test_lower = FALSE,
+#' info_scale = "h0_h1_info") |> to_integer()
+#'
+#' # calculate conditional power
+#' # case 1: currently at IA1, compute conditional power at IA2
+#' gs_cp_npe2(# IA1 and IA2's theta
+#' theta = x$analysis$theta[1:2],
+#' # IA1 and IA2's information fraction
+#' t = x$analysis$info_frac[1:2],
+#' # IA1 and IA2's statistical information
+#' info = x$analysis$info[1:2],
+#' # IA2's futility bound
+#' a = -Inf,
+#' # IA2's efficacy bound
+#' b = x$bound$z[x$bound$bound == "upper" & x$bound$analysis == 2],
+#' # IA1's Z-score
+#' 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:4],
+#' # IA1, IA2, IA3 and FA's information fraction
+#' t = x$analysis$info_frac[1:4],
+#' # IA1, IA2, IA3 and FA's statistical information
+#' info = x$analysis$info[1:4],
+#' # IA2, IA3 and FA's futility bound
+#' a = c(-Inf, -Inf, -Inf),
+#' # IA2, IA3 and FA's efficacy bound
+#' b = x$bound$z[x$bound$bound == "upper" & x$bound$analysis %in% c(2, 3, 4)],
+#' # IA1's Z-score
+#' 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 checking
+ # ------------------------------ #
+ if(length(theta) != length(t))
+ stop("The input of theta should have the same length of the input of t.")
+
+ if(length(theta) != length(info))
+ stop("The input of theta should have the same length of the input of info.")
+
+ if(length(theta) - length(b) != 1)
+ stop("The length of theta should equal to length of b + 1. ")
+
+ if(length(b) != length(a))
+ stop("The length of b should equal to length of a. ")
+ # ------------------------------ #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
+ n_future_analysis <- length(b) # = j-i
+ prob_alpha <- rep(0, n_future_analysis)
+ prob_alpha_plus <- rep(0, n_future_analysis)
+ prob_beta <- rep(0, n_future_analysis)
+
+ for(x in 1:n_future_analysis){
+ # x ranges from 1 to j-i, represents cases for alpha_{i,i+1}, ..., {alpha_i,j-1}, alpha_{i,j}
+ # x is the increment from i
+
+ # ------------------------------ #
+ # Build the asymptotic
+ # mean of B_x - B_i
+ # vector of length x
+ # ------------------------------ #
+
+ mu <- sapply(seq_len(x), function(k){
+ idx <- k + 1 # i.e., start from i+1
+ theta[idx] * sqrt(t[idx] * info[idx]) - theta[1] * sqrt(t[1] * info[1]) #first element of `theta` is the treatment effect of IA i.
+ })
+
+ # ---------------------------------- #
+ # Build the asymptotic
+ # covariance matrix of B_x - B_i
+ # ---------------------------------- #
+ cov_matrix <- matrix(0, nrow = x, ncol = x)
+
+ for(k in seq_len(x)) {
+ for(l in seq_len(x)) {
+ cov_matrix[k, l] <- t[1 + min(k, l)] - t[1]
+ }
+ }
+
+ # ----------------------------------------------------------------------------------------- #
+ # Calculate Lower/upper bounds for alpha
+ # first cross efficacy bound at analysis x given no efficacy/futility bound crossing before
+ # ----------------------------------------------------------------------------------------- #
+ # Integration limits: D_m = B_m - B_i
+ # for D_{i+1},...,D_{j-1} use [a, b); for D_j use [b_j, +Inf)
+ # lower bound
+ lower_alpha <- rep(0, x)
+ if(x == 1){
+ lower_alpha[x] <- b[x] * sqrt(t[x + 1]) - c * sqrt(t[1])
+ }else{
+ for (m in 1:(x - 1)) {
+ 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])
+ }
+
+ # upper bound
+ upper_alpha <- rep(0, x)
+ if(x == 1){
+ upper_alpha[x] = Inf
+ }else{
+ for(m in 1:(x - 1)){
+ upper_alpha[m] <- b[m] * sqrt(t[m + 1]) - 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 = cov_matrix)[1]
+
+ # --------------------------------------------------------------------------------- #
+ # Calculate Lower/upper bounds for alpha_plus
+ # first cross efficacy bound at analysis x given no efficacy bound crossing before
+ # --------------------------------------------------------------------------------- #
+ # Integration limits: D_m = B_m - B_i
+ # for D_{i+1},...,D_{j-1} use (-Inf, bm); for D_j use [b_j, +Inf)
+ # lower bound
+ lower_alpha_plus <- rep(0, x)
+ if(x == 1){
+ lower_alpha_plus[x] = b[x] * sqrt(t[x + 1]) - 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)
+ 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
+ }
+
+ prob_alpha_plus[x] <- mvtnorm::pmvnorm(
+ lower = lower_alpha_plus,
+ upper = upper_alpha_plus,
+ mean = mu,
+ sigma = cov_matrix)[1]
+
+ # ---------------------------------------------------------------------------------------------- #
+ # Calculate Lower/upper bounds for beta
+ # First crossing the futility bound at analysis x given no efficacy/futility bound crossing before
+ # ---------------------------------------------------------------------------------------------- #
+ # Integration limits: D_m = B_m - B_i
+ # for D_{i+1},...,D_{j-1} use [a, b); for D_j use [-Inf, a_j)
+ # lower bound
+ lower_beta <- rep(0, x)
+ if(x == 1){
+ lower_beta[x] <- -Inf
+ }else{
+ for(m in 1:(x - 1)){
+ lower_beta[m] <- a[m] * sqrt(t[m + 1]) - c * sqrt(t[1])
+ }
+ lower_beta[x] <- -Inf
+ }
+
+ # upper bound
+ upper_beta <- rep(0, x)
+ if(x == 1){
+ upper_beta[x] <- b[x] * sqrt(t[x + 1]) - c * sqrt(t[1])
+ }else{
+ for(m in 1:x){
+ 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 = cov_matrix)[1]
+
+ }
+
+
+ return(prob = list(prob_alpha = prob_alpha,
+ prob_alpha_plus = prob_alpha_plus,
+ prob_beta = prob_beta))
+
+}
diff --git a/R/gs_update_ahr.R b/R/gs_update_ahr.R
index 023bc07f..c4d8d6d7 100644
--- a/R/gs_update_ahr.R
+++ b/R/gs_update_ahr.R
@@ -302,7 +302,7 @@ gs_update_ahr <- function(
# Tidy outputs #
# ----------------------------------- #
ans <- list()
-
+
ans$design <- x$design
ans$enroll_rate <- x$enroll_rate
diff --git a/_pkgdown.yml b/_pkgdown.yml
index 59f7f563..7f9dd0d3 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -78,6 +78,7 @@ reference:
Functions for conditional power.
contents:
- gs_cp_npe
+ - gs_cp_npe2
- title: "Input definition"
desc: >
Helper functions to define inputs for study design.
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
new file mode 100644
index 00000000..92949eda
--- /dev/null
+++ b/man/gs_cp_npe2.Rd
@@ -0,0 +1,122 @@
+% 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, 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.}
+
+\item{c}{Interim z-value at analysis i (scalar).}
+}
+\value{
+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
+}
+\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)
+library(dplyr)
+library(mvtnorm)
+# Example 1 ----
+# Calculate conditional power under arbitrary theta, info and lower/upper bound
+# In practice, the value of theta and info commonly comes from a design.
+# More examples are available at the pkgdown vignettes.
+gs_cp_npe2(theta = c(0.1, 0.2, 0.3),
+ t = c(0.15, 0.35, 0.6),
+ info = c(15, 35, 60),
+ a = c(-0.5, -0.5),
+ b = c(1.8, 2.1),
+ 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.6, 0.8, 1), analysis_time = 30,
+ binding = FALSE,
+ upper = gs_spending_bound,
+ upar = list(sf = sfLDOF, total_spend = 0.025, param = NULL),
+ lower = gs_b,
+ lpar = c(-Inf, -Inf, -Inf),
+ h1_spending = TRUE,
+ test_lower = FALSE,
+ info_scale = "h0_h1_info") |> to_integer()
+
+# calculate conditional power
+# case 1: currently at IA1, compute conditional power at IA2
+gs_cp_npe2(# IA1 and IA2's theta
+ theta = x$analysis$theta[1:2],
+ # IA1 and IA2's information fraction
+ t = x$analysis$info_frac[1:2],
+ # IA1 and IA2's statistical information
+ info = x$analysis$info[1:2],
+ # IA2's futility bound
+ a = -Inf,
+ # IA2's efficacy bound
+ b = x$bound$z[x$bound$bound == "upper" & x$bound$analysis == 2],
+ # IA1's Z-score
+ 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:4],
+ # IA1, IA2, IA3 and FA's information fraction
+ t = x$analysis$info_frac[1:4],
+ # IA1, IA2, IA3 and FA's statistical information
+ info = x$analysis$info[1:4],
+ # IA2, IA3 and FA's futility bound
+ a = c(-Inf, -Inf, -Inf),
+ # IA2, IA3 and FA's efficacy bound
+ b = x$bound$z[x$bound$bound == "upper" & x$bound$analysis \%in\% c(2, 3, 4)],
+ # IA1's Z-score
+ c = -gsDesign::hrn2z(hr = 0.8, n = 150+180, ratio = 1))
+}
diff --git a/tests/testthat/_snaps/independent_as_gt.md b/tests/testthat/_snaps/independent_as_gt.md
index 5980018c..1e83b9d5 100644
--- a/tests/testthat/_snaps/independent_as_gt.md
+++ b/tests/testthat/_snaps/independent_as_gt.md
@@ -2,9 +2,9 @@
\begin{table}[t]
\caption*{
- {\large Fixed Design under AHR Method\textsuperscript{\textit{1}}}
+ {\fontsize{20}{25}\selectfont Fixed Design under AHR Method\textsuperscript{\textit{1}}\fontsize{12}{15}\selectfont }
}
- \fontsize{12.0pt}{14.4pt}\selectfont
+ \fontsize{12.0pt}{14.0pt}\selectfont
\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}}lrrrrrrr}
\toprule
Design & N & Events & Time & AHR & Bound & alpha & Power \\
@@ -21,9 +21,9 @@
\begin{table}[t]
\caption*{
- {\large Custom Title\textsuperscript{\textit{1}}}
+ {\fontsize{20}{25}\selectfont Custom Title\textsuperscript{\textit{1}}\fontsize{12}{15}\selectfont }
}
- \fontsize{12.0pt}{14.4pt}\selectfont
+ \fontsize{12.0pt}{14.0pt}\selectfont
\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}}lrrrrrrr}
\toprule
Design & N & Events & Time & AHR & Bound & alpha & Power \\
@@ -40,9 +40,9 @@
\begin{table}[t]
\caption*{
- {\large Fixed Design under Fleming-Harrington Method\textsuperscript{\textit{1}}}
+ {\fontsize{20}{25}\selectfont Fixed Design under Fleming-Harrington Method\textsuperscript{\textit{1}}\fontsize{12}{15}\selectfont }
}
- \fontsize{12.0pt}{14.4pt}\selectfont
+ \fontsize{12.0pt}{14.0pt}\selectfont
\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}}lrrrrrrr}
\toprule
Design & N & Events & Time & AHR & Bound & alpha & Power \\
@@ -59,10 +59,10 @@
\begin{table}[t]
\caption*{
- {\large Bound summary for AHR design} \\
- {\small AHR approximations of \textasciitilde{}HR at bound}
+ {\fontsize{20}{25}\selectfont Bound summary for AHR design\fontsize{12}{15}\selectfont } \\
+ {\fontsize{14}{17}\selectfont AHR approximations of \textasciitilde{}HR at bound\fontsize{12}{15}\selectfont }
}
- \fontsize{12.0pt}{14.4pt}\selectfont
+ \fontsize{12.0pt}{14.0pt}\selectfont
\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}}lrrrrr}
\toprule
& & & & \multicolumn{2}{c}{Cumulative boundary crossing probability} \\
@@ -84,10 +84,10 @@
\begin{table}[t]
\caption*{
- {\large Bound summary for AHR design} \\
- {\small AHR approximations of \textasciitilde{}HR at bound}
+ {\fontsize{20}{25}\selectfont Bound summary for AHR design\fontsize{12}{15}\selectfont } \\
+ {\fontsize{14}{17}\selectfont AHR approximations of \textasciitilde{}HR at bound\fontsize{12}{15}\selectfont }
}
- \fontsize{12.0pt}{14.4pt}\selectfont
+ \fontsize{12.0pt}{14.0pt}\selectfont
\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}}lrrrrr}
\toprule
& & & & \multicolumn{2}{c}{Cumulative boundary crossing probability} \\
@@ -120,10 +120,10 @@
\begin{table}[t]
\caption*{
- {\large Bound summary for WLR design} \\
- {\small WLR approximation of \textasciitilde{}wHR at bound}
+ {\fontsize{20}{25}\selectfont Bound summary for WLR design\fontsize{12}{15}\selectfont } \\
+ {\fontsize{14}{17}\selectfont WLR approximation of \textasciitilde{}wHR at bound\fontsize{12}{15}\selectfont }
}
- \fontsize{12.0pt}{14.4pt}\selectfont
+ \fontsize{12.0pt}{14.0pt}\selectfont
\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}}lrrrrr}
\toprule
& & & & \multicolumn{2}{c}{Cumulative boundary crossing probability} \\
@@ -146,10 +146,10 @@
\begin{table}[t]
\caption*{
- {\large Bound summary for WLR design\textsuperscript{\textit{1}}} \\
- {\small WLR approximation of \textasciitilde{}wHR at bound}
+ {\fontsize{20}{25}\selectfont Bound summary for WLR design\textsuperscript{\textit{1}}\fontsize{12}{15}\selectfont } \\
+ {\fontsize{14}{17}\selectfont WLR approximation of \textasciitilde{}wHR at bound\fontsize{12}{15}\selectfont }
}
- \fontsize{12.0pt}{14.4pt}\selectfont
+ \fontsize{12.0pt}{14.0pt}\selectfont
\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}}lrrrrr}
\toprule
& & & & \multicolumn{2}{c}{Cumulative boundary crossing probability\textsuperscript{\textit{2}}} \\
@@ -184,10 +184,10 @@
\begin{table}[t]
\caption*{
- {\large Bound summary for MaxCombo design} \\
- {\small MaxCombo approximation}
+ {\fontsize{20}{25}\selectfont Bound summary for MaxCombo design\fontsize{12}{15}\selectfont } \\
+ {\fontsize{14}{17}\selectfont MaxCombo approximation\fontsize{12}{15}\selectfont }
}
- \fontsize{12.0pt}{14.4pt}\selectfont
+ \fontsize{12.0pt}{14.0pt}\selectfont
\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}}lrrrr}
\toprule
& & & \multicolumn{2}{c}{Cumulative boundary crossing probability} \\
@@ -220,10 +220,10 @@
\begin{table}[t]
\caption*{
- {\large Bound summary of Binary Endpoint} \\
- {\small measured by risk difference}
+ {\fontsize{20}{25}\selectfont Bound summary of Binary Endpoint\fontsize{12}{15}\selectfont } \\
+ {\fontsize{14}{17}\selectfont measured by risk difference\fontsize{12}{15}\selectfont }
}
- \fontsize{12.0pt}{14.4pt}\selectfont
+ \fontsize{12.0pt}{14.0pt}\selectfont
\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}}lrrrrr}
\toprule
& & & & \multicolumn{2}{c}{Cumulative boundary crossing probability} \\
@@ -244,10 +244,10 @@
\begin{table}[t]
\caption*{
- {\large Bound summary of Binary Endpoint} \\
- {\small measured by risk difference}
+ {\fontsize{20}{25}\selectfont Bound summary of Binary Endpoint\fontsize{12}{15}\selectfont } \\
+ {\fontsize{14}{17}\selectfont measured by risk difference\fontsize{12}{15}\selectfont }
}
- \fontsize{12.0pt}{14.4pt}\selectfont
+ \fontsize{12.0pt}{14.0pt}\selectfont
\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}}lrrrrr}
\toprule
& & & & \multicolumn{2}{c}{Cumulative boundary crossing probability} \\
@@ -278,10 +278,10 @@
\begin{table}[t]
\caption*{
- {\large Bound Summary} \\
- {\small from gs\_power\_wlr}
+ {\fontsize{20}{25}\selectfont Bound Summary\fontsize{12}{15}\selectfont } \\
+ {\fontsize{14}{17}\selectfont from gs\_power\_wlr\fontsize{12}{15}\selectfont }
}
- \fontsize{12.0pt}{14.4pt}\selectfont
+ \fontsize{12.0pt}{14.0pt}\selectfont
\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}}lrrrrr}
\toprule
& & & & \multicolumn{2}{c}{Cumulative boundary crossing probability} \\
@@ -315,10 +315,10 @@
\begin{table}[t]
\caption*{
- {\large Bound summary for WLR design} \\
- {\small WLR approximation of \textasciitilde{}wHR at bound}
+ {\fontsize{20}{25}\selectfont Bound summary for WLR design\fontsize{12}{15}\selectfont } \\
+ {\fontsize{14}{17}\selectfont WLR approximation of \textasciitilde{}wHR at bound\fontsize{12}{15}\selectfont }
}
- \fontsize{12.0pt}{14.4pt}\selectfont
+ \fontsize{12.0pt}{14.0pt}\selectfont
\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}}lrrrrr}
\toprule
& & & & \multicolumn{2}{c}{Cumulative probability to cross boundaries} \\
@@ -352,10 +352,10 @@
\begin{table}[t]
\caption*{
- {\large Bound summary for WLR design\textsuperscript{\textit{1}}} \\
- {\small WLR approximation of \textasciitilde{}wHR at bound}
+ {\fontsize{20}{25}\selectfont Bound summary for WLR design\textsuperscript{\textit{1}}\fontsize{12}{15}\selectfont } \\
+ {\fontsize{14}{17}\selectfont WLR approximation of \textasciitilde{}wHR at bound\fontsize{12}{15}\selectfont }
}
- \fontsize{12.0pt}{14.4pt}\selectfont
+ \fontsize{12.0pt}{14.0pt}\selectfont
\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}}lrrrrr}
\toprule
& & & & \multicolumn{2}{c}{Cumulative boundary crossing probability\textsuperscript{\textit{2}}} \\
@@ -390,10 +390,10 @@
\begin{table}[t]
\caption*{
- {\large Bound summary for WLR design} \\
- {\small WLR approximation of \textasciitilde{}wHR at bound}
+ {\fontsize{20}{25}\selectfont Bound summary for WLR design\fontsize{12}{15}\selectfont } \\
+ {\fontsize{14}{17}\selectfont WLR approximation of \textasciitilde{}wHR at bound\fontsize{12}{15}\selectfont }
}
- \fontsize{12.0pt}{14.4pt}\selectfont
+ \fontsize{12.0pt}{14.0pt}\selectfont
\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}}lrrrrr}
\toprule
& & & & \multicolumn{2}{c}{Cumulative boundary crossing probability} \\
@@ -424,10 +424,10 @@
\begin{table}[t]
\caption*{
- {\large Bound summary for WLR design} \\
- {\small WLR approximation of \textasciitilde{}wHR at bound}
+ {\fontsize{20}{25}\selectfont Bound summary for WLR design\fontsize{12}{15}\selectfont } \\
+ {\fontsize{14}{17}\selectfont WLR approximation of \textasciitilde{}wHR at bound\fontsize{12}{15}\selectfont }
}
- \fontsize{12.0pt}{14.4pt}\selectfont
+ \fontsize{12.0pt}{14.0pt}\selectfont
\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}}lrrrr}
\toprule
& & & \multicolumn{2}{c}{Cumulative boundary crossing probability} \\
diff --git a/tests/testthat/test-developer-gs_cp_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
new file mode 100644
index 00000000..0b75b4f2
--- /dev/null
+++ b/tests/testthat/test-developer-gs_cp_npe2.R
@@ -0,0 +1,106 @@
+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]
+ )
+
+ gsDesign2_simple_cp <- gs_cp_npe1(theta = c(-log(0.8), -log(0.8)),
+ info = x_gsd$n.I[1:2] / 4,
+ c = -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])
+
+ # --------------------------------------------------- #
+ # 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])
+})