[R] fitting models with gnls

[Douglas Bates]

"Antonio Olinto" <aolinto at bignet.com.br> writes:

> Dear R-list members,

> Some months ago I wrote a message on the usage of gnls (nlme
> library) and here I come again.

Sorry about the delay.  As some people have noticed, I am seriously
over-committed right now and it is difficult to find time to look at
problems with the R version of nlme.

As Peter Dalgaard said, the problem is in an assignment of the form
  fact[] <- value
where fact is a factor.  The particular line that is giving the
problem is in the gnls function itself
        grpShrunk[] <- grpShrunk[revOrderShrunk]

If you remove the first set of brackets, so it reads
        grpShrunk <- grpShrunk[revOrderShrunk]
you can fit your model in R.

> growth.gnls <- gnls(lt ~ Linf*(1-exp(-K*(age-t0))),
+   data=growth.dat, params= Linf +K + t0 ~ 1, start=list(Linf=500,K=0.2,t0=0),
+   control = list(returnObject = T), corr = corAR1(form=~fish|age))
> summary(growth.gnls)
Generalized nonlinear least squares fit
  Model: lt ~ Linf * (1 - exp(-K * (age - t0))) 
  Data: growth.dat 
       AIC      BIC    logLik
  757.3853 770.4112 -373.6927

Correlation Structure: AR(1)
 Formula: ~fish | age 
 Parameter estimate(s):
       Phi 
0.04672911 

Coefficients:
        Value Std.Error  t-value p-value
Linf 509.7454  8.278899 61.57164  <.0001
K      0.1988  0.007966 24.95591  <.0001
t0     0.2987  0.046531  6.41996  <.0001

 Correlation: 
   Linf   K     
K  -0.963       
t0 -0.647  0.783

Standardized residuals:
       Min         Q1        Med         Q3        Max 
-1.7863143 -0.5340824 -0.1925364  0.3771429  2.5162105 

Residual standard error: 10.32122 
Degrees of freedom: 100 total; 97 residual

I agree with Peter that this is best fixed by changing [<-.factor in
R.  The temporary fix is to create a private copy of gnls and make the
change described above.  I would prefer not to make a permanent change
like that in the nlme package.

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._