Bioequivalence and Bioavailability Forum • Example

Example [Bioanalytics]

posted by Helmut – Vienna, Austria, 2019-03-17 02:56 (2261 d 02:23 ago) – Posting: # 20043
Views: 7,179

Hi ElMaestro,

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

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/x^1.3355 is the winner! Duno why the ICs of 1/s_y² 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

Complete thread:

Handling BLOQ values (Fisher Info etc.) mittyri 2019-03-16 14:41
- Handling BLOQ values (Fisher Info etc.) Helmut 2019-03-16 15:39
  - The smartest solution ElMaestro 2019-03-16 22:35
    - An old solution Helmut 2019-03-16 23:05
      - old, obvious? ElMaestro 2019-03-19 08:26
    - ExampleHelmut 2019-03-17 01:56
  - ADA example mittyri 2019-03-18 22:20
    - My beloved Ada Helmut 2019-03-19 01:34
- Handling BLOQ values (Fisher Info etc.) Ohlbe 2019-03-17 14:32
  - Handling BLOQ values (Fisher Info etc.) nobody 2019-03-18 08:14
    - Don’t weight by 1/s² Helmut 2019-03-18 11:19