inflection point for 3-parameter weibull in R?

Question

The drc package in R contains self-starting functions to fit a 3-parameter Weibull model to data. The package gives 2 parameterizations of the 3-parameter weibull (https://cran.r-project.org/web/packages/drc/drc.pdf):

f(x) = 0 + (d − 0) exp(− exp(b(log(x) − log(e))))
f(x) = 0 + (d − 0)(1 − exp(− exp(b(log(x) − log(e)))))

Documentation for the package suggests that d represents the upper asymptote, b is the rate parameter, and that the curve has an inflection point at dose e. However, plotting these functions and using R to solve for when the 2nd derivative = 0 gives a value different than e for the inflection point. For example:

dd=1
bb = -5
ee = 30
curve(dd*exp(-exp(bb * (log(x) - log(ee))) ) , xlim=c(0,100))
abline(v=ee)

# 1st deriv
g <- function(x) {}
body(g) <- D( expression(dd*exp(-exp(bb * (log(x) - log(ee))) )), "x")
curve(g, xlim=c(0, 100))
abline(v=ee)

# 2nd deriv
g <- function(x) {}
body(g) <- D(D( expression(dd*exp(-exp(bb * (log(x) - log(ee))) )), "x"), "x")
curve(g, xlim=c(0, 100))
abline(v=ee)

uniroot(g, c(20,50))  # should be 30, but is not?

Is the documentation incorrect? Or does this have to do with the particular logged parameterization of the weibull (that constrains it to be positive)?

Answer 1

I think the documentation is wrong.

exp(-exp(x)) has an inflection point at 0:

DD <- function(expr, name, order = 1) {
  if(order < 1) stop("'order' must be >= 1")
  if(order == 1) D(expr, name)
  else DD(D(expr, name), name, order - 1)
}

DD( expression(exp(-exp(x))) , "x" , 2)
#-(exp(-exp(x)) * exp(x) - exp(-exp(x)) * exp(x) * exp(x))

exp(0) equals 1. Thus, we have: -(1/e - 1*e) = 0

I can see how one could then assume that exp(-exp(log(x) - log(c))) would also have an inflection point when the argument of the inner exponential gets zero. Bbut that's just not the case because the log changes the curvature and the second derivative:

DD( expression(exp(-exp(log(x) - log(c)))) , "x" , 2)
#-(exp(-exp(log(x) - log(c))) * (exp(log(x) - log(c)) * (1/x) * 
#    (1/x) - exp(log(x) - log(c)) * (1/x^2)) - exp(-exp(log(x) - 
#    log(c))) * (exp(log(x) - log(c)) * (1/x)) * (exp(log(x) - 
#    log(c)) * (1/x)))

If we put x = c, we have:

-(1/e * (1 * 1/c * 1/c - 1 * 1/c^2) - 1/e * 1 * 1/c * 1 * 1/c)

= -(0 - 1/e * 1/c^2)

= 1/e * 1/c^2

And that is not 0.

Just to check my maths:

g <- function(x) {}
body(g) <- DD( expression(exp(-exp(log(x) - log(c)))) , "x" , 2)
c <- 2
g(2)
#[1] 0.09196986
1/exp(1) * 1/2^2
#[1] 0.09196986

(For extra credit you should do the second derivative by hand. I don't have that much time.)