[R] Question about fitting power

[Ruben Roa]

> -----Original Message-----
> From:	r-help-bounces at stat.math.ethz.ch [SMTP:r-help-bounces at stat.math.ethz.ch] On Behalf Of Ana Quitério
> Sent:	Wednesday, January 25, 2006 3:33 PM
> To:	r-help at stat.math.ethz.ch
> Subject:	[R]  Question about fitting power
> 
>From: Ana Quitério
> > 
> > Hi R users
> > 
> > I'm trying to fit a model y=ax^b.
> > 
> > I know if I made ln(y)=ln(a)+bln(x)  this is a linear regression.
> > 
> > But I obtain differente results with nls() and lm()
> > 
> > My commands are:
> nls(CV ~a*Est^b, data=limiares, start =list(a=100,b=0), trace = TRUE)  
> for nonlinear regression, and :  
> lm(ln_CV~ln_Est, data=limiares) 
> for linear regression
> > 
> > Nonlinear regression model:   a=738.2238151  and    b=-0.3951013 
> > Linear regression:   Coefficients:
> >                   Estimate    Std. Error    t value Pr(>|t|)    
> > (Intercept)  7.8570224  0.0103680   757.8   <2e-16 ***
> ln_Est      -0.5279412  0.0008658  -609.8   <2e-16 ***
> > 
> > I think it should be   a=exp("(Intercept)  ") = 
> > exp(7.8570224) = 2583.815 and b=ln_Est
> > 
> > Probably I'm wrong, but why??
> > 
> > Thanks in advance.
> 
> Ana Quiterio

In addition to what Andy said, maybe your nls() model is not natural because a power 
statistical model should produce multiplicative outcomes. An alternative is to try the 
nonlinear power model with multiplicative deviates assuming lognormality (which is
what the linear model does but without the need for back-transformation), with nlm(), 
for instance:

fn<-function(p){
+  CV_mod=p[1]*Est^p[2];
+  squdiff=(log(limiares$CV)-log(CV_mod))^2;
+  lik=(length(limiares$CV)/2)*log(sum(squdiff)/length(limiares$CV))
+  }
CV.lik<-nlm(fn,p=c(100,0),hessian=TRUE)
fec.lik
covmat<-solve(CV.lik$hessian)
covmat

Ruben