Required Packages

> library(FSA)     # growmodel.sim

Preparing Data

The working directory is set, the data file (i.e., box5_4.txt) is read, and the structure is observed with,

> setwd("c://aaaWork//web//fishR//bookex//AIFFD//Box5_4")
> d <- read.table("box5_4.txt", header = TRUE)
> str(d)
'data.frame':   95 obs. of  2 variables:
 $ tl : int  388 418 438 428 539 432 444 421 438 419 ...
 $ age: int  4 4 4 5 10 4 7 4 4 4 ...

Getting Initial Values for the Model Parameters

As noted in Box 5.4, the computer algorithms for fitting non-linear models are iterative and require starting values for the parameters. An alternative to what was described in Box 5.4 for identifying starting values for the traditional von Bertalanffy equation is to plot the data and use interactive graphics to visually "fit" the von Bertalanffy model. The parameters from this process can then serve as the starting values for the more objective non-linear modeling algorithm. The growmodel.sim() function in the FSA package provides a mechanism for visually fitting the von Bertalanffy model. For this purpose, growmodel.sim() has three arguments,

Alternatively, if the data frame that contains the length and age variables is attached then the vector of observed ages can be sent to the x= argument and the vector of lengths to the y= argument.

The starting values for the spotted sucker data can be obtained by manipulating the slider bars on the graphic produced with,

> growmodel.sim("vb", tl ~ age, data = d)
Demonstration of Finding Von Bertalanffy Starting Values

From this process, I (your fitting may have resulted in different values) obtained starting values of 490 for Linfinity, 0.4 for K, and 0 for t0. For ease in the non-linear model fitting procedure described next these values should be entered into a saved list object,

> sv <- list(Linf = 490, K = 0.4, t0 = 0)

Fitting the Von Bertalanffy Model

The von Bertalanffy model can be placed into a formula object as follows,

> vbmdl <- tl ~ Linf * (1 - exp(-K * (age - t0)))

where you must make sure that tl and age are the exact names of the variables in your data frame and Linf, K, and t0 are the exact names of the parameters that you stored in your starting values list.

The non-linear model fitting procedure in R is implemented with nls(), which requires the model formula, the list of starting values, and the data frame as arguments. The results of the model fit should be stored into an object. The von Bertalanffy model is fit to the spotted sucker data with,

> ssvb <- nls(vbmdl, start = sv, data = d)

The coefficient estimates and the correlations among coefficient estimates are extracted with,

> summary(ssvb, correlation = TRUE)
Formula: tl ~ Linf * (1 - exp(-K * (age - t0)))

Parameters:
     Estimate Std. Error t value Pr(>|t|)
Linf 516.2333    48.7987  10.579   <2e-16
K      0.1900     0.1194   1.591    0.115
t0    -4.5422     3.4083  -1.333    0.186

Residual standard error: 18.11 on 92 degrees of freedom

Correlation of Parameter Estimates:
   Linf  K
K  -0.99
t0 -0.96  0.99

Number of iterations to convergence: 31
Achieved convergence tolerance: 9.94e-06

which closely match the results in Box 5.4

Warning One should pay very close attention to the correlations among parameter estimates shown in the results above. These values illustrate VERY strong correlations among the parameters which hinders interpretation as many other parameter triplets would provide nearly the same model fit (this is illustrated with the lack of change in SS for the last 15 or so iterations of the model fitting as shown in Box 5.4).

Finally, one can plot the model fit with,

> plot(tl ~ age, data = d, pch = 16, xlab = "Age", ylab = "Total Length (mm)", xlim = c(0, 11), ylim = c(0, 550))
> ages <- seq(0, 11, 0.1)
> preds <- predict(ssvb, data.frame(age = ages))
> lines(preds ~ ages, lwd = 2, col = "red")
Box5_4a-009.png
Warning The predicted length for age-0 fish suggests issues with the data or the use of this model with these data.

Comparison Among Groups

The authors of Chapter 5 hint at comparisons of von Bertalanffy model parameters between groups but do not provide an example of such an analysis. An example of this analysis in R can be found in my Von Bertalanffy vignette.