Bioequivalence and Bioavailability Forum

Dear all,

following a question at LinkedIn’s non-public user-group “BA/BE professionals”, code which combines some of PowerTOST’s functions to

• calculate the CVs of several studies by the function CVfromCI(),
  
• calculate 95% CI of the CVs by the function CVCL(),
  
• pool the CVs and calculate the 80% (upper/lower) CL of the pooled CV by the function CVpooled(),
  
• check whether the studies were balanced, and
  
• estimate the sample size of a planned study (given its design, expected T/R-ratio, desired power, and regulatory body) by the functions sampleN.TOST() – crossover, paired, parallel; sampleN.scABEL() – for the EMA or Health Canada; sampleN.RSABE() – for the FDA. The FDA’s reference-scaling method for NTIDs by sampleN.NTIDFDA() is not implemented yet.

###################
# your data below #
###################
design <- c("3x6x3", "2x2x2", "2x2x4")
N      <- c(17, 34, 33)             # (total) sample sizes
lower  <- c(0.8666, 1.0211, 0.9750) # lower CL
upper  <- c(1.2845, 1.2665, 1.1390) # upper CL
CI     <- c(95, 90, 90)             # level of the CI in %
theta0 <- 1.10                      # expected T/R-ratio
target <- 0.80                      # target power
des    <- "2x2x4"                   # planned design
RSABE  <- FALSE                     # reference-scaling desired?
reg    <- "EMA"                     # any of "EMA", "FDA", "HC"
adjust <- FALSE                     # iteratively adjust alpha to
                                    # prevent inflation of the TIE?
##########################################################
# don't change below unless you know what your are doing #
##########################################################
library(PowerTOST)
alpha  <- (1-CI/100)/2
digits <- max(nchar(as.character(c(lower, upper))))-2
if (RSABE) RSABE.des <- "yes" else RSABE.des <- "no"
PE     <- round(sqrt(lower*upper), digits)
CL     <- df <- calc <- bal <- CVs <- CLlo <- CLhi <-vector()
# calculate the CVs and their confidence intervals
for (j in seq_along(N)) {
  CL[j]   <- paste0(CI[j], "%")
  n       <- N[j]
  df[j]   <- eval(parse(text=known.designs()[which(known.designs()[["design"]] == design[j]), "df"],
                        srcfile=NULL))
  calc[j] <- known.designs()[which(known.designs()[["design"]] == design[j]), "df"]
  steps   <- known.designs()[which(known.designs()[["design"]] == design[j]), "steps"]
  if (n %% steps == 0) bal[j] <- "yes" else bal[j] <- "no"
  CVs[j]  <- CVfromCI(upper=lower[j], lower=upper[j], n=N[j], design=design[j], alpha=alpha[j])
  CLlo[j] <- CVCL(CV=CVs[j], df=df[j], side="2-sided", alpha=0.05)[[1]]
  CLhi[j] <- CVCL(CV=CVs[j], df=df[j], side="2-sided", alpha=0.05)[[2]]
}
CVdata <- data.frame(PE, CL, lower, upper, CVs,  CLlo, CLhi, N, design,
                     paste("study", 1:length(N)), bal, calc, paste("=", df))
names(CVdata) <- c("PE", "CI", "lo", "hi", "CV", "lower CL", "upper CL",
                   "n", "design", "source", "balanced", "df", "")
