TwoStageInf <- function(p0, p1, n1, R1, nT, RT, x1, m2=nT-n1, x2, conf.level=0.9, two.sided=TRUE, ...){
                        ## Last update: 2021/8/5
## ----------------------------------------------- ##
## Proper Inference from Simon's Two-Stage Designs ##
## ----------------------------------------------- ##
## Please cite
## Koyama T and Chen H. ``Proper inference from Simon's two-stage designs''
## Statistics in Medicine 27(16), 2008. PMID: 17960777. PMCID: PMC6047527.
## --------------------- ##

################
## User input ##
################
# p0 ... P[Success under H0]
# p1 ... P[Success under H1]
# n1 ... Stage 1 sample size
# R1 ... Stage 1 critical value
# nT ... Final planned sample size
# RT ... Final planned critical value
# Convention is to use r1 (= R1 - 1) and rT (= RT - 1)

# x1 ... Stage 1 data
# m2 ... New stage 2 sample size (if changed)
# x2 ... Stage 2 data

#######################################
## Design and design characteristics ##
#######################################
TwoStageInf.prob <- function(p0, p1, n1, R1, nT, RT){
    ## Continue if X1 >= R1.
    ## Reject if Xt >= RT.
    ## Convention is to use r1 = R1-1 and rT = RT-1
    n2 <- nT-n1
    x1v <- 0:n1
    px0 <- dbinom(x1v, n1, p0)
    px1 <- dbinom(x1v, n1, p1)

    cont <- x1v >= R1
    R2 <- pmax(RT - x1v, 0)
    cp0 <- replace(1 - pbinom(R2 - 1, n2, p0), x1v < R1, 0)
    cp1 <- replace(1 - pbinom(R2 - 1, n2, p1), x1v < R1, 0)

    PROB <- data.frame(x1v, px0, px1, cp0, cp1)
        ## Stage 1 pmf and conditional power
    POWE <- c( sum(px0 * cp0), sum(px1 * cp1) )
        ## Unconditional type I erro rate and power
    PET  <- c( pbinom(R1-1, n1, p0), pbinom(R1-1, n1, p1) )
        ## Probability of early termination 
    EN   <- n1 + (nT-n1) * (1-PET)
        ## Expected sample size

    DESIGN <- data.frame(p0=p0, p1=p1, n1=n1, R1=R1, nT=nT, RT=RT)
    CHARA <- data.frame(POWE, PET, EN)
        row.names(CHARA) <- c('NULL','ALT')
    list(DESIGN, CHARA)
}

#############
## p-value ##
#############
pval.st2 <- function(p, n1, R1, nT, RT, m2, x1, x2){
    ## p is p0 (because this is for p.value) unless computing confidence interval.
    ## m2 is the actual sample size for stage 2. If there is no change, use m2 = nT-n1.
    n2 <- nT - n1
    ## Continue if X1 >= R1.
    ## Convention is to use r1 = R1-1 and rT = RT-1
    x1c <- R1:n1
    x2c <- RT - x1c - 1
    R2 <- RT - x1
        d1m <- dbinom(x1c, n1, p)
    ## pi.star is such that A(x1,n2,pi.star) = conditional p.value
    cond.pval <- pbinom(x2-1, m2, p, lower.tail=FALSE)
        ## pi.star <- qbeta(cond.pval, R2, m2-R2+1) ## Gone back to original 6/16/2020
        ## pi.star <- qbeta(cond.pval, R2, n2-R2+1) ## Fixed 12/29/2019
        pi.star <- qbeta(cond.pval, max(R2,0), n2-R2+1) 
        ## p2m <- pbinom(x2c, m2, pi.star, lower.tail=FALSE)  ## Gone back to original 6/16/2020
        ## p2m <- pbinom(x2c, n2, pi.star, lower.tail=FALSE) ## Fixed 12/29/2019
        p2m <- pbinom(x2c, n2, pi.star, lower.tail=FALSE) ## Fixed 12/29/2019
    sum(d1m*p2m)
}

#########################
## Confidence interval ##
#########################
conf.int.st2 <- function(n1, R1, nT, RT, m2, x1, x2, conf.level=0.9, two.sided=TRUE){
    alp <- (1 - conf.level) / (1 + as.numeric(two.sided))
        f1 <- function(x,a,n2,R1,nT,RT,m2,x1,x2) pval.st2(x, n1, R1, nT, RT, m2, x1, x2) - a

    lower.bound <- uniroot(f1, interval=c(0,1), a=alp, n2=n2,R1=R1,nT=nT,RT=RT,m2=m2,x1=x1,x2=x2)$root
    median.estimate <- uniroot(f1, interval=c(0,1), a=0.5, n2=n2,R1=R1,nT=nT,RT=RT,m2=m2,x1=x1,x2=x2)$root
    upper.bound <- 1
    if (two.sided){
    upper.bound <- uniroot(f1, interval=c(0,1), a=1-alp, n2=n2,R1=R1,nT=nT,RT=RT,m2=m2,x1=x1,x2=x2)$root
    }

    out <- list()
    out$Conf.int <- c(lower.bound, upper.bound)
    out$median.est <- median.estimate
    out
}

##########################
## Putting all together ##
##########################
    result <- list()
        param <- TwoStageInf.prob(p0, p1, n1, R1, nT, RT)
    result$parameter <- param[[1]]
    result$character <- param[[2]]
    result$data <- data.frame(x1=x1, n2=nT-n1, m2=m2, x2=x2)
    if (x1 < R1)
    {
        result$phat <- x1/n1
        result$p.value <- 1-pbinom(x1-1, n1,p)
    }
    else
    {
        result$p.value <- pval.st2(p0,n1,R1,nT,RT,m2,x1,x2)
        conf <- conf.int.st2(n1,R1,nT,RT,m2,x1,x2, ...)
        result$conf <- conf[[1]]
        result$med.esti <- conf[[2]]
    }
    result
}

## --- ##
## END ##
## --- ##

## EXAMPLE in Koyama and Chen (2008)
# In Section 3.4 
# OPTIMAL design for p0=0.1, p1=0.3, alpha=0.05, beta=0.20
    # p0 <- 0.1 ; p1 <- 0.3
    # n1 <- 10 ; R1 <- 2 ; nT <- 29 ; RT <- 6
    # x1 <- 2 ; x2 <- 4
    # m2 <- nT-n1 ## No sample size change
TwoStageInf(p0=0.1, p1=0.3, n1=10, R1=2, nT=29, RT=6, x1=2, x2=4)

# Section 4.1 and Section 4.2
# MINIMAX design for p0=0.3, p1=0.5, alpha=0.05, beta=0.20
    # p0 <- 0.3 ; p1 <- 0.5
    # n1 <- 19 ; R1 <- 7 ; nT <- 39 ; RT <- 17
    # x1 <- 7 ; x2 <- 10
    # sample size change for stage 2
    # m2 <- 23
# With sample size change n2=20, m2=23
TwoStageInf(p0=0.3, p1=0.5, n1=19, R1=7, nT=39, RT=17, x1=7, m2=23, x2=10)

## Another example 
# p0 <- 0.05 ; p1 <- 0.25
# n1 <- 13 ; R1 <- 1 ; nT <- 20 ; RT <- 3
# x1 <- 5
# m2 <- nT - n1 ## No sample size change
# x2 <- 2
TwoStageInf(p0=0.05, p1=0.25, n1=13, R1=1, nT=20, RT=3, x1=5, m2=7, x2=2)