Introduction

Fish have indeterminate growth, which means they don’t stop growing throughout their lives, but their growth tapers off with age towards what fisheries biologist call their asymptotic size \(L_\infty\). To calculate this size, one must fit growth models to data of length to age relationships. It is important, to make the right predictions of \(L_\infty\), that the right curve is chosen. Today we will work with two growth curvres:

The Von Bertalanffy curve, which assumes a linear decrease in growth with age.

\[ L(x) = L_\infty \cdot (1-e^{-k\cdot(x-L_0)}) \] where \(L(x)\) is the length at age \(x\), \(L_\infty\) is the assymptotic growth rate, \(k\) is the growth rate, and \(L_0\) is the initial size at sage 0.

The Logistic curve, which assumes an increase and then a decrease of growth with age.

\[ L(x) = L_\infty \cdot (1+e^{-k\cdot(x-flip)})^{-1} \] where the only parameter that is different here is \(flip\), the inflexion point of the growth curve (when growth rate switches from increasing to decreasing with age)

Before we can do any model fitting, we need the data:

data <- read.table('../data/rockfish.txt', header = TRUE)
head(data)
##   age len
## 1   2  18
## 2   3  15
## 3   3  17
## 4   3  14
## 5   3  15
## 6   3  26

It’s a simple two-column dataset with age ans length of the captured fish

We can plot it to have a better feel for it.

plot(len ~ age, data, xlab = 'Age', ylab = 'Length (cm)', pch = 16, cex = 0.7)

Building likelihood functions for growth curves

Let’s start with Von Bertalanffy and make a function that takes any set of paraleters + some data and yields the likelihood of those parameters.

vonbert <- function(Loo, L0, k, sigma){
  
  L <- data$len
  x <- data$age
  
  L.pred <- Loo*(1-exp(-k*(x-L0)))
  
  lik <- dnorm(L, L.pred, sigma,  log = TRUE)
  return(-sum(lik))
}

Let’s try it out

vonbert(Loo = 50, L0 = 10, k = 0.2, sigma = 5)
## [1] 27005.95

That’s the negative log-likelikood for those parameters!

Finding the Maximum Likelihood Estimate

Now that we have a likelihood function that works, we need to find what parameters maximize that likelihood. Package stats4 has a nice function mle to do that.

library(stats4)

Let’s try with the vonBertalanffy growth.

fit.vb <- mle(vonbert, 
              start = list(Loo=50, L0=10,k=0.2,sigma=0.5),
              method="L-BFGS-B",
              upper = list(Loo=100, L0=30,k=2,sigma=10), 
              lower = list(Loo=1, L0=-10,k=-1,sigma=0.1))
         

fit.vb
## 
## Call:
## mle(minuslogl = vonbert, start = list(Loo = 50, L0 = 10, k = 0.2, 
##     sigma = 0.5), method = "L-BFGS-B", upper = list(Loo = 100, 
##     L0 = 30, k = 2, sigma = 10), lower = list(Loo = 1, L0 = -10, 
##     k = -1, sigma = 0.1))
## 
## Coefficients:
##         Loo          L0           k       sigma 
## 53.57623278 -3.80550672  0.05503021  4.42318387