‘Two-at-a-Time’ approach; ≥3 treatments: new 🇷 script [Design Issues]
❝ 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)
library(randomizeBE)
RL4(n = 20, seqs = c("ADCBEF", "BFDEAC", "CAEDFB", "DBAFCE", "ECFABD", "FEBCDA"),
seed = 820697)
Warning message:
Unbalanced design! # of subj. in sequences ADCBEF BFDEAC CAEDFB DBAFCE ECFABD FEBCDA: 3 4 4 3 3 3
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
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)
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)
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 dropout-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)))))
}
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:
- ‘Two-at-a-Time’ approach with >3 treatments Helmut 2022-06-12 11:37 [Design Issues]
- ‘Two-at-a-Time’ approach with >3 treatments Shuanghe 2022-06-12 19:02
- ‘Two-at-a-Time’ approach with >3 treatments Helmut 2022-06-12 20:10
- ‘Two-at-a-Time’ approach with >3 treatments Shuanghe 2022-07-10 00:05
- ‘Two-at-a-Time’ approach; ≥3 treatments: new 🇷 scriptHelmut 2022-07-10 16:02
- ‘Two-at-a-Time’ approach; ≥3 treatments: new 🇷 script Shuanghe 2022-07-13 15:35
- You’re a NE🇷D! Helmut 2022-07-13 16:46
- ‘Two-at-a-Time’ approach; ≥3 treatments: new 🇷 script Shuanghe 2022-07-13 15:35
- ‘Two-at-a-Time’ approach; ≥3 treatments: new 🇷 scriptHelmut 2022-07-10 16:02
- ‘Two-at-a-Time’ approach with >3 treatments Shuanghe 2022-07-10 00:05
- ‘Two-at-a-Time’ approach with >3 treatments Helmut 2022-06-12 20:10
- ‘Two-at-a-Time’ approach with >3 treatments Shuanghe 2022-06-12 19:02