I am wanting to compare a mono- and bi-exponential fit on the following data. The mono-exponential function fits fine, but I can't find starting values for the bi-exponential function (I have also tried entering some manually based on plausible values).
Can someone explain why a bi-exp function won't fit these data and the cause? Is there another package that will potentially work or is this a limitation of the data?
Thanks
# Dummy data
dat <- structure(list(x = c(5, 10, 30, 60, 120, 180, 240), y = c(0.000215875066736858,
0.000213143286462547, 0.000184642176655826, 0.000154069098478106,
0.000100308278926715, 6.56265655274e-05, 4.80349446283348e-05
)), class = c("tbl_df", "tbl", "data.frame"), row.names = c(NA,
-7L))
# Mono-exponential model
nls_mod_mono <- nls(y ~ statforbiology::NLS.expoDecay(x, C0, lambda1), dat)
# Predict and plot
xNew <- seq(0, 240, length.out = 100) # new grid of times
yNew <- predict(nls_mod_mono, list(x = xNew)) # predict activity
dfNew <- data.frame(x = xNew, y = yNew)
ggplot(data = dfNew, aes(x = x, y = y*100)) +
geom_line() +
geom_point(data = dat, aes(x = x, y = y*100), size = 3)
# Bi-exponential model
params <- getInitial(y ~ SSbiexp(x, intercept_initial, log_slope_initial, intercept_terminal, log_slope_terminal), data = dat)
#...failing to find initial values
nls_mod_bi <- nls(y ~ C0*exp(-lambda1*x) + C1*exp(-lambda2*x), start = list(C0 = params[1], lambda1 = exp(params[2]), C1 = params[3], lambda2 = exp(params[4])), data = dat)
I dug around with this for a bit but don't have a definitive answer for why it doesn't work, with the obvious "why it should be hard" reasons:
You can debug your way through the initialization machinery via debug(attr(SSbiexp, "initial"))
, or trace its progress by adding trace = TRUE
to your getInitial
call ... the internal machinery uses the "plinear" (partially linear) algorithm, which should be more robust, but fails anyway ...
I managed to get a decent fit, eventually, by brute-forcing with Nelder-Mead optimization via the bbmle
package ... unfortunately if you have a lot of data sets like this and you want to fit them in an unsupervised way, without lots of tweaking of settings and initial parameters, this will be quite a headache ...
library(bbmle)
## fit monoexponential
fit0 <- bbmle::mle2(y ~ dnorm(mean = A1*exp(-exp(lrc1)*x),
sd = exp(logsd)),
start = list(A1 = 1e-4, lrc1 = -5, logsd = -3),
data = dat, method = "Nelder-Mead",
trace = TRUE)
## tweak monoexponential fit slightly for starting values
fit <- bbmle::mle2(y ~ dnorm(mean = A1*exp(-exp(lrc1)*x) + A2*exp(-exp(lrc2)*x),
sd = exp(logsd)),
start = list(A1 = 1e-4, A2 = 1e-4,
lrc1 = -5, lrc2 = -3, logsd = -13),
data = dat, method = "Nelder-Mead",
trace = TRUE)
The biexponential fit only improves the log-likelihood by logLik(fit) - logLik(fit0)
= 0.18 units, which is quite small ...
plot(y~x, data = dat)
xvec <- 0:250
lines(xvec, predict(fit0, newdata = data.frame(x = xvec)), col = 4)
lines(xvec, predict(fit, newdata = data.frame(x = xvec)), col = 2)