Bioequivalence and Bioavailability Forum

❝ I would look at the original formula

❝      \(s_\text{wR}^2=\Large\frac{\sum\limits_{i=1}^{m}\sum_\limits{j=1}^{n_i}\left(D_{ij}-\overline{D}_{i\,\cdot}\right)^2}{2\left(n-m\right)}\)

❝ So it is impossible to get D for `TRT` sequence. Could we just omit one sequence here?

• The 4-period full replicate
  swR       = 0.45168 (>= 0.294, RSABE acceptable)
CVwR      =  47.57%
PE        = 115.46% (pass)
90% CI    = 106.39% - 125.31% (informative only)
critbound = -0.09454 (<= 0, pass)
aggregate = pass
  If I would apply the workaround though not necessary:
  swR       = 0.45874
CVwR      =  48.40%
  In Phoenix/WinNonin (cross-validated with SAS):
  swR       = 0.44645 (>= 0.294, RSABE acceptable)
CVwR      =  46.96%
PE        = 115.46% (pass)
90% CI    = 106.39% - 125.31% (informative only)
critbound = -0.09208 (<= 0, pass)
aggregate = pass
  Why? BTW, my code for ABE is seemingly completely wrong.
  
  
• The 3-period full replicate
  swR       = 0.53372 (>= 0.294, RSABE acceptable)
CVwR      =  57.41%
PE        = 124.52% (pass)
90% CI    = 113.72% - 136.34% (informative only)
critbound = -0.096683 (<= 0, pass)
aggregate = pass
Note: Based on workaround; questionable!

# for Detlew’s original code see
# https://forum.bebac.at/forum_entry.php?id=15490
# libraries for the mixed-effects model of ABE
suppressMessages(library(lmerTest)) # requires also 'lme4' and 'Matrix'
alpha      <- 0.05 # level of the tests
ds         <- "per4full" # EMA data set I (TRTR|RTRT)
ds         <- "per3full" # period 3 of I removed (TRT|RTR)
if (ds == "per4full") {
  f <- "https://bebac.at/downloads/ds01.csv"
} else {
  f <- "https://bebac.at/downloads/ds03.csv"
}
data        <- read.csv(file = f, stringsAsFactors = FALSE)
names(data) <- c("subj", "per", "seq", "treat", "pk")
# 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]
# now T-R analysis, ilat in the SAS code
# modified for TRT|RTR and TTR|RRT
data_r$TR <- NA
if (ds == "per4full") {
  data_r$TR <- (data_r$logpk.T1 + data_r$logpk.T2 -
                data_r$logpk.R1 - data_r$logpk.R2) / 2
} else {
  for (i in 1:nrow(data_r)) { # clumsy
    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]
  }
}
# now 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
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
# now get rid of subjects where R is not repeated (RR == NA)
data_r    <- data_r[complete.cases(data_r$RR), ]
if (ds == "per4full") {
  # now the equivalent of SAS code for R-R, dlat analysis
  m2      <- lm(RR ~ seq, data = data_r)
  dfRR    <- m2$df.residual
  s2wR    <- summary(m2)$sigma^2 / 2
} else {
  # workaround because lm() with only one sequence not possible
  # make a column with intra-subject variances
  for (i in 1:nrow(data_r)) {
    data_r$varR[i] <- var(c(data_r$logpk.R1[i],
                            data_r$logpk.R2[i]), na.rm = TRUE)
  }
  dfRR    <- nrow(data_r) - 2
  s2wR    <- mean(data_r$varR)
}
CVwR      <- sqrt(exp(s2wR) - 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, subjects random: evaluation by lmer() of package lme4
# Satterthwaite df by package lmerTest
fac       <- c("subj", "per", "seq", "treat")
data[fac] <- lapply(data[fac], factor)
data$treat <- relevel(data$treat, "T")
m3        <- lmer(log(pk) ~ seq + per + treat + (1 | subj),
                            na.action = na.omit, data = data)
sumry     <- summary(m3, ddf = "Satterthwaite")
# direct access of "treatT" does not always work (sometimes "treat1")
r         <- 1 + length(unique(data$seq)) - 1 +
             length(unique(data$per)) - 1 +
             length(unique(data$treat)) - 1
pe2       <- exp(sumry$coef[r, "Estimate"])
ci2       <- exp(confint(m3, level = 1 - 2 * alpha, method = "Wald")[r + 2, ])
CVw       <- sqrt(exp(m3@sigma^2) - 1)
options(oc) # restore options
RSABE     <- ABE <- "fail" # the nulls
if (ds == "per3full") {
  note <- "\nNote: Based on workaround; questionable!\n"
} else {
  note <- "\n"
}
# assess RSABE or ABE dependent on swR
if (swR >= 0.294) { # assess RSABE (critbound <=0 and pe within 0.8–1.25)
  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 <- "(pass)" else txt2 <- "(fail)"
  if (critbound <= 0) {
    txt3 <- paste0("(", substr("\\<= 0", 2, 9), ", pass)")
  } else {
    txt3 <- "(> 0, fail)"
  }
  cat("\nswR       =", sprintf("%.5g", swR), txt1,
      "\nCVwR      =", sprintf("%6.2f%%", 100 * CVwR),
      "\nPE        =", sprintf("%6.2f%%", 100 * pe1), txt2,
      "\n90% CI    =", sprintf("%6.2f%% \u002D %.2f%%", 100 * ci1[1], 100 * ci1[2]),
      "(informative only)",
      "\ncritbound =", sprintf("%+.5g", critbound), txt3,
      "\naggregate =", RSABE, note)
} else {            # assess ABE (ci within 0.8–1.25)
  if (round(ci2[1], 4) >= 0.8 & round(ci2[2], 4) <= 1.25) ABE <- "pass"
  if (swR >= 0.294) {
    txt1 <- paste0("(", substr("\\>= 0.294", 2, 9), ", apply RSABE instead)")
  } else {
    txt1 <- "(< 0.294, ABE applied)"
  }
  cat("\nswR    =", sprintf("%.5g", swR), txt1,
      "\nCVw    =", sprintf("%6.2f%%", 100 * CVw),
      "\nPE     =", sprintf("%.2f%%", 100 * pe2),
      "\n90% CI =", sprintf("%.2f%% \u002D %.2f%%", 100 * ci2[1], 100 * ci2[2]),
      "\nABE    =", ABE, "\n")
}