29
$\begingroup$

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.

$\endgroup$
7
  • $\begingroup$ I would suggest starting your deciphering by searching our site for the error message. $\endgroup$ – whuber Jul 8 '15 at 20:02
  • 3
    $\begingroup$ Actually, I did that and my search for the full error turned up a half baked question with three data points and no answer. But your more specific search does get some results. Possibly because you have more experience here and know which terms stand out as relevant. $\endgroup$ – Amanda Jul 8 '15 at 20:05
  • 1
    $\begingroup$ One thing I have found about software errors is that a search for the specific error message (usually in quotation marks) is the surest way to find out whether it has been discussed before. (This holds Internet-wide, not just on SE sites.) As our "on hold" message says, if your additional research does not resolve your issue, then please do come back and push back on us a little: this question is at the intersection of statistics and computing and might expose some issues of great interest here. $\endgroup$ – whuber Jul 8 '15 at 20:11
  • 1
    $\begingroup$ The fit for your start values is very far from the data; compare exp(50) and exp(95) to the y-values at x=50 and x=95. If you set c=0 and take log of y (making a linear relationship), you can use regression to get initial estimates for log($a$) and $b$ that will suffice for your data (or if you fit a line through the origin, you can leave $a$ at 1 and just use the estimate for $b$; that also suffices for your data). If $b$ is much outside a fairly narrow interval around those two values, you'll run into some problems. [Alternatively try a different algorithm] $\endgroup$ – Glen_b Jul 8 '15 at 23:09
  • 1
    $\begingroup$ Thanks @Glen_b. I was hoping I could use R in lieu of a graphing calculator to work through a stats intro textbook (and leapfrog the course itself) so I'm starting with only the barest statistical insight, but lots of experience doing other slicing and dicing in R. $\endgroup$ – Amanda Jul 9 '15 at 2:55
47
$\begingroup$

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)

Figure

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.

$\endgroup$
1
  • $\begingroup$ Sorry for my question, I am applying this to a similar model but my point start from 0, so when I start estimating the parameters c.0 <- min(q24$cost.per.car) * 0.5, in my case is 0 since my data starts in (0,0) and then increasing. Could you recommend something for this case? I understand that if the curve that you can see in the point is similar to the one that you got in this example, the I am working with a positive exponencial, otherwise it would be a negative, but in my case I think it is positive. Could you really give me some advice, please.Thank you. $\endgroup$ – Rachel Feb 8 at 23:52
2
$\begingroup$

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))
$\endgroup$
2
  • $\begingroup$ That function seems to be called nls.lm now. $\endgroup$ – Matt May 23 '16 at 14:53
  • $\begingroup$ The above comment is sort of right sort of wrong. nlsLM() and nls.lm() are two different functions that essentially do the same thing. nlsLM() is just nicer for people who already know how to use the nls() function because they use similar formats. $\endgroup$ – obewanjacobi Oct 7 '20 at 21:12
-1
$\begingroup$

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")
$\endgroup$
2
  • 2
    $\begingroup$ This doesn't answer the question you asked--it changes it to something different (and rather less interesting, IMHO). $\endgroup$ – whuber Jul 14 '15 at 21:36
  • 6
    $\begingroup$ Although, it may solve the problem of fitting a function to represent the data, your accepted answers is not the expect for your question. Mr. @whuber provided you an excellent explanation and deserve the accepted answer. $\endgroup$ – Lourenco Nov 10 '15 at 17:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.