# 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))) }