THX! Another attempt… [RSABE / ABEL]

posted by Helmut Homepage – Vienna, Austria, 2023-07-08 14:15 (143 d 18:13 ago) – Posting: # 23664
Views: 1,147

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
22,811 posts in 4,783 threads, 1,642 registered users;
16 visitors (0 registered, 16 guests [including 6 identified bots]).
Forum time: 07:28 CET (Europe/Vienna)

Every man gets a narrower and narrower field of knowledge
in which he must be an expert in order to compete with other people.
The specialist knows more and more about less and less
and finally knows everything about nothing.    Konrad Lorenz

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