Example [Bioanalytics]

posted by Helmut Homepage – Vienna, Austria, 2019-03-17 01:56  – Posting: # 20043
Views: 672

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)

Cheers,
Helmut Schütz
[image]

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

Complete thread:

Activity
 Mix view
Bioequivalence and Bioavailability Forum |  Admin contact
19,489 posts in 4,135 threads, 1,336 registered users;
online 16 (1 registered, 15 guests [including 8 identified bots]).
Forum time (Europe/Vienna): 11:48 CEST

No rational argument will have a rational effect on a man
who does not want to adopt a rational attitude.    Karl R. Popper

The BIOEQUIVALENCE / BIOAVAILABILITY FORUM is hosted by
BEBAC Ing. Helmut Schütz
HTML5