Search code examples
rregressionnlsnon-linear

Why the singular gradient error for this bi-exponential model


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)

enter image description here


Solution

  • 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:

    • this is a small data set (7 total points, two of which are close together and don't provide distinct information, so effectively about 6 df, to estimate 4 nonlinear parameters),
    • the pattern is hardly distinguishable from a monoexponential, so there's not much to work with.

    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)
    

    plot of data showing an exponential-like curve, of 7 data points (x range 0-250, y range 5e-5 to 2.5e-3). A red curve (biexponential fit) and a blue curve (monoexponential curve) are nearly indistinguishable.