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

posted by d_labes  – Berlin, Germany, 2013-01-09 15:23  – Posting: # 9803
Views: 6,694

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 :cool::
# 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

Regards,

Detlew

Complete thread:

Activity
 Admin contact
20,797 posts in 4,354 threads, 1,445 registered users;
online 12 (0 registered, 12 guests [including 3 identified bots]).
Forum time: 13:52 CEST (Europe/Vienna)

The function of genius is not to give new answers, but to pose
new questions which time and mediocrity can resolve.    H.R. Trevor-Roper

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