gs_spending_bound bug when IA is close to FA

[LittleBeannie]

library(survival)
library(dplyr)
library(simtrial)
library(gsDesign)
library(gsDesign2)
devtools::load_all()

# set parameters
analysisTimes <- c(44, 56, 66)
duration <- 44
sampleSize <- 385
enrollRates <- tibble(Stratum = "All", duration = duration, rate = sampleSize / duration)
T  <-  200 		        # Total duration 
eta <- -log(0.95)/12	# Dropoff rate
sfu <- sfLDOF         #Spending function
sfupar <- c(0)        # similar to OF

tswitch1 <- 2
tswitch2 <- 10
tswitch3 <- 24
tswitch4 <- 38
t <- 100
lambda1 <- -log(0.97) / 2 
lambda2 <- -log(0.4948454) / 8 #48/97;10-2
lambda3 <- -log(0.2916667) / 28 #14/48;38-10
lambda4 <- lambda3
lambda5 <- -log(10/14) / (72 - 38)

hi <- rep(c(lambda1, lambda2, lambda3, lambda4, lambda5),
          times = c(tswitch1,tswitch2-tswitch1,tswitch3-tswitch2,tswitch4-tswitch3,t-tswitch4))
S <- rep(1, times = t-1)
hr1 <- 1
hr2 <- 0.65
delay <- 2
ratio <- 1

failRates <- tibble(Stratum = "All", 
                    duration = rep(c(1),each=length(hi)), 
                    failRate = hi, 
                    hr = rep(c(hr1,hr2),times = c(delay,length(hi)-delay)), 
                    dropoutRate = eta)

ahr <- AHR(enrollRates = enrollRates,
           failRates = failRates,
           totalDuration = analysisTimes,
           ratio = 1)
events <-ahr$Events
alpha  <-  0.025 # 1-sided Type I error 
# alpha = 0.025
beta  <-  0.105 # Type II error (1 - power)
k <- 3

# run `gs_design_ahr` with a bug
gs_design_ahr(enrollRates = enrollRates,
              failRates = failRates,
              ratio = 1, alpha = alpha, beta = beta,
              analysisTimes = analysisTimes,
              upper = gs_spending_bound,
              upar = list(sf = gsDesign::sfLDOF, total_spend = alpha))

# the above bug lies in `gs_design_npe` -> `gs_power_npe` -> `gs_spending_bound`
y <- gs_info_ahr(enrollRates, failRates, ratio = ratio, events = NULL, analysisTimes=analysisTimes)
gs_design_npe(theta = y$theta,
              theta1 = y$theta,
              info = y$info,
              info0 = y$info0,
              info1 = y$info0,
              alpha = alpha,
              beta = beta,
              binding = FALSE,
              upper = gs_spending_bound,
              lower = gs_b,
              upar = list(sf = gsDesign::sfLDOF, total_spend = alpha),
              lpar = c(qnorm(.1), -Inf, -Inf),
              test_upper = TRUE,
              test_lower = TRUE,
              r = 18,
              tol = 1e-6)


gs_power_npe(theta = y$theta, 
             theta1 = y$theta,
             info = y$info,
             info0 = y$info0,
             info1 = y$info0,
             binding = FALSE,
             upper = gs_spending_bound,
             lower = gs_b,
             upar = list(sf = gsDesign::sfLDOF, total_spend = alpha),
             lpar = c(qnorm(.1), -Inf, -Inf),
             test_upper = TRUE,
             test_lower = TRUE,
             r = 18,
             tol = 1e-6)

# debug with `gs_power_npe`
# it works for 1st IA, and 2nd IA, but fails at FA
# initialize
theta = y$theta
theta1 = y$theta
info = y$info
info0 = y$info0
info1 = y$info0
binding = FALSE
upper = gs_spending_bound
lower = gs_b
upar = list(sf = gsDesign::sfLDOF, total_spend = alpha)
lpar = c(qnorm(.1), -Inf, -Inf)
test_upper = TRUE
test_lower = TRUE
r = 18
tol = 1e-6

