Bioequivalence and Bioavailability Forum

Hi Shuanghe,

❝ I'll be in Amsterdam this September for the 5th global BE harmonisation conference, so hopefully we'll have time to have a drink together

Great! BioBridges (22/23 September) as well?

❝ By the way, the Article section is new, am I right? (Don't tell me I missed this section all those years...)

<1½ years.

❝ I am supprised that no one mentioned it in the forum.

Well, I do. Actually I was bored of answering similar questions over and over again.
Now I reply only with: “, for details see article XYZ.”

❝ I don't have time to read the lengthy ones yet but I really enjoy the reading among the few short ones that I've picked.

THX!

❝ What I understood about balanced IBD here is that, e.g., the number for A vs. F comparison should be same as number of F vs. A. So there should be 10 for each comparison for 20 in total.

❝ If that is a correct understanding, …

Of course, it is.

❝ … it seems that for B, C and D, yes, but for A and E, not... There're 11 AF (FE) but 9 FA (EF). Shouldn't it be 10 and 10 like B/C/D vs F?

THX for discovering such a nasty bug!

If debugging is the process of removing bugs,
then programming must be the process of putting them in.     Edsger W. Dijkstra

New (lengthy) script at the end. More flexible (theta0, theta1, theta2, target, do.rate, and a vector of references ref can be specified).
Call the function with e.g., make.ibds(ntmt = 6, CV = 0.2):

 subject seqno sequence  IBD 1  IBD 2  IBD 3  IBD 4  IBD 5
       1     5   ECFABD ••FA•• ••F•B• •CF••• ••F••D E•F•••
       2     6   FEBCDA F••••A F•B••• F••C•• F•••D• FE••••
       3     5   ECFABD ••FA•• ••F•B• •CF••• ••F••D E•F•••
       4     6   FEBCDA F••••A F•B••• F••C•• F•••D• FE••••
       5     3   CAEDFB •A••F• ••••FB C•••F• •••DF• ••E•F•
       6     1   ADCBEF A••••F •••B•F ••C••F •D•••F ••••EF
       7     2   BFDEAC •F••A• BF•••• •F•••C •FD••• •F•E••
       8     4   DBAFCE ••AF•• •B•F•• •••FC• D••F•• •••F•E
       9     3   CAEDFB •A••F• ••••FB C•••F• •••DF• ••E•F•
      10     2   BFDEAC •F••A• BF•••• •F•••C •FD••• •F•E••
      11     1   ADCBEF A••••F •••B•F ••C••F •D•••F ••••EF
      12     4   DBAFCE ••AF•• •B•F•• •••FC• D••F•• •••F•E
      13     4   DBAFCE ••AF•• •B•F•• •••FC• D••F•• •••F•E
      14     3   CAEDFB •A••F• ••••FB C•••F• •••DF• ••E•F•
      15     1   ADCBEF A••••F •••B•F ••C••F •D•••F ••••EF
      16     2   BFDEAC •F••A• BF•••• •F•••C •FD••• •F•E••
      17     5   ECFABD ••FA•• ••F•B• •CF••• ••F••D E•F•••
      18     2   BFDEAC •F••A• BF•••• •F•••C •FD••• •F•E••
      19     6   FEBCDA F••••A F•B••• F••C•• F•••D• FE••••
      20     3   CAEDFB •A••F• ••••FB C•••F• •••DF• ••E•F•

Reference             : F
Tests                 : A, B, C, D, E
Sequences             : ADCBEF, BFDEAC, CAEDFB, DBAFCE, ECFABD, FEBCDA
Subjects per sequence : 3 | 4 | 4 | 3 | 3 | 3 (imbalanced)
Estimated sample size : 20
Achieved power        : 0.8347
Randomized            : 2022-07-13 22:55:28 CEST
Seed                  : 820697
Runtime               : 30 ms (1 call for balance of IBDs)

Since we are working in a regulated environment, perhaps you will be asked by an inspector how you performed the randomization. You need to provide the sample size, all sequences and the seed to reproduce it. Try:

library(randomizeBE)
RL4(n = 20, seqs = c("ADCBEF", "BFDEAC", "CAEDFB", "DBAFCE", "ECFABD", "FEBCDA"),
    seed = 820697)

Note that in this case you get at the end:

Warning message:
Unbalanced design! # of subj. in sequences  ADCBEF  BFDEAC  CAEDFB  DBAFCE  ECFABD  FEBCDA:  3  4  4  3  3  3

Is that important? Although we are interested in balanced IBDs – which they are – they are not balanced for period effects. Let’s check how often the treatments are administered in each period:

  p1 p2 p3 p4 p5 p6
A  3  4  3  3  4  3
B  4  3  3  3  3  4
C  4  3  3  3  3  4
D  3  3  4  4  3  3
E  3  3  4  4  3  3
F  3  4  3  3  4  3

If you want period-balance (which is a desirable property), use bal = TRUE in the call:

tmp <- make.ibds(ntmt = 6, CV = 0.2, print = FALSE, bal = TRUE)
cat(tmp$txt, "\n")
Reference                : F
Tests                    : A, B, C, D, E
Sequences                : AEBDFC, BAFECD, CFDBEA, DCEFAB, EDACBF, FBCADE
Subjects per sequence    : 4 | 4 | 4 | 4 | 4 | 4 (balanced)
Estimated sample size    : 20
Achieved power           : 0.8347
Adjustment to obtain period-balance
 Adjusted sample size    : 24
 Achieved power          : 0.8960
Randomized               : 2022-07-15 11:13:50 CEST
Seed                     : 9503537
Runtime                  : 20 ms (1 call for balance of IBDs)

There is no free lunch; we need 24 subjects instead of 20.

If you are interested in the performance:

calls <- et <- integer(1000)
for (j in seq_along(calls)) {
  tmp      <- make.ibds(ntmt = 6, CV = 0.2, print = FALSE)
  calls[j] <- as.integer(tmp$calls)
  et[j]    <- as.integer(tmp$et)
}
probs <- c(0, 0.025, 0.25, 0.5, 0.75, 0.975, 1)
perf  <- matrix(data = c(quantile(calls, probs = probs),
                         round(quantile(et, probs = probs))),
                         nrow = 2, ncol = length(probs), byrow = TRUE,
                dimnames = list(c("calls", "runtimes (ms)"),
                                paste0(100 * probs, "%")))
print(perf)

On my aged machine:

              0% 2.5% 25% 50% 75% 97.5% 100%
calls          1    1   3   7  14    35   59
runtimes (ms)  0   10  30  70 140   360  610

More arguments (NTID, T/R-ratio 0.975, power ≥90%, adjusting for anticipated dro­pout-rate 5%, runtime suppressed):

make.ibds(ntmt = 4, CV = 0.1, theta0 = 0.975, theta1 = 0.9, target = 0.9,
          do.rate = 0.05, sep = "–", details = FALSE)
 subject seqno sequence IBD 1 IBD 2 IBD 3
       1     2     BDAC  –DA–  BD––  –D–C
       2     3     CADB  –AD–  ––DB  C–D–
       3     1     ABCD  A––D  –B–D  ––CD
       4     4     DCBA  D––A  D–B–  DC––
       5     3     CADB  –AD–  ––DB  C–D–
       6     4     DCBA  D––A  D–B–  DC––
       7     1     ABCD  A––D  –B–D  ––CD
       8     2     BDAC  –DA–  BD––  –D–C
       9     1     ABCD  A––D  –B–D  ––CD
      10     4     DCBA  D––A  D–B–  DC––
      11     3     CADB  –AD–  ––DB  C–D–
      12     3     CADB  –AD–  ––DB  C–D–
      13     4     DCBA  D––A  D–B–  DC––
      14     1     ABCD  A––D  –B–D  ––CD
      15     2     BDAC  –DA–  BD––  –D–C
      16     2     BDAC  –DA–  BD––  –D–C
      17     2     BDAC  –DA–  BD––  –D–C
      18     2     BDAC  –DA–  BD––  –D–C
      19     4     DCBA  D––A  D–B–  DC––
      20     1     ABCD  A––D  –B–D  ––CD
      21     4     DCBA  D––A  D–B–  DC––
      22     1     ABCD  A––D  –B–D  ––CD
      23     3     CADB  –AD–  ––DB  C–D–
      24     3     CADB  –AD–  ––DB  C–D–
      25     3     CADB  –AD–  ––DB  C–D–
      26     1     ABCD  A––D  –B–D  ––CD
      27     2     BDAC  –DA–  BD––  –D–C
      28     1     ABCD  A––D  –B–D  ––CD
      29     4     DCBA  D––A  D–B–  DC––
      30     4     DCBA  D––A  D–B–  DC––
      31     3     CADB  –AD–  ––DB  C–D–
      32     2     BDAC  –DA–  BD––  –D–C

Reference                : D
Tests                    : A, B, C
Sequences                : ABCD, BDAC, CADB, DCBA
Subjects per sequence    : 8 | 8 | 8 | 8 (balanced)
Estimated sample size    : 30
Achieved power           : 0.9167
 Anticipated dropout-rate: 5%
 Adjusted sample size    : 32
 Achieved power          : 0.9318
Randomized               : 2022-07-13 23:02:13 CEST
Seed                     : 874461

You can also assign the stuff to a list-variable for checking. No counting any more. 

tmp <- make.ibds(ntmt = 6, CV = 0.2)
tmp[2:6]
$n
[1] 20

$sequences
[1] "ACFDEB" "BEDFCA" "CDABFE" "DBCEAF" "EFBADC" "FAECBD"

$n.seq
[1] 4 4 3 3 3 3

$info
[1] "imbalanced"

$ibds
      seq  n seq  n
IBD 1  AF 10  FA 10
IBD 2  BF 10  FB 10
IBD 3  CF 10  FC 10
IBD 4  DF 10  FD 10
IBD 5  EF 10  FE 10

Comparison of one test with two references:

make.ibds(ntmt = 3, CV = 0.15, ref = c("B", "C"), sep = "–")
 subject seqno sequence IBD 1 IBD 2
       1     3      BAC   BA–   –AC
       2     3      BAC   BA–   –AC
       3     1      ABC   AB–   A–C
       4     6      CBA   –BA   C–A
       5     1      ABC   AB–   A–C
       6     2      ACB   A–B   AC–
       7     6      CBA   –BA   C–A
       8     5      CAB   –AB   CA–
       9     4      BCA   B–A   –CA
      10     4      BCA   B–A   –CA
      11     5      CAB   –AB   CA–
      12     2      ACB   A–B   AC–

References            : B, C
Test                  : A
Sequences             : ABC, ACB, BAC, BCA, CAB, CBA
Subjects per sequence : 2 | 2 | 2 | 2 | 2 | 2 (balanced)
Estimated sample size : 12
Achieved power        : 0.8305
Randomized            : 2022-07-13 23:04:24 CEST
Seed                  : 9664508
Runtime               : 20 ms (1 call for balance of IBDs)

If all pairwise comparisons are required (i.e., additionally between the tests), the number rises quickly:
$$N_\text{comp}=\frac{n!}{2\,(n-2)!}$$ For six treatments we have already 15 pairwise comparisons.

---

make.ibds <- function(ntmt = 3, CV, alpha = 0.05, theta0 = 0.95,
                      theta1, theta2, target = 0.80, ref, do.rate = 0,
                      bal = FALSE, sep = "•", print = TRUE, details = TRUE) {
  ##################################################################
  # Estimates the sample size of a Higher-Order (Williams’) design #
  # to be evaluated by the ‘Two at a Time’ approach, i.e., with    #
  # extracted Incomplete Block Designs. All IBDs will be balanced. #
  # Note: Balanced IBDs cannot be extracted from Latin Squares.    #
  # Optionally the estimated sample size can be adjusted for an    #
  # anticipated dropout-rate and/or uprounded to obtain balanced   #
  # sequences.                                                     #
  # -------------------------------------------------------------- #
  # TODO: All pairwise comparisons, i.e., additionally between the #
  #       tests.                                                   #
  ##################################################################
  if (ntmt < 3) stop("ntmt must be at least 3.")
  if (missing(CV)) stop("CV must be given.")
  if (missing(theta1) & missing(theta2)) theta1 <- 0.80
  if (missing(theta1)) theta1 <- 1/theta2
  if (missing(theta2)) theta2 <- 1/theta1
  if (theta0 < theta1 | theta0 > theta2) stop("theta0 must lie within BE limits.")
  if (target <= 0.5) stop("target ", target, " does not make sense.")
  if (alpha > 0.5) stop("alpha ", alpha, " does not make sense.")
  make.equal <- function(n, ns) {
    return(as.integer(ns * (n %/% ns + as.logical(n %% ns))))
  }
  nadj <- function(n, do.rate, ns) {
    return(as.integer(make.equal(n / (1 - do.rate), ns)))
  }
  suppressMessages(require(PowerTOST))
  suppressMessages(require(randomizeBE))
  st     <- proc.time()[[3]]
  # Sample size for evaluation by ‘Two at a Time’ approach.
  # Note that we are using the default design = "2x2".
  smpl   <- sampleN.TOST(alpha = alpha, CV = CV, theta0 = theta0,
                         theta1 = theta1, theta2 = theta2,
                         targetpower = target, print = FALSE)
  n.orig <- smpl[["Sample size"]]
  p.orig <- smpl[["Achieved power"]]
  # Get sequences of a Williams’ design.
  seqs   <- williams(ntmt)
  ns     <- 2 # because 2x2!
  # Adjust for dropout-rate.
  n      <- nadj(n.orig, do.rate, ns)
  # Eventual uprounding if balance is requested.
  if (bal) n <- make.equal(n, length(seqs))
  pwr    <- power.TOST(alpha = alpha, CV = CV, theta0 = theta0,
                       theta1 = theta1, theta2 = theta2, n = n)
  # Extract treatment codes.
  trts   <- sort(
              unique(
                unlist(
                  strsplit(Reduce(function(x, y) paste0(x, y), seqs), ""))))
  if (missing(ref)) { # The last treatment in lexical order is the reference.
    s   <- Reduce(function(x, y) paste0(x, y), trts)
    ref <- substr(s, nchar(s), nchar(s))
  }
  if (sum(ref %in% trts) == 0)
    stop("ref \'", ref, "\'not in trts \'", trts, "\'")
  # Extract test(s).
  tests  <- trts[!trts %in% ref]
  calls  <- 0
  repeat { # Call randomizations until all IBDs are balanced.
    # Get a randomization list.
    calls <- calls + 1
    tmp   <- suppressWarnings( # Not relevant for TaT.
               RL4(nsubj = n, seqs = seqs, randctrl = FALSE))
    rand  <- tmp$rl
    # Generate columns for the IBDs.
    n.ibd <-  length(tests) * length(ref)
    for (j in 1:n.ibd) {
      rand[[paste0("IBD.", j)]] <- NA
    }
    # Build the IBDs.
    for (j in 1:nrow(rand)) {
      c <- 3
      for (k in seq_along(ref)) {
        for (m in seq_along(tests)) {
          c          <- c  + 1
          excl       <- tests[!tests == tests[m]]
          excl       <- c(excl, ref[!ref == ref[k]])
          excl       <- paste0("[", paste(excl, collapse = ", "), "]")
          rand[j, c] <- gsub(excl, sep, rand$sequence[j])
        }
      }
    }
    # Check the IBDs for balance.
    checks  <- NA
    ibd.seq <- as.data.frame(
                 matrix(data = NA,
                        nrow = n.ibd,
                        ncol = 4, byrow = TRUE,
                        dimnames = list(names(rand)[4:ncol(rand)],
                                        paste0(rep(c("seq.", "n."), 2),
                                               rep(1:2, each = 2)))))
    for (j in 1:n.ibd) {
      ibd                 <- gsub("[^[:alnum:], ]", "", rand[3 + j])
      ibd                 <- gsub("c", "", ibd)
      ibd                 <- gsub(" ", "", unlist(strsplit(ibd, ",")))
      ibd.seq[j, c(1, 3)] <- sort(unique(ibd))
      ibd.seq[j, c(2, 4)] <- c(length(ibd[ibd == sort(unique(ibd))[1]]),
                               length(ibd[ibd == sort(unique(ibd))[2]]))
      checks[j]           <- ibd.seq[j, 2] == ibd.seq[j, 4]
    }
    if (sum(checks) == length(trts) - 1) break # Success!
  }
  ifelse (do.rate == 0, f <- "\n%-22s: ", f <- "\n%-25s: ") # Cosmetics.
  if (bal) f <- "\n%-25s: "
  ifelse (length(ref) == 1,
    txt <- paste0(sprintf(f, "Reference"), ref),
    txt <- paste0(sprintf(f, "References"), paste(ref, collapse = ", ")))
  ifelse (length(tests) == 1,
    txt <- paste0(txt, sprintf(f, "Test"), tests),
    txt <- paste0(txt, sprintf(f, "Tests"), paste(tests, collapse = ", ")))
  txt    <- paste0(txt, sprintf(f, "Sequences"), paste(seqs, collapse = ", "))
  txt    <- paste0(txt, sprintf(f, "Subjects per sequence"),
                   paste(tmp$ninseqs, collapse = " | "))
  ifelse (length(unique(tmp$ninseqs)) == 1,
    info <- "balanced", info <- "imbalanced")
  txt    <- paste0(txt, " (", info, ")")
  txt    <- paste0(txt, sprintf(f, "Estimated sample size"), n.orig,
                   sprintf(f, "Achieved power"), sprintf("%.4f", p.orig))
  if (!n == n.orig) {
    if (bal) {
      txt <- paste0(txt, "\nAdjustment to obtain period-balance")
    }
    if (do.rate > 0) {
      txt <- paste0(txt, sprintf("\n %-24s: ", "Anticipated dropout-rate"),
                    sprintf("%.3g%%", 100 * do.rate))
    }
    txt <- paste0(txt, sprintf("\n %-24s: ", "Adjusted sample size"), n,
                  sprintf("\n %-24s: ", "Achieved power"), sprintf("%.4f", pwr))
  }
  txt    <- paste0(txt, sprintf(f, "Randomized"), format(tmp$date, usetz = TRUE),
                   sprintf(f, "Seed"), tmp$seed)
  ibd    <- names(rand)[4:ncol(rand)]
  names(rand)[4:ncol(rand)] <- gsub("\\.", " ", ibd)
  row.names(ibd.seq)        <- names(rand)[4:ncol(rand)]
  names(ibd.seq)            <- rep(c("seq", "n"), 2)
  et     <- signif(1000 * (proc.time()[[3]] - st), 3)
  if (details) {          # Only if requested.
    ifelse (et < 10, rt <- "<10 ms", rt <- paste(et, "ms"))
    txt <- paste0(txt, sprintf(f, "Runtime"), rt, " (", calls, " call")
    if (calls > 1) txt <- paste0(txt, "s")
    txt <- paste(txt, "for balance of IBDs)")
  }
  if (print) {            # To the console.
    print(rand, row.names = FALSE)
    txt <- paste0(txt, "\n")
    cat(txt)
  }
  # Results as a list.
  return(invisible(list(rand = rand, n = n, sequences = seqs,
                        n.seq = tmp$ninseqs, info = info,
                        ibds = ibd.seq, seed = tmp$seed,
                        calls = calls, et = et,
                        txt = substr(txt, 2, nchar(txt)))))
}