Why is nls() giving me "singular gradient matrix at initial parameter estimates" errors?

Question

I have some basic data on emission reductions and cost per car:

q24 <- read.table(text = "reductions  cost.per.car
    50  45
    55  55
    60  62
    65  70
    70  80
    75  90
    80  100
    85  200
    90  375
    95  600
    ",header = TRUE, sep = "")

I know that this is an exponential function, so I expect to be able to find a model that fits with:

    model <- nls(cost.per.car ~ a * exp(b * reductions) + c, 
         data = q24, 
         start = list(a=1, b=1, c=0))

but I'm getting an error:

Error in nlsModel(formula, mf, start, wts) : 
  singular gradient matrix at initial parameter estimates

I've read through a ton of questions on the error I'm seeing and I'm gathering that the problem is probably that I need better/different start values (the initial parameter estimates makes a little more sense) but I'm not sure, given the data that I have, how I'd go about estimating better parameters.

Answer 1

Automatically finding good starting values for a nonlinear model is an art. (It's relatively easy for one-off datasets when you can just plot the data and make some good guesses visually.) One approach is to linearize the model and use least squares estimates.

In this case, the model has the form

$$\mathbb{E}(Y) = a \exp(b x) + c$$

for unknown parameters $a,b,c$. The presence of the exponential encourages us to use logarithms--but the addition of $c$ makes it difficult to do that. Notice, though, that if $a$ is positive then $c$ will be less than the smallest expected value of $Y$--and therefore might be a little less than the smallest observed value of $Y$. (If $a$ could be negative you will also have to consider a value of $c$ that is a little greater than the largest observed value of $Y$.)

Let us, then, take care of $c$ by using as initial estimate $c_0$ something like half the minimum of the observations $y_i$. The model can now be rewritten without that thorny additive term as

$$\mathbb{E}(Y) - c_0 \approx a \exp(b x).$$

That we can take the log of:

$$\log(\mathbb{E}(Y) - c_0) \approx \log(a) + b x.$$

That is a linear approximation to the model. Both $\log(a)$ and $b$ can be estimated with least squares.

Here is the revised code:

c.0 <- min(q24$cost.per.car) * 0.5
model.0 <- lm(log(cost.per.car - c.0) ~ reductions, data=q24)
start <- list(a=exp(coef(model.0)[1]), b=coef(model.0)[2], c=c.0)
model <- nls(cost.per.car ~ a * exp(b * reductions) + c, data = q24, start = start)

Its output (for the example data) is

Nonlinear regression model
  model: cost.per.car ~ a * exp(b * reductions) + c
   data: q24
        a         b         c 
 0.003289  0.126805 48.487386 
 residual sum-of-squares: 2243

Number of iterations to convergence: 38 
Achieved convergence tolerance: 1.374e-06

The convergence looks good. Let's plot it:

plot(q24)
p <- coef(model)
curve(p["a"] * exp(p["b"] * x) + p["c"], lwd=2, col="Red", add=TRUE)

It worked well!

When automating this, you might perform some quick analyses of the residuals, such as comparing their extremes to the spread in the ($y$) data. You might also need analogous code to deal with the possibility $a\lt 0$; I leave that as an exercise.

---

Another method to estimate initial values relies on understanding what they mean, which can be based on experience, physical theory, etc. An extended example of a (moderately difficult) nonlinear fit whose initial values can be determined in this way is described in my answer at https://stats.stackexchange.com/a/15769.

Visual analysis of a scatterplot (to determine initial parameter estimates) is described and illustrated at https://stats.stackexchange.com/a/32832.

In some circumstances, a sequence of nonlinear fits is made where you can expect the solutions to change slowly. In that case it's often convenient (and fast) to use the previous solutions as initial estimates for the next ones. I recall using this technique (without comment) at https://stats.stackexchange.com/a/63169.

Answer 2

This library was able to resolve my problem with nls's singular gradient: http://www.r-bloggers.com/a-better-nls/ An example:

library(minpack.lm)
nlsLM(function, start=list(variable=2,variable2=12))

Answer 3

So ... I think I mis-read this as an exponential function. All I needed was poly()

model <- lm(cost.per.car ~ poly(reductions, 3), data=q24)
new.data <- data.frame(reductions = c(91,92,93,94))
predict(model, new.data)

plot(q24)
lines(q24$reductions, predict(model, list(reductions = q24$reductions)))

Or, using lattice:

xyplot(cost.per.car ~ reductions, data = q24,
       panel = function(x, y) {
         panel.xyplot(x, y)
         panel.lines(x, predict(model,list(reductions = x) ))
       }, 
       xlab = "Reductions", 
       ylab = "Cost per car")