Getting the right starting values for an nls model in R

Question

I'm trying to fit a simple power law model to a data set that is as follows:

mydf:

rev     weeks
17906.4 1
5303.72 2
2700.58 3
1696.77 4
947.53  5
362.03  6

The goal being to pass the power line through and use it to predict rev vlaues for future weeks. A bunch of research has led me to the nls function, which I implemented as follows.

newMod <- nls(rev ~ a*weeks^b, data=modeldf, start = list(a=1,b=1))
predict(newMod, newdata = data.frame(weeks=c(1,2,3,4,5,6,7,8,9,10)))

While this works for an lm model, I get a singular gradient error, which I understand has to do with my starting values a and b. I tried different values, even going so far as to plot this in Excel, pass a lone, get an equation, then use the values from the equation, but I still get the error. I looked at a bunch of answers like this one and tried the second answer (couldn't understand the first), but to no result.

I could really use some help here on how to find the right starting values. Or alternatively, what other function I can use instead of nls.

In case you want to recreate mydf with ease:

mydf <- data.frame(rev=c(17906.4, 5303.72, 2700.58 ,1696.77 ,947.53 ,362.03), weeks=c(1,2,3,4,5,6)) 

Answer 1

This is a common problem with nonlinear least squares models; if your start values are very far from the optimum the algorithm may not converge, even though it may be well behaved near the optimum.

If you start by taking logs of both sides and fit a linear model, you get estimates of $\log(a)$ and $b$ as the slope and intercept ( 9.947 and -2.011 ) (edit: that's natural log)

If you use those to guide the starting values for $a$ and $b$ everything seems to work okay:

 newMod <- nls(rev ~ a*weeks^b, data=mydf, start = list(a=exp(9.947),b=-2.011))
 predict(newMod, newdata = data.frame(weeks=c(1,2,3,4,5,6,7,8,9,10)))
 [1] 17919.2138  5280.7001  2584.0109  1556.1951  1050.1230   761.4947   580.3091   458.6027
 [9]   372.6231   309.4658

Answer 2

Try

 newMod <- nls(rev ~ a*weeks^b, data=mydf, startlist(a=17919.2127344,b=-1.76270557120))

I've been asked to expand this answer a bit. This problem is so simple I'm kind of surprised that nls fails at it. The real problem however is with the entire R approach and philosophy of nonlinear model fitting. In the real world one would scale x to lie between -1 and 1 and y and y to lie between 0 an 1 (y=ax^b). That would probably be enough to get nls to converge. Of course as Glen points out you can fit the corresponding log-linear model. That relies on the fact that there exists a simple transformation which linearizes the model. That is often not the case. The problem with R routines like nls is that they do not offer support for reparameterizing the model. In this case the reparameterization is simple, just rescale/recentre x and y. However having fit the model the user will have different parameters a and b from the original ones. While it is simple to calculate the original ones from these, the other difficulty is that it is not so simple in general to get the estimated standard deviations for these parameter estimates. This is done by the delta method which involves the Hessian of the log-likelihood and some derivatives. Nonlinear parameter estimation software should supply these calculations automatically, so that reparameterization of the model is easily supported. Another thing which software should support is the notion of phases. You can think of first fitting the model with Glen's version as phase 1. The "real" model is fit in stage 2.

I fit your model with AD Model Builder which supports phases in a natural way. In the first phase only a was estimated. This gets your model into the ballpark. In the second phase a and b are estimated to get the solution. AD Model Builder automatically calculates the standard deviations for any function of the model parameters via the delta method so that it encourages stable reparameterization of the model.

Answer 3

The Levenberg-Marquardt algorithm can help:

modeldf <- data.frame(rev=c(17906.4, 5303.72, 2700.58 ,1696.77 ,947.53 ,362.03), weeks=c(1,2,3,4,5,6))

require(minpack.lm)
fit <- nlsLM(rev ~ a*weeks^b, data=modeldf, start = list(a=1,b=1))

require(broom)
fit_data <- augment(fit)

plot(.fitted~rev, data=fit_data)

Answer 4

In my experience a good way of finding starting values for parameters of NLR models is to use an evolutionary algorithm. From an initial population (100) of random estimates (parents) in a search space choose the best 20 (offspring) and use these to help define a search in a succeeding population. Repeat until convergence. No need for gradients or hessians, just SSE evaluations. If you are not too greedy this very often works. The problems that people often have is that they are using a local search (Newton-Raphson) to perform the work of a global search. As always it is a matter of using the correct tool for the job at hand. It makes more sense to use an EA global search to find starting values for the Newton local search, and then let this run down to the minimum. But, as with all things, the devil is in the detail.