[R] nls function

[peter dalgaard]

On Apr 10, 2012, at 22:03 , nerak13 wrote:

> Hi,
> 
> I've got the following data:
> 
> x<-c(1,3,5,7)
> y<-c(37.98,11.68,3.65,3.93)
> penetrationks28<-dataframe(x=x,y=y)
> 
> now I need to fit a non linear function so I did:
> 
> fit <- nls(y ~ I(a+b*exp(1)^(-c * x)), data = penetrationks28, start =
> list(a=0,b = 1,c=1), trace = T)


That's the second time you have been posting code with that weird sort of formula. 

Please: 
y ~ a + b * exp(-c * x)

I() is only needed for linear models where arithmetic operators have special meaning, and exp(1)^x for exp(x) just hurts our eyes.

Apart from that, you just need better starting values, it seems:

> fit <- nls(y ~ a + b * exp(-c * x),data = penetrationks28, start =
+ list(a=3,b = 15,c=1), trace = T)
932.0747 :   3 15  1 
310.2288 :   2.9756409 25.5127564  0.3610715 
137.811 :   5.6501434 42.5276986  0.6921207 
2.187083 :   2.6843590 71.7142051  0.7180404 
1.9709 :   2.6657517 71.6213474  0.7058412 
1.970713 :   2.6726912 71.6839169  0.7069178 
1.970712 :   2.6719547 71.6797621  0.7068432 
1.970712 :   2.6720049 71.6800584  0.7068484 
> summary(fit)

Formula: y ~ a + b * exp(-c * x)

Parameters:
  Estimate Std. Error t value Pr(>|t|)  
a   2.6720     1.4758   1.811   0.3212  
b  71.6801     7.7101   9.297   0.0682 .
c   0.7068     0.1207   5.856   0.1077  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 1.404 on 1 degrees of freedom

Number of iterations to convergence: 7 
Achieved convergence tolerance: 3.067e-06 

> fitted(fit)
[1] 38.024329 11.271188  4.763691  3.180792
attr(,"label")
[1] "Fitted values"
 
..........
As David points out, three parameters to 4 data values is a bit absurd. Although not technically infeasible, I would not trust those 1 df t-tests at all.

Concerning the starting values, finding them is a bit of an art, but you clearly need something better than grabbing values out of thin air. Various techniques for finding starting values can be extracted by looking at the code for the self-start routines that come with nls() (the ones with names starting with SS). 

In the present case, you might try guessing a as the asymptote for large x so a value somewhat lower than the last value of y should do, say a=3. Subtract a from both sides and you get y-a = b*exp(-c*x) which you can make linear as log(y-a) = log(b) - c*x. So fit this using

> lm(log(y-3)~x,data=penetrationks28)

Call:
lm(formula = log(y - 3) ~ x, data = penetrationks28)

Coefficients:
(Intercept)            x  
     3.9979      -0.6737  

which suggests that a=3, b=exp(4), c=.67 might work, and it does, even though the starting values aren't exactly spot on:

> fit <- nls(y ~ a + b * exp(-c * x),data = penetrationks28, start =
+ list(a=3,b = exp(4),c=.67), trace = T)
53.23186 :   3.00000 54.59815  0.67000 
2.188135 :   2.7131685 71.6770898  0.7187123 
1.970939 :   2.6654559 71.6169034  0.7057594 
...


-- 
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Email: pd.mes at cbs.dk  Priv: PDalgd at gmail.com