# calculate the pooled CV and the 80% upper one-sided CL
CVpooled <- CVpooled(CVdata[, c(5, 8:10)], alpha=0.2)
CVdata[, 5:7] <- round(CVdata[, 5:7], digits) # precision like input
# estimate sample size and expected power for the pooled CV
# and its upper CL
alpha.1  <- rep(0.05, 2)
if (RSABE) { # reference-scaling
  CL.lo <- CVCL(CV=CVpooled[["CV"]], df=CVpooled[["df"]],
                side="lower", alpha=0.2)[["lower CL"]]
  CL.hi <- CVCL(CV=CVpooled[["CV"]], df=CVpooled[["df"]],
                side="upper", alpha=0.2)[["upper CL"]]
  txt   <- paste("Pooled CV =", round(CVpooled[["CV"]], digits),
                 "with", CVpooled[["df"]], "degrees of freedom")
  if (reg == "FDA") { # the FDA's RSABE
    if (des == "2x3x3" | des == "2x2x4" | des == "2x2x3") {
      tmp <- sampleN.RSABE(CV=CVpooled[["CV"]], theta0=theta0,
                           targetpower=target, design=des,
                           details=FALSE, print=FALSE)
      n.est1   <- tmp[["Sample size"]]
      pwr.est1 <- round(tmp[["Achieved power"]], digits)
      L.1      <- scABEL(CV=CVpooled[["CV"]], regulator=reg)[["lower"]]
      U.1      <- scABEL(CV=CVpooled[["CV"]], regulator=reg)[["upper"]]
      tmp1 <- sampleN.RSABE(CV=CL.lo, theta0=theta0,
                            targetpower=target, design=des,
                            details=FALSE, print=FALSE)
      tmp2 <- sampleN.RSABE(CV=CL.hi, theta0=theta0,
                            targetpower=target, design=des,
                            details=FALSE, print=FALSE)
      if (tmp2[["Sample size"]] >= tmp1[["Sample size"]]) { # critical: upper CL
        n.est2   <- tmp2[["Sample size"]]
        pwr.est2 <- round(tmp2[["Achieved power"]], digits)
        L.2      <- scABEL(CV=CL.hi, regulator=reg)[["lower"]]
        U.2      <- scABEL(CV=CL.hi, regulator=reg)[["upper"]]
        txt      <- paste(txt, "\nUpper 80% confidence limit of the CV =",
                          round(CL.hi, digits), "\n")
      } else { # critical: lower CL
        n.est2   <- tmp1[["Sample size"]]
        pwr.est2 <- round(tmp1[["Achieved power"]], digits)
        L.2      <- scABEL(CV=CL.lo, regulator=reg)[["lower"]]
        U.2      <- scABEL(CV=CL.lo, regulator=reg)[["upper"]]
        txt      <- paste(txt, "\nUpper 80% confidence limit of the CV =",
                          round(CL.lo, digits), "\n")
      }
    } else {
      stop("Design not implemented.")
    }
  } else { # the EMA's and Health Canada's ABEL
    if (des == "2x3x3" | des == "2x2x4" | des == "2x2x3") {
      if (!adjust) { # according to the GL: alpha 0.05
        tmp <- sampleN.scABEL(CV=CVpooled[["CV"]], theta0=theta0,
                              targetpower=target, design=des,
                              regulator=reg, details=FALSE, print=FALSE)
        n.est1   <- tmp[["Sample size"]]
        pwr.est1 <- round(tmp[["Achieved power"]], digits)
        L.1      <- scABEL(CV=CVpooled[["CV"]], regulator=reg)[["lower"]]
        U.1      <- scABEL(CV=CVpooled[["CV"]], regulator=reg)[["upper"]]
        tmp1 <- sampleN.scABEL(CV=CL.lo, theta0=theta0,
                               targetpower=target, design=des,
                               regulator=reg, details=FALSE, print=FALSE)
        tmp2 <- sampleN.scABEL(CV=CL.hi, theta0=theta0,
                               targetpower=target, design=des,
                               regulator=reg, details=FALSE, print=FALSE)
        if (tmp2[["Sample size"]] >= tmp1[["Sample size"]]) { # critical: upper CL
          n.est2   <- tmp2[["Sample size"]]
          pwr.est2 <- round(tmp2[["Achieved power"]], digits)
          L.2      <- scABEL(CV=CL.hi, regulator=reg)[["lower"]]
          U.2      <- scABEL(CV=CL.hi, regulator=reg)[["upper"]]
          txt      <- paste(txt, "\nUpper 80% confidence limit of the CV =",
                            round(CL.hi, digits), "\n")
        } else { # critical: lower CL
          n.est2   <- tmp1[["Sample size"]]
          pwr.est2 <- round(tmp1[["Achieved power"]], digits)
          L.2      <- scABEL(CV=CL.lo, regulator=reg)[["lower"]]
          U.2      <- scABEL(CV=CL.lo, regulator=reg)[["upper"]]
          txt      <- paste(txt, "\nUpper 80% confidence limit of the CV =",
                            round(CL.lo, digits), "\n")
        }
      } else { # iteratively adjusted alpha
        tmp <- sampleN.scABEL.ad(CV=CVpooled[["CV"]], theta0=theta0,
                                 targetpower=target, design=des,
                                 regulator=reg, details=FALSE, print=FALSE)
        n.est1   <- tmp[["Sample size"]]
        pwr.est1 <- round(tmp[["Achieved power"]], digits)
        L.1      <- scABEL(CV=CVpooled[["CV"]], regulator=reg)[["lower"]]
        U.1      <- scABEL(CV=CVpooled[["CV"]], regulator=reg)[["upper"]]
        if (!is.na(tmp[["adj. alpha"]])) alpha.1[1] <- tmp[["adj. alpha"]]
        tmp1 <- sampleN.scABEL.ad(CV=CL.lo, theta0=theta0,
                                 targetpower=target, design=des,
                                 regulator=reg, details=FALSE, print=FALSE)
        tmp1 <- sampleN.scABEL.ad(CV=CL.hi, theta0=theta0,
                                 targetpower=target, design=des,
                                 regulator=reg, details=FALSE, print=FALSE)
        if (tmp2[["Sample size"]] >= tmp1[["Sample size"]]) { # critical: upper CL
          n.est2   <- tmp[["Sample size"]]
          pwr.est2 <- round(tmp[["Achieved power"]], digits)
          L.2      <- scABEL(CV=CL.hi, regulator=reg)[["lower"]]
          U.2      <- scABEL(CV=CL.hi, regulator=reg)[["upper"]]
        } else { # critical: lower CL
          n.est2   <- tmp1[["Sample size"]]
          pwr.est2 <- round(tmp1[["Achieved power"]], digits)
          L.2      <- scABEL(CV=CL.lo, regulator=reg)[["lower"]]
          U.2      <- scABEL(CV=CL.lo, regulator=reg)[["upper"]]
        }
        if (!is.na(tmp[["adj. alpha"]])) alpha.1[2] <- tmp[["adj. alpha"]]
      }
    } else {
      stop("Design not implemented.")
    }
  }
} else { # conventional ABE
  tmp <- sampleN.TOST(CV=CVpooled[["CV"]], theta0=theta0,
                      targetpower=target, design=des, print=FALSE)
  n.est1   <- tmp[["Sample size"]]
  pwr.est1 <- round(tmp[["Achieved power"]], digits)
  tmp <- sampleN.TOST(CV=CVpooled[["CVupper"]], theta0=theta0,
                      targetpower=target, design=des, print=FALSE)
  n.est2   <- tmp[["Sample size"]]
  pwr.est2 <- round(tmp[["Achieved power"]], digits)
  L.2 <- L.1 <- 0.8
  U.2 <- U.1 <- 1.25
  txt <- paste("Pooled CV =", round(CVpooled[["CV"]], digits), "with",
               CVpooled[["df"]], "degrees of freedom",
               "\nUpper 80% confidence limit of the CV =",
               round(CVpooled[["CVupper"]], digits), "\n")
}
Nest <- data.frame(rep(des, 2), rep(RSABE.des, 2), rep(reg, 2),
                   rep(sprintf("%.4f", theta0), 2),
                   sprintf("%.4f", c(L.1, L.2)),
                   sprintf("%.4f", c(U.1, U.2)), alpha.1,
                   round(c(CVpooled[["CV"]], CVpooled[["CVupper"]]), digits),
                   c(n.est1, n.est2), c(pwr.est1, pwr.est2))