K <- length(info)
if (length(theta) == 1 && K > 1) theta <- rep(theta, K)
if (is.null(theta1)){theta1 <- theta}else if (length(theta1)==1) theta1 <- rep(theta1,K)
if (length(test_upper) == 1 && K > 1) test_upper <- rep(test_upper, K)
if (length(test_lower) == 1 && K > 1) test_lower <- rep(test_lower, K)
a <- rep(-Inf, K)
b <- rep(Inf, K)
hgm1_0 <- NULL
hgm1_1 <- NULL
hgm1 <- NULL
upperProb <- rep(NA, K)
lowerProb <- rep(NA, K)

# k = 1
k = 1
a[k] <- lower(k = k, par = lpar, hgm1 = hgm1_1, theta = theta1, info = info1, r = r, tol = tol, test_bound = test_lower,
              efficacy = FALSE)
b[k] <- upper(k = k, par = upar, hgm1 = hgm1_0, info = info0, r = r, tol = tol, test_bound = test_upper)

upperProb[1] <- if(b[1] < Inf) {pnorm(b[1], mean = sqrt(info[1]) * theta[1], lower.tail = FALSE)}else{0}
lowerProb[1] <- if(a[1] > -Inf){pnorm(a[1], mean = sqrt(info[1]) * theta[1])}else{0}

hgm1_0 <- h1(r = r, theta = 0,         I = info0[1], a = if(binding){a[1]}else{-Inf}, b = b[1])
hgm1_1 <- h1(r = r, theta = theta1[1], I = info1[1], a = a[1], b = b[1])
hgm1   <- h1(r = r, theta = theta[1],  I = info[1],  a = a[1], b = b[1])

a
b
upperProb
lowerProb
# k = 2
k = 2
a[k] <- lower(k = k, par = lpar, hgm1 = hgm1_1, theta = theta1, info = info1, r = r, tol = tol, test_bound = test_lower,
              efficacy = FALSE)
b[k] <- upper(k = k, par = upar, hgm1 = hgm1_0, info = info0, r = r, tol = tol, test_bound = test_upper)

upperProb[k] <- if(b[k]< Inf){
  sum(hupdate(r = r, theta = theta[k], I = info[k], a = b[k], b = Inf,
              thetam1 = theta[k - 1], Im1 = info[k - 1], gm1 = hgm1)$h)
}else{0}

lowerProb[k] <- if(a[k] > -Inf){
  sum(hupdate(r = r, theta = theta[k], I = info[k], a = -Inf, b = a[k],
              thetam1 = theta[k - 1], Im1 = info[k - 1], gm1 = hgm1)$h)}else{0}

hgm1_0 <- hupdate(r = r, theta = 0,         I = info0[k], a = if(binding){a[k]}else{-Inf}, b = b[k],
                  thetam1 = 0,           Im1 = info0[k-1], gm1 = hgm1_0)

hgm1_1 <- hupdate(r = r, theta = theta1[k], I = info1[k], a = a[k], b = b[k],
                  thetam1 = theta1[k-1], Im1 = info1[k-1], gm1 = hgm1_1)

hgm1   <- hupdate(r = r, theta = theta[k],  I = info[k],  a = a[k], b = b[k],
                  thetam1 = theta[k-1],  Im1 = info[k-1],  gm1 = hgm1)

a
b
upperProb
lowerProb

# k = 3
k = 3
a[k] <- lower(k = k, par = lpar, hgm1 = hgm1_1, theta = theta1, info = info1, r = r, tol = tol, test_bound = test_lower,
              efficacy = FALSE)
b[k] <- upper(k = k, par = upar, hgm1 = hgm1_0, info = info0, r = r, tol = tol, test_bound = test_upper)

debug(upper)
upper(k = k, par = upar, hgm1 = hgm1_0, info = info0, r = r, tol = tol, test_bound = test_upper)
undebug()