#####################################################################
### Optimal Allocation Sample size calculation for 3 arm NI trial ###
##  "HT design" with survival endpoints                            ##
##                                                                 ##
##-- Purpose                                                     --##
##   This R code is a sample size caluculation according to the    ##
##   equation (11) and (A3.2)* for 3 arm NI trial with binary      ##
##   endpoints proposed by Hida E and Tango T.                     ##
##   *: Design and analysis of a three-arm non-inferiority trial   ##
##    with a prespecified margin for the hazard ratio.             ##
##                                                                 ##
##-- Parameters setting                                          --##
##    TruePara (survival.rate, survival.rateP): survival rate of   ##
##          an Experimental(E) = a Reference(R), and a Placebo(P)  ##
##    HR: hazard ratio of the reference treatment to the placebo   ##
##    time: follow-up time of 3NI trial                            ##
##    Delta: non-inferiority margin based on the hazard ratio      ##
##    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) )              ##
##    accur <- 0.25  # approximation accuracy for CP >= 0.25       ##
##                                                                 ##
##  For the setting of TruePara                                    ##
##    lambdaR, lambdaP: hazard rates of the reference and placebo  ##
##                      ( HR=( lambdaR / lambdaP )^nu )            ##
##    nu: Shape parameter of the Weibull model                     ##
##    survival.rate     <- exp(-lambdaR*time^nu)                   ##
##    survival.rateP    <- exp(-lambdaP*time^nu)                   ##
##                                                                 ##
##-- 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 ( Table 3 )                              ##
##                                                                 ##
##    > HR    <- 1 / 2                                             ##
##    > time  <- 1.2                                               ##
##    > Delta <- 0.25                                              ##
##    > alpha <- 0.05                                              ##
##    > power <- 0.80                                              ##
##    > CR    <- 1.00                                              ##
##    > accur <- 0.01                                              ##
##                                                                 ##
##    > nu               <- 1                                      ##
##    > lambdaR          <- HR                                     ##
##    > survival.rate    <- exp(-lambdaR*time^nu)                  ##
##    > survival.rateP   <- exp(-lambdaR/HR*time^nu)               ##
##    > TruePara         <- c(survival.rate, survival.rateP)       ##
##                                                                 ##
##    > SSC3NI_SurvOA( TruePara, CR, CP, Delta, alpha, power )     ##
##                                                                 ##
##      $Optimal.Allocation.Min.N                                  ##
##         CP        rho  nE  nR  nP   N     Power                 ##
##      1 0.22 -0.3564705 717 717 158 1592 0.800056                ##
##                                                                 ##
##      $Around.OA.Min.N #( within Min.N+2 )                       ##
##          CP        rho  nE  nR  nP   N     Power                ##
##     1 0.22 -0.3564705 717 717 158 1592 0.8000560                ##
##     2 0.20 -0.3438769 724 724 145 1593 0.8000022                ##
##     3 0.21 -0.3503039 721 721 152 1594 0.8005181                ##
##                                                                 ##
##      $All.Result                                                ##
##           CP         rho   nE   nR   nP   N     Power           ##
##     1   0.01 -0.08732667 4593 4593  46 9232 0.8000569           ##
##     2   0.02 -0.12256740 2345 2345  47 4737 0.8000979           ##
##     3   0.03 -0.14899878 1633 1633  49 3315 0.8000768           ##
##     4   0.04 -0.17078973 1317 1317  53 2687 0.8001280           ##
##     5   0.05 -0.18957132 1146 1146  58 2350 0.8003745           ##
##     6   0.06 -0.20618830 1039 1039  63 2141 0.8004092           ##
##     7   0.07 -0.22114737  966  966  68 2000 0.8004223           ##
##      .....                                                      ##
#####################################################################

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

#--- Definition of the power function ---#
SSC3NI_SurvOA <- function( TruePara, CR, CP, Delta, alpha, power ){
CP          <- seq( accur, 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 ){
prop.Surv       <- ( 1 - TruePara[2] ) / ( 1 - TruePara[1] )
rho         <- -sqrt( prop.Surv * CP[ i ] ) / sqrt( ( 1 + CR ) * ( CR + prop.Surv * CP[ i ] ) )

for( nE in 10:10000 ){
# the rejection region #
mu1     <- ( 1 + CR ) * sqrt( 1 + Delta ) / ( 1 + Delta + CR)
mu2     <- ( CR * HR + CP[ i ] ) / ( CR + ( 1+ Delta ) * CP[ i ] ) * sqrt( ( 1 + Delta ) / HR )
cons1   <- Delta / ( 1 + Delta + CR ) * sqrt( CR * ( 1 + CR ) * ( 1 - TruePara[1]) ) 
cons2   <- ( 1 - ( 1 + Delta ) * HR ) / ( CR + ( 1+ Delta ) * CP[ i ] ) * sqrt( CR * CP[ i ] * ( 1 - TruePara[1] ) * ( CR + prop.Surv * CP[ i ] ) / HR )

l1      <- -mu1 * z.alpha.par2 + cons1 * sqrt( nE )
l2      <- -mu2 * z.alpha.par2 + cons2 * sqrt( nE )

mu      <- -c( 0, 0 )
sigma   <- matrix( c( 1, rep( rho, 2 ), 1 ), ncol = 2 )
PowerD  <- -1 + pnorm( l1, mean = 0, sd = 1) + pnorm( l2, mean = 0, sd = 1) + 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 )
 }