##################################################################### ### Optimal Allocation 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) ## ## 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. 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 ## ## ## ## > SSC3NI_BinomOA( TruePara, CR, CP, Delta, alpha, power ) ## ## ## ## $Optimal.Allocation.Min.N ## ## CP rho nE nR nP N Power ## ## 1 0.78 -0.5393972 460 460 359 1279 0.8001454 ## ## ## ## $Around.OA.Min.N #( within Min.N+2 ) ## ## CP rho nE nR nP N Power ## ## 1 0.78 -0.5393972 460 460 359 1279 0.8001454 ## ## 2 0.81 -0.5436227 455 455 369 1279 0.8000386 ## ## 3 0.77 -0.5379387 462 462 356 1280 0.8004877 ## ## 4 0.80 -0.5422385 457 457 366 1280 0.8005461 ## ## 5 0.83 -0.5463211 452 452 376 1280 0.8001351 ## ## ..... ## ## ## ## $All.Result ## ## CP rho nE nR nP N Power ## ## 1 0.01 -0.09362249 10000 10000 100 20100 0.6036026 ## ## 2 0.02 -0.13125670 7389 7389 148 14926 0.8000133 ## ## 3 0.03 -0.15938884 4983 4983 150 10116 0.8000037 ## ## 4 0.04 -0.18250735 3781 3781 152 7714 0.8001220 ## ## 5 0.05 -0.20237118 3059 3059 153 6271 0.8001198 ## ## 6 0.06 -0.21989274 2577 2577 155 5309 0.8000019 ## ## 7 0.07 -0.23562008 2234 2234 157 4625 0.8001570 ## ## ..... ## ##################################################################### #-- Load library --# library( mvtnorm ) #--- Definition of the power function ---# SSC3NI_BinomOA <- 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 ){ pq1 <- TruePara[1] * ( 1 - TruePara[1] ) pq2 <- TruePara[2] * ( 1 - TruePara[2] ) pq3 <- TruePara[3] * ( 1 - TruePara[3] ) rho <- -pq2 * sqrt( CP[ i ]/( ( pq1 * CR + pq2 ) * ( pq2 * CP[ i ] + pq3 * CR) ) ) for( nE in 10:10000 ){ # restricted maximum likelihood estimator of piR # a1 <- 1 + CR + CP[ i ] b1 <- -( 1 + CR + CP[ i ] + TruePara[1] + CR * TruePara[2] + CP[ i ] * TruePara[3] + Delta * ( 1 + 2 * CR + CP[ i ] ) ) c1 <- CR * Delta^2 + ( 1 + CR + CP[ i ] + 2 * CR * TruePara[2] ) * Delta + TruePara[1] + CR * TruePara[2] + CP[ i ] * 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 ) / sqrt( pq1 * CR + pq2 ) tau2 <- sqrt( ( ( rmleR - Delta) * ( 1 - rmleR + Delta ) * CR + rmleR * ( 1 - rmleR ) * CP[ i ]) / (pq2 * CP[ i ] + pq3 * CR ) ) cons2 <- ( TruePara[2] - ( TruePara[3] + Delta )) * sqrt( CR * CP[ i ] ) / sqrt( pq2 * CP[ i ] + pq3 * CR ) 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 ){ 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 ) }