Bioequivalence and Bioavailability Forum

Dear Daniel!

❝ […] my basic query

Well, most bioanalysts simply ignore this issue. No chromatography data system I know calculates the confidence interval of back-calculated concentrations.

❝ but I tried to find out the same values and I couldn't. Maybe I did something wrong.

Let’s see (coded in for simplicity):
x   <- c(1, 2, 4, 8, 16, 32)
y   <- c(0.032, 0.062, 0.128, 0.280, 0.520, 1.020)
S.y <- c(0.0064, 0.0091, 0.0173, 0.0392, 0.0640, 0.1224)
LR  <- lm(y ~ x, weights = 1/S.y^2)
summary(LR)

Should give:
Call:
lm(formula = y ~ x, weights = 1/S.y^2)
Weighted Residuals:
       1        2        3        4        5        6
 0.10895 -0.21644 -0.07534  0.51095 -0.02040 -0.19588
Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.0013642  0.0019957  -0.684    0.532
x            0.0326669  0.0007375  44.291 1.55e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.3018 on 4 degrees of freedom
Multiple R-squared: 0.998,      Adjusted R-squared: 0.9975
F-statistic:  1962 on 1 and 4 DF,  p-value: 1.554e-06

Now let’s set up a data frame containing the x-values where we want to predict the confidence interval of y. I’ll ask R for 100 values in order to get a ‘smooth’ curve later on.
newdata <- data.frame(x=seq(min(x), max(x), length.out=100))

… and calculate the 95% CI at the 100 locations (new) based on the fitted model (LR).
newdata <- data.frame(x = seq(min(x), max(x), length.out = 100))
pred    <- predict(LR, newdata, interval = "confidence", level = 0.95)

Now we can plot this stuff:
plot(x, y, xlim = c(0, 32), ylim = c(0, 1.1), pch = 19, cex = 1.35,
     col = "blue", las = 1)
grid()
arrows(x0 = x, y0 = y-S.y, x1 = x, y1 = y+S.y, length = 0.025,
       angle = 90, code = 3, col = "blue")
abline(LR, lwd = 2, col = "blue")
lines(newdata$x, pred[, 2], col = "red")
lines(newdata$x, pred[, 3], col = "red")

I’m a little bit too short in time in order to set up the bisection method to extract the back-calculated CI of estimated concentrations.* The bisection method would fiddle around with the prediction range in such a way that the CI equals the observed y. We would start from a given response (y) and try to find the two intersections with the CI (to the ‘left’ with the upper CL and to the ‘right’ with the lower CL). Once we have these two values we can calculate the CI of the concentration as usual. Hint: Type pred in R. You see three columns – denoted fit, lwr, and upr.

A quick estimate of the upper value close to the LLOQ in my example (y=0.045):
lo2   <- data.frame(x = seq(1.2, 1.6, by = 0.01))
pred2 <- as.data.frame(predict(LR, lo2, interval = "confidence", level = 0.95))
loc   <- c(which.min(abs(pred2$upr - 0.045)),
           which.min(abs(pred2$fit - 0.045)),
           which.min(abs(pred2$lwr - 0.045)))

By calling print(pred2[loc, ], digits = 5) we get:
        fit      lwr      upr
10 0.040776 0.036501 0.045051
23 0.045023 0.040804 0.049241
36 0.049269 0.045091 0.053448

The 10^th and 36^th rows of pred2 are pretty close to the intersection we want (at 0.045). Therefore:
print(lo2[range(loc), ])
[1] 1.29 1.55

The 10^th and 36^th values of lo2 are 1.29 and 1.55 – a first approximation of the CI.

Homework:
For the CI of x=1.42 (from y=0.045) I get by the exact method* in Gnumeric 1.288422087 – 1.547226655. What do you get?

---

• See also the end of this post. The intersection of the CI with a given y-value can be solved analytically – but only for a linear model (irrespective of the weighting scheme). The formula is quite lengthy and it took me a lot of paper/brain to come to an end. Since quadratic calibration is also common in analytics, I would only set up a numeric method anyway.

---

“After-the-dinner-quickshot” for the CI of x with a response of 0.045 (requires the fitted model LR from above):
prec  <- 1e-11                      # target precision
y.obs <- 0.045                      # response
x.est <- (y.obs - coef(LR)[[1]])/coef(LR)[[2]]
for (j in 1:2) {                    # 1: upper CL, 2: lower CL
  if (j == 1) {
    lower <- min(x); upper <- x.est # 1: from LLOQ to x
  } else {
    lower <- x.est; upper <- max(x) # 2: from x to ULOQ
  }
  x.range <- data.frame(x = c(lower, upper))
  CL.est <- predict(LR, x.range, interval="confidence", level=0.95)
  if (j == 1) {
    f.lo <- CL.est[1, 3]; f.hi <- CL.est[2, 3]
  } else {
    f.lo <- CL.est[1, 2]; f.hi <- CL.est[2, 2]
  }
  delta <- upper - lower            # interval
  while (abs(delta) > prec) {       # continue until precision reached
    x.new <- data.frame(x = (lower+upper)/2) # bisection
    if (j == 1) {
      f.new <- predict(LR, x.new, interval="confidence", level=0.95)[1, 3]
    } else {
      f.new <- predict(LR, x.new, interval="confidence", level=0.95)[1, 2]
    }
    if (abs(f.new - y.obs) <= prec) break # precision reached
    if (f.new > y.obs) upper <- x.new     # move limit downwards
    if (f.new < y.obs) lower <- x.new     # move limit upwards
    delta <- upper - lower                # new interval
  }
  ifelse (j == 1, CL.lo <- as.numeric(x.new), CL.hi <- as.numeric(x.new))
}
cat(paste0("response (y): ", y.obs,
           "\nconcent. (x): ", signif(x.est, 5),
           " (95% CI: ", signif(CL.lo, 5), " \u2013 ", signif(CL.hi, 5),
           ")\n"))

Which gives on my machine:
response (y): 0.045
concent. (x): 1.4193 (95% CI: 1.2884 – 1.5472)]

Or alternatively:
x.range <- data.frame(x = c(CL.lo, x.est, CL.hi))
CL.est  <- predict(LR, x.range, interval="confidence", level=0.95)
result  <- cbind(x.range, CL.est)              # add x-values
result  <- result[,names(result)[c(1,4,2,3)]]  # nicer order
print(result, row.names=FALSE)                 # Q.E.D.

Which gives:
      x      upr      fit      lwr
 1.2884 0.045000 0.040725 0.036449
 1.4193 0.049219 0.045000 0.040781
 1.5472 0.053358 0.049179 0.045000

If you still have the plot open:
lines(x=c(-pretty(x)[2], rep(x.est, 2)),
      y=c(rep(y.obs, 2), -pretty(y)[2]), col="blue")
lines(x=rep(CL.lo, 2),
      y=c(coef(LR)[[1]]+coef(LR)[[2]]*CL.lo, -pretty(x)[2]),
      col="red")
lines(x=rep(CL.hi, 2),
      y=c(coef(LR)[[1]]+coef(LR)[[2]]*CL.hi, -pretty(x)[2]),
      col="red")