THX! Another attempt… [RSABE / ABEL]

posted by Helmut Homepage – Vienna, Austria, 2023-07-08 14:15 (373 d 16:37 ago) – Posting: # 23664
Views: 2,184

Hi mittyri,

❝ I suggest to substitute the workaround to dinosaur method:

❝ […]

❝ at least it gives exactly the same results as classical one.


THX, that’s what’s given in the formula of the guidance. Avoing its overhead, more than 10times faster than lm():-D

❝ Yes, ABE is still under fire, the best approach you've found is here


Of course, this model is not exactly the one of the FDA. Seemingly (always?) a bit more conservative. From a regulatory perspective this should be acceptable. But what if the FDA’s passes and the other fails? Sponsors will not like that.

New [image]-script at the end. Allows also testing with local files. Below always RSABE followed by ABE. Results are similar but not the same.

❝ Sorry, right now don't have a clue regarding differences with Winnonlin in RSABE part


So do I… We cross-validated it, right?

It’s funny that the FDA’s draft guidance of December 2022 states:

For example, a three-period design [TRT|RTR], could be used. A fully replicated design can estimate the subject-by-formulation interaction variance components.

True, but not with the FDA’s code.

PS @Achievwin: You keep us busy. ;-)


get.data <- function(file) { # import data and attach libraries
  require(nlme)    # for the mixed-effects model of ABE
  require(emmeans) # for its PE and CI
  # must be a flat file with column separator ',' and decimal separator '.'
  # arbitrary column names with the order below; untransformed PK-response

  data        <- read.csv(file = file,
                          stringsAsFactors = FALSE) # for R <4.0.0
  names(data) <- c("subj", "per", "seq", "treat", "pk")
  nseq        <- length(unique(data$seq))
  nper        <- length(unique(data$per))
  if (nper == 2 & nseq == 2) stop ("Simple crossover not supported!")
  if (nper == 3 & nseq == 3) stop ("Partial replicate design not supported!")
  if (nper < 3 & !nseq == 2) stop ("Design not supported!")
  if (nper == 4) ds <- "repl4"
  if (nper == 3) ds <- "repl3"
  res <- list(data = data, ds = ds)
  return(res)
}
###########################################################
# The FDA’s RSABE and ABE in full replicate designs with  #
# two sequences (a desperate attempt)                     #
# ------------------------------------------------------- #
# Detlew’s original code                                  #
# https://forum.bebac.at/forum_entry.php?id=15490         #
# Michael’s direct implementation of dlat                 #
# https://forum.bebac.at/forum_entry.php?id=23663         #
# ABE with the mixed-effects model with unequal variances #
# https://forum.bebac.at/forum_entry.php?id=20503         #
###########################################################
# select one of the examples or a local file

file      <- "https://bebac.at/downloads/ds01.csv" # EMA data set I (TRTR|RTRT)
file      <- "https://bebac.at/downloads/ds03.csv" # period 3 of I removed (TRT|RTR)
res       <- get.data(file)
data      <- res$data
ds        <- res$ds
# sort by subject, treatment and period to get replicate number
# this works in contrast to the SAS code also for
# other sequences than TRTR|RTRT

data      <- data[order(data$subj, data$treat, data$per), ]
data$repl <- 1
for (i in c(2:nrow(data))) {
  if (data$subj[i] != data$subj[i-1] | data$treat[i] != data$treat[i-1]) {
    data$repl[i] <- 1
  } else {
    data$repl[i] <- data$repl[i-1] + 1
  }
}
# make a column with treatment and repl, i.e T1,T2, R1,R2
data$code2 <- paste0(data$treat, data$repl)
data       <- data[order(data$subj, data$per), ]
# reshape to wide with colums subj, seq, pk.R1, pk.R2, pk.T1, pk.T2
# and log-transform pk

data_r     <- reshape(data = data, direction = "wide", v.names = "pk",
                      idvar = c("subj", "seq"), timevar = "code2",
                      drop = c("per", "treat", "repl"))
data_r[, 3:6]      <- log(data_r[, 3:6])
names(data_r)[3:6] <- sub("pk.", "logpk.", names(data_r))[3:6]
# T-R analysis, ilat in the SAS code
# modified for TRT|RTR and TTR|RRT

