Example [Bioanalytics]

posted by Helmut Homepage – Vienna, Austria, 2019-03-17 02:56 (1718 d 15:45 ago) – Posting: # 20043
Views: 4,971

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 🖖🏼 Довге життя Україна! [image]
Helmut Schütz
[image]

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

Complete thread:

UA Flag
Activity
 Admin contact
22,811 posts in 4,783 threads, 1,639 registered users;
21 visitors (0 registered, 21 guests [including 8 identified bots]).
Forum time: 18:42 CET (Europe/Vienna)

Every man gets a narrower and narrower field of knowledge
in which he must be an expert in order to compete with other people.
The specialist knows more and more about less and less
and finally knows everything about nothing.    Konrad Lorenz

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