# Getting the right starting values for an nls model in R

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

• Although stated in terms of R (it really has to be stated in some language), how to find appropriate starting values for a non-linear model fit is sufficiently statistical to be on-topic here, IMO. It isn't really a programming Q, eg. Nov 26 '15 at 20:47

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

• That's extremely helpful, thank you so much! I do have a question about how you got your "a" value here though. I tried running lm(log10(rev) ~ log10(weeks)) and then using the "summary" function, and while I get the same "b" value, my "a" value comes out to 4.3201. What did you do differently to arrive at a=9.947? Nov 26 '15 at 23:57
• Notice that I used exp to take it back to unlogged values, which is a clue indicating that I used the plain log function. As long as you are consistent with which log and antilog you use, you'll get the same answer for the start value. So you can do base 10 and I can do base $e$ and everything's the same. Nov 27 '15 at 1:04
• Ah, you're totally right. Amateur mistake on my part. Kept thinking of mathematical notation expecting "log" to mean log base 10 and "ln" for natural log. Thanks for the clarification. Nov 27 '15 at 7:23
• To many mathematicians (and many statisticians), an unadorned "log" is the natural log, much as an unadorned argument to a sin function is in radians. [Clashing conventions can lead to confusion, unfortunately, but when I started using R, for example, I didn't think twice about the use of the log function since R and I share the same convention.] Nov 27 '15 at 8:38

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.

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)


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.