Bioequivalence and Bioavailability Forum

Dear all,

after some preliminary sim’s I think the most challenging part is not noise, but coping with lee’s greedy behavior – fitting two lines more often than necessary. If we feed lee’s greed we would loose many data points in monophasic profiles (compared to TTT and even more to R²_adj). So the problem is: Dampen lee’s sensitivity but keep it high enough in order to pick true multiphasic profiles. New code (as a function for easier testing):

require(PK)
lambdaz <- function(time, conc, damp=NA, hl=NA, print=FALSE) {
  data0<- data.frame(time=time, conc=conc)
  TTT  <- time[min(seq_along(data0$time)[data0$time>=data0$time[which.max(data0$conc)]*2])]
                 # Ace’s ingenious one-liner (THX!)
                 # see http://forum.bebac.at/forum_entry.php?id=2548
  data1<- subset(data0, data0$time >= TTT)
                 # Post-absorption data only
  ET   <- max(time) # t-last
  meth <- c("ols", "lad", "hub", "npr")
                 # Methods of lee()
  ow   <- options("warn"); options(warn=-1)
  BP   <- NULL; BPT <- NULL; lz <- NULL; t12 <- NULL
  cat("\nFunction lee(), method:\n")
  for(j in 1:4) {# Run all methods
    res <- lee(data1$conc, data1$time, method=meth[j])
    t12.ratio <- res$parms[1, 2]/res$parms[1, 1]
    if(print) cat("Half live ratio (slow/fast):", round(t12.ratio, 2), "\n")
    ifelse(!is.na(damp) & t12.ratio < damp, esc <- "too close", esc <- "")
    if(is.na(res$chgpt) | esc == "too close") {
      # no distinct or too similar phases: fallback / force to TTT
      fit <- lm(log(data1$conc[data1$time >= TTT & data1$time <= ET])
                ~ data1$time[data1$time >= TTT & data1$time <= ET])
      BP <- c(BP, NA); BPT <- c(BPT, NA) # for completeness
      lz  <- c(lz, abs(as.numeric(fit$coeff[2])))
      t12 <- c(t12, log(2)/lz[j])
      if(is.na(res$chgpt)) {
        cat(sprintf("%s%s%s ", "\u201c", meth[j], "\u201d monophasic         (=> TTT):"))
      } else {
        cat(sprintf("%s%s%s ", "\u201c", meth[j], "\u201d two similar phases (=> TTT):"))
      }
      cat(sprintf("%5.2f%s%5.2f %s%2i%s %s %6.5f %s %5.2f",
        TTT, "\u2013", ET,
        "(n=", length(time[time >= TTT]), ")",
        "\u03bbz:", lz[j], "(t\u00bd", t12[j]))
    } else {
      BP  <- c(BP, as.numeric(res$chgpt))     # extract breakpoint
      BPT <- c(BPT, min(time[time >= BP[j]])) # time point which is at least BP
      fit <- lm(log(data1$conc[data1$time >= BPT[j] & data1$time <= ET])
                ~ data1$time[data1$time >= BPT[j] & data1$time <= ET])
      lz  <- c(lz, abs(as.numeric(fit$coeff[2])))
      t12 <- c(t12, log(2)/lz[j])
      cat(sprintf("%s%s%s %5.2f, %s %5.2f%s%5.2f %s%2i%s %s %6.5f %s %5.2f",
        "\u201c", meth[j], "\u201d breakpoint:", BP[j],
        "proposed:", BPT[j], "\u2013", ET,
        "(n=", length(time[time >= BPT[j]]), ")",
        "\u03bbz:", lz[j], "(t\u00bd", t12[j]))
    }
    if(is.na(hl)) {
      cat(sprintf("%s", ")\n"))
    } else {
      cat(sprintf(", %s %+4.2f%%%s", "RE", 100*(t12[j]-hl)/hl, ")\n"))
    }
  }
  options(ow) # reset warnings
}

Parameters:

• damp defaults to NA. No dampening (also with damp=1) – allows lee() to fit two lines ad libitum. Higher values increase the threshold (i.e., ratios of half lives < damp will be forced to the simple TTT-method). Values between 1.5 and 2.5 are some reasonable starters.
  
• hl defaults to NA. If you support a known half life (e.g., of simulated data), bias (as %RE) will be calculated.
  
• print defaults to FALSE. Set to true to see half lives of the segments.

For the one-compartment data set of my last post I get with
lambdaz(time, conc, damp=1.5, hl=24, print=F)
Function lee(), method:
“ols” two similar phases (=> TTT):  7.00–72.00 (n=10) λz: 0.02680 (t½ 25.87, RE +7.78%)
“lad” two similar phases (=> TTT):  7.00–72.00 (n=10) λz: 0.02680 (t½ 25.87, RE +7.78%)
“hub” monophasic         (=> TTT):  7.00–72.00 (n=10) λz: 0.02680 (t½ 25.87, RE +7.78%)
“npr” two similar phases (=> TTT):  7.00–72.00 (n=10) λz: 0.02680 (t½ 25.87, RE +7.78%)
The ratios of half lives (with the exception of “hub”, which detected only one phase) are 1.17, 1.20, 1.22. We know from modeling that such small differences are simply impossible to detect with the full arsenal of PK software. In other words: lee’s ols/lad/npr-methods are crap. So a threshold of 1.5 was already very optimistic. Now for the good news: With the same threshold lee() would still detect two phases in the other two data sets.

Of course for a user this is not the way to go. In simulations we should try to find out a value which reliably distinguishes between mono- and multiphasic profiles. This could be the default in your future ‘semi-automatic’ method.

If you want to give it a try: Run lambdaz() and use the code from this post to simulate profiles. Change the relevant lines to:
n      <- 1       # no of sims
AErr   <- 0.15    # analytical error CV 15%, lognormal
curve(F.d*D/V.d*k01.d/(k01.d-k10.d)*(exp(-k10.d*(x-tlag.d))-exp(-k01.d*(x-tlag.d))),
      from=tlag.d, to=24, n=(24-tlag.d)*5+1, type="l", lwd=2, col="red", xlim=c(0,24),
      ylim=c(1,350), log="y", xlab="time", ylab="concentration", add=FALSE)
  for (j in 1:length(t)){
    if (j==1)
      AErr1 <- rtruncnorm(n=1, a=-3*abs(C[j]*AErr), b=+3*abs(C[j]*AErr),
                          mean=0, sd=abs(C[j]*AErr))
    else
      AErr1 <- c(AErr1, rtruncnorm(n=1, a=-3*abs(C[j]*AErr), b=+3*abs(C[j]*AErr),
                                   mean=0, sd=abs(C[j]*AErr)))
  }
Call lambdaz(t, C, damp=2, hl=log(2)/k10, print=FALSE)
Worked nicely with low bias until I increased the analytical error to 25%. Quite often I was kicked by:
Error in lee(data1$conc, data1$time, method = meth[j]) :
  not enough points for inital phase
Ah, the lag-time! Sometimes t_max is at 2.5 and TTT therefore at 5. Leaves only five data point for two phases (or just even four if the LLOQ hits). Since absorption starts at the lagtime, in such cases the subset should not start at 2×t_max, but rather at 2×t_max-t_lag. Maybe tomorrow…