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

posted by Helmut Homepage – Vienna, Austria, 2017-08-20 18:58  – Posting: # 17727
Views: 23,514

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

Cheers,
Helmut Schütz
[image]

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

Complete thread:

Activity
 Admin contact
20,250 posts in 4,262 threads, 1,398 registered users;
online 21 (0 registered, 21 guests [including 3 identified bots]).
Forum time (Europe/Vienna): 18:59 CET

No written law has ever been more binding than
unwritten custom supported by popular opinion.    Carrie Chapman Catt

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