[R] newbie question: fitting power law

[Douglas Bates]

Thomas Hoffmann <hoffmann at giub.uni-bonn.de> writes:

> Dear R-helpers, dear Andrew
> 
> Andrew C. Ward schrieb:
> 
> >Dear Thomas,
> >
> > I wonder if you have any NAs in xy$x or xy$y.
> 
> The Dataframe looks like the following:
>  > xy
>    x          y
> 1 3.0  121.50951
> 2 2.0   30.61258
> 3 4.0  323.14837
> 4 8.2 3709.92471
> For testing reasons I calculated y by:
>  > y <- 2.9 x^3.4

Did you read the help page for nls that emphasizes

     *Do not use 'nls' on artificial "zero-residual" data.*

and explains why?

In any case, the best way to fit such a model to these data is to use
the plinear algorithm for partially linear models.  Either

> nls(y ~ x^exp(logpow), data=xy, start = c(logpow=0.1), alg='plinear',trace = TRUE)
2661838 :   0.1000 281.0589 
87270.38 :  0.9550582 15.4392528 
1573.086 : 1.181199 3.903272 
1.645442 : 1.222351 2.929585 
2.301169e-06 : 1.223774 2.900035 
1.483770e-11 : 1.223775 2.900000 
1.483769e-11 : 1.223775 2.900000 
Nonlinear regression model
  model:  y ~ x^exp(logpow) 
   data:  xy 
  logpow     .lin 
1.223775 2.900000 
 residual sum-of-squares:  1.483769e-11 

or

> nls(y ~ exp(log(x) * pow), data=xy, start = c(pow=1.1), alg = 'plinear', trace = TRUE)
2684530 :   1.1000 283.5010 
408828.2 :  2.045211 47.930710 
33333.73 : 2.848788 9.185464 
857.465 : 3.293871 3.622750 
1.307911 : 3.395685 2.926367 
3.723389e-06 : 3.399993 2.900044 
1.483772e-11 : 3.4 2.9 
1.483769e-11 : 3.4 2.9 
Nonlinear regression model
  model:  y ~ exp(log(x) * pow) 
   data:  xy 
 pow .lin 
 3.4  2.9 
 residual sum-of-squares:  1.483769e-11 

> >Also, I think
> >you could take logs of your equation and end up with a
> >linear expression?
> >
> if I use:
> 
> plot (x,y)
> abline lm ( log(y)~log(x),xy)
> 
> the line does not seem to fit the plotted datapoints at all. While the
> fitted exponent seems to be okay, the obtained intercept value is
> wrong.

That's because your plot is in the x,y scale and the line is in the
log(x), log(y) scale.  You need to convert back to the original scale
to put the line on the plot.

-- 
Douglas Bates                            bates at stat.wisc.edu
Statistics Department                    608/262-2598
University of Wisconsin - Madison        http://www.stat.wisc.edu/~bates/