thr3ads.net - R help - [R] fitting models with gnls [Sep 2001]

If this information is useful, please help other people find it:
Share via:

Antonio Olinto

2001-Sep-07 14:26 UTC

[R] fitting models with gnls

Dear R-list members,

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

Let me give an example:

I have a 10 year length-at-age data set of 10 fishes (see growth.dat at the end
of this message) and I want to fit a von Bertalanffy growth model, Li=
Linf*(1-exp(-k*(ti-t0))) where Li = length at age i, Linf= asymptotic length, k=
curvature parameter, ti= age i, and t0= age of length 0.

I'd solved the problem using nls but the points are not independent. The
length-at-age data are correlated for each fish because I'm following the
growth history of each specimen.

Trying to use gnls I wrote to R-list and to Mr. Jose Pinheiro who kindly advised
me on this matter. The problem was that using S-Plus he could fit the model and
I, using R with the same commands, got an error message:
> 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))


Error in "[<-.factor"(*tmp*, , value = grpShrunk[revOrderShrunk]) :
        Argument "i" is missing, with no default

If I use nls function I can estimate the parameters:
> growth.nls <- nls(lt ~ Linf*(1-exp(-K*(age-t0))),+ data=growth.dat, start=list(Linf=500,K=0.2,t0=0))
> summary(growth.nls)
Formula: lt ~ Linf * (1 - exp(-K * (age - t0)))

Parameters:
      Estimate Std. Error t value Pr(>|t|)    
Linf 509.72841    7.93394  64.247  < 2e-16 ***
K      0.19888    0.00764  26.031  < 2e-16 ***
t0     0.29873    0.04462   6.696  1.4e-09 ***
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` '
1

Residual standard error: 10.32 on 97 degrees of freedom

Correlation of Parameter Estimates:
      Linf      K
K  -0.9631       
t0 -0.6472 0.7827

Well, I would be grateful to receive any help.

Best regards,

Antonio Olinto
Fisheries Institute
Sao Paulo - BRASIL

-------------
> growth.dat
    fish age    lt
1      1   1  65.3
2      1   2 144.1
3      1   3 208.6
4      1   4 261.4
5      1   5 304.7
6      1   6 340.1
7      1   7 369.1
8      1   8 392.8
9      1   9 412.2
10     1  10 428.1
11     2   1  68.5
12     2   2 149.8
13     2   3 215.0
14     2   4 267.3
15     2   5 309.3
16     2   6 343.0
17     2   7 370.1
18     2   8 391.8
19     2   9 409.2
20     2  10 423.2
21     3   1  60.4
22     3   2 134.4
23     3   3 196.3
24     3   4 248.0
25     3   5 291.1
26     3   6 327.2
27     3   7 357.3
28     3   8 382.5
29     3   9 403.5
30     3  10 421.0
31     4   1  64.1
32     4   2 142.2
33     4   3 206.7
34     4   4 260.0
35     4   5 304.1
36     4   6 340.6
37     4   7 370.8
38     4   8 395.8
39     4   9 416.4
40     4  10 433.5
41     5   1  71.4
42     5   2 156.0
43     5   3 223.9
44     5   4 278.5
45     5   5 322.2
46     5   6 357.3
47     5   7 385.5
48     5   8 408.1
49     5   9 426.3
50     5  10 440.8
51     6   1  66.6
52     6   2 147.0
53     6   3 212.8
54     6   4 266.7
55     6   5 310.8
56     6   6 346.9
57     6   7 376.5
58     6   8 400.7
59     6   9 420.5
60     6  10 436.7
61     7   1  64.5
62     7   2 143.0
63     7   3 207.9
64     7   4 261.5
65     7   5 305.9
66     7   6 342.6
67     7   7 373.0
68     7   8 398.1
69     7   9 418.8
70     7  10 436.0
71     8   1  66.0
72     8   2 146.3
73     8   3 212.7
74     8   4 267.6
75     8   5 313.0
76     8   6 350.6
77     8   7 381.6
78     8   8 407.3
79     8   9 428.5
80     8  10 446.1
81     9   1  74.2
82     9   2 162.3
83     9   3 232.9
84     9   4 289.6
85     9   5 335.1
86     9   6 371.6
87     9   7 400.9
88     9   8 424.4
89     9   9 443.3
90     9  10 458.5
91    10   1  62.2
92    10   2 138.4
93    10   3 202.1
94    10   4 255.3
95    10   5 299.7
96    10   6 336.8
97    10   7 367.8
98    10   8 393.7
99    10   9 415.3
100   10  10 433.4
-------------- next part --------------
An HTML attachment was scrubbed...
URL:
https://stat.ethz.ch/pipermail/r-help/attachments/20010907/d001a160/attachment.html

Peter Dalgaard BSA

2001-Sep-07 15:40 UTC

head link

[R] fitting models with gnls

"Antonio Olinto" <aolinto at bignet.com.br> writes:
> Error in "[<-.factor"(*tmp*, , value =
grpShrunk[revOrderShrunk]) :
>         Argument "i" is missing, with no default

Hmm. That probably comes from an assignment f[]<-something where f is
a factor. The same thing is seen with

f<-factor(1:3)
g<-factor(3:1)
f[]<-g

This might count as a bug in R. You could try fixing [<-.factor
yourself. The offending line would be

  x[i] <- m

and might be happier as

  if (missing(i)) x[] <- m else x[i] <- m

or

  x[if (missing(i)) TRUE else i] <- m

-- 
   O__  ---- Peter Dalgaard             Blegdamsvej 3  
  c/ /'_ --- Dept. of Biostatistics     2200 Cph. N   
 (*) \(*) -- University of Copenhagen   Denmark      Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk)             FAX: (+45) 35327907
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

Douglas Bates

2001-Sep-07 17:17 UTC

head link

[R] fitting models with gnls

"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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

Antonio Olinto

2001-Sep-07 19:24 UTC

head link

[R] fitting models with gnls

Dear Bates and Dalgaard,

Thanks for your attention. Removing the brackets I could fit the model.

Sincerely,

Antonio Olinto

----- Original Message -----
From: "Douglas Bates" <bates at stat.wisc.edu>
To: "Antonio Olinto" <aolinto at bignet.com.br>
Cc: "R-help" <r-help at stat.math.ethz.ch>
Sent: Friday, September 07, 2001 2:17 PM
Subject: Re: [R] fitting models with gnls

> "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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._

Seemingly Similar Threads

Search for more seemingly similar threads

R help - Sep 2001 - fitting models with gnls

[R] fitting models with gnls

[R] fitting models with gnls

[R] fitting models with gnls

[R] fitting models with gnls

Seemingly Similar Threads