Bioequivalence and Bioavailability Forum

Dear D. Labes!

❝ Save your time for nicer things, f.i. NLYW .

Oh, dear! I'm not sure whether my girlfriend would like that - especially the Y in NLYW.

❝ Set the call to minimizing the weighted sum of squares in the function to:

❝ fit <- optimize(minfun, c(0.1, 10), x=x, y=y, xv=xstd^2, yv=ystd^2, tol=1.e-8) and look what happens .

Great! For the archive:
Call:
deming(x = x, y = y, xstd = c(1, 0), ystd = c(1, 0), jackknife = TRUE)

               Coef    se(coef)            z         p
Intercept 0.1857505 0.140423178 123.97018229 0.0000000
Slope     0.9891061 0.005457953  -0.01334144 0.4893563

Intercept (95% CI): -0.1586471 0.5301481
Slope     (95% CI):  0.9758092 1.002403

Code for the plot with the fit and its CI:

df  <- length(x)-2
t   <- qt(1-0.05/2, df=df)
x1  <- seq(min(x), max(x), length.out=250)
CI  <- matrix(nrow=250, ncol=2, byrow=TRUE,
        dimnames=list(NULL, c("CL.lo", "CL.hi")))
for (j in 1:length(x1)) {
  CI[j, 1] <- (fit$coefficient[[1]] - t*sqrt(fit$variance[1, 1])) +
              (fit$coefficient[[2]] - t*sqrt(fit$variance[2, 2]))*x1[j]
  CI[j, 2] <- (fit$coefficient[[1]] + t*sqrt(fit$variance[1, 1])) +
              (fit$coefficient[[2]] + t*sqrt(fit$variance[2, 2]))*x1[j]
}
plot(x, y, main="method comparison", sub="calibrators", xlab="raw",
  ylab="smooth", cex=2, col="red", cex.sub=0.9)
lines(x=range(x), y=fit$coefficient[[1]]+fit$coefficient[[2]]*range(x),
  col="blue", lwd=2)
lines(x=x1, y=CI[, 1], col="blue")
lines(x=x1, y=CI[, 2], col="blue")

You could also ask whether the smoothed results (y) are within the 95% CI of the fit:
CI <- CI[1:length(x), ] for (j in 1:length(x)) {
  CI[j, 1] <- (fit$coefficient[[1]]-t*sqrt(fit$variance[1,1])) +
              (fit$coefficient[[2]]-t*sqrt(fit$variance[2,2]))*x[j]   CI[j, 2] <- (fit$coefficient[[1]]+t*sqrt(fit$variance[1,1])) +
              (fit$coefficient[[2]]+t*sqrt(fit$variance[2,2]))*x[j] }
CI <- cbind(CI, y)
CI <- as.data.frame(round(CI, 2))
CI <- cbind(CI, CI[, 3] >= CI[, 1] & CI[, 3] <= CI[, 2])
names(CI)[4]<- "within CI?"
print(CI, row.names=FALSE)

  CL.lo  CL.hi      y within CI?
 101.46 104.92 102.97       TRUE
  87.69  90.78  88.58       TRUE
  72.03  74.69  74.37       TRUE
  48.54  50.56  49.79       TRUE
  19.82  21.05  20.13       TRUE
   9.00   9.94   9.58       TRUE
   3.76   4.56   4.21       TRUE
   1.81   2.56   1.96       TRUE