Power2Stage / power.2stage() [Power / Sample Size]

posted by Helmut Homepage – Vienna, Austria, 2018-03-24 15:49 (1309 d 12:41 ago) – Posting: # 18593
Views: 6,246

Hi Pjs,

extending Detlew’s post. Try this code:

library(Power2Stage)
CV      <- 0.4
n1      <- 12
method  <- "C"
pmethod <- c("shifted", "nct", "exact")
st      <- proc.time()[[3]]
res     <- data.frame(method=method, CV=CV, n1=n1, pmethod=pmethod,
                      power=NA, TIE=NA, n.tot=NA, pct05=NA, pct50=NA,
                      pct95=NA, pct.stg2=NA)
for (j in seq_along(pmethod)) {
  x <- power.2stage(alpha=rep(0.0294, 2), method=method, CV=CV, n1=n1,
                    GMR=0.95, theta0=0.95, pmethod=pmethod[j],
                    npct=c(0.05, 0.5, 0.95))
  y <- power.2stage(alpha=rep(0.0294, 2), method=method, CV=CV, n1=n1,
                    GMR=0.95, theta0=1.25, pmethod=pmethod[j])
  res[j, "power"]    <- round(x$pBE, 4)
  res[j, "TIE"]      <- round(y$pBE, 4)
  res[j, "n.tot"]    <- round(x$nmean, 1)
  res[j, 8:10]       <- x$nperc
  res[j, "pct.stg2"] <- round(x$pct_s2, 1)
}
et      <- proc.time()[[3]] - st
cat("run time:", signif(et/60, 1), "minutes\n"); print(res, row.names=FALSE)

Gives:

run time: 8.9 minutes
 method  CV n1 pmethod  power    TIE n.tot pct05 pct50 pct95 pct.stg2
      C 0.4 12 shifted 0.7509 0.0342  79.0    34    74   142     99.0
      C 0.4 12     nct 0.7506 0.0345  78.9    32    74   142     98.9
      C 0.4 12   exact 0.7506 0.0345  78.9    32    74   142     98.9

Columns: power overall power, TIE empiric Type I Error (for true ratio 1.25), n.tot expected average total sample size; pct05, pct50, pct95 5%, 50%, 95% percentiles of expected total sample size, pct.stg2 percent of studies expected to proceed to the 2nd stage.
Reported by Potvin et al. in Tables I & II:

method  CV n1 pmethod  power    TIE n.tot pct05 pct50 pct95 pct.stg2
      C 0.4 12 shifted
0.7505 0.0346  79.0    32    74   142     99.0

Potvin et al. employed FORTRAN-code (hence, we don’t know which pseudorandom number generator was used). In later papers (Montague et al., Xu et al.) R was used as well where the Mersenne-Twister is the default PRNG. However, we don’t know which seed was used. Therefore, slight differences in the simulations’ results are normal. BTW, Fuglsang employed also the Mersenne-Twister (C-code) in his papers.
Try this code with a different seed in each run:

runs <- 25
res  <- data.frame(power=rep(NA, runs), TIE=rep(NA, runs))
st   <- proc.time()[[3]]
for (j in 1:runs) {
  res[j, "power"]  <- power.2stage(alpha=rep(0.0294, 2), method=method,
                                   CV=CV, n1=n1, GMR=0.95, theta0=0.95,
                                   pmethod="shifted", setseed=FALSE)$pBE
  res[j, "TIE"]    <- power.2stage(alpha=rep(0.0294, 2), method=method,
                                   CV=CV, n1=n1, GMR=0.95, theta0=1.25,
                                   pmethod="shifted", setseed=FALSE)$pBE
}
et   <- proc.time()[[3]] - st
cat("run time:", signif(et/60, 2), "minutes\n"); summary(res)

Gives:

run time: 6.4 minutes
     power             TIE         
 Min.   :0.7482   Min.   :0.03403 
 1st Qu.:0.7503   1st Qu.:0.03443 
 Median :0.7511   Median :0.03455 
 Mean   :0.7511   Mean   :0.03451 
 3rd Qu.:0.7522   3rd Qu.:0.03460 
 Max.   :0.7537   Max.   :0.03485


If you want to estimate the sample size of the 2nd stage I recommend the function sampleN2.TOST() of Power2Stage instead of sampleN.TOST() of PowerTOST and subtract n1. In the final analysis we have one degree of freedom less than in the conventional model (due to the additional stage-term). Comparison:

x <- sampleN2.TOST(alpha=0.0294, CV=CV, n1=n1, theta0=0.95)
y <- sampleN.TOST(alpha=0.0294, CV=CV, theta0=0.95, details=FALSE, print=FALSE)
print(x, row.names=FALSE) # n2 directly with correct df
print(y, row.names=FALSE) # total n with incorrect df

Gives:

 Design  alpha  CV theta0 theta1 theta2 n1 Sample size Achieved power Target power
    2x2 0.0294 0.4   0.95    0.8   1.25 12          68      0.8107469          0.8
 Design  alpha  CV theta0 theta1 theta2 Sample size Achieved power Target power
    2x2 0.0294 0.4   0.95    0.8   1.25          80      0.8105103          0.8


Let’s see whether we can reproduce Potvin’s examples for ‘Method C’:

