THX! Another attempt… [RSABE / ABEL]

posted by Helmut Homepage – Vienna, Austria, 2023-07-08 14:15 (563 d 03:38 ago) – Posting: # 23664
Views: 3,416

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,363 posts in 4,906 threads, 1,675 registered users;
174 visitors (0 registered, 174 guests [including 12 identified bots]).
Forum time: 16:53 CET (Europe/Vienna)

The combination of some data and an aching desire
for an answer does not ensure that a reasonable answer
can be extracted from a given body of data.    John W. Tukey

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