[R] Question about fitting power
Ruben Roa
RRoa at fisheries.gov.fk
Wed Jan 25 19:05:52 CET 2006
> -----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
More information about the R-help
mailing list