#   In this file we calculate the operating characteristics for the design
#      with a lag time.  Source this file to see the output.

#### INPUT PARAMETERS :

postCut <- 0.96            # Posterior probability threshold at final analysis
hypothesisP <- 0.50         # The null hypothesis value (0.50 in example)
nCuts <- c(50,75,100)  # The sample size cut points for the analysis
predSuccessCut <- 0.90      # Stop for SUCCESS
predFutilCut <- 0.05        # Stop for futility (Pred prob)
nLag <- 20                  # The lag time in observing subjects

nsims <- 100000           # Number of simulated trials
simulateP <- 0.65      # The simulated probability of success

###################################################

### These are calculated values

k <- length(nCuts)     # The number of 'looks' including the 'final'

## Here we find the cuts for winning at each nCuts value
winCuts <- numeric(k)
for (i in 1:k){
   quants <- qbeta(1-postCut,1 + 0:nCuts[i], 1 + nCuts[i] - (0:nCuts[i]))
   winCuts[i] <- min((0:nCuts[i])[quants > hypothesisP],nCuts[i]+1)
}

## The dbetabinom() needed
dbetabinom <- function(x,n=10,a=1,b=1){
   num <- lgamma(a+b)+lgamma(n+1)+lgamma(x+a)+lgamma(n-x+b)
   den <- lgamma(a)+lgamma(b)+lgamma(x+1)+lgamma(n-x+1)+lgamma(n+a+b)
   prob <- exp(num - den)
   prob
}

## Result trackers...
win <- logical(nsims)
ss <- numeric(nsims)

## Simulate a trial

for (i in 1:nsims){
   x <- rbinom(nCuts[k],1,simulateP)      # The k subject results, in order
   iresult <- 0                           # Tracks interim result win (=1) or loss (=2)
   result <- 0                            # Tracks final result win (=1) or loss (=2)
   look <- 0
   while (result == 0){                   # This loops through looks 
      look = look + 1
      if (look == k){
         result <- ifelse(sum(x) >= winCuts[k],1,2)
      }
      else
      {
         nCurrent <- max(0,nCuts[look] - nLag)
         xCurrent <- sum(x[1:nCurrent])
         predCurrent <- numeric(nCuts[look]+1)      #All possible at look 
         predCurrent[(xCurrent+1):(xCurrent+nLag+1)] <- dbetabinom(0:nLag,nLag,1+xCurrent,1+nCurrent-xCurrent)
         pCurrent <- sum(predCurrent[(winCuts[look]+1):(nCuts[look]+1)])
         if (pCurrent >=  predSuccessCut) iresult <- 1

         predCap <- numeric(nCuts[k]+1)             #All possible at CAP
         nMore <- nCuts[k] - nCurrent
         predCap[(xCurrent+1):(xCurrent+nMore+1)] <- dbetabinom(0:nMore,nMore,1+xCurrent,1+nCurrent-xCurrent)
         pCap <- sum(predCap[(winCuts[k]+1):(nCuts[k]+1)])
         if (pCap <  predFutilCut) iresult <- 2

         if (iresult > 0){       #Trial has stopped, look at final data
            result <- ifelse(sum(x[1:nCuts[look]]) >= winCuts[look],1,2)
         }
      }
   }

# Probability of Success
   win[i] <- (result == 1)
# Sample size
   ss[i] <- nCuts[look]
}

## Print out values
out1 <- c(length(win[win])/nsims,mean(ss),sqrt(var(ss)))
out2 <- c(length(ss[ss == 50 & win]),length(ss[ss == 75 & win]),length(ss[ss == 100 & win]))
out3 <- c(length(ss[ss == 50 & win==FALSE]),length(ss[ss == 75 &
win==FALSE]),length(ss[ss == 100 & win==FALSE]))
out2 <- out2/nsims
out3 <- out3/nsims
out2 <- rbind(out2,out3)

names(out1) <- c('Pr(win)','Mean SS','SD SS')
colnames(out2) <- as.character(nCuts)
rownames(out2) <- c('Win','Lose')

print(out1)
print(out2)