Fitting von Bertalanffy Growth Function

The von Bertalanffy Growth Function describes the body growth in length of most fish, i.e.

\[ L_t = L_\inf (1 - \exp(-K(t-t_0))) \]

where \( L_t \) is the body length at age \( t \), \( L_\inf \) is the asymptotic length, \( K \) is the growth rate coefficient, and \( t_0 \) is the age at \( L=0 \). All you need to estimate the three parameters of the model are observed data-pairs of length and age, i.e.,

DATA_SECTION

age <- 1:12
lt <- c(10, 19, 24, 28, 32, 34, 38, 40, 41, 42, 42, 43)
dat <- data.frame(age, lt)
dat
##    age lt
## 1    1 10
## 2    2 19
## 3    3 24
## 4    4 28
## 5    5 32
## 6    6 34
## 7    7 38
## 8    8 40
## 9    9 41
## 10  10 42
## 11  11 42
## 12  12 43

EXPLORATORY_DATA_SECTION

To explore the data, you need the ggplot2 library

library(ggplot2)

A plot of the length against age is describing the observed growth for a fish population, i.e.

qplot(age, lt, data = dat, xlab = "Age (yrs)", ylab = "Length (cm)", xlim = c(0, 
    13), ylim = c(0, 50))

plot of chunk ggplot2ex

PARAMETERS_SECTION

You need to guess initial parameters, which can be in a vector called \( theta \), i.e.,

theta <- c(45, 0.2, -0.5)

PROCEDURE_SECTION

You need a function to estimate the parameters with \( optim \), i.e.,

SSQ <- function(theta, x) {
    Linf <- theta[1]
    K <- theta[2]
    t0 <- theta[3]
    epsilon <- rep(0, length(age))
    lpred <- rep(0, length(age))
    for (i in 1:length(age)) {
        lpred[i] <- Linf * (1 - exp(-K * (age[i] - t0)))
        epsilon[i] <- (lt[i] - lpred[i])^2
    }
    ssq <- sum(epsilon)
    return(ssq)
}

Once ran the previous function, then you can use \( optim \) to solve for the parameters, i.e.,

out <- optim(theta, fn = SSQ, method = "BFGS", x = age, hessian = TRUE)
out$V <- solve(out$hessian)  #solve the hessian
out$S <- sqrt(diag(out$V))  #Standard Error
out$R <- out$V/(out$S %o% out$S)  #Correlation

REPORT_SECTION

To known the parameters than have been estimated, you can list the following

out$par
## [1] 45.78855  0.23959 -0.08253
out$S
## [1] 0.82927 0.01577 0.11253

Finally, you can add the following to see the fitted curve to observed data:

lp <- out$par[1] * (1 - exp(-out$par[2] * (age - out$par[3])))
plot(age, lt, xlab = "Age (yrs)", ylab = "Length (cm)", xlim = c(0, 13), ylim = c(0, 
    50))
lines(sort(age), sort(lp), col = 3, lwd = 2)

plot of chunk unnamed-chunk-5

SUGGESTIONS

Who Am I?

Luis A. Cubillos

Fisheries Biologist