[R] fitting weibull curve to data using nls

[John.Morrongiello at csiro.au]

Thanks very much for identifying this error Michael- of course I don't want to exponentiate an exponential! I had copied the provided formula a little too literally and hadn't done my QA/QCs

John

-----Original Message-----
From: R. Michael Weylandt [mailto:michael.weylandt at gmail.com] 
Sent: Thursday, 1 November 2012 6:38 PM
To: Morrongiello, John (CMAR, Hobart)
Cc: r-help at r-project.org
Subject: Re: [R] fitting weibull curve to data using nls

On Thu, Nov 1, 2012 at 2:44 AM,  <John.Morrongiello at csiro.au> wrote:
> Hi
> I'd like to fit an asymmetrical curve function to some physiological data. I've been told a weibull curve is a good place to start, but I'm having trouble specifying and fitting the function with nls and was wondering if someone could help.
>
> After some reading, I think the function specification I want is
> y=c*(x/a)^(b-1)*e^(-(x/a)^b)
>
> (from 
> http://www.mathworks.se/products/statistics/demos.html?file=/products/
> demos/shipping/stats/cfitdfitdemo.html)
> Where *a* scales curve along the horizontal axis, *b* defines the 
> curve shape and *c* scales curve on y axis
>
> However, when I try to fit this function using nls I get the following error: 'Error in exp^(-(temp/a)^b) : non-numeric argument to binary operator'
>
> ##example dataset
> temp<-c(7,9,10,13,13.5,13.5,14,20,22,22,25,25,30,30,35,35,40,40,42,42.
> 5,45,45.6,46,46,47,47,47.5,48,48)
> response<-c(0.16,0.17,0.18,0.2,0.18,0.25,0.2,0.344,0.4,0.45,0.628,0.55
> ,0.783,0.7,0.9,0.85,1,1.05,1,0.95,0.982,0.97,0.85,0.9,0.8,0.85,0.7,0.6
> ,0.65)
> data<-data.frame(cbind(temp,response))

Don't do this. It doesn't actually hurt here, but the
data.frame(cbind(...)) idiom generally will screw things up. Rather just data.frame(temp, response).

> plot(response~temp,data)
>
> ##fit function with nls
> m1<-nls(response~(c*((temp/a)^(b-1)))*exp^(-(temp/a)^b),start=list(a=4
> 0,b=3,c=1.5),data=data)
> ##returns following error:
> Error in exp^(-(temp/a)^b) : non-numeric argument to binary operator 
> ### Is this because my starting values are wrong or something amiss 
> with the function specification?
>

Yes, something's wrong with the function specification.

You want

exp(-(temp/a)^b)

[Note the removed ^]

exp() is a function which does exponentiation, you don't want to try to exponentiate it.

> nls(response ~ c*((temp/a)^(b-1))*exp(-(temp/a)^b), start = list(a = 
> 40, b = 3,c = 1.5), data = dd)
Nonlinear regression model
  model:  response ~ c * ((temp/a)^(b - 1)) * exp(-(temp/a)^b)
   data:  dd
     a      b      c
43.128  3.365  2.504
 residual sum-of-squares: 0.1678

Number of iterations to convergence: 9
Achieved convergence tolerance: 6.384e-06

Cheers,
Michael Weylandt