potvin.examples <- function(alpha0=0.05, alpha=0.0294, n1,
                            log.GMR1, S.sq1, log.GMR2, S.sq2) {
  require(PowerTOST)
  require(Power2Stage)
  method   <- "C"
  ltheta1  <- log(0.80)
  ltheta2  <- log(1.25)
  ltheta0  <- log(0.95)
  pass.1   <- pass.2 <- FALSE
  pwr.1    <- power.TOST(alpha=alpha0, CV=mse2CV(S.sq1), n=n1,
                         method="shifted")
  if (pwr.1 >= 0.80) {
    CI.stg1 <- round(100*CI.BE(alpha=alpha0, pe=exp(log.GMR1),
                               CV=mse2CV(S.sq1), n=n1), 2)
  } else {
    CI.stg1 <- round(100*CI.BE(alpha=alpha, pe=exp(log.GMR1),
                               CV=mse2CV(S.sq1), n=n1), 2)
  }
  if (CI.stg1[["lower"]] < 0.80 | CI.stg1[["upper"]] > 1.25 & pwr.1 < 0.80) {
    n2 <- sampleN2.TOST(alpha=alpha, CV=mse2CV(S.sq1), n1=n1,
                        theta0=0.95, method="shifted")[["Sample size"]]
  } else {
    pass.1 <- TRUE
  }
  if (!pass.1) {
    df2     <- n1 + n2 - 3
    sem     <- sqrt(2*S.sq2/(n1+n2))
    CI.stg2 <- round(100*exp(log.GMR2 + c(-1, +1) *
                             qt(1-alpha, df2)*sem), 2) # balanced!
    names(CI.stg2) <- c("lower", "upper")
    if (CI.stg2[["lower"]] >= 0.80 | CI.stg2[["upper"]] <= 1.25) pass.2 <- TRUE
  }
  txt <- paste0("Power in the interim = ", round(100*pwr.1, 1), "%.",
                "\nSince power is ")
  if (pwr.1 >= 0.80) {
    txt <- paste0(txt, "\u226580%, BE assessed with alpha ", alpha0, ".")
  } else {
    txt <- paste0(txt, "<80%, BE assessed with alpha ", alpha, ".")
  }
  if (pass.1) {
    txt <- paste0(txt, "\nStudy passed in the 1st stage: ")
    if (pwr.1 >= 0.80) {
      txt <- paste0(txt, "90% CI: ")
    } else {
      txt <- paste0(txt, 100*(1-2*alpha), "% CI = ")
    }
    txt <- paste0(txt, CI.stg1[["lower"]], "% \u2013 ", CI.stg1[["upper"]], "%.")
  } else {
    txt <- paste0(txt, "\nStudy failed in the 1st stage: ")
    if (pwr.1 >= 0.80) {
      txt <- paste0(txt, "90% CI = ")
    } else {
      txt <- paste0(txt, 100*(1-2*alpha), "% CI = ")
    }
    txt <- paste0(txt, CI.stg1[["lower"]], "% \u2013 ", CI.stg1[["upper"]], "%.")
    if (pwr.1 < 0.80)
      txt <- paste0(txt, "\n2nd stage with ", n2, " subjects initiated.")
  }
  if (!pass.1 & pwr.1 < 0.80) {
    if (pass.2) {
      txt <- paste0(txt, "\nStudy passed in the 2nd stage: ")
    } else {
      txt <- paste0(txt, "\nStudy failed in the 2nd stage: ")
    }
  }
  if (!pass.1 & pwr.1 < 0.80) {
    pwr.21  <- Power2Stage:::.calc.power(alpha=alpha, ltheta1=ltheta1,
                                         ltheta2=ltheta2, diffm=ltheta0,
                                         sem=sem, df=df2, method="shifted")
    pwr.22  <- Power2Stage:::.calc.power(alpha=alpha, ltheta1=ltheta1,
                                         ltheta2=ltheta2, diffm=log.GMR2,
                                         sem=sem, df=df2, method="shifted")
    txt <- paste0(txt, 100*(1-2*alpha), "% CI = ",
                  CI.stg2[["lower"]], "% \u2013 ", CI.stg2[["upper"]], "%.")
    txt <- paste0(txt, "\n(Irrelevant) post hoc powers:",
                  "\nFor GMR 0.95 = ", round(100*pwr.21, 1), "%, ",
                  "for GMR ", round(exp(log.GMR2), 4), " = ",
                  round(100*pwr.22, 1), "%.")
  }
  cat(txt, "\n")
}

Example 1:

potvin.examples(n1=12, log.GMR1=0.16785, S.sq1=0.020977,
                log.GMR2=0.14401, S.sq2=0.0212240)
Power in the interim = 84.1%.
Since power is ≥80%, BE assessed with alpha 0.05.
Study failed in the 1st stage:
90% CI = 106.26% – 131.66%.

Example 2:

potvin.examples(n1=12, log.GMR1=0.08396, S.sq1=0.032634,
                log.GMR2=0.014439, S.sq2=0.045896)
Power in the interim = 64.9%.
Since power is <80%, BE assessed with alpha 0.0294.
Study failed in the 1st stage:
94.12% CI = 92.93% – 127.28%.
2nd stage with
8 subjects initiated.
Study passed in the 2nd stage:
94.12% CI = 88.45% – 116.38%.
(Irrelevant) post hoc powers:
For GMR 0.95 =
66.3%, for GMR 1.0145 = 76.9%.


Dif-tor heh smusma 🖖
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
21,753 posts in 4,548 threads, 1,544 registered users;
online 3 (0 registered, 3 guests [including 3 identified bots]).
Forum time: Sunday 05:31 CEST (Europe/Vienna)

They were “so intent of making everything numerical”
that they frequently missed seeing
what was there to be seen.    Barbara McClintock

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