Bioequivalence and Bioavailability Forum

❝ I was always wondering: why the limit is ≤0.0001?

limit ≤ 1e-04
    lambda.z           bias (%)          lz.start          lz.n
 Min.   :0.003178   Min.   :-97.249   Min.   : 3.25   Min.   : 3.000
 1st Qu.:0.107161   1st Qu.: -7.239   1st Qu.: 5.25   1st Qu.: 3.000
 Median :0.118384   Median :  2.475   Median :14.50   Median : 4.000
 3rd Qu.:0.136400   3rd Qu.: 18.070   3rd Qu.:17.25   3rd Qu.:10.000
 Max.   :0.266735   Max.   :130.890   Max.   :17.25   Max.   :13.000
 NA's   :11         NA's   :11        NA's   :11      NA's   :11

limit ≤ 0.01
    lambda.z           bias (%)          lz.start          lz.n
 Min.   :0.006491   Min.   :-94.381   Min.   : 3.75   Min.   : 3.000
 1st Qu.:0.106485   1st Qu.: -7.825   1st Qu.:10.25   1st Qu.: 3.000
 Median :0.119452   Median :  3.400   Median :14.50   Median : 4.000
 3rd Qu.:0.136873   3rd Qu.: 18.479   3rd Qu.:17.25   3rd Qu.: 6.000
 Max.   :0.278439   Max.   :141.022   Max.   :17.25   Max.   :12.000
 NA's   :8          NA's   :8         NA's   :8       NA's   :8

limit ≤ 0
    lambda.z           bias (%)          lz.start          lz.n
 Min.   :0.001576   Min.   :-98.636   Min.   : 3.25   Min.   : 3.000
 1st Qu.:0.107295   1st Qu.: -7.124   1st Qu.: 5.25   1st Qu.: 3.000
 Median :0.118507   Median :  2.582   Median :14.50   Median : 4.000
 3rd Qu.:0.136410   3rd Qu.: 18.078   3rd Qu.:17.25   3rd Qu.:10.000
 Max.   :0.264027   Max.   :128.546   Max.   :17.25   Max.   :13.000
 NA's   :10         NA's   :10        NA's   :10      NA's   :10

❝ Not simple <0? Or maybe ≤0.01?

sim.el <- function(D, f, V, t12.a, t12.e, tlag, t,
                   CV0, limit = 0.0001) {
  one.comp <- function(f, D, V, k01, k10, tlag, t) {
    # one-compartment model, first order absorption
    # and elimination; optional lag time
    if (!isTRUE(all.equal(k01, k10))) { # common: k01 != k10
      C    <- f * D * k01 / (V * (k01 - k10)) *
              (exp(-k10 * (t - tlag)) - exp(-k01 * (t - tlag)))
      tmax <- log(k01 / k10) / (k01 - k10) + tlag
      Cmax <- f * D * k01 / (V * (k01 - k10)) *
              (exp(-k10 * tmax) - exp(-k01 * tmax))
    } else {                            # flip-flop
      k    <- k10
      C    <- f * D / V * k * (t - tlag) * exp(-k * (t - tlag))
      tmax <- 1 / k
      Cmax <- f * D / V * k * tmax * exp(-k * tmax)
    }
    C[C <= 0] <- 0                     # correct negatives due to lag-time
    res <- list(C = C, Cmax = Cmax, tmax = tmax)
    return(res)
  }
  k01    <- log(2) / t12.a   # absorption rate constant
  k10    <- log(2) / t12.e   # elimination rate constant
  C0     <- one.comp(f, D, V, k01, k10, tlag, t)$C # model wo error
  CV     <- CV0 - C0 * 0.005 # noise increases with decreasing C
  varlog <- log(CV^2 + 1)
  C      <- numeric()
  for (j in 1:length(C0)) {
    C[j] <- rlnorm(1, meanlog = log(C0[j]) - 0.5 * varlog[j],
                   sdlog = sqrt(varlog[j]))
  }
  data   <- data.frame(t = t, C = C)
  data   <- data[complete.cases(data), ]    # discard NAs
  data   <- data[data$t > t[C == max(C)], ] # discard tmax and earlier
  lz.end <- tail(data$t, 1)
  tmp    <- tail(data, 3)
  r2     <- a <- b <- numeric()
  m      <- lm(log(C) ~ t, data = tmp)
  a[1]   <- coef(m)[[1]]
  b[1]   <- coef(m)[[2]]
  r2[1]  <- summary(m)$adj.r.squared
  k      <- 1
  for (j in 4:nrow(data)) {
    k         <- k + 1
    tmp       <- tail(data, j)
    m         <- lm(log(C) ~ t, data = tmp)
    a[k]      <- coef(m)[[1]]
    b[k]      <- coef(m)[[2]]
    r2[k]     <- summary(m)$adj.r.squared
    if (r2[k] < r2[k-1] | abs(r2[k] - r2[k-1]) <= limit) break
  }
  loc <- which(r2 == max(r2))
  if (b[loc] >= 0) { # positive slope not meaningful
    intcpt <- lambda.z <- lz.n <- lz.start <- lz.end <- NA
  } else {
    intcpt   <- a[loc]
    lambda.z <- -b[loc]
    lz.start <- tmp$t[2]
    lz.n     <- nrow(tmp) - 1
  }
  res <- data.frame(limit = limit, intcpt = intcpt, lambda.z = lambda.z,
                    lz.start = lz.start, lz.end = lz.end, lz.n = lz.n)
  return(res)
}
sum.simple <- function(x, digits = 4) {
  # nonparametric summary:
  # remove arithmetic mean whilst keeping eventual NAs
  res <- summary(x)
  if (nrow(res) == 6) {
    res <- res[c(1:3, 5:6), ]
  } else {
    res <- res[c(1:3, 5:7), ]
  }
  return(res)
}
D      <- 200L   # dose
f      <- 2/3    # fraction absorbed (BA)
V      <- 3      # volume of distribution
t12.a  <- 0.75   # absorption half life
t12.e  <- 6      # elimination half life
tlag   <- 0.5    # lag time
t      <- c(0, 0.75, 1.25, 2, 2.5, 3.25, 3.75, 4.5, 5.25, 6.25,
            7.25, 8.75, 10.25, 12.25, 14.5, 17.25, 20.25, 24)
CV0    <- 0.20   # maximum CV (at low concentrations)
limit  <- 0.0001 # stopping rule for R2adj
nsims  <- 1e5L   # number of simulations
aggr   <- data.frame()
pb     <- txtProgressBar(0, 1, 0, char = "\u2588", width = NA, style = 3)
for (j in 1:nsims) {
  aggr <- rbind(aggr, sim.el(D, f, V, t12.a, t12.e, tlag, t, CV0, limit))
  setTxtProgressBar(pb, j / nsims)
}
close(pb)
aggr$bias <- 100 * (aggr$lambda.z - log(2) / t12.e) / (log(2) / t12.e)
names(aggr)[7] <- "bias (%)"
cat("limit \u2264", limit, "\n"); sum.simple(aggr[, c(3, 7, 4, 6)])