### ESETT TRIAL
### Program Copyright Berry Consultants LLC, 2011, 2012
### This code may not be altered and used without permission.
### Primary author : Jason Connor
### 317-877-1084   or   jason@ConfluenceStat.com



rm(list=ls())

## All times in months
## This sets up a scenario
library(VGAM)
v = list(
  S3 = c(## There are success rates for the four groups.
         ## CHANGE THE NEXT THREE LINES TO STEP THROUGH SCENARIOS
         ## GOES **VERY** slow in scenarios where the 3 rates are equal.
     0.50,      # fPHT
     0.50,      # LVT
     0.50       # VPA
  ),
    
  # Maximum sample size & max sample size for Stage 1
  MaxN = 720,
  # Priors
  a = rep(1, 3),
  b = rep(1, 3),
  # First look and look every
  firstlook = 300,
  firststop = 400
  lookevery = 100,
  # Min to randomized
  minpr = 0.05,
  # simulations
  nsims = 100000,
  badlim = 0.25, 
  # critv to (a) for 'best'
  #          (b) for 'worst'
  #          (c) to stop for futility (i.e Pred prob a winner or loser id'd)
  #          (d) for worse than 25%
 critv = c(.975, .975, 0.05, 0.05) 
)


## This function simulates the trials, produces one big output matrix
simtrials <- function(v){
  
  co <- ppcutoffs(v$critv[3])
  
  #out.mat
  # (1) N
  # (2-4) N per group
  # (5-7) Rank as 1, 2, 3 (according to prob best)
  # (8)   Sig best (1 2 or 3 or 0 if none)
  # (9)   Sig worst (1 2 or 3 or 0 if none)
  # (10) Final conclusion
  #           1 = overall futility stop, 
  #           2 = stop early for winner
  #           3 = stop early for winner & loser
  #           4 = stop early for loser and futility (not possible in ours)
  #           5 = max overall futility
  #           6 = max and loser
  #           7 = max and winner
  #           8 = max & winner & loser
  #  (11-13) Final Pr(best)
  #  (14-16) Final Pr(2nd)
  #  (17-19) Final Pr(worst)
  #  (20-22) Successes per group
  #  (23-25) Ever drop arm? (rand goes to 0 at any pt)
  
  out.mat <- matrix(NA, nrow=v$nsims, ncol=25)
  for(s in 1:v$nsims){
    ad <- c(1,1,1)
    ## Get in time for each patient 
    ## Rand assignment for first FirstLook pts
    group <- rep(NA, v$MaxN)
    group[1:v$firstlook] <- rand.new(v$firstlook, c(1,1,1))
    y <- rep(NA, v$MaxN)
    y[1:v$firstlook] <- sim.endpoint(group[1:v$firstlook], v$S3)
    look1 <- interim(v$firstlook, y, group, v, co)
#    print(round(look1,3))
    # Track if arm every dropped
    ad <- ad * as.numeric(look1[12:14]>0)
    n.now <- v$firstlook
    print(c(s,n.now))
 ## Now loop through Stage 1
    while(look1[1]==1){
      new <- min(v$MaxN-n.now, v$lookevery)
      group[(n.now+1):(n.now+new)] <- rand.new(new, look1[12:14])
      y[(n.now+1):(n.now+new)] <- sim.endpoint(group[(n.now+1):(n.now+new)], v$S3)
      look1 <- interim(n.now+new, y, group, v, co)
#      print(round(look1,3))
      ad <- ad * as.numeric(look1[12:14]>0)
      n.now <- n.now+new
      print(c(s,n.now))
    }
    mx <- look1[3:5];    mn <- look1[6:8]
    winner <- ifelse(max(mx) > v$critv[1], (1:3)[mx==max(mx)], 0)
    loser <-  ifelse(max(mn) > v$critv[2], (1:3)[mn==max(mn)], 0)
   if(look1[2]==1){
       whystop <- 1    ## futility
     }else if(look1[2]==3){
       if(loser>0){
         whystop <- 3
       }else{
         whystop <- 2
       }
     }else if(look1[2]==2){
       if(winner==0 & loser==0){
         whystop <- 5
       }else if(winner>0 & loser>0){
         whystop <- 8
       }else if(winner>0){
         whystop <- 7
       }else if(loser>0){
         whystop <- 6
       }else{
         print("error why stop at max?")
       }
     }else{
      print("error, why did trial stop?")
     }
    out.mat[s,1:25] <- c(n.now, look1[18:20], order(mx), winner, loser,
                 whystop,look1[c(3,4,5,9,10,11,6,7,8,15,16,17)],1-ad)
}
    out.mat <- data.frame(out.mat)
    names(out.mat) <- c("N","N1","N2","N3","R1","R2","R3","Best","Worst","Dec",
                      "PrBest1","PrBest2","PrBest3","PrMid1","PrMid2","PrMid3",
                      "PrWorst1","PrWorst2","PrWorst3","X1","X2","X3","AD1","AD2","AD3")
     return(out.mat)
}   


