Example [Bioanalytics]

posted by Helmut Homepage – Vienna, Austria, 2019-03-17 01:56 (766 d 14:07 ago) – Posting: # 20043
Views: 3,609

Hi ElMaestro,

played with an example of a study I have on my desk. Chiral GC/MS, quadratic model, w=1/x2.

ObjF1 <- function(x) {
  w <- 1/Conc^x
  M <- lm(Ratio ~ Conc + I(Conc^2), weights=w)
  return(sum(abs(resid(M)/Conc)))
}
ObjF2 <- function(x) {
  w <- 1/Ratio^x
  M <- lm(Ratio ~ Conc + I(Conc^2), weights=w)
  return(sum(abs(resid(M)/Conc)))
}
IC <- function(m, n) {
  return(list(AIC=signif(extractAIC(m, k=2)[2],5),
              BIC=signif(extractAIC(m, k=log(n))[2]),5))
}
Acc <- function(m, x, y) {
  if (coef(m)[[3]] == 0) stop("panic!")
  if (coef(m)[[3]] < 0 {
    return(100*(-(coef(m)[[2]]/2/coef(m)[[3]] +
                  sqrt((coef(m)[[2]]/2/coef(m)[[3]])^2-
                       (coef(m)[[1]]-y)/coef(m)[[3]])))/x)
  } else {
    return(100*(-(coef(m)[[2]]/2/coef(m)[[3]] -
                  sqrt((coef(m)[[2]]/2/coef(m)[[3]])^2-
                       (coef(m)[[1]]-y)/coef(m)[[3]])))/x)
  }
}
Conc  <- c(0.1, 0.1, 0.3, 0.3, 0.9, 0.9, 2, 2, 6, 6, 12, 12, 24, 24)
Ratio <- c(0.022, 0.024, 0.073, 0.068, 0.193, 0.204, 0.438, 0.433,
           1.374, 1.376, 2.762, 2.732, 5.616, 5.477)
n     <- length(Conc)
w.x1  <- 1/Conc
w.x2  <- 1/Conc^2
x.opt <- optimize(ObjF1,  c(0, 10))$minimum
w.xo  <- 1/Conc^x.opt
w.y1  <- 1/Ratio
w.y2  <- 1/Ratio^2
y.opt <- optimize(ObjF2,  c(0, 10))$minimum
w.yo  <- 1/Ratio^x.opt
dupl  <- sum(duplicated(Conc))
var   <- n/2
for (j in 1:dupl) {
  var[j] <- var(c(Ratio[j], Ratio[j+1]))
}
w.var <- 1/rep(var, each=2)
m.1   <- lm(Ratio ~ Conc + I(Conc^2))
m.2   <- lm(Ratio ~ Conc + I(Conc^2), weights=w.x1)
m.3   <- lm(Ratio ~ Conc + I(Conc^2), weights=w.x2)
m.4   <- lm(Ratio ~ Conc + I(Conc^2), weights=w.xo)
m.5   <- lm(Ratio ~ Conc + I(Conc^2), weights=w.y1)
m.6   <- lm(Ratio ~ Conc + I(Conc^2), weights=w.y2)
m.7   <- lm(Ratio ~ Conc + I(Conc^2), weights=w.yo)
m.8   <- lm(Ratio ~ Conc + I(Conc^2), weights=w.var)
mods  <- c("w=1", "w=1/x", "w=1/x^2", "w=1/x^opt",
           "w=1/y", "w=1/y^2", "w=1/y^opt", "w=1/sd.y^2")
AIC   <- c(IC(m.1, n=n)$AIC, IC(m.2, n=n)$AIC, IC(m.3, n=n)$AIC, IC(m.4, n=n)$AIC,
           IC(m.5, n=n)$AIC, IC(m.6, n=n)$AIC, IC(m.7, n=n)$AIC, IC(m.8, n=n)$AIC)
BIC   <- c(IC(m.1, n=n)$BIC, IC(m.2, n=n)$BIC, IC(m.3, n=n)$BIC, IC(m.4, n=n)$BIC,
           IC(m.5, n=n)$BIC, IC(m.6, n=n)$BIC, IC(m.7, n=n)$BIC, IC(m.8, n=n)$BIC)
res1  <- data.frame(model=mods, exp=signif(c(0:2, x.opt, 1:2, y.opt, NA),5),
                    AIC=signif(AIC,5), BIC=signif(BIC,5))
res2  <- data.frame(Conc=Conc,
                    Acc(m=m.1, x=Conc, y=Ratio), Acc(m=m.2, x=Conc, y=Ratio),
                    Acc(m=m.3, x=Conc, y=Ratio), Acc(m=m.4, x=Conc, y=Ratio),
                    Acc(m=m.5, x=Conc, y=Ratio), Acc(m=m.6, x=Conc, y=Ratio),
                    Acc(m=m.7, x=Conc, y=Ratio), Acc(m=m.8, x=Conc, y=Ratio))
names(res2) <- c("Conc", mods)
cat("\nAkaike & Bayesian Information Critera (smaller is better)\n");print(res1);cat("\nAccuracy (%)\n");print(round(res2, 2), row.names=F)


I got:

Akaike & Bayesian Information Critera (smaller is better)
       model    exp      AIC      BIC
1        w=1 0.0000  -94.099  -92.181
2      w=1/x 1.0000 -127.480 -125.560
3    w=1/x^2 2.0000 -131.720 -129.800

4  w=1/x^opt 1.3355 -132.920 -131.010
5      w=1/y 1.0000 -106.670 -104.750
6    w=1/y^2 2.0000  -90.571  -88.654
7  w=1/y^opt 2.5220 -105.150 -103.230

8 w=1/sd.y^2     NA   62.387   64.304

Accuracy (%)
 Conc    w=1  w=1/x w=1/x^2 w=1/x^opt  w=1/y w=1/y^2 w=1/y^opt w=1/sd.y^2
  0.1 115.66  96.07   94.53     94.63  96.48   94.95     95.02      99.04
  0.1 124.45 104.96  103.49    103.56 105.37  103.93    103.96     107.83
  0.3 113.24 107.57  107.64    107.46 107.74  107.97    107.65     107.71
  0.3 105.92 100.17  100.18    100.02 100.33  100.49    100.21     100.39
  0.9  96.30  95.06   95.52     95.29  95.14   95.77     95.41      94.46
  0.9 101.66 100.48  100.98    100.74 100.56  101.24    100.86      99.83
  2.0  97.07  97.06   97.63     97.40  97.12   97.85     97.49      96.25
  2.0  95.97  95.96   96.51     96.29  96.01   96.74     96.38      95.15
  6.0 100.54 101.06  101.53    101.37 101.10  101.71    101.43     100.29
  6.0 100.69 101.21  101.68    101.51 101.24  101.85    101.58     100.44
 12.0 100.48 100.86  101.08    101.03 100.89  101.18    101.07     100.38
 12.0  99.40  99.78  100.01     99.95  99.81  100.10     99.99      99.30
 24.0 101.23 101.10  100.81    100.98 101.10  100.75    100.98     101.23
 24.0  98.77  98.66   98.40     98.56  98.67   98.35     98.56      98.77


Hey, yours with w=1/x1.3355 is the winner! Duno why the ICs of 1/sy² are that bad. Coding error? The accuracy looks fine. Try a plot:

plot(Conc, Ratio, type="n", log="xy", las=1)
points(Conc, Ratio, pch=21, cex=1.5, col="blue", bg="#CCCCFF80")
curve(coef(m.4)[[1]]+coef(m.4)[[2]]*x+coef(m.4)[[3]]*x^2, range(Conc),
      lwd=2, col="darkgreen", add=TRUE)
curve(coef(m.8)[[1]]+coef(m.8)[[2]]*x+coef(m.8)[[3]]*x^2, range(Conc),
      lwd=2, col="red", add=TRUE)

Dif-tor heh smusma 🖖
Helmut Schütz
[image]

The quality of responses received is directly proportional to the quality of the question asked. 🚮
Science Quotes

Complete thread:

Activity
 Admin contact
21,419 posts in 4,475 threads, 1,510 registered users;
online 17 (0 registered, 17 guests [including 7 identified bots]).
Forum time: Wednesday 17:04 CEST (Europe/Vienna)

In the Middles Ages the lingua franca of science was Latin.
Nowadays the language of science is bad English.    Anonymous

The Bioequivalence and Bioavailability Forum is hosted by
BEBAC Ing. Helmut Schütz
HTML5