Bioequivalence and Bioavailability Forum

R-code for all ABE designs [Power / Sample Size]

posted by Helmut  – Vienna, Austria, 2014-08-12 19:21 (4326 d 15:23 ago) – Posting: # 13370
Views: 22,886

Dear Detlew and all,

below the code for all designs covered in PowerTOST(), conventional (unscaled) ABE. You have to paste the functions only once to the R-console and call the function Sens() with values defined in the end or call it directly, e.g.,
Sens(CV=0.24, GMR=1.07, pwr.target=0.8, pwr.min=0.7, des="3x6x3")



Note that the code tries to keep the degree of imbalance as close as possible between sequences. In real life other combinations may occur. As already noted by Detlew in this post the more im­bal­anced a study is, the lower the power.*
Example (3 drop-outs from planned 18 ⇒ 15)
power2.TOST(CV=0.2, theta0=0.95, n=c(3, 3, 3, 2, 2, 2), design="3x6x3")
[1] 0.7089454
power2.TOST(CV=0.2, theta0=0.95, n=c(3, 3, 3, 3, 2, 1), design="3x6x3")
[1] 0.6464662

---

require(PowerTOST)
##################################################################
# Don't modify the functions unless you know what you are doing! #
# Input area area the end of the code.                           #
##################################################################
Sens <- function(CV, GMR, pwr.target, pwr.min, des) {
  pwrCV  <- function(x, ...) {
              power.TOST(CV=x, n=n.est, theta0=GMR,
              design=des)-pwr.min }
  pwrGMR <- function(x, ...) {
              power.TOST(CV=CV, n=n.est, theta0=x,
              design=des)-pwr.min }
  designs <- c(
    "2x2", "2x2x2", "2x2x3", "2x2x4",
    "3x3", "2x3x3",
    "4x4", "2x4x4", "2x4x2",
    "3x6x3",
    "paired", "parallel")
  res     <- sampleN.TOST(CV=CV, theta0=GMR, targetpower=pwr.target, design=des,
    print=F, details=F)
  n.est   <- res[1, "Sample size"   ]
  pwr.est <- res[1, "Achieved power"]
  seg     <- 50; segs <- seq(seg-1)
  ########################################
  # max. CV for minimum acceptable power #
  ########################################
  CV.max <- uniroot(pwrCV,  c(CV, 10*CV), tol=1e-7, n=n.est, theta0=GMR)$root
  CVs    <- seq(CV, CV.max, length.out=seg)
  pBECV  <- NULL
  for(j in 1:length(CVs)) {
    pBECV <- c(pBECV, 100*power.TOST(CV=CVs[j], n=n.est, theta0=GMR, design=des))
  }
  ######################################
  # min. GMR for minimum accept. power #
  ######################################
  ifelse(GMR <= 1, int <- c(0.8, 1), int <- c(GMR, 1.25))
  GMR.min <- uniroot(pwrGMR, int, tol=1e-7, n=n.est, CV=CV, design=des)$root
  GMRs    <- seq(GMR.min, GMR, length.out=seg)
  pBEGMR  <- NULL
  for(j in 1:length(GMRs)) {
    pBEGMR <- c(pBEGMR, 100*power.TOST(CV=CV, n=n.est, theta0=GMRs[j], design=des))
  }
  ####################################
  # min. n for minimum accept. power #
  # workaround, since uniroot() does #
  # not accept two vectors as limits #
  ####################################
  Ns <- seq(n.est, 12, by=-1) # don't drop below 12 subjects
  j <- 0
  pwrN <- pwr.est
  n.min <- NULL; pBEn <- NULL
  if(des == "paired") {
    while(pwrN >= pwr.min) {
      j <- j+1
      pwrN <- power.TOST(CV=CV, n=Ns[j], theta0=GMR, design=des)
      if(pwrN >= pwr.min) {
        n.min <- c(n.min, Ns[j])
        pBEn <- c(pBEn, 100*pwrN)
      } else {
        break
      }
    }
  }
  if(des == "3x3" | des == "2x3x3"){ # three sequences
    while(pwrN >= pwr.min) {
      j <- j+1
      n1 <- floor(Ns[j]/3)
      n2 <- floor(Ns[j]*2/3-n1)
      n3 <- Ns[j]-(n1+n2)
      n <- c(n1, n2, n3)
      pwrN <- power2.TOST(CV=CV, n=n, theta0=GMR, design=des)
      if(pwrN >= pwr.min) {
        n.min <- c(n.min, sum(n))
        pBEn <- c(pBEn, 100*pwrN)
      } else {
        break
      }
    }
  }
  if(des == "4x4" | des == "2x4x4" | des == "2x4x2"){ # four seq's
    while(pwrN >= pwr.min) {
      j <- j+1
      n1 <- floor(Ns[j]/4)
      n2 <- floor(Ns[j]*2/4-n1)
      n3 <- floor(Ns[j]*3/4-(n1+n2))
      n4 <- Ns[j]-(n1+n2+n3)
      n <- c(n1, n2, n3, n4)
      pwrN <- power2.TOST(CV=CV, n=n, theta0=GMR, design=des)
      if(pwrN >= pwr.min) {
        n.min <- c(n.min, sum(n))
        pBEn <- c(pBEn, 100*pwrN)
      } else {
        break
      }
    }
  }
  if(des == "3x6x3"){
    while(pwrN >= pwr.min) { # six seq's
      j <- j+1
      n1 <- floor(Ns[j]/6)
      n2 <- floor(Ns[j]*2/6-n1)
      n3 <- floor(Ns[j]*3/6-(n1+n2))
      n4 <- floor(Ns[j]*4/6-(n1+n2+n3))
      n5 <- floor(Ns[j]*5/6-(n1+n2+n3+n4))
      n6 <- Ns[j]-(n1+n2+n3+n4+n5)
      n <- c(n1, n2, n3, n4, n5, n6)
      pwrN <- power2.TOST(CV=CV, n=n, theta0=GMR, design=des)
      if(pwrN >= pwr.min) {
        n.min <- c(n.min, sum(n))
        pBEn <- c(pBEn, 100*pwrN)
      } else {
        break
      }
    }
  } else { # parallel and all two sequence designs
    while(pwrN >= pwr.min) {
      j <- j+1
      n1 <- floor(Ns[j]/2)
      n2 <- Ns[j]-n1
      n <- c(n1, n2)
      pwrN <- power2.TOST(CV=CV, n=n, theta0=GMR, design=des)
      if(pwrN >= pwr.min) {
        n.min <- c(n.min, sum(n))
        pBEn <- c(pBEn, 100*pwrN)
      } else {
        break
      }
    }
  }
  op <- par()
  par(mar=c(c(4, 4, 2.5, 0.5))+0.1) # default for B, L, T, R: c(5, 4, 4, 2) + 0.1
  clr <- colorRampPalette(c("blue", "red"))(seg)
  split.screen(c(2, 2))
  screen(1) ### Sensitivity of CV (GMR and n constant) ###
  plot(100*CVs, pBECV, type="n", main=paste0("Higher variability\n",
    "GMR = ", GMR, ", n = ", n.est), lwd=2, xlab="CV%", ylab="% power", las=1,
    cex.main=1)
  abline(h=c(100*pwr.target, 80, 100*pwr.min), lty=3, col="grey50")
  segments(100*CVs[segs], pBECV[segs], 100*CVs[segs+1], pBECV[segs+1], lwd=2, col=clr[segs])
  points(100*CVs[1], pBECV[1], col=clr[1], pch=16, cex=1.25)
  points(100*CVs[seg], pBECV[seg], col=clr[seg], pch=16, cex=1.25)
  text(100*CV, 100*(pwr.min+(pwr.est-pwr.min)*0.1), labels=paste0("CV = ",
    signif(100*CV.max, 4), "% (", round(100*pwr.min, 2), "%)"), cex=0.9, pos=4)
  screen(2) ### Sensitivity of GMT (CV and n constant) ###
  clr <- rev(clr)
  plot(GMRs, pBEGMR, type="n", main=paste0("Larger deviation from 1\n",
    "CV = ", 100*CV, "%, n = ", n.est), lwd=2, xlim=c(GMR, GMR.min),
    xlab="GMR", ylab="% power", las=1, cex.main=1)
  abline(h=c(100*pwr.target, 80, 100*pwr.min), lty=3, col="grey50")
  segments(GMRs[segs], pBEGMR[segs], GMRs[segs+1], pBEGMR[segs+1], lwd=2, col=clr[segs])
  points(GMRs[1], pBEGMR[1], col=clr[1], pch=16, cex=1.25)
  points(GMRs[seg], pBEGMR[seg], col=clr[seg], pch=16, cex=1.25)
  text(GMR, 100*(pwr.min+(pwr.est-pwr.min)*0.1), labels=paste0("GMR = ",
    signif(GMR.min, 4), " (", round(100*pwr.min, 2), "%)"), cex=0.9, pos=4)
  screen(3) ### Sensitivity of n (GMR and CV constant) ###
  plot(n.min, pBEn, type="n", main=paste0("Drop-outs\n",
    "GMR = ", GMR, ", CV = ", 100*CV, "%"), lwd=2, xlim=c(max(n.min), min(n.min)),
    ylim=c(100*pwr.min, 100*pwr.est), xlab="n", ylab="% power", las=1, cex.main=1)
  abline(h=c(100*pwr.target, 80, 100*pwr.min), lty=3, col="grey50")
  clr <- colorRampPalette(c("blue", "red"))(length(n.min))
  points(n.min, pBEn, pch=16, cex=0.8, col=clr)
  points(n.min[1], pBEn[1], col=clr[1], pch=16, cex=1.25)
  points(n.min[length(n.min)], pBEn[length(n.min)], col=clr[length(n.min)],
    pch=16, cex=1.25)
  text(max(n.min), 100*(pwr.min+(pwr.est-pwr.min)*0.1), labels=paste0("n = ",
    min(n.min), " (", signif(min(pBEn), 4), "%)"), cex=0.9, pos=4)
  screen(4) ### Some basic information ###
  plot(1, type="n", axes=F, xlab="", ylab="")
  legend("topleft", legend=c(paste0(des, " design", "; assumed:"),
    sprintf("  %s %1.0f%%%s %5.4f", "CV =", 100*CV, ", GMR =", GMR), "power:",
    sprintf("  %s %2.0f%%", "target =", 100*pwr.target),
    sprintf("  %s %5.2f%% %s %i%s", "expected =", 100*pwr.est, "(n =", n.est, ")"),
    sprintf("  %s %2.0f%%", "minimum acceptable =", 100*pwr.min), "acceptable rel. deviations:",
    sprintf("  %s %+5.1f%%", "CV =", 100*(CV.max-CV)/CV),
    sprintf("  %s %+5.2f%%", "GMR =", 100*(GMR.min-GMR)/GMR),
    sprintf("  %s %+5.1f%%", "n =", 100*(min(n.min)-n.est)/n.est)), bty="n", cex=1)
  close.screen(all=T)
  par <- op
  cat(sprintf("%s %s%1.0f%%%s%5.4f%s%i%s%5.2f%% %s %2.0f%%%s", des,
    "design:\nCV=", 100*CV, ", GMR=", GMR, ", n=", n.est, ", expected power=",
    100*pwr.est, "(>target of", 100*pwr.target, ")\n\n"))
}
#########################################################
# You need to paste the functions above only once to R. #
# Define values below (as fractions, not percent):      #
# #######################################################
CV     <- 0.20  # intra-subject for cross-overs & paired, total for parallel.
GMR    <- 0.95  # [sic]
target <- 0.90  # The power aimed at in designing the study.
min    <- 0.70  # The lowest power which is acceptable for the actual study.
des    <- "2x2" # See on top of the function for the implemented strings.
                # Call known.designs() to get information from PowerTOST
                # (dfs, design constants, ...)
