How to avoid "error:singular gradient matrix at initial parameter estimates" using nlsLM?

Question

This is my first post, so I hope that I include everything correctly and thank you for any help that you may be able to give.

I am trying to fit a curve, to a set of raw data which I have collected. The data is sigmoid and should be explained by a Shae-Ackers model (commonly used in biology). I initially tried to do this using the 'nls' function, after looking over other examples on this site. However I received the following error:

"error:singular gradient matrix at initial parameter estimates".

Subsequently, I tried to use the 'nlsLM' function. As I understand from other questions on this site, this should use an alternative algorithm to fit my model which avoids the error. However, I still get the same error:

"error:singular gradient matrix at initial parameter estimates".

Here is the code which I am currently trying to get to work:

library("minpack.lm") 
library("nls2")

L = c(0, 0.0001, 0.0005, 0.001, 0.005, 
0.01, 0.05, 0.1, 0.5, 1, 5, 10, 50) # X- values, inducer concentration.
Fluo = c(23.10263117,
     21.5111054,
     25.01080225,
     32.63906667,
     82.6671047,
     287.1788694,
     812.2339928,
     1308.71973,
     2822.260637,
     3675.085354,
     4634,
     4399.131349,
     4096.759224) # Y-values, fluorescence measurement.

model <- nlsLM(Fluo ~ Fmax*((k1 + k3C0*((L^n)/(kd^n + L^n)))/(1 + k1 + k2C0*((L^n)/(kd^n + L^n)) + k3C0*((L^n)/(kd^n + L^n)))),
           start = c(k1 = 0.009,k2C0 = 37.5,k3C0 = 3.4,kd = 0.09,n = 2.8,Fmax = 7650), algorithm = "LM", trace = TRUE)

plot(log10(L), log10(Fluo), main = "data")
lines(log10(L), log10(fitted(model)), col = "red", lty = 2)  

The starting estimates which I have provided should be very close to the values I am looking for, as these have been previously reported in literature for the system I am trying to fit.

So my overall questions are this:
1. Is there a way that I can find alternative starting estimates which will avoid this error?
2. Which other functions other than 'nls' or 'nlsLM' may I be able to use to fit my model? (I have also looked at 'nls.lm' but am unsure how to implement this correctly).

Thank you very much for any help and I would be happy to provide any other information if needed!

Answer 1

The problem isn't which optimizer to use. The problem is that the model parameters are not identifiable.

Let the right hand side be:

rhs <- function(k1, k2C0, k3C0, kd, n, Fmax, L) {
  r <- L^n / (kd^n + L^n)
  Fmax*((k1 + k3C0 * r)/(1 + k1 + k2C0*r + k3C0*r))
}

Then note that these two calls to rhs give the same values (except for 0 which is degenerate):

st <- c(k1 = 0.009, k2C0 = 37.5, k3C0 = 3.4, 
  kd = 0.09, n = 2.8, Fmax = 7650)

rhs1 <- with(as.list(st), rhs(k1, k2C0, k3C0, kd, n, Fmax, L))

r <- with(as.list(st), L^n/(kd^n+L^n))
rhs2 <- with(as.list(st), rhs(2*k1, (k2C0 * r - k1 - r * k3C0)/r, 2*k3C0, 
 kd, n, Fmax/2, L))

all.equal(rhs1[-1], rhs2[-1])
## [1] TRUE

Singular values

We can look at the singular values and we see that the smallest two singular values are very small and even the third smallest is several orders of magnitude smaller than the largest so we are going to have to fix 2 or 3 of the parameters in order to apply any optimization algorithm.

library(nls2)
fm2 <- nls2(Fluo ~ rhs(k1, k2C0, k3C0, kd, n, Fmax, L), st = st, alg = "brute")
svd(fm2$m$Rmat())[-2]

giving:

$`d`
[1] 1.815654e+04 2.063566e+03 4.100935e+02 4.992959e+01 1.234634e-06
[6] 1.526028e-09

$v
              [,1]          [,2]          [,3]          [,4]          [,5]
[1,] -9.993013e-01  0.0373287823  0.0014810192  1.127961e-03 -2.744777e-09
[2,]  9.184235e-05  0.0058555553 -0.0879560839  3.069373e-03 -9.960997e-01
[3,] -1.388630e-03 -0.0755201221  0.9926417526 -3.431623e-02 -8.820121e-02
[4,]  3.726825e-02  0.9958422212  0.0744652781 -3.692667e-02 -8.316237e-04
[5,]  2.458570e-03  0.0341651105  0.0371291870  9.987232e-01 -2.370720e-09
[6,] -1.845720e-06 -0.0000361664  0.0004803576 -1.661631e-05  2.309941e-03
              [,6]
[1,] -1.187056e-06
[2,]  2.343452e-03
[3,] -2.763882e-04
[4,]  1.622277e-06
[5,] -1.799465e-11
[6,]  9.999972e-01

Fixing parameters

For example, if we fix parameters 2, 5 and 6 (k2C0, n, Fmax) at the starting values given in the question then we get the following (also shown as the blue fit in the graph at the end).

fm <- with(as.list(st), nls(Fluo ~ rhs(k1, k2C0, k3C0, kd, n, Fmax, L), 
  start = st[-c(2, 5:6)]))

giving:

Nonlinear regression model
  model: Fluo ~ rhs(k1, k2C0, k3C0, kd, n, Fmax, L)
   data: parent.frame()
      k1     k3C0       kd 
 0.04214 49.21873  2.02506 
 residual sum-of-squares: 1737050

Number of iterations to convergence: 8 
Achieved convergence tolerance: 1.947e-06

Note

If you are willing to use a different formula then this one only has 4 parameters (vs. 6 in the question) and provides a reasonable fit. It is shown as the red fit in the graph at the end.

By the way, in the question log10(Fluo) is plotted against log10(L) but minimizing the sum of the squares of log10(Flou) - log10(rhs(...)) is not the same as minimizing the sum of squares of Flou - rhs(...) although the answer may be close.

fm2 <- nls(Fluo ~ L^d / (a + b * L^c), 
 start = c(a = 10e-5, b = 10e-5, c = 1, d = 1))

plot(Fluo ~ L)
lines(fitted(fm) ~ L, col = "blue")
lines(fitted(fm2) ~ L, col = "red")
legend("bottomright", legend = c("fm", "fm2"),
  col = c("blue", "red"), lty = 1, lwd = 2)