names(Nest) <- c("Design", "RSABE desired", "regulator", "theta0",
                 "theta1", "theta2", "alpha", "CV", "n", "Achieved power")
if (length(unique(design)) > 1 & "parallel" %in% design)
  warning("Mixing crossover and parallel studies = nonsense!")
print(as.data.frame(CVdata), row.names=F);cat(txt);print(as.data.frame(Nest), row.names=FALSE)

Which gives:
     PE  CI     lo     hi     CV lower CL upper CL  n design  source balanced    df     
 1.0551 95% 0.8666 1.2845 0.2832   0.2247   0.3844 17  3x6x3 study 1       no 2*n-4 = 30
 1.1372 90% 1.0211 1.2665 0.2667   0.2132   0.3574 34  2x2x2 study 2      yes   n-2 = 32
 1.0538 90% 0.9750 1.1390 0.2736   0.2386   0.3210 33  2x2x4 study 3       no 3*n-4 = 95
Pooled CV = 0.2741 with 157 degrees of freedom
Upper 80% confidence limit of CV = 0.2888
 Design RSABE desired regulator theta0 theta1 theta2 alpha     CV  n Achieved power
  2x2x4            no       EMA 1.1000 0.8000 1.2500  0.05 0.2741 28         0.8015
  2x2x4            no       EMA 1.1000 0.8000 1.2500  0.05 0.2888 32         0.8135

Doing the math is just the first step. Before you pool CVs I would suggest to inspect whether the confidence intervals of the CVs overlap (columns 6–7). If not, try to find out why (different CROs, populations, bioanalytical methods, …). Use common sense to decide which CVs are reliable enough to pool.

Another question at LinkedIn was whether one should use the highest CV of several studies for conservatism. Maybe, maybe not. If the sample sizes were very close and the CVs not too different, why not. Otherwise (or if different designs were used) I prefer the upper CL of the pooled CV.

Note:

• Since in the example two studies are unbalanced, the function CVfromCI() throws these messages:
    Unbalanced 3x6x3 design. n(i)= 3/3/3/3/3/2 assumed.
  Unbalanced 2x2x4 design. n(i)= 17/16 assumed.
CVfromCI() tries to keep sequences as balanced as possible – which gives a conservative estimate of the CV if the study was more unbalanced. If you know the subjects / sequences calculate the CV and use CVpooled() “manually” (see help(CVfromCI) and help(CVpooled) how to do that).
  
• The CV calculated from the CI of a replicate design is always pooled from CV_wR and CV_wT. Since they are not necessarily equal, be aware.
  
• If you opt for reference-scaling (RSABE <- TRUE) the upper CL of the CV might not be the most conservative solution. It is possible that a lower CV results in less power (since less scaling is applicable) and therefore, requires a higher sample size. The script checks for both possibilities and shows the worst case (lower/upper CL). I suggest to explore how power behaves. After executing this script run:
  plot(pa.scABE(CV=round(CVpooled[["CV"]], digits), theta0=theta0,
              targetpower=target, design=des, regulator=reg),
              pct=FALSE)
  And:
  if (tmp2[["Sample size"]] >= tmp1[["Sample size"]]) {
  plot(pa.scABE(CV=round(CL.hi, digits), theta0=theta0,
                targetpower=target, design=des, regulator=reg),
                pct=FALSE)
} else {
  plot(pa.scABE(CV=round(CL.lo, digits), theta0=theta0,
                targetpower=target, design=des, regulator=reg),
                pct=FALSE)
}