Bioequivalence and Bioavailability Forum • The ultimate crackpot iteration!

The ultimate crackpot iteration! [Two-Stage / GS Designs]

posted by Helmut – Vienna, Austria, 2017-08-20 20:58 (2434 d 10:16 ago) – Posting: # 17727
Views: 27,455

Hi ElMaestro,

❝ You mean if you do the array ini as described here once before each run, and then use GetNps (or equivalent) then that will be 500x slower than using Zhang + iteration??

❝ That would surprise me a lot. Please elaborate.

Had to dive deep into the guts of PowerTOST’s internal functions. Hope I got it right. Detlew?

library(PowerTOST)

library(microbenchmark)

sampleN.Stg2 <- function(alpha, CV, theta1=0.8, theta2=1.25,

                         theta0=0.95, targetpower=0.8,

                         method="shifted", imax=100) {

  d.no  <- PowerTOST:::.design.no("2x2")

  ades  <- PowerTOST:::.design.props(d.no)

  dfe   <- PowerTOST:::.design.df(ades, robust=FALSE)

  steps <- ades$steps/2

  bk    <- ades$bk

  df    <- n <- 0

  while (df < 1) {

    n  <- n + 1

    df <- eval(dfe) - 1

  }

  nmin <- as.integer(steps*trunc(n/steps))

  nmin <- nmin + steps*(nmin < n)

  ltheta1 <- log(theta1)

  ltheta2 <- log(theta2)

  diffm   <- log(theta0)

  se      <- CV2se(CV)

  n       <- PowerTOST:::.sampleN0_3(alpha, targetpower, ltheta1,

                                     ltheta2, diffm, se, steps, bk)

  if (n < nmin) n <- nmin

  df  <- eval(dfe) - 1

  pow <- PowerTOST:::.calc.power(alpha, ltheta1, ltheta2, diffm,

                                 sem=se*sqrt(bk/n), df, method)

  iter <- 0

  down <- FALSE; up <- FALSE

  while (pow > targetpower) {

    if (n <= nmin) break

    down <- TRUE

    n    <- n-steps

    iter <- iter + 1

    df   <- eval(dfe) - 1

    pow  <- PowerTOST:::.calc.power(alpha, ltheta1, ltheta2, diffm,

                                    sem=se*sqrt(bk/n), df, method)

    if (iter > imax) break

  }

  while (pow<targetpower) {

    up   <- TRUE; down <- FALSE

    n    <- n+steps

    iter <- iter + 1

    df   <- eval(dfe)-1

    pow <- PowerTOST:::.calc.power(alpha, ltheta1, ltheta2, diffm,

                                   sem=se*sqrt(bk/n), df, method)

    if (iter > imax) break

  }

  return(n)

}



EM <- function(CV, alpha1=0.05, alpha2=0.0294, theta1=0.8,

               theta2=1.25, GMR=0.95, targetpower=0.8) {

  vecQT1 <- qt(p=1-alpha1, df=1:2000)

  vecQT2 <- qt(p=1-alpha2, df=1:2000)

  Power.calc <- function(Nps, GMR, CV, Is.St2=FALSE) {

    s    <- sqrt(log(CV*CV+1))

    nom1 <- log(theta2/GMR)

    nom2 <- log(theta1/GMR)

    Nps2 <- 2*Nps

    den  <- s*sqrt(2/Nps2)

    df   <- Nps2-2

    if (!Is.St2) {

      q <- vecQT1[df]


    } else {

      df <- df-1

      q  <- vecQT2[df]


    }

    pw <- pt((nom1/den)-q, df) - pt((nom2/den)+q, df)

    if (pw < 0) pw <- 0.01

    return(pw)

  }

  FindCritCV <- function(Pwr.T, GMR, Nps, Is.St2=TRUE, toler=0.000001) {

    CV1 <- 0.0123

    CV2 <- 8.7654

    while (abs(Power.calc(Nps, GMR, CV1, Is.St2)) - 

           abs(Power.calc(Nps, GMR, CV2, Is.St2)) > toler) {

      CV3 <- 0.5*(CV1+CV2)

      tmp <- Pwr.T-Power.calc(Nps, GMR, CV3, Is.St2)

      ifelse (tmp > 0, CV2 <- CV3, CV1 <- CV3)

    }

    return(0.5*(CV2+CV1))

  }

  vecCritCVs <- numeric()

  for (Nps in 2:1000)

    vecCritCVs[Nps] <- FindCritCV(targetpower, GMR, Nps, Is.St2=TRUE, 1e-12)

  GetNps <- function(CV) {

    Nps <- 2

    while (vecCritCVs[Nps] <= CV) Nps <- Nps+1

    return(Nps)

  }

  N <- GetNps(CV)

  return(N)

}

