## ‘Two-at-a-Time’ approach; ≥3 treat­ments: new 🇷 script [Design Issues]

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  20 $sequences  "ACFDEB" "BEDFCA" "CDABFE" "DBCEAF" "EFBADC" "FAECBD"$n.seq  4 4 3 3 3 3 $info  "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()[]   # 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))]),                                length(ibd[ibd == sort(unique(ibd))]))       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()[] - 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))))) }

Dif-tor heh smusma 🖖🏼 Довге життя Україна! Helmut Schütz The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes  Ing. Helmut Schütz 