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

## Common inputs are

     # DES contains the design parameter (n1, nt, r1, rt)
     # 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
     # xt ... Total number of responses
     # alp ... type I error rate of a design (for computing a confidence interval)
     # p ... any p

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


## Some main functions
     # all(DES, p0)
     # For a given design (DES) under p0, this computes 
     # pvalues (correct and naive) and confidence interval for each outcome (nt).  

     # find.p(DES)
     # For a given design (DES), this computes four estimates of p for each outcome (nt).

     # find.Epall(DES, p)
     # For a given design (DES) and the assumed true p, this computes four estimates of p and their bias, variance and MSE.  

## -- ## -- ## -- ##
px <- function(DES, p, xt){
   n1 <- DES[1] ; nt <- DES[2] ; r1 <- DES[3] ; rt <- DES[4] ; n2 <- nt - n1
   temp1 <- dbinom(xt, n1, p)
   pos.x1 <- seq(r1+1, xt)
   cor.x2 <- xt - pos.x1
   temp2 <- sum(dbinom(pos.x1, n1, p) * dbinom(cor.x2, n2, p))
   
   c(temp1, temp2)[1+(xt > r1)]}   
## -- ## -- ## -- ##

## -- ## -- ## -- ##
px.all <- function(DES, p){
   n1 <- DES[1] ; nt <- DES[2] ; r1 <- DES[3] ; rt <- DES[4] ; n2 <- nt - n1
   xtv <- pxt <- seq(0, nt) ; L <- length(xtv)
   
   for(i in 1:L){pxt[i] <- px(DES, p, xtv[i])}
   
   out <- data.frame(cbind(xtv, pxt))
   names(out) <- c('Total X', 'P[X = x]')
   out}
## -- ## -- ## -- ##

## -- ## -- ## -- ##
pval <- function(DES, p0, xt){
   temp <- px.all(DES, p0)
   
   1-sum(temp[0:xt, 2])}
## -- ## -- ## -- ##

## -- ## -- ## -- ##
pval.all <- function(DES, p0){
   
   n1 <- DES[1] ; nt <- DES[2] ; r1 <- DES[3] ; rt <- DES[4] ; n2 <- nt - n1
   xtv <- pv <- seq(0, nt) ; L <- length(xtv)
   
   for(i in 1:L){pv[i] <- pval(DES, p0, xtv[i])}
   
   out <- data.frame(cbind(xtv, pv))
   names(out) <- c('Total X', 'pvalue')
   out}
## -- ## -- ## -- ##

## -- ## -- ## -- ##
pval.naive <- function(DES, p0){

   n1 <- DES[1] ; nt <- DES[2] ; r1 <- DES[3] ; rt <- DES[4] ; n2 <- nt - n1
   xtv <- seq(0, nt)
   
   pvn <- 1 - pbinom(xtv-1, n1 + n2* (xtv > r1), p0)
   cbind(xtv, pvn)}
## -- ## -- ## -- ##

## -- ## -- ## -- ##
pval.both <- function(DES, p0, truncate=F){
   temp1 <- pval.all(DES, p0)
   temp2 <- pval.naive(DES, p0)
   x <- temp1[,1] ; pv <- temp1[,2] ; pvn <- temp2[,2]
   
   ok <- rep(T, length(x))
   if(truncate){
   ok1 <- pmin(pv, pvn) > 1e-5
   ok2 <- (1 - pmax(pv, pvn)) > 1e-5
   ok <- ok1 & ok2}
   
   out <- cbind(x, formatC(cbind(pv, pvn), format='f'))[ok, ]
   data.frame(out)}
## -- ## -- ## -- ##