Let’s see:

CV <- seq(0.1, 0.8, 0.05)

N2 <- data.frame(GMR=0.95, target=0.8, CV, EM=NA, PT=NA)

for (j in seq_along(CV)) {

  N2[j, "EM"] <- EM(CV=CV[j])

  N2[j, "PT"] <- ceiling(sampleN.Stg2(alpha=0.0294, CV=CV[j])/2)

}

print(N2, row.names=FALSE)



  GMR target   CV  EM  PT

 0.95    0.8 0.10   4   4

 0.95    0.8 0.15   7   7

 0.95    0.8 0.20  12  12

 0.95    0.8 0.25  17  17

 0.95    0.8 0.30  24  24

 0.95    0.8 0.35  31  31

 0.95    0.8 0.40  40  40

 0.95    0.8 0.45  49  49

 0.95    0.8 0.50  59  59

 0.95    0.8 0.55  69  69

 0.95    0.8 0.60  80  80

 0.95    0.8 0.65  92  92

 0.95    0.8 0.70 104 104

 0.95    0.8 0.75 116 116

 0.95    0.8 0.80 128 128

So far, so good. Speed?

# wrappers for nice output

EM.f <- function() EM(CV=0.165)

PT.f <- function() ceiling(sampleN.Stg2(alpha=0.0294, CV=0.165)/2)

Check first:

EM.f()

[1] 9

PT.f()

[1] 9

Looking good. Now:

res <- microbenchmark(EM.f(), PT.f(), times=50L,

                      control=list("random", warmup=10))

options(microbenchmark.unit="relative")

print(res)



Unit: relative

   expr     min       lq     mean   median       uq      max neval cld

 EM.f() 1936.08 1765.027 1430.441 1304.101 1245.526 1141.647    50   b

 PT.f()    1.00    1.000    1.000    1.000    1.000    1.000    50  a

—
Dif-tor heh smusma 🖖🏼 Довге життя Україна!
Helmut Schütz

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes

Complete thread:

Initial sample size guess for the Potvin methods ElMaestro 2017-08-19 15:04 [Two-Stage / GS Designs]
- Initial sample size guess for the Potvin methods Helmut 2017-08-19 16:06
  - Initial sample size guess for the Potvin methods ElMaestro 2017-08-19 16:14
    - Initial sample size guess for the Potvin methods Helmut 2017-08-19 17:12
      - Initial sample size guess for the Potvin methods ElMaestro 2017-08-19 17:17
        
        Confuse-a-Cat Helmut 2017-08-19 17:33
        
        Confuse-a-Cat ElMaestro 2017-08-19 17:56
    - loop ↔ vectorized ↔ direct Helmut 2017-08-20 14:40
      - loop ↔ vectorized ↔ direct ElMaestro 2017-08-20 15:22
        
        loop ↔ vectorized ↔ direct Helmut 2017-08-20 16:23
        
        loop ↔ vectorized ↔ direct ElMaestro 2017-08-20 17:22
- The n ext crackpot iteration ElMaestro 2017-08-19 20:04
  - The n ext crackpot iteration Helmut 2017-08-20 02:20
- The ultimate crackpot iteration! ElMaestro 2017-08-20 14:35
  - The ultimate crackpot iteration! Helmut 2017-08-20 15:11
    - The ultimate crackpot iteration! ElMaestro 2017-08-20 15:28
      - The ultimate crackpot iteration! Helmut 2017-08-20 16:06
        
        The ultimate crackpot iteration! ElMaestro 2017-08-20 16:15
        
        The ultimate crackpot iteration!Helmut 2017-08-20 18:58
        
        The ultimate crackpot iteration! ElMaestro 2017-08-20 19:32
        
        Suggested code ElMaestro 2017-08-21 18:13
        
        Nitpicker! Helmut 2017-08-22 13:33
        
        Nitpicker! ElMaestro 2017-08-22 17:27
        
        Nitpicker! Helmut 2017-08-22 17:49
        
        Nitpicker! ElMaestro 2017-08-22 17:59
        
        Nitpicker! Helmut 2017-08-22 19:15
        
        Benchmark code ElMaestro 2017-08-22 22:29
        
        Benchmark code Helmut 2017-08-23 01:48