Solved – Singular gradient error in nls with correct starting values

nlsnonlinear regressionr

I'm trying to fit a line+exponential curve to some data. As a start, I tried to do this on some artificial data. The function is:
$$y=a+b\cdot r^{(x-m)}+c\cdot x$$
It is effectively an exponential curve with a linear section, as well as an additional horizontal shift parameter (m). However, when I use R's nls() function I get the dreaded "singular gradient matrix at initial parameter estimates" error, even if I use the same parameters that I used to generate the data in the first place.
I've tried the different algorithms, different starting values and tried to use optim to minimise the residual sum of squares, all to no avail. I've read that a possible reason for this could be an over-parametrisation of the formula, but I don't think it is (is it?)
Does anyone have a suggestion for this problem? Or is this just an awkward model?

A short example:

#parameters used to generate the data
reala=-3
realb=5
realc=0.5
realr=0.7
realm=1
x=1:11 #x values - I have 11 timepoint data
#linear+exponential function
y=reala + realb*realr^(x-realm) + realc*x
#add a bit of noise to avoid zero-residual data
jitter_y = jitter(y,amount=0.2)
testdat=data.frame(x,jitter_y)

#try the regression with similar starting values to the the real parameters
linexp=nls(jitter_y~a+b*r^(x-m)+c*x, data=testdat, start=list(a=-3, b=5, c=0.5, r=0.7, m=1), trace=T)

Thanks!

Best Answer

I've got bitten by this recently. My intentions were the same, make some artificial model and test it. The main reason is the one given by @whuber and @marco. Such model is not identified. To see that, remember that NLS minimizes the function:

$$\sum_{i=1}^n(y_i-a-br^{x_i-m}-cx_i)^2$$

Say it is minimized by the set of parameters $(a,b,m,r,c)$. It is not hard to see that the the set of parameters $(a,br^{-m},0,r,c)$ will give the same value of the function to be minimized. Hence the model is not identified, i.e. there is no unique solution.

It is also not hard to see why the gradient is singular. Denote

$$f(a,b,r,m,c,x)=a+br^{x-m}+cx$$

Then

$$\frac{\partial f}{\partial b}=r^{x-m}$$

$$\frac{\partial f}{\partial m}=-b\ln rr^{x-m}$$

and we get that for all $x$

$$b\ln r\frac{\partial f}{\partial b}+\frac{\partial f}{\partial m}=0.$$

Hence the matrix

\begin{align} \begin{pmatrix} \nabla f(x_1)\\\\ \vdots\\\\ \nabla f(x_n) \end{pmatrix} \end{align}

will not be of full rank and this is why nls will give the the singular gradient message.

I've spent over a week looking for bugs in my code elsewhere till I noticed that the main bug was in the model :)

Related Question