Bioequivalence and Bioavailability Forum

Hi PVRC,

❝ Anybody have the reference datasets for four treatment, four sequence, four period design and the R code.

I don’t. See this article first.

An -script for the 4×4 and 6×3 Williams’s designs, evaluation by the ‘Two at at Time’ approach (incomplete block designs preferred by the EMA). No period effects, no carryover, balanced sequences (sorry, 200+ lines of code):
make.williams <- function(treatments = 3,
                          sequences = 6,
                          seqs, n) {
  # returns a data.frame of 3- or 4-treatment
  # Williams' designs
  # seq: vector of sequences (optional)
  # n  : number of subjects
  if (treatments < 3 | treatments > 4)
    stop("Only 3 or 4 treatments implemented.")
  if (!sequences %in% c(4, 6))
    stop("Only 4 or 6 sequences implemented.")
  if (treatments == 3 & !sequences == 6)
    stop("6 sequences required for 3 treatments.")
  if (treatments == 4 & !sequences == 4)
    stop("4 sequences required for 4 treatments.")
  if (n %% sequences != 0) {
    msg <- paste0("\n", n, " is not a multiple of ",
                  sequences, ".",
                 "\nOnly balanced sequences implemented.")
    stop(msg)
  }
  # defaults if seqs not given
  if (is.null(seqs)) {
    if (treatments == 3) {
      seqs <- c("ABC", "ACB", "BAC",
                "BCA", "CAB", "CBA")
    } else {
      seqs <- c("ADBC", "BACD",
                "CBDA", "DCAB")
    }
  }
  if (!length(seqs) == sequences)
    stop("number of elements of 'seqs' have to equal 'sequences'.")
  if (!typeof(seqs) == "character")
    stop("'seqs' has to be a character vector.")
  subject  <- rep(1:n, each = treatments)
  sequence <- rep(seqs, each = n / length(seqs) * treatments)
  data     <- data.frame(subject  = subject,
                         period   = 1:treatments,
                         sequence = sequence)
  for (j in 1:nrow(data)) {
    data$treatment[j] <- strsplit(data$sequence[j],
                                  "")[[1]][data$period[j]]
  }
  return(data)
}

