Singular gradient error in nls with correct starting values

Question

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!

Answer 1

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

Answer 2

The answers above are, of course, correct. For what its worth, in addition to the explanations given, if you are trying this on an artificial data set, according to the nls help page found at: http://stat.ethz.ch/R-manual/R-patched/library/stats/html/nls.html

R's nls wont be able to handle it. The help page specifically states:

Warning

Do not use nls on artificial "zero-residual" data.

The nls function uses a relative-offset convergence criterion that compares the numerical imprecision at the current parameter estimates to the residual sum-of-squares. This performs well on data of the form

y = f(x, θ) + eps

(with var(eps) > 0). It fails to indicate convergence on data of the form

y = f(x, θ)

because the criterion amounts to comparing two components of the round-off error. If you wish to test nls on artificial data please add a noise component, as shown in the example below.

So, no noise==no good for R's nls.