[R] NLS equation self starting non linear

[Peter Ehlers]

On 2010-09-02 22:32, Gabor Grothendieck wrote:
> On Thu, Sep 2, 2010 at 7:39 PM, Marlin Keith Cox<marlinkcox at gmail.com>  wrote:
>> This data are kilojoules of energy that are consumed in starving fish over a
>> time period (Days).  The KJ reach a lower asymptote and level off and I
>> would like to use a non-linear plot to show this leveling off.  The data are
>> noisy and the sample sizes not the largest.  I have tried selfstarting
>> weibull curves and tried the following, both end with errors.
>>
>> Days<-c(12, 12, 12, 12, 22, 22, 22, 22, 32, 32, 32, 32, 37, 38, 38, 39, 39,
>> 39, 39, 39, 39, 39,  1,  1,  1,  1)
>> Joules<-c(8.782010,  8.540524,  8.507526, 11.296904,  4.232690, 13.026588,
>> 10.213342,  4.771482,  4.560359,  6.146684, 9.651727,  8.064329,  9.419335,
>> 7.129264,  6.079012,  7.095888, 17.996794,  7.028287,  8.028352,  5.497564,
>> 11.607090,  9.987215, 11.065307, 21.433094, 20.366385, 22.475717)
>> X11()
>> par(cex=1.6)
>> plot(Joules~Days,xlab="Days after fasting was initiated",ylab="Mean energy
>> per fish (KJ)")
>> model<-nls(joules~a+b*exp(-c*Days),start=list(a=8,b=9,c=-.229),
>> control=list(minFactor=1e-12),trace=TRUE)
>> summary(model)
>
> Note that Joules is defined above but joules is used as well.  We have
> assumed they are the same.
>
> Also note that the formula is linear in log(joules) if a=0 so lets
> assume a=0 and use lm to solve.  Also switch to the "plinear"
> algorithm which only requires starting values for the nonlinear
> parameters -- in this case only c (which we get from the lm).
>
>> joules<- Joules
>> coef(lm(log(joules) ~ Days))
> (Intercept)        Days
>   2.61442015 -0.01614845
>>
>> # so use 0.016 as starting value of c
>> fm<- nls(joules~cbind(1, exp(-c*Days)), start = list(c = 0.016), alg = "plinear"); fm
> Nonlinear regression model
>    model:  joules ~ cbind(1, exp(-c * Days))
>     data:  parent.frame()
>        c   .lin1   .lin2
>   0.2314  8.3630 13.2039
>   residual sum-of-squares: 290.3
>

Keith,

You would get Gabor's solution from your nls() call if you did:

1. replace 'Joules' with 'joules'
2. replace 'c=-.229' with 'c=.229' in your start vector. You
already have the 'minus' in your formula.

I find it useful to add the curve you get with your start values
to the scatterplot to see how reasonable the values are:

plot(Joules~Days)
curve(8 + 9 * exp(-.229 * x), add=TRUE)

After you get a convergent solution, you can add a curve with
the model values.

   -Peter Ehlers