data_r$TR <- NA
if (ds == "repl4") {          # 4-period full replicate
  data_r$TR <- (data_r$logpk.T1 + data_r$logpk.T2 -
                data_r$logpk.R1 - data_r$logpk.R2) / 2
} else {                      # 3-period full replicate
  for (i in 1:nrow(data_r)) {
    if (data_r$seq[i] == "RTR" | data_r$seq[i] == "RRT")
      data_r$TR[i] <- data_r$logpk.T1[i] -
                      (data_r$logpk.R1[i] + data_r$logpk.R2[i]) / 2
    if (data_r$seq[i] == "TRT" | data_r$seq[i] == "TTR")
      data_r$TR[i] <- (data_r$logpk.T1[i] + data_r$logpk.T2[i]) / 2 -
                       data_r$logpk.R1[i]   }
}
# get rid of subjects with missing periods (TR == NA)
data_r    <- data_r[complete.cases(data_r$TR), ]
# ilat analysis, ANOVA with seq as effect
data_r$seq <- as.factor(data_r$seq)
# with standard contr.treatment we get not the wanted intercept!
# with that the intercept is intercept+seq1

oc        <- options("contrasts") # save options
options(contrasts = c("contr.sum", "contr.poly"))
m1        <- lm(TR ~ seq, data = data_r)
# intercept is estimate of µt-µR
alpha     <- 0.05 # level of the tests
est       <- coef(m1)[[1]]
pe1       <- exp(est)
ci1       <- as.numeric(confint(m1, "(Intercept)", level = 1 - 2 * alpha))
dfTR      <- m1$df.residual
# stderr, "the unknown x" for unbiased estimate of (µT-µR)^2
stderr    <- summary(m1)$coefficients["(Intercept)", "Std. Error"]
# linearized scABE criterion component x
x         <- est^2 - stderr^2
boundx    <- max(abs(ci1))^2
ci1       <- exp(ci1)
data_r$RR <- data_r$logpk.R1 - data_r$logpk.R2
data_r$TT <- data_r$logpk.T1 - data_r$logpk.T2
# get rid of subjects where R or T is not repeated (RR == NA or TT = NA)
data_rr   <- data_r[complete.cases(data_r$RR), c(1:2, 8)]
data_tt   <- data_r[complete.cases(data_r$TT), c(1:2, 9)]
# dlat analysis acc. to the formula in the guidance; THX to Michael!
# lm() with only one sequence in 3-period full replicates not possible
# gives the same result as lm() for 4-period full replicate
# more than 10times faster than lm(RR ~ seq, data = data_rr)

seqs      <- unique(as.character(data_rr$seq))
n         <- nrow(data_rr)
m         <- length(seqs)
denom     <- 2 * (n - m)
numer     <- 0
for (i in seqs) {
  D       <- data_rr$RR[data_rr$seq == i]
  numer   <- numer + sum((D - mean(D, na.rm = TRUE))^2)
}
dfRR      <- n - m
s2wR      <- numer / denom
# the same for T-T (informative!)
seqs      <- unique(as.character(data_tt$seq))
n         <- nrow(data_tt)
m         <- length(seqs)
denom     <- 2 * (n - m)
numer     <- 0
for (i in seqs) {
  D       <- data_tt$TT[data_tt$seq == i]
  numer   <- numer + sum((D - mean(D, na.rm = TRUE))^2)
}
dfTT      <- n - m
s2wT      <- numer / denom
CVwR1     <- sqrt(exp(s2wR) - 1)
CVwT1     <- sqrt(exp(s2wT) - 1)
# regulatory setting for HVD(P)s
theta     <- ((log(1.25)) / 0.25)^2
# linearized scABE criterion component y
y         <- -theta * s2wR
boundy    <- y * dfRR / qchisq(1 - alpha, dfRR)
swR       <- sqrt(s2wR)
# linearized scABE criterion, 95% upper bound
critbound <- (x +y ) + sqrt(((boundx - x)^2) + ((boundy - y)^2))
# ABE: evaluation by lme() of 'nlme' and confint() of 'emmeans'
# Satterthwaite df

fac       <- c("subj", "per", "seq", "treat")
data[fac] <- lapply(data[fac], factor)
m2        <- lme(log(pk) ~ seq + per + treat,
                           random = ~ treat | subj,
                           data = data, weights = varIdent(~ treat),
                           method = "REML", na.action = na.exclude,
                           control = list(opt = "optim"))
res       <- confint(emmeans(m2, pairwise ~ treat, mode = "satterth"),
                     level = 1 - 2 * alpha)