## This function takes the output matrix from simtrials, and produces prettier summary output
sumtrial <- function(outmat){
  mat <- matrix(nrow=4, ncol=9)
  out <- table(factor(outmat[,10], levels=1:8))
#              Ntotal SDN phat Rank1  Rank2  Rank3  SigBest SigWorst Drop
#      fPHT
#      LVT
#      VPA             --
#      Total
  mat[1:3,1] <- apply(outmat[,2:4], 2, mean)
  mat[1:3,2] <- apply(outmat[,2:4], 2, sd)
  mat[1:3,3] <- c(mean(outmat[,20]/outmat[,2]), mean(outmat[,21]/outmat[,3]), mean(outmat[,22]/outmat[,4]))
  mat[1,4:6] <- table(factor(outmat[,5], levels=3:1))/dim(outmat)[1]
  mat[2,4:6] <- table(factor(outmat[,6], levels=3:1))/dim(outmat)[1]
  mat[3,4:6] <- table(factor(outmat[,7], levels=3:1))/dim(outmat)[1]
  mat[1:3,7] <- table(factor(outmat[,8], levels=1:3))/dim(outmat)[1]
  mat[1:3,8] <- table(factor(outmat[,9], levels=1:3))/dim(outmat)[1]
  mat[1:3,9] <- apply(outmat[,23:25], 2, mean)
  mat[4,1] <- mean(outmat[,1])
  mat[4,2] <- sd(outmat[2])
  mat[4,3] <- mean(rowSums(outmat[,20:22]) / rowSums(outmat[2:4]))
  mat[4,4:6] <- NA
  mat[4,7] <- sum(mat[1:3,7])
  mat[4,8] <- sum(mat[1:3,8])
  mat[4,9] <- NA
  mat <- data.frame(mat)
  names(mat) <- c("N","SD","Phat","Best","Mid","Worst","SigBest","SigWorst","Drop")
  dimnames(mat)[[1]] <- c("fPHT","LVT","VPA","Total")
  return(list(out, mat))
}

  
  