sim.williams <- function(alpha = 0.05, treatments = 3, sequences = 6, n,
                         T, R, theta0, theta1 = 0.80, theta2 = 1.25,
                         CV, seqs, nested = FALSE, show.data = FALSE,
                         details = FALSE, setseed = TRUE) {
  # 3-trt Williams' design (comparisons = 2)
  #   A vs B and A vs C: (A = T , B = R1, C = R2)
  #   A vs C and B vs C: (A = T1, B = T2, C = R)
  # 4-trt Williams' design (comparisons = 2)
  #   A vs B and C vs D: (A = Tfed   , B = Rfed
  #                       C = Tfasted, D = Rfasted)
  # 4-trt Williams' design (comparisons = 3)
  # dose proportionality (A = reference dose level)
  #   B vs A and C vs A and D vs A
  # T: vector of T-codes (character, e.g., c("A", "C")
  # R: vector of R-codes (character, e.g., c("B", "D")
  # theta0: vector of T/R-ratios
  # CV: vector of CVs of the T/R-comparisons
  if (treatments < 3)
    stop("At least 3 treatments required.")
  if (treatments > 4)
    stop("> 4 treatments not implemented.")
  comp <- paste0(T, R)
  comparisons <- length(comp)
  if (treatments == 3 & comparisons > 2)
    stop("Only 2 comparisons for 3 treatments implemented.")
  if (!length(theta0) == comparisons)
    stop("theta0 has to be a vector with ", comparisons,
         " elements.")
  if (!length(CV) == comparisons)
    stop("CV has to be a vector with ", comparisons,
         " elements.")
  # vector of ordered treatment codes
  trts <- c(T, R)[order(c(T, R))]
  heading.aovIII <- c(paste("Type III ANOVA:",
                            "Two at a Time approach"),
                             "sequence vs")
  # data.frame with columns subject, period, sequence, treatment
  data <- make.williams(treatments = treatments,
                        sequences = sequences,
                        seq = NULL, n = n)
  if (setseed) set.seed(1234567)
  sd   <- sqrt(log(CV^2 + 1))
  Tsim <- data.frame(subject = 1:n,
                     treatment = rep(T, each = n),
                     PK = NA)
  cuts <- seq(n, nrow(Tsim), n)
  k    <- 1
  for (j in seq_along(cuts)) {
    Tsim$PK[k:cuts[j]] <- exp(rnorm(n = n,
                              mean = log(theta0[j]),
                              sd = sd[j]))
    k <- k + n
  }
  Rsim <- data.frame(subject = 1:n,
                     treatment = rep(R, each = n),
                     PK = NA)
  cuts <- seq(n, nrow(Rsim), n)
  k    <- 1
  for (j in seq_along(cuts)) {
    Rsim$PK[k:cuts[j]] <- exp(rnorm(n = n, mean = 0,
                                    sd = sd[j]))
    k <- k + n
  }
  TRsim      <- rbind(Tsim, Rsim)
  data       <- merge(data, TRsim,
                      by.x = c("subject", "treatment"),
                      by.y = c("subject", "treatment"),
                      sort = FALSE)
  cols       <- c("subject", "period",
                  "sequence", "treatment")
  data[cols] <- lapply(data[cols], factor)
  seqs       <- unique(data$sequence)
  if (show.data) print(data, row.names = FALSE)
  # named vectors of geometric means
  per.mean <- setNames(rep(NA, treatments), 1:treatments)
  trt.mean <- setNames(rep(NA, treatments), trts)
  seq.mean <- setNames(rep(NA, length(seqs)), seqs)
  for (j in 1:treatments) {
    per.mean[j] <- exp(mean(log(data$PK[data$period == j])))
    trt.mean[j] <- exp(mean(log(data$PK[data$treatment == trts[j]])))
  }
  for (j in 1:sequences) {
    seq.mean[j] <- exp(mean(log(data$PK[data$sequence == seqs[j]])))
  }
  txt <- "Geometric means of PK response\n  period"
  for (j in seq_along(per.mean)) {
    txt <- paste0(txt, "\n    ", j,
                  paste(sprintf(": %.4f", per.mean[j])))
  }
  txt <- paste(txt, "\n  sequence")
  for (j in seq_along(seq.mean)) {
    txt <- paste0(txt, "\n    ", seqs[j],
                  paste(sprintf(": %.4f", seq.mean[j])))
  }
  txt <- paste(txt, "\n  treatment")
  for (j in seq_along(trt.mean)) {
    txt <- paste0(txt, "\n    ", trts[j],
                  paste(sprintf(": %.4f", trt.mean[j])))
  }
  if (details) cat(txt, "\n")
  for (j in 1:comparisons) { # IBD comparisons
    comp.trt      <- strsplit(comp[j], "")[[1]]
    TaT           <- data[data$treatment %in% comp.trt, ]
    TaT$treatment <- as.character(TaT$treatment)
    TaT$treatment[TaT$treatment == comp.trt[1]] <- "T"
    TaT$treatment[TaT$treatment == comp.trt[2]] <- "R"
    TaT$treatment <- as.factor(TaT$treatment)
    if (nested) { # bogus
      model <- lm(log(PK) ~ sequence + subject%in%sequence +
                            period + treatment,
                            data = TaT)
    } else {      # simple for unique subject codes
      model <- lm(log(PK) ~ sequence + subject +
                            period + treatment,
                            data = TaT)
    }
    TR.est <- exp(coef(model)[["treatmentT"]])
    TR.CI  <- as.numeric(exp(confint(model, "treatmentT",
                                     level = 1 - 2 * alpha)))
    m.form <- toString(model$call)
    m.form <- paste0(substr(m.form, 5, nchar(m.form)),
                     " (", paste(comp.trt, collapse=""), ")")
    aovIII <- anova(model)
    if (nested) {
      MSdenom <- aovIII["sequence:subject", "Mean Sq"]
      df2     <- aovIII["sequence:subject", "Df"]
    } else {
      MSdenom <- aovIII["subject", "Mean Sq"]
      df2     <- aovIII["subject", "Df"]
    }
    fvalue <- aovIII["sequence", "Mean Sq"] / MSdenom
    df1    <- aovIII["sequence", "Df"]
    aovIII["sequence", 4]      <- fvalue
    aovIII["sequence", 5]      <- pf(fvalue, df1, df2,
                                     lower.tail = FALSE)
    attr(aovIII, "heading")[1]   <- m.form
    attr(aovIII, "heading")[2:3] <- heading.aovIII
    if (nested) {
      attr(aovIII, "heading")[3] <- paste(heading.aovIII[2],
                                          "sequence:subject")
    } else {
      attr(aovIII, "heading")[3] <- paste(heading.aovIII[2],
                                          "subject")
    }
    CV.TR.est <- sqrt(exp(aovIII["Residuals", "Mean Sq"])-1)
    if (details) {
      info <- paste0("\nIncomplete blocks extracted from ",
                     length(seqs), "\u00D7", treatments,
                     " Williams\u2019 design\n")
      cat(info); print(aovIII, digits = 4, signif.legend = FALSE)
    }
    txt <- paste0("\nPE ",
                  comp.trt[1], "/", comp.trt[2],
                  "     ",
                  sprintf("%6.4f", TR.est),"\n",
                  sprintf("%.f%% %s", 100*(1-2*alpha),
                  "CI     "),
                   paste(sprintf("%6.4f", TR.CI),
                         collapse = " \u2013 "))
    if (round(TR.CI[1], 4) >= theta1 &
        round(TR.CI[2], 4) <= theta2) {
      txt <- paste(txt, "(pass)")
    } else {
      txt <- paste(txt, "(fail)")
    }
    txt <- paste(txt, "\nCV (intra)",
                 sprintf("%6.4f", CV.TR.est), "\n")
    cat(txt)
  }
}

