#####################################################################
###    Sample size calculation for 3 arm NI trial "HT design"     ###
##     with binary endpoints                                       ##
##                                                                 ##
##-- Purpose                                                     --##
##   This R code is a sample size caluculation according to the    ##
##   equation (13) and (15)* for 3 arm NI trial with binary        ##
##   endpoints proposed by Hida E and Tango T.                     ##
##   *: Three-Arm Noninferiority Trials with a Prespecified        ##
##   Margin for Inference of the Difference in the Proportions of  ##
##   Binary Endpoints. Journal of Biopharmaceutical Statistics,    ##
##   23(4). 774-789. 2013.                                         ##
##                                                                 ##
##-- Parameters setting                                          --##
##    TruePara (piE, piR, piP): response probabilities of          ##
##      an Experimental(E), a Reference(R), a Placebo(P) treatment ##
##    Delta: non-inferiority margin                                ##
##    alpha: two-tails alpha error                                 ##
##    power: power (1-beta error)                                  ##
##    CR, CP: allocation ratio of a Reference(R) and a Placebo(P)  ##
##                                                                 ##
##-- Output                                                      --##
##    rho: correlation coefficient of between two hypotheses       ##
##    nE, nR, nP: sample size of an Experimental, a Reference,     ##
##                a Placebo treatment                              ##
##    N: total sample size                                         ##
##    Power: actual power                                          ##
##                                                                 ##
##-----------------------------------------------------------------##
##  Examples of execution (3. Application; 3.2. Sample size )      ##
##                                                                 ##
## > TruePara  <- c( 0.35, 0.35, 0.15 ) # c( piE, piR, piP )       ##
## > Delta     <- 0.1                                              ##
## > alpha     <- 0.05                                             ##
## > power     <- 0.80                                             ##
## > CR        <- 1.00                                             ##
## > CP        <- 0.78                                             ##
##                                                                 ##
##    > SSC3NI_Binom( TruePara, CR, CP, Delta, alpha, power )      ##
##                                                                 ##
##       CR   CP        rho  nE  nR  nP    N     Power             ##
##     1  1 0.78 -0.5393972 460 460 359 1279 0.8001454             ##
##                                                                 ##
#####################################################################

#-- Load library --#
library( mvtnorm )

#--- Definition of the power function ---#
SSC3NI_Binom <- function( TruePara, CR, CP, Delta, alpha, power ){
z.alpha.par2      <- qnorm( 1-alpha/2, 0, 1 )

pq1     <- TruePara[1] * ( 1 - TruePara[1] )
pq2     <- TruePara[2] * ( 1 - TruePara[2] )
pq3     <- TruePara[3] * ( 1 - TruePara[3] )

rho     <- -pq2 * sqrt( CP/( ( pq1 * CR + pq2 ) * ( pq2 * CP + pq3 * CR) ) )

# restricted maximum likelihood estimator of piR #
a1      <- 1 + CR + CP
b1      <- -( 1 + CR + CP + TruePara[1] + CR * TruePara[2] + CP * TruePara[3] + Delta * ( 1 + 2 * CR + CP ) )
c1      <- CR * Delta^2 + ( 1 + CR + CP + 2 * CR * TruePara[2] ) * Delta + TruePara[1] + CR * TruePara[2] + CP * TruePara[3] 
d1      <- - CR * TruePara[2] * Delta * ( 1 + Delta )
v1      <- b1^3 / ( 27 * a1^3 ) - ( b1 * c1 )/( 6 * a1^2 ) + d1 / ( 2 * a1 )
u1      <- sign( v1 ) * sqrt( b1^2 / ( 9 * a1^2 ) - c1 / (3 * a1 ))
w1      <-( pi + acos( v1 / u1^3 )) / 3

rmleR       <- 2 * u1 * cos(w1) - b1 / ( 3 * a1 )

# the rejection region #
tau1       <- sqrt( ( ( rmleR - Delta) * ( 1 - rmleR + Delta ) * CR + rmleR * ( 1 - rmleR )) / (pq1 * CR + pq2 ) )
cons1      <- ( TruePara[1] - ( TruePara[2] - Delta )) * sqrt( CR /( pq1 * CR + pq2 ) )
tau2       <- sqrt( ( ( rmleR - Delta) * ( 1 - rmleR + Delta ) * CR + rmleR * ( 1 - rmleR ) * CP) / (pq2 * CP + pq3 * CR ) )
cons2      <- ( TruePara[2] - ( TruePara[3] + Delta )) * sqrt( CR * CP / ( pq2 * CP + pq3 * CR ) )

Output0   <- data.frame( 
CR, 
CP, 
rho, 
nE    = 0, 
nR    = 0, 
nP    = 0,
N     = 0,
Power = 0
 )

for( nE in 10:10000 ){
l1         <- tau1 * z.alpha.par2 - cons1 *  sqrt( nE )
l2         <- tau2 * z.alpha.par2 - cons2 *  sqrt( nE )

mu      <- -c( 0, 0 )
sigma   <- matrix( c( 1, rep( rho, 2 ), 1 ), ncol = 2 )
PowerD  <- pmvnorm( lower = c( l1, l2 ), upper = Inf, mean = mu, sigma = sigma )

if( PowerD >= power ){
nR      <- ceiling( nE*CR )
nP      <- ceiling( nE*CP )
N       <- nE + ceiling( nE*CR ) + ceiling( nE*CP )

Output0  <- data.frame( CR, CP, rho, nE, nR, nP, N, PowerD[1] )
 break
 }
if( nE == 10000 ){
nR      <- ceiling( nE * CR )
nP      <- ceiling( nE * CP )
N       <- nE + ceiling( nE*CR ) + ceiling( nE*CP )

Output0  <- data.frame( CR, CP, rho, nE, nR, nP, N, PowerD[1] )
 }
  }

names(Output0)  <- c( "CR", "CP", "rho", "nE", "nR", "nP", "N", "Power" )
return( Output0 )
   }