# note: order reversed because none of relevel(), reverse = TRUE,
#       or ref = "R" works as desired (gives always contrast R - T)

pe2       <- exp(-res$contrasts$estimate)
ci2       <- exp(-c(res$contrasts$upper.CL, res$contrasts$lower.CL))
CVwR2     <- sqrt(exp(res$emmeans$SE[1] * 2) - 1)
CVwT2     <- sqrt(exp(res$emmeans$SE[2] * 2) - 1)
options(oc) # restore options
# assess RSABE or ABE dependent on swR
if (swR >= 0.294) { # assess RSABE (critbound <=0 and pe within 0.8–1.25)
  IL      <- setNames(exp(c(-1, +1) * log(1.25) / 0.25 *swR),
                      c("lower", "upper"))
  RSABE   <- "fail" # the null
  if (critbound <= 0 & (pe1 >= 0.8 & pe1 <= 1.25)) RSABE <- "pass"
  if (swR >= 0.294) {
    txt1 <- paste0("(", substr("\\>= 0.294", 2, 9), ": RSABE acceptable)")
  } else {
    txt1 <- "(< 0.294, apply ABE instead)"
  }
  if (pe1 >= 0.8 & pe1 <= 1.25) {
    txt2 <- "  (within 80% \u002D 125%: pass)"
  } else {
    txt2 <- "  (outside 80% \u002D 125%: fail)"
  }
  if (critbound <= 0) {
    txt3 <- " (nonpositive: pass)"
  } else {
    txt3 <- " (positive: fail)"
  }
  cat("\nswR           :", sprintf("%9.5g", swR), txt1,
      "\nCVwR          :", sprintf("%6.2f%%", 100 * CVwR1),
      "\nCVwT          :", sprintf("%6.2f%%", 100 * CVwT1),
      "  (informative only)",
      "\nimplied limits:", sprintf("%6.2f%% \u002D %.2f%%",
                                   100 * IL[["lower"]], 100 * IL[["upper"]]),
      "\nPE            :", sprintf("%6.2f%%", 100 * pe1), txt2,
      "\n90% CI        :", sprintf("%6.2f%% \u002D %.2f%%",
                                   100 * ci1[1], 100 * ci1[2]),
      "(informative only)",
      "\ncritbound     :", sprintf("%+8.5g", critbound), txt3,
      "\naggregate     :", RSABE, "\n")
} else {            # assess ABE (ci within 0.8–1.25)
  fails <- 0
  if (round(ci2[1], 4) < 0.8) fails <- fails + 1
  if (round(ci2[2], 4) > 1.25) fails <- fails + 1
  if (fails == 1) ABE <- "fail      (one CL outside 80.00% \u002D 125.00%)"
  if (fails == 2) ABE <- "fail      (CI not within 80.00% \u002D 125.00%)"
  if (round(ci2[1], 4) > 1.25 | round(ci2[1], 4) < 0.8)
    ABE <- "inequivalent (CI entirely outside 80.00% \u002D 125.00%)"
  if (round(ci2[1], 4) >= 0.8 & round(ci2[2], 4) <= 1.25)
    ABE <- "pass     (CI within 80.00% \u002D 125.00%)"
  if (swR >= 0.294) {
    txt1 <- paste0("(", substr("\\>= 0.294", 2, 9),
                   ": should apply RSABE instead)")
  } else {
    txt1 <- "(< 0.294, ABE applied)"
  }
  cat("\nswR    :", sprintf("%9.5g", swR), txt1,
      "\nCVwR   :", sprintf("%6.2f%%", 100 * CVwR2), "  (informative only)",
      "\nCVwT   :", sprintf("%6.2f%%", 100 * CVwT2), "  (informative only)",
      "\nPE     :", sprintf("%6.2f%%", 100 * pe2),
      "\n90% CI :", sprintf("%6.2f%% \u002D %.2f%%",
                            100 * ci2[1], 100 * ci2[2]),
      "\nABE    :", ABE, "\n")
}


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
23,101 posts in 4,857 threads, 1,644 registered users;
70 visitors (0 registered, 70 guests [including 7 identified bots]).
Forum time: 06:52 CEST (Europe/Vienna)

One can show the following: given any rule, however “fundamental”
or “necessary” for science, there are always circumstances
when it is advisable not only to ignore the rule,
but to adopt its opposite.    Paul Feyerabend

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