## -- ## -- ## -- ##
conf <- function(DES, xt, alp=.05){
   
   LB <- UB <- 0
   if(xt!=0){
   n1 <- DES[1] ; nt <- DES[2] ; r1 <- DES[3] ; rt <- DES[4] ; n2 <- nt - n1

   f.L <- function(p0t, DESt, xtt, alpt){pval(DESt, p0t, xtt) - alpt}
   f.U <- function(p0t, DESt, xtt, alpt){pval(DESt, p0t, xtt) - (1-alpt)}
   
   LB <- uniroot(f.L, c(0, 1), DESt=DES, xtt=xt, alpt=alp)$root   
   UB <- uniroot(f.U, c(0, 1), DESt=DES, xtt=xt, alpt=alp)$root
   }
   cbind('lower'=LB, 'upper'=UB)}
## -- ## -- ## -- ##

## -- ## -- ## -- ##
conf.all <- function(DES, alp=.05){

   p0 <- p01[1] ; p1 <- p01[2]
   n1 <- DES[1] ; nt <- DES[2] ; r1 <- DES[3] ; rt <- DES[4] ; n2 <- nt - n1

   xt <- seq(0, nt) ; L <- nt + 1
   
   conf.m <- matrix(0, ncol=2, nrow=L)
   for(i in 1:L){conf.m[i,] <- conf(DES, xt[i], alp)}
   

   out <- data.frame(conf.m)
   names(out) <- c('Lower', 'Upper')
   out}
## -- ## -- ## -- ##

## -- ## -- ## -- ##   
all <- function(DES, p0, alp=.05){

   px <- formatC(px.all(DES, p0)[,2], format='f')
   temp1 <- pval.both(DES, p0)
   temp2 <- round(conf.all(DES, alp), 5)
   
   temp3 <- data.frame(cbind(temp1[,1], px, temp1[,c(2,3)], temp2))
   names(temp3) <- c('xt', 'px', 'correct', 'naive', 'lower', 'upper')
   temp3}
## -- ## -- ## -- ##

## -- ## -- ## -- ##
find.phat <- function(DES){

   n1 <- DES[1] ; nt <- DES[2] ; r1 <- DES[3] ; rt <- DES[4] ; n2 <- nt - n1

   xtv <- seq(0, nt)
   phat <- xtv / (n1 + n2* (xtv > r1))
   cbind(xtv, phat)}   
## -- ## -- ## -- ##

## -- ## -- ## -- ##
find.phat1 <- function(DES, xt){
   temp <- find.phat(DES)
   temp[,2][temp[,1] == xt]}
## -- ## -- ## -- ##

## -- ## -- ## -- ##
find.Ep <- function(DES, p){
   temp <- find.phat(DES)
   pxx <- px.all(DES, p)
   sum(temp[,2] * pxx[,2])}   
## -- ## -- ## -- ##

## -- ## -- ## -- ##
find.pcup1 <- function(DES, xt){
   
   insidepcup <- function(pt, DESt, xtt){
      one <- find.Ep(DES, pt)
      two <- find.phat(DES)[xtt+1, 2]
      one - two}
   
   n1 <- DES[1] ; nt <- DES[2] ; r1 <- DES[3] ; rt <- DES[4] ; n2 <- nt - n1

   out <- 0
   if( (xt!=0) & (xt<nt)) {   
   out <- uniroot(insidepcup, c(0, 1), DESt = DES, xtt = xt)$root}
   
   if( xt>=nt) { out <- 1}
   
   out}
## -- ## -- ## -- ##

## -- ## -- ## -- ##
find.pcup <- function(DES){
   
   n1 <- DES[1] ; nt <- DES[2] ; r1 <- DES[3] ; rt <- DES[4] ; n2 <- nt - n1

   xtv <- pcup <- seq(0, nt) ; L <- length(pcup)
   for(i in 1:L){pcup[i] <- find.pcup1(DES, xtv[i])}
   
   cbind(xtv, pcup)}
## -- ## -- ## -- ##

