Estimating non-monotonic bi-exponential curve fit

Question

I am doing some pharmacokinetic analyses and am fine with non-compartmental methods. But I trying to also learn some non-linear curve-fitting techniques.

If we have the following data:

df <-  data.frame(x = c(0,2,4,8,12,24,48,72), y = c(0.000,0.05,0.0671,0.068,0.05,0.0250,0.0103,0.0043))
plot(df$x, df$y)

and I want to fit a non-monotonic bi-exponential curve to it as described on page 579 of this book, using the equation highlighted below:

I am fine with the standard bi-exponential form and finding starting values using Sbiexp, but do not know how to approach finding starting values for this equation. Would someone be able to give me some pointers?

> nls_mod_bi <-  nls(y ~ k * (exp(-b1*x) - exp(-b2*x)), start = c(k = 1, b1 = 1, b2 = 1), data = df)
Error in nlsModel(formula, mf, start, wts, scaleOffset = scOff, nDcentral = nDcntr) : 
  singular gradient matrix at initial parameter estimates

Answer 1

Option 1: nls + SSbiexp

I think you can use SSbiexp in nls, and you don't need to manually setup the initial starting point

fit <- nls(y ~ SSbiexp(x, k1, b1, k2, b2), data = df)

f <- function(x, p = coef(fit)) {
  p["k1"] * exp(-exp(p["b1"]) * x) + p["k2"] * exp(-exp(p["b2"]) * x)
}

plot(y ~ x, df, type = "p", col = "blue")
curve(f, range(df$x), col = "red", add = TRUE)

Option 2: constrOptim

x <- df$x
y <- df$y
f <- function(p) {
  norm(p[1] * (exp(-p[2] * x) - exp(-p[3] * x)) - y, "2")
}
sol <- constrOptim(c(1, 0.1, 0.1), f, grad = NULL, ui = diag(3), ci = rep(0,3))
sol$par

which gives

> sol$par # (k, b1, b2)
[1] 0.10780726 0.05793686 0.43720768

which is really close to the nls approach outcome

> fit <- nls(y ~ SSbiexp(x, k1, b1, k2, b2), data = df)

> p <- coef(fit)

> c(p["k1"], exp(p["b1"]), p["k2"], exp(p["b2"]))
         k1          b1          k2          b2
-0.10802463  0.43855265  0.10774247  0.05791062