######################################################
# Once you have defined values, call the function... #
######################################################
Sens(CV, GMR, target, min, des)
########################################
# ... or call it directly with values: #
########################################
# Sens(CV=0.24, GMR=1.07, pwr.target=0.8, pwr.min=0.7, des="3x6x3")

@R-Freaks: The construction of imbalanced sequences is both lengthy and clumsy. Any suggestions are welcome.

Too busy to develop the code for RSABE and ABEL at the moment.

---

• If a sequence is completely lost, power2.TOST() – correctly? – gives power=0. See these two examples:
  power2.TOST(CV=0.2, n=c(1, 2, 3, 3, 3, 3), design="3x6x3")
power2.TOST(CV=0.2, n=c(0, 3, 3, 3, 3, 3), design="3x6x3")

—
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:

• Deviating from assumptions Helmut 2014-08-08 15:44 
  • Some Nitpicking d_labes 2014-08-12 09:12
    • Some Nitpicking Helmut 2014-08-12 10:49
    • R-code for all ABE designsHelmut 2014-08-12 17:21
      • Sensitivity analysis for all ABE designs d_labes 2014-08-13 09:30
        • Sensitivity analysis for all ABE designs Helmut 2014-08-13 14:49
      • R-code shortening d_labes 2014-08-13 16:32
        • Suggestions / Sneak Preview Helmut 2014-08-13 16:42
          • Suggestions / Sneak Preview d_labes 2014-08-15 09:01
            • Suggestions / Sneak Preview Helmut 2014-08-15 12:02
          • Mehl returned! d_labes 2014-08-15 11:27
            • Mehl returned! Helmut 2014-08-15 11:42