Example [Bioanalytics]

posted by Helmut Homepage – Vienna, Austria, 2019-03-17 01:56 (713 d 17:12 ago) – Posting: # 20043
Views: 3,467

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,357 posts in 4,459 threads, 1,492 registered users;
online 6 (0 registered, 6 guests [including 2 identified bots]).
Forum time: Saturday 19:09 CET (Europe/Vienna)

The way to do research is to attack the facts
at the point of greatest astonishment.    Celia Green

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