# setwd("~/research_drive/AdaptiveTrialDesign/simulations/binaryOutcomes/")
options(width = 120, scipen = 6)
library(VGAM) # Has function for beta-binomial distribution.
############################################################################
## bionomialOneArm is created from the specs on pages example 5.2, 197-203 in the Berry book.
## It has stopping for success based on posterior probability.
##
## Includes stopping for futility based on predictive probability
## The timing of the interim looks is based on the number of patients accrued/ number of observations.
#################################################################################
binomialOneArm <- function(nreps = 100, # number of complete trials simulated
nmax = 100, # maximum number of patients accrued, which occurs when the trial does not stop early
futilitybound = 0.05,
postcut = 0.97, # The posterior probability cut off for the final analysis.
checkpoints = c(50, 75), # The times of the interim looks are based on the number of patients enrolled.
betapriorparam1 = 1,
betapriorparam2 = 1,
p_true = 0.5,
p_0 = 0.5,
seed = 50){
set.seed(seed)
####################################################################
## For each complete trial, record final sample size, reason for stopping, and final posterior probability.
finalsamplesize <- rep(NA, nreps)
finalPostProb <- rep(NA, nreps)
stopearly<- factor(rep("No", nreps), levels = c("No", "Efficacy", "Futility"))
####################################################################
# This loop finds the number of successes needed at Nmax to have a successful trial that concludes efficacy.
for(j in seq.int(nmax, 0)){
if(pbeta(q = p_0, shape1 = betapriorparam1 + j, shape2 = betapriorparam2 + nmax - j, lower.tail = FALSE) < postcut){
targetatnmax <- j + 1
break
}
}
checkpoints <- sort(checkpoints)
for(reps in seq_len(nreps)){
X <- 0 # Initialize X as the number of successes observed.
for(i in seq_len(length(checkpoints))){
currentn <- checkpoints[i] # currentn is the current number ENROLLED
X <- rbinom(n = 1L, size = currentn - ifelse(i != 1, checkpoints[i - 1], 0), prob = p_true) + X # Find the current number of observed successes.
# Decide if need to stop for efficacy based on posterior prob.
posteriorProbSuperiority <- pbeta(q = p_0, shape1 = betapriorparam1 + X, shape2 = betapriorparam2 + currentn - X, lower.tail = FALSE)
if(posteriorProbSuperiority > postcut){
stopearly[reps] <- "Efficacy"
break
}
# Decide if need to stop for futility based on predicted probability of success.
# Calculate predicted probability of success at final analysis with nmax.
# This is based on the calculation on page 143.
# If the number of successes we need for a success at Nmax is larger than the remaining trials (m), then PP = 0.
remainingSuccessesNeeded <- targetatnmax - X
m <- nmax - currentn
if(remainingSuccessesNeeded <= m){
vals <- seq.int(remainingSuccessesNeeded, m)
PP <- sum(dbetabinom.ab(x = vals, size = m, shape1 = betapriorparam1 + X, shape2 = betapriorparam2 + currentn - X))
} else{PP <- 0}
if(PP < futilitybound){ #*
#as.numeric(pbeta(q = p_0, shape1 = betapriorparam1 + X + vals, shape2 = betapriorparam2 + nmax - X - vals, lower.tail = FALSE) > postcut)) < futilitybound){
stopearly[reps] <- "Futility"
break
}
} # End of checkpoints loop
cat(paste(reps, "\n"))
# Observe remaining data until nmax, if nmax is not in checkpoints.
if(stopearly[reps] == "No" & !(nmax %in% checkpoints)){
X <- rbinom(n = 1L, size = nmax - checkpoints, prob = p_true) + X # Find the current number of observed successes.
finalsamplesize[reps] <- nmax
finalPostProb[reps] <- pbeta(q = p_0, shape1 = betapriorparam1 + X, shape2 = betapriorparam2 + nmax - X, lower.tail = FALSE)
} else {
finalPostProb[reps] <- posteriorProbSuperiority
finalsamplesize[reps] <- currentn}
} # End of reps loop
output <- c("p_true" = p_true, "Pr(win)" = mean(finalPostProb > postcut), "E(N)" = mean(finalsamplesize), "sd(N)" = sd(finalsamplesize),
"PET Futility" = mean(stopearly == "Futility"), "PET Efficacy" = mean(stopearly == "Efficacy"))
return(list(output = output, stopreasons = table(stopearly), stopbyss = prop.table(table(stopearly, finalsamplesize)), all = data.frame(stopearly, finalsamplesize, finalPostProb)))
}