[R] fitting models with gnls

Fri Sep 7 16:26:41 CEST 2001

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