#####################################################################
### Optimal Allocation 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)                                  ##
##    CP: allocation ratio of a Placebo(P)                         ##
##    CR = 1 ( optimal allocation of a Reference(R) )              ##
##                                                                 ##
##-- 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 # (fixed value for the optimal allocation)     ##
##                                                                 ##
##    > SSC3NI_NormOA( TruePara, CR, CP, Delta, alpha, power )     ##
##                                                                 ##
##      $Optimal.Allocation.Min.N                                  ##
##          CP        rho  nE  nR  nP   N    Power                 ##
##      1 0.8 -0.4714045 151 151 121 423 0.800577                  ##
##                                                                 ##
##      $Around.OA.Min.N #( within Min.N+2 )                       ##
##          CP        rho  nE  nR  nP   N     Power                ##
##      1  0.80 -0.4714045 151 151 121 423 0.8005770               ##
##      2  0.75 -0.4629100 154 154 116 424 0.8003556               ##
##      3  0.77 -0.4663841 153 153 118 424 0.8013564               ##
##      4  0.82 -0.4746311 150 150 124 424 0.8009924               ##
##      5  0.84 -0.4777665 149 149 126 424 0.8012003               ##
##      .....                                                      ##
##                                                                 ##
##      $All.Result                                                ##
##            CP         rho   nE   nR  nP     N     Power         ##
##      1   0.01 -0.07035975 5359 5359  54 10772 0.8000072         ##
##      2   0.02 -0.09901475 2706 2706  55  5467 0.8000029         ##
##      3   0.03 -0.12067770 1822 1822  55  3699 0.8000705         ##
##      4   0.04 -0.13867505 1380 1380  56  2816 0.8001367         ##
##      5   0.05 -0.15430335 1115 1115  56  2286 0.8002719         ##
##      6   0.06 -0.16823165  938  938  57  1933 0.8002653         ##
##      7   0.07 -0.18085984  812  812  57  1681 0.8004657         ##
##      .....                                                      ##
#####################################################################

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

#--- Definition of the power function ---#
SSC3NI_NormOA <- function( TruePara, CR, CP, Delta, alpha, power ){
CP                <- seq( 0.01, 1.00, 0.01 ); lenCP <- length( CP )
z.alpha.par2      <- qnorm( 1-alpha/2, 0, 1 )

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

for( i in 1:lenCP ){
rho         <- -sqrt( CP[ i ]/( (1+CR)*(CR+CP[ i ]) ) )

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[i]*nE/( CR+CP[i] ) )

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 ){
nP      <- ceiling( nE * CP[i] )
N       <- 2 * nE + ceiling( nE * CP[i] )

Output0[ i, ]  <- c( CP[ i ], rho, nE, nE, nP, N, PowerD[1] )
 break
 }
if( nE == 10000 ){
nP      <- ceiling( nE * CP[i] )
N       <- 2 * nE + ceiling( nE * CP[i] )

Output0[ i, ]  <- c( CP[ i ], rho, nE, nE, nP, N, PowerD[1] )
   }
  }
 }

Output1.d      <- Output0[ which( Output0$N == min( Output0$N ) ), ]
if( nrow( Output1.d ) == 1 ){ 
Output1        <- Output1.d
  }else{
Output1        <- Output1.d[ which( Output0$nP == min( Output0$nP ) ), ]
  }
rownames( Output1 ) <- NULL

Output2.d      <- Output0[ which( Output0$N <= min( Output0$N )+2 ), ]
Output2        <- Output2.d[ order( Output2.d$N ), ]
rownames( Output2 ) <- NULL

Output99       <- list( Output1, Output2, Output0 )
names( Output99 ) <- c( "Optimal.Allocation.Min.N", "Around.OA.Min.N", "All.Result" )
return( Output99 )
 }