## Inference procedure from a Simon's two-stage design.

## Common inputs are
     # p0 ... P[success] under H0
     # p1 ... P[success] under H1
     # n1 ... Stage I sample size
     # r1 ... Stage I critical value
     # nt ... Total sample size
     # rt ... Critical value at the end of stage II

## Some functions additionally require
     # x1 ... Number of responses in stage I
     # xt ... Total number of responses
     # alp ... type I error rate of a design (for computing a confidence interval)

## In Stage I, a sample of size n1 is collected.  
## If there are r1 or fewer responses, then the trial is terminated for futility.  
## Otherwise an additional sample of size (nt-n1) is collected in stage II.  
## If there are total of rt or fewer responses, then the futility is concluded at the end of stage II.
## Otherwise efficacy is concluded.  


## Explanation of the functions:

# simon.design
     # Given sample sizes and critical values and response rate under null and alternative, this 
     # gives design characteristics such as conditional power for each stage I observation.  
     # This also computes type I error rate, power and expected sample sizes under null and alternative.  
     # Note that this function is NOT necessary to run any of the functions that follow.  

# p.val.simon
     # This computes p-value from a simon's two-stage design.  
     # By construction, p-value < alpha is equivalent to rejecting H0.  

# p.val.naive
     # This computes an incorrect p-value from a simon's two-stage design.
     # Computation is carried out assuming that all the data are collected in 
     # a conventional single stage design.  

# p.val.both
    # This simply desplays the design parameters and correct and incorrect pvalues.

# conf.simon
     # This computes a (1-alpha) * 100 % confidence interval based on the data collected in
     # a two-stage design.  
     # Note that the confidence interval is one sided, i.e., the output of this function is the
     # lower bound of a confidence interval.  



## -- ## -- ## -- ##
simon.design <- function(p0, p1, n1, r1, nt, rt){

   x1 <- seq(0, n1)
   p01 <- round(dbinom(x1, n1, p0), 5)
   p11 <- round(dbinom(x1, n1, p1), 5)

   both0 <- (p01 + p11) == 0
   
   continue <- x1 > r1
   st2.cr <- replace(rt - x1, !continue, NA)

   cond.pow0 <- 1- round(pbinom(st2.cr, nt-n1, p0), 4)
   cond.pow1 <- 1- round(pbinom(st2.cr, nt-n1, p1), 4)

   n2 <- nt-n1
   out0 <- data.frame(p0, p1, n1, r1, nt, rt, n2)
   out1 <- data.frame(x1, p01, p11, continue, st2.cr, 
      cond.pow0, cond.pow1)[!both0, ]

   futility <- round(c(sum(p01[!continue]), sum(p11[!continue])), 4)
   p.contin <- round(c(sum(p01[continue]), sum(p11[continue])), 4)
   power <- round(c(sum(p01*replace(cond.pow0, is.na(cond.pow0), 0)), 
              sum(p11*replace(cond.pow1, is.na(cond.pow1), 0))), 4)
   EN <- round(c(n1 + n2 * p.contin[1], n1 + n2 * p.contin[2]), 2)

   out2 <- data.frame(futility, p.contin, power, EN,
      row.names=c('null', 'alt'))

   list(out0, out1, out2)}
## -- ## -- ## -- ##

## -- ## -- ## -- ##
p.val.simon <- function(p0, p1, n1, r1, nt, rt, x1, x2){
   
   # Correct p value from Simon's two-stage designs.

   p.val1 <- function(p0, p1, n1, r1, nt, rt, x1, x2){
      # Correct p value from Simon's two-stage designs.
      # Only for terminated in stage I (i.e., x1 <= r1).
      round(1- pbinom(x1-1, n1, p0), 4)}

   p.val2 <- function(p0, p1, n1, r1, nt, rt, x1, x2){
      # Correct p value from Simon's two-stage designs.
      # Only for terminated in stage II (i.e., x1 > r1).

      n2 <- nt - n1
      tot.x <- x1 + x2
      possible.n1 <- seq(r1+1, n1) ; lpn <- length(possible.n1)

      temp.cond <- rep(0, lpn)

      for(i in 1:lpn){temp.cond[i] <- 1-pbinom(tot.x-possible.n1[i]-1, n2, p0)}
      temp.unco <- temp.cond * dbinom(possible.n1, n1, p0)   

      round(sum(temp.unco), 4)}

   which.f <- list(p.val1, p.val2)[[1 + (x1 > r1)]]
   which.f(p0, p1, n1, r1, nt, rt, x1, x2)}
## -- ## -- ## -- ##

## -- ## -- ## -- ##
p.val.naive <- function(p0, p1, n1, r1, nt, rt, x1, x2){

   # Incorrect naive p value from Simon's two-stage designs.

   xt <- x1 + x2 * (x1 > r1)
   n2 <- nt - n1
   nt <- n1 + n2 * (x1 > r1)

   ptilda <- xt / nt
   1 - round(pbinom(xt-1, nt, p0), 4)}
## -- ## -- ## -- ##

## -- ## -- ## -- ##
p.val.both <- function(p0, p1, n1, r1, nt, rt, x1, x2){
   pvs <- p.val.simon(p0, p1, n1, r1, nt, rt, x1, x2)
   pvn <- p.val.naive(p0, p1, n1, r1, nt, rt, x1, x2)

   xt <- x1 + x2 * (x1 > r1)
   n2 <- nt - n1
   nt <- n1 + n2 * (x1 > r1)

   cbind(n1, r1, nt, rt, x1, xt, n1, nt, pvs, pvn)}
## -- ## -- ## -- ##

## -- ## -- ## -- ##
conf.simon <- function(p0, p1, n1, r1, nt, rt, x1, x2, alp){
   
   if(x1 == 0) {stop('\n', "can't compute a confidence interval when x1=0", '\n', '\n')}
   if(x1 > n1) {stop('\n', "x1 > n1", '\n', '\n')}
   if( (x1 > r1) & (x2 > nt-n1) ) {stop('\n', "x2 > n2", '\n', '\n')}

   cont <- x1 > r1 ; n2 <- nt - n1
   n2 <- replace(n2, !cont, 0) 
   x2 <- replace(x2, !cont, 0)
   nt <- n1 + n2

   f <- function(p0, p1, n1, r1, nt, rt, x1, x2, alp){
      p.val.simon(p0, p1, n1, r1, nt, rt, x1, x2) - alp}

   uniroot(f, c(0, 1), p1=p1, n1=n1, r1=r1, nt=nt, rt=rt, x1=x1, x2=x2, alp=alp)$root}
## -- ## -- ## -- ##