Bioequivalence and Bioavailability Forum

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

posted by d_labes  – Berlin, Germany, 2013-01-09 16:23 (4117 d 08:19 ago) – Posting: # 9803
Views: 7,588

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)
[1] 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
[1] 0.033578
set.seed(123)
power.sim2x2(CV=0.3, n=12, GMR=1.25, roundit=T)
[1] 0.033605

Complete thread:

• Imbalanced Xovers: α-inflation? Helmut 2012-12-27 17:21 [General Sta­tis­tics] 
  • Imbalanced Xovers: α-inflation? ElMaestro 2012-12-28 09:30
    • Imbalanced Xovers: α-inflation? Helmut 2012-12-28 14:37
      • please explain ElMaestro 2012-12-29 08:44
        • I’ll try Helmut 2012-12-29 15:42
  • Imbalanced Xovers: α-inflation? d_labes 2012-12-31 13:50
  • One Mio sims of 2x2x2 in R in under 1 secd_labes 2013-01-09 15:23
    • Wow! Helmut 2013-01-09 16:10
      • Unbalanced d_labes 2013-01-09 16:56
        • α-inflation ≈ urban myth Helmut 2013-01-09 18:07
          • Beer d_labes 2013-01-10 11:07
            • Beer Helmut 2013-01-10 13:45