# Does the interim analysis of N patients based on outcome y and group assignment 'group. 
interim <- function(N, y, group, v, co){
   ## Runs trial returns:
   #  (1) go (0=stop, 1=keep going)
   #  (2) why stop (1=3-way fut, 2=max n, 3=1 winner)
   #  (3-5) Pr each is best
   #  (6-8) Pr each is worst
   #  (9-14) x/N for each group
   #  (15-17) rand probs
   ns <- table(factor(group[1:N], levels=1:3))
   tab <- table(factor(group[1:N],levels=1:3), factor(y[1:N], levels=0:1))
   post1 <- rbeta(10000, v$a[1]+tab[1,2], v$b[1]+tab[1,1])
   post2 <- rbeta(10000, v$a[2]+tab[2,2], v$b[2]+tab[2,1])
   post3 <- rbeta(10000, v$a[3]+tab[3,2], v$b[3]+tab[3,1])
   vr <- as.numeric(( (v$a+tab[,2])*(v$b+tab[,1])) / ((v$a+v$b+ns)^2 * (v$a+v$b+ns+1)))
   top <- apply(cbind(post1,post2,post3), 1, max)
   bot <- apply(cbind(post1,post2,post3), 1, min)
 
   best <- c(mean(post1==top), mean(post2==top), mean(post3==top))
   worst <- c(mean(post1==bot), mean(post2==bot), mean(post3==bot))
   middle <- 1-best-worst

   toobad <- 1-c(pbeta(v$badlim, v$a[1]+tab[1,2], v$b[1]+tab[1,1]),
               pbeta(v$badlim, v$a[2]+tab[2,2], v$b[2]+tab[2,1]),
               pbeta(v$badlim, v$a[3]+tab[3,2], v$b[3]+tab[3,1]))

   wt <- sqrt(best * vr / as.numeric(ns))
   wt <- wt/sum(wt)
   wt[wt < v$minpr] <- 0
   wt[toobad < v$critv[4]] <- 0
   if(sum(wt) > 0){
   wt <- wt/sum(wt)
   }
    
   ##### NEED TO ADD PRED PROBS; only do if all 3 arms left
   if((N >= v$firststop) & (N < v$MaxN) & (prod(wt>0)> 0)){
     drop <- 0
     left <- v$MaxN - N
     left <- ceiling(rep(left/3, 3))
     ns.total <- ns+left
     winlose <- 0
     counter <- 1
     
     while((winlose < co[counter,1]) & (winlose >= co[counter,2]) & (counter < 1000)){
       y.end <- tab[,2] + rbetabin.ab(3, left, v$a+tab[,2], v$b+tab[,1])
       post1f <- rbeta(10000, v$a[1]+y.end[1], v$b[1]+ns.total[1]-y.end[1])
       post2f <- rbeta(10000, v$a[2]+y.end[2], v$b[2]+ns.total[2]-y.end[2])
       post3f <- rbeta(10000, v$a[3]+y.end[3], v$b[3]+ns.total[3]-y.end[3])
       topf <- apply(cbind(post1f,post2f,post3f), 1, max)
       botf <- apply(cbind(post1f,post2f,post3f), 1, min)
       bestf <- c(mean(post1f==topf), mean(post2f==topf), mean(post3f==topf))
       worstf <- c(mean(post1f==botf), mean(post2f==botf), mean(post3f==botf))
       winlose <- winlose + ifelse((max(bestf)>v$critv[1]) | (max(worstf)>v$critv[2]), 1, 0)
       counter <- counter + 1
#       print(c(winlose/counter, counter))
     }
     ppwin <- winlose/counter
   }else{
     drop <- 1
     ppwin <- v$critv[3]+1 #  If missing just make bigger than the crit value.
   }

   
   ## Stopping:
   if(N < v$firststop){
     go <- 1
     whystop <- NA
   }else if(N >= v$MaxN){
     go <- 0
     whystop <- 2
   }else if(max(best) > v$critv[1]){
     go <- 0
     whystop <- 3
   }else if(ppwin < v$critv[3]){
     go <- 0
     whystop <- 1
   }else if(wt[1]==0 & wt[2]==0 & wt[3]==0){
     go <- 0
     whystop <- 1
   }else{
     go <- 1
     whystop <- NA
   }
     
  return(as.numeric(c(go, whystop, best, worst, middle, wt, tab[,2], ns, ppwin, drop)))
}
     

## Simulates data for new patients using inputs group assignement and success rate (length 3)
sim.endpoint <- function(group, successrate){
  out <- rbinom(length(group), 1, successrate[group])
}

rand.new <- function(N, p, minp){
  ### Returns randomization codes (1:3) for N patients
  ### requires prob vector, p, of lenght 3.
  ### If 
  if(prod(p ==c(1,1,1))==1){
     out <- rep(sample(1:3,3), ceiling(N/3))
     out <- out[1:N]
   }else{
     out <- rep(sample(1:3, N, prob=p, replace=T))
   }
    return(out)  
} 

  
### Creates a lookup matrix to make the predictive probability stopping algorithm go faster.
### Creates a 99.9% confidence interval, then basically sees if its' highly likely that the stop rate is less than the cutoff
ppcutoffs <- function(critv){
    whenstop <- cbind(rep(0,1000),rep(0,1000))
    for(i in 50:1000){
    x <- ceiling(critv*i)
    while(as.numeric(binom.test(x,i,conf.level=0.999)$conf.int[1])<critv){
    x <- x+1 
    }
    whenstop[i,1] <- x
    }
    whenstop[1:49,1] <- whenstop[50,1]
    
    for(i in 50:1000){
    x <- ceiling(critv*i)
    while(as.numeric(binom.test(x,i,conf.level=0.999)$conf.int[1])>=critv){
    x <- x-1 
    }
    whenstop[i,2] <- x
    }
    return(whenstop)
}

  
  
  
##### This executes the code
##### Change response rates in "V" list to step through scenarios.

# Does the sims
# See "V" at the top which sets up a scenario.  
scene1 <- simtrials(v)
# Produces a summary
print(sumtrial(scene1))
# Saves 
save.image("1-early.Rdata")

#### RETURN to top and change  "S3" in object "v" to run additional scenarios.