## -- ## -- ## -- ##
find.pube1 <- function(DES, xt){
   n1 <- DES[1] ; nt <- DES[2] ; r1 <- DES[3] ; rt <- DES[4] ; n2 <- nt - n1

   x1 <- seq( max(r1+1, xt-n2), min(n1, xt))
   num <- sum( choose(n1-1, x1-1) * choose(n2, xt - x1) ) 
   den <- sum( choose(n1, x1)     * choose(n2, xt - x1) )
   num/den}
## -- ## -- ## -- ##

## -- ## -- ## -- ##
find.pube <- function(DES){
   n1 <- DES[1] ; nt <- DES[2] ; r1 <- DES[3] ; rt <- DES[4] ; n2 <- nt - n1
   xtv <- pube <- seq(0, nt) ; L <- length(pube)
   for(i in 1:L){pube[i] <- find.pube1(DES, xtv[i])}
   
   cbind(xtv, pube)}
## -- ## -- ## -- ##

## -- ## -- ## -- ##
find.pchu1 <- function(DES, xt){
   phatv <- find.phat(DES)
   phat <- phatv[,2][phatv[,1] == xt]
   fphat <- find.Ep(DES, phat)

   2 * phat - fphat}
## -- ## -- ## -- ##

## -- ## -- ## -- ##
find.pchu <- function(DES){

   n1 <- DES[1] ; nt <- DES[2] ; r1 <- DES[3] ; rt <- DES[4] ; n2 <- nt - n1
   xtv <- pKoyama <- seq(0, nt) ; L <- length(pKoyama)
   for(i in 1:L){pKoyama[i] <- find.pchu1(DES, xtv[i])}
   
   cbind(xtv, pKoyama)}
## -- ## -- ## -- ##

## -- ## -- ## -- ##
find.p <- function(DES){
   phat <- find.phat(DES)
      xt <- phat[,1] ; phat <- phat[,2]
   pcup <- find.pcup(DES)[,2]
   pube <- find.pube(DES)[,2]
   pchu <- find.pchu(DES)[,2]

   data.frame(xt, phat, pcup, pube, pchu)}
## -- ## -- ## -- ##

## -- ## -- ## -- ##
find.Epall <- function(DES, p){
   phat <- find.phat(DES)[,2]
   pcup <- find.pcup(DES)[,2]
   pube <- find.pube(DES)[,2]
   pchu <- find.pchu(DES)[,2]
   
   pxx <- px.all(DES, p)[,2]
   
   e.phat <- sum( pxx * phat )
   e.pcup <- sum( pxx * pcup )
   e.pube <- sum( pxx * pube )
   e.pchu <- sum( pxx * pchu )
   
   esv <- c(e.phat, e.pcup, e.pube, e.pchu)
   biasv <- p - esv
   pbv <- biasv / p * 100

   mse.phat <- sum((phat - p)^2 * pxx)
   mse.pcup <- sum((pcup - p)^2 * pxx)
   mse.pube <- sum((pube - p)^2 * pxx)
   mse.pchu <- sum((pchu - p)^2 * pxx)

   msev <- c(mse.phat, mse.pcup, mse.pube, mse.pchu)

   var.phat <- mse.phat - biasv[1]^2
   var.pcup <- mse.pcup - biasv[2]^2
   var.pube <- mse.pube - biasv[3]^2
   var.pchu <- mse.pchu - biasv[4]^2

   varv <- c(var.phat, var.pcup, var.pube, var.pchu)

   out <- data.frame(rbind(esv, biasv, pbv, msev, varv), 
      row.names=c('Estimate', 'Bias', 'Bias%', 'MSE', 'VAR'))
   names(out)=c('mle', 'wh', 'ub', 'Koyama')
   out}
## -- ## -- ## -- ##

## -- ## -- ## -- ##   

p0 <- .2 ; p1 <- .4
n1 <- 24 ; r1 <- 5 ; nt <- 45 ; rt <- 13
x1 <- 8 ; x2 <- 5 ; xt <- x1 + x2

p01 <- c(p0, p1)
DES <- c(n1, nt, r1, rt)