#####################################################################
###    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)                                  ##
##    CR, CP: allocation ratio of a Reference(R) and a Placebo(P)  ##
##                                                                 ##
##  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 (3. Application; 3.2. Sample size )      ##
##                                                                 ##
##    > HR    <- 1 / 2                                             ##
##    > time  <- 12                                                ##
##    > Delta <- 0.33                                              ##
##    > alpha <- 0.05                                              ##
##    > power <- 0.80                                              ##
##    > CR    <- 1.00                                              ##
##    > CP    <- 0.40                                              ##
##                                                                 ##
##    > nu               <- 1                                      ##
##    > lambdaR          <- - log(0.5) /10^nu                      ##
##    > survival.rate    <- exp(-lambdaR*time^nu)                  ##
##    > survival.rateP   <- exp(-lambdaR/HR*time^nu)               ##
##    > TruePara         <- c(survival.rate, survival.rateP)       ##
##                                                                 ##
##    > SSC3NI_Surv( TruePara, CR, CP, Delta, alpha, power )       ##
##                                                                 ##
##       CR  CP        rho  nE  nR  nP   N     Power               ##
##     1  1 0.4 -0.4270365 382 382 153 917 0.8003667               ##
##                                                                 ##
#####################################################################

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

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

prop.Surv   <- ( 1 - TruePara[2] ) / ( 1 - TruePara[1] )
rho         <- -sqrt( prop.Surv * CP ) / sqrt( ( 1 + CR ) * ( CR + prop.Surv * CP ) )

# the rejection region #
mu1     <- ( 1 + CR ) * sqrt( 1 + Delta ) / ( 1 + Delta + CR)
mu2     <- ( CR * HR + CP ) / ( CR + ( 1+ Delta ) * CP ) * 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 ) * sqrt( CR * CP * ( 1 - TruePara[1] ) * ( CR + prop.Surv * CP ) / HR )

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

for( nE in 10:10000 ){
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 ){
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 )
   }