Log in
RegisterTHX! Another attempt… [RSABE / ABEL]
posted by Helmut 
– Vienna, Austria, 2023-07-08 14:15 (556 d 15:10 ago) – Posting: # 23664
Views: 3,347
❝ 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(). 
❝ 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]](https://static.bebac.at/img/Rlogo_15_12.svg)
-script at the end. Allows also testing with local files. Below always RSABE followed by ABE. Results are similar but not the same.- The 4-period full replicate
swR : 0.45168 (>= 0.294: RSABE acceptable)
CVwR : 47.57%
CVwT : 35.49% (informative only)
implied limits: 66.82% - 149.65%
PE : 115.46% (within 80% - 125%: pass)
90% CI : 106.39% - 125.31% (informative only)
critbound : -0.09454 (nonpositive: pass)
aggregate : pass
Phoenix/WinNonlin
swR : 0.44645
CVwR : 46.96%
PE : 115.46%
90% CI : 106.39% - 125.31%
critbound : -0.092076
swR : 0.45168 (>= 0.294: should apply RSABE instead)
CVwR : 48.09% (informative only)
CVwT : 46.72% (informative only)
PE : 115.80%
90% CI : 107.13% - 125.17%
ABE : fail (one CL outside 80.00% - 125.00%)
Phoenix/WinNonlin
CVwR : 48.09%
CVwT : 35.29%
PE : 115.66%
90% CI : 107.10% - 124.89% (pass!)
- The 3-period full replicate
swR : 0.54127 (>= 0.294: RSABE acceptable)
CVwR : 58.34%
CVwT : 30.67% (informative only)
implied limits: 61.69% - 162.11%
PE : 124.52% (within 80% - 125%: pass)
90% CI : 113.72% - 136.34% (informative only)
critbound : -0.1022 (nonpositive: pass)
aggregate : pass
Phoenix/WinNonlin not possible
swR : 0.54127 (>= 0.294: should apply RSABE instead)
CVwR : 49.63% (informative only)
CVwT : 46.67% (informative only)
PE : 124.38%
90% CI : 113.10% - 136.78%
ABE : fail (one CL outside 80.00% - 125.00%)
Phoenix/WinNonlin
swR : 0.46667
CVwR : 49.33%
CVwT : 31.50%
PE : 124.28%
90% CI : 113.17% - 136.49%
❝ 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]](https://static.bebac.at/pics/Blue_and_yellow_ribbon_UA.png)
Helmut Schütz
![[image]](https://static.bebac.at/img/CC by.png)
The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes
Complete thread:
- any SAS or R-macro for 3 period two sequence RSABEIs Achievwin 2023-07-05 18:38 [RSABE / ABEL]
![[Mix]](https://static.bebac.at/img/mix.png "open in Mix view")
- Use the one for TRTR | RTRT Helmut 2023-07-06 08:46
- RSABE R coded by Detlew mittyri 2023-07-06 14:07
- ‘dlat’ for 3-period full replicate? Helmut 2023-07-07 14:41
- ‘dlat’ for 3-period full replicate? mittyri 2023-07-07 16:04
- Likely wrong workaround? Helmut 2023-07-07 20:06
- a bit corrected workaround mittyri 2023-07-07 22:52
- THX! Another attempt…Helmut 2023-07-08 12:15
- complete.cases for dlat and ilat? mittyri 2023-07-09 21:45
- a bit corrected workaround Achievwin 2023-07-10 20:43
- THX! Another attempt…Helmut 2023-07-08 12:15
- a bit corrected workaround mittyri 2023-07-07 22:52
- Likely wrong workaround? Helmut 2023-07-07 20:06
- ‘dlat’ for 3-period full replicate? mittyri 2023-07-07 16:04
- ‘dlat’ for 3-period full replicate? Helmut 2023-07-07 14:41
- RSABE R coded by Detlew mittyri 2023-07-06 14:07
- Use the one for TRTR | RTRT Helmut 2023-07-06 08:46
Ing. Helmut Schütz