Say you have a 4×4 design and are interested in the comparisons A/B and C/D. We simulate a data set with different T/R-ratios and CVs.

sim.williams(treatments = 4, sequences = 4, n = 28,
             T = c("A", "C"), R = c("B", "D"),
             theta0 = c(0.95, 0.90),
             CV = c(0.25, 0.20), details = TRUE)

Gives
Geometric means of PK response
  period
    1: 0.9153
    2: 0.9901
    3: 0.9372
    4: 0.9407
  sequence
    ADBC: 0.9375
    BACD: 0.8875
    CBDA: 0.9705
    DCAB: 0.9894
  treatment
    A: 0.9423
    B: 0.9740
    C: 0.8859
    D: 0.9825

Incomplete blocks extracted from 4×4 Williams’ design
log(PK) ~ sequence + subject + period + treatment, TaT (AB)
Type III ANOVA: Two at a Time approach
sequence vs subject
          Df Sum Sq Mean Sq F value Pr(>F)
sequence   3 0.1263 0.04209   1.212  0.327
subject   24 0.8338 0.03474   0.450  0.972
period     3 0.2614 0.08714   1.129  0.357
treatment  1 0.0153 0.01532   0.198  0.660
Residuals 24 1.8531 0.07721               

PE A/B     0.9675
90% CI     0.8520 – 1.0985 (pass)
CV (intra) 0.2833

Incomplete blocks extracted from 4×4 Williams’ design
log(PK) ~ sequence + subject + period + treatment, TaT (CD)
Type III ANOVA: Two at a Time approach
sequence vs subject
          Df Sum Sq Mean Sq F value Pr(>F) 
sequence   3 0.1644 0.05480   1.696 0.1945 
subject   24 0.7754 0.03231   0.920 0.5799 
period     3 0.0544 0.01812   0.516 0.6751 
treatment  1 0.1500 0.15000   4.272 0.0497 *
Residuals 24 0.8427 0.03511                 

PE C/D     0.9017
90% CI     0.8276 – 0.9823 (pass)
CV (intra) 0.1890

If you want only the results use details = FALSE. To retrieve the simulated data, use show.data = TRUE.

This is just one of the six possible 4×4 Williams’ designs (ADBC|BACD|CBDA|DCAB, ADCB|BCDA|CABD|DBAC, ACDB|BDCA|CBAD|DABC, ACBD|BADC|CDAB|DBCA, ABDC|BCAD|CDBA|DACB, ABCD|BDAC|CADB|DCBA). If you are interested in another one, specify it as a character vector in the argument seq.
If you insist in the bogus nested structure (why?), use nested = TRUE.
If you want other simulations use setseed = FALSE.

Dose proportionality (4 levels, dose-normalized PK-responses A/D, B/D, C/D).
sim.williams(treatments = 4, sequences = 4, n = 28,
             T = c("A", "B", "C"), R = "D",
             theta0 = c(0.93, 0.97, 0.99),
             CV = c(0.20, 0.22, 0.25))

PE A/D     0.9404
90% CI     0.8695 – 1.0171 (pass)
CV (intra) 0.1727

PE B/D     0.9703
90% CI     0.8888 – 1.0593 (pass)
CV (intra) 0.1940

PE C/D     0.9815
90% CI     0.8870 – 1.0860 (pass)
CV (intra) 0.2241

6×3 Williams’ design, one test (A) compared to references from different regions (B, C).
sim.williams(treatments = 3, sequences = 6, n = 24,
             T = c("A"), R = c("B", "C"),
             theta0 = c(0.93, 0.95),
             CV = c(0.25, 0.20))

PE A/B     0.9793
90% CI     0.8866 – 1.0817 (pass)
CV (intra) 0.2022

PE A/C     0.9351
90% CI     0.8457 – 1.0340 (pass)
CV (intra) 0.2045

Results of the three examples confirmed in Phoenix/WinNonlin v8.1.