Simulations (lengthy 🇷 script) [🇷 for BE/BA]

posted by Helmut Homepage – Vienna, Austria, 2021-04-23 14:37 (647 d 18:15 ago) – Posting: # 22320
Views: 1,223

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 [image]-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.

Dif-tor heh smusma 🖖🏼 Довге життя Україна! [image]
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes

Complete thread:

UA Flag
Activity
 Admin contact
22,477 posts in 4,708 threads, 1,603 registered users;
10 visitors (0 registered, 10 guests [including 3 identified bots]).
Forum time: 07:52 CET (Europe/Vienna)

The mediocre teacher tells.
The good teacher explains.
The superior teacher demonstrates.
The great teacher inspires.    William Arthur Ward

The Bioequivalence and Bioavailability Forum is hosted by
BEBAC Ing. Helmut Schütz
HTML5