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

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

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,135 posts in 4,245 threads, 1,386 registered users;
online 6 (0 registered, 6 guests [including 1 identified bots]).
Forum time (Europe/Vienna): 23:52 CET

A Camel is a Horse designed by committee.    Anonymous

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