Re: [R] fitting "power model" in nls()

From: Rolf Turner <r.turner_at_auckland.ac.nz>
Date: Mon, 3 Dec 2007 08:50:26 +1300

Rule number 1: Read the help for the function you are using. You must supply starting values for the fit --- which the code you gave doesn't do.

Rule number 2: Don't use nls()! Endless grief results.

Instead try:

foo <- function(par,x,y){
  Const <- par[1]
  B <- par[2]
  A <- par[3]
  sum((y-(Const+B*(x^A)))^2)
  }

fit <- optim(c
(1,2.5,0.2),foo,x=area,y=richness,method="BFGS",control=list (maxit=1000))

plot(area,richness)
ccc <- fit$par
curve(ccc[1]+ccc[2]*(x^ccc[3]),from=0,to=2000,add=TRUE,col="red") # Fit doesn't look too bad.

HTH.                 cheers,

                        Rolf Turner

On 3/12/2007, at 8:08 AM, Milton Cezar Ribeiro wrote:

> Dear all,
> I am still fighting against my "power model".
> I tryed several times to use nls() but I canīt run it.
> I am sending my variables and also the model which I would like to
> fit.
> As you can see, this "power model" is not the best model to be fit,
> but I really need also to fit it.
>
> The model which I would like to fit is Richness = B*(Area^A)
>
> richness<-c
> (44,36,31,39,38,26,37,33,34,48,25,22,44,5,9,13,17,15,21,10,16,22,13,20
> ,9,15,14,21,23,23,32,29,20,
> 26,31,4,20,25,24,32,23,33,34,23,28,30,10,29,40,10,8,12,13,14,56,47,44,
> 37,27,17,32,31,26,23,31,34,
> 37,32,26,37,28,38,35,27,34,35,32,27,22,23,13,28,13,22,45,33,46,37,21,2
> 8,38,21,18,21,18,24,18,23,22,
> 38,40,52,31,38,15,21)
> area<-c
> (26.22,20.45,128.68,117.24,19.61,295.21,31.83,30.36,13.57,60.47,205.30
> ,40.21,
> 7.99,1.18,5.40,13.37,4.51,36.61,7.56,10.30,7.29,9.54,6.93,12.60,
> 2.43,18.89,15.03,14.49,28.46,36.03,38.52,45.16,58.27,67.13,92.33,1.17,
> 29.52,84.38,87.57,109.08,72.28,66.15,142.27,76.41,105.76,73.47,1.71,30
> 5.75,
> 325.78,3.71,6.48,19.26,3.69,6.27,1689.67,95.23,13.47,8.60,96.00,436.97
> ,
> 472.78,441.01,467.24,1169.11,1309.10,1905.16,135.92,438.25,526.68,88.8
> 8,31.43,21.22,
> 640.88,14.09,28.91,103.38,178.99,120.76,161.15,137.38,158.31,179.36,21
> 4.36,187.05,
> 140.92,258.42,85.86,47.70,44.09,18.04,127.84,1694.32,34.27,75.19,54.39
> ,79.88,
> 63.84,82.24,88.23,202.66,148.93,641.76,20.45,145.31,27.52,30.70)
> plot(richness~area)
>
> I tryed to fit the following model:
>
> m1<-nls(richness ~ Const+B*(area^A))
>
> Thanks a lot,
> miltinho
> Brazil.
>
>
>
> para armazenamento!
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> R-help_at_r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-
> guide.html
> and provide commented, minimal, self-contained, reproducible code.

######################################################################
Attention:
This e-mail message is privileged and confidential. If you are not the intended recipient please delete the message and notify the sender. Any views or opinions presented are solely those of the author.

This e-mail has been scanned and cleared by MailMarshal www.marshalsoftware.com

######################################################################

______________________________________________
R-help_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code. Received on Sun 02 Dec 2007 - 19:54:28 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Sun 02 Dec 2007 - 20:30:17 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.