## Example [Bioanalytics]

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

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