#####################################################################
###    Sample size calculation for 3 arm NI trial "HT design"     ###
##     with continuous endpoints                                   ##
##                                                                 ##
##-- Purpose                                                     --##
##   This R code is a sample size caluculation according to the    ##
##   equation (15) and (A1)* for 3 arm NI trial with continuous    ##
##   endpoints proposed by Hida E and Tango T.                     ##
##   *: On the three‐arm non‐inferiority trial including a placebo ##
##      with a prespecified margin. Statistics in Medicine, 30(3). ##
##      224-231. 2011.                                             ##
##                                                                 ##
##-- Parameters setting                                          --##
##    TruePara (muE, muR, muP, sigma): mean of an Experimental(E), ##
##          a Reference(R), a Placebo(P) treatment and c common SE ##
##    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. Example; 3.2. Sample size )          ##
##                                                                 ##
## > TruePara  <- c( 10, 10, 5, 6.5 ) # c( muE, muR, muP, sigma )  ##
## > Delta     <- 2.5                                              ##
## > alpha     <- 0.05                                             ##
## > power     <- 0.80                                             ##
## > CR        <- 1.00                                             ##
## > CP        <- 0.80                                             ##
##                                                                 ##
##    > SSC3NI_Norm( TruePara, CR, CP, Delta, alpha, power )       ##
##                                                                 ##
##       CR  CP        rho  nE  nR  nP   N    Power                ##
##     1  1 0.8 -0.4714045 151 151 121 423 0.800577                ##
##                                                                 ##
#####################################################################

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

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

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

for( nE in 10:10000 ){
l1  <- z.alpha.par2-( TruePara[1]-TruePara[2]+Delta )/TruePara[4] * sqrt( CR*nE/( 1+CR ) )
l2  <- z.alpha.par2-( TruePara[2]-TruePara[3]-Delta )/TruePara[4] * sqrt( CR*CP*nE/( CR+CP ) )

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 )
   }