## One Mio sims of 2x2x2 in R in under 1 sec [General Sta­tis­tics]

Dear Helmut,

from this post I know that you are working on simulations for the 2x2x2 crossover, at least on balanced design.

Do you believe the subject line?
If not, have a look on the following code :
# simulate 2x2x2 crossover studies with log-normal data # Author: dlabes ############################################################################### # We are simulating via the way described in the MATOST paper # namely via the distribution properties of the stats that are necessary # to calculate the 90% CI. Assuming a log-normal distribution of the PK metrics # it follows that the sample mean of ln(T)-ln(R) is normal distributed with # mean=log(GMR) and sd=sqrt(2*mse/n) in case of a balanced 2x2 crossover. # The sample mse is idependently distributet via a chi-squared distribution # of df*mses/mse with df=n-2. mse and GMR are the assumed true vales. # # For unbalanced crossover designs the question is: We are using LSMeans. # Are the distribution properties the same except that 2/n has to be replaced # by 0.5*(1/n1+1/n2)? I think so but are not quite sure. # 'commercial' rounding: 0-4 down and 5-9 up cround <- function(x, digits=0) sign(x)*trunc(abs(x)*10^digits + 0.5)/10^digits power.sim2x2 <- function(CV, GMR, n, nsims=1E6, alpha=0.05, lBEL=0.8, uBEL=1.25,                          roundit=FALSE, details=FALSE) {   ptm <- proc.time()   # n is total if given as simple number   # to be correct n must then be even!    if (length(n)==1) {     nsum <- n     fact <- 2/n   }   if (length(n)==2) {     nsum <- sum(n)     fact <- 0.5*sum(1/n)   }   mse    <- log(1.0 + CV^2)   df     <- nsum-2   tval   <- qt(1-alpha,df)   # Attention! With nsims=1E8 memory of my machine (4 GB) is too low   # Thus work in chunks of 1E6 if nsims>1E6.   chunks <- 1   ns     <- nsims   if (nsims>1e6) {     chunks <- trunc(nsims/1E6)     ns     <- 1E6     if (chunks*1e6!=nsims){       nsims <- chunks*1e6       warning("nsims truncated to", nsims)     }   }   BEcount <- 0   sdm     <- sqrt(fact*mse)   mlog    <- log(GMR)   for (i in 1:chunks)   {      # simulate sample mean via its normal distribution     means  <- rnorm(ns, mean=mlog, sd=sdm)     # simulate sample mse via chi-square distribution of df*mses/mse     mses   <- mse*rchisq(ns,df)/df     hw     <- tval*sqrt(fact*mses)     lCL <- means - hw     uCL <- means + hw     # point  <- exp(means)     lCL    <- exp(lCL)     uCL    <- exp(uCL)     if (roundit){       # use our commercial rounding       lCL  <- cround(lCL,4)       uCL  <- cround(uCL,4)       # in case of rounding we use also the EMA rule for BE       BE     <- (lBEL<=lCL & uCL<=uBEL)     } else {       BE     <- (lBEL<lCL & uCL<uBEL)     }     BEcount <- BEcount + sum(BE)   }   if (details) {     cat(nsims,"sims. Time elapsed (sec):\n")     print(proc.time()-ptm)   }   BEcount/nsims }

Using these functions:
require(PowerTOST) power.TOST(CV=0.3, n=12, theta0=1.25)  0.03374229 set.seed(123) power.sim2x2(CV=0.3, n=12, GMR=1.25, details=T) 1e+06 sims. Time elapsed (sec):    user  system elapsed    0.50    0.05    0.55  0.033578 set.seed(123) power.sim2x2(CV=0.3, n=12, GMR=1.25, roundit=T)  0.033605

Regards,

Detlew

### Complete thread:

Admin contact
21,419 posts in 4,475 threads, 1,509 registered users;
online 9 (0 registered, 9 guests [including 3 identified bots]).
Forum time: Sunday 23:28 UTC (Europe/Vienna)

There is one certainty in drug development
and statistics that one can depend on:
the data are always late.    Scott Patterson and Byron Jones

The Bioequivalence and Bioavailability Forum is hosted by Ing. Helmut Schütz 