vonbert: Von Bertalanffy Growth Model

View source: R/vonbert.R

vonbertR Documentation

Von Bertalanffy Growth Model

Description

Fit a von Bertalanffy growth model to otoliths and/or tags, using the Schnute-Fournier parametrization.

Usage

vonbert(par, data, silent = TRUE, ...)

vonbert_curve(t, L1, L2, k, t1, t2)

vonbert_objfun(par, data)

Arguments

par

a parameter list.

data

a data list.

silent

passed to MakeADFun.

...

passed to MakeADFun.

t

age (vector).

L1

predicted length at age t1.

L2

predicted length at age t2.

k

growth coefficient.

t1

age where predicted length is L1.

t2

age where predicted length is L2.

Details

The main function vonbert creates a model object, ready for parameter estimation. The auxiliary functions vonbert_curve and vonbert_objfun are called by the main function to calculate the regression curve and objective function value. The user can also call the auxiliary functions directly for plotting and model exploration.

The par list contains the following elements:

  • log_L1, predicted length at age t1

  • log_L2, predicted length at age t2

  • log_k, growth coefficient

  • log_sigma_min, growth variability at the shortest observed length in the data

  • log_sigma_max (*), growth variability at the longest observed length in the data

  • log_age (*), age at release of tagged individuals (vector)

*: The parameter log_sigma_max can be omitted to estimate growth variability that does not vary with length. The parameter vector log_age can be omitted to fit to otoliths only.

The data list contains the following elements:

  • Aoto (*), age from otoliths (vector)

  • Loto (*), length from otoliths (vector)

  • Lrel (*), length at release of tagged individuals (vector)

  • Lrec (*), length at recapture of tagged individuals (vector)

  • liberty (*), time at liberty of tagged individuals in years (vector)

  • t1, age where predicted length is L1

  • t2, age where predicted length is L2

*: The data vectors Aoto and Loto can be omitted to fit to tagging data only. The data vectors Lrel, Lrec, and liberty can be omitted to fit to otoliths only.

Value

The vonbert function returns a TMB model object, produced by MakeADFun.

The vonbert_curve function returns a numeric vector of predicted length at age.

The vonbert_objfun function returns the negative log-likelihood as a single number, describing the goodness of fit of par to the data.

Note

The Schnute-Fournier parametrization used in vonbert reduces parameter correlation and improves convergence reliability compared to the traditional parametrization used in vonberto. Therefore, the vonbert parametrization can be recommended for general usage, as both parametrizations produce the same growth curve. However, there can be some use cases where the traditional parametrization (Linf, k, t0) is preferred over the Schnute-Fournier parametrization (L1, L2, k).

The von Bertalanffy (1938) growth model, as parametrized by Schnute and Fournier (1980), predicts length at age as:

\hat L_t ~=~ L_1 \;+\; (L_2-L_1)\, \frac{1\,-\,e^{-k(t-t_1)}}{\,1\,-\,e^{-k(t_2-t_1)}\,}

The variability of length at age increases linearly with length,

\sigma_L ~=~ \alpha \,+\, \beta \hat L

where the slope is \beta=(\sigma_{\max}-\sigma_{\min}) / (L_{\max}-L_{\min}), the intercept is \alpha=\sigma_{\min} - \beta L_{\min}, and L_{\min} and L_{\max} are the shortest and longest observed lengths in the data. Alternatively, growth variability can be modelled as a constant \sigma_L=\sigma_{\min} that does not vary with length, by omitting log_sigma_max from the parameter list (see above).

The negative log-likelihood is calculated by comparing the observed and predicted lengths:

  nll_Loto <- -dnorm(Loto, Loto_hat, sigma_Loto, TRUE)
  nll_Lrel <- -dnorm(Lrel, Lrel_hat, sigma_Lrel, TRUE)
  nll_Lrec <- -dnorm(Lrec, Lrec_hat, sigma_Lrec, TRUE)
  nll <- sum(nll_Loto) + sum(nll_Lrel) + sum(nll_Lrec)

References

von Bertalanffy, L. (1938). A quantitative theory of organic growth. Human Biology, 10, 181-213. https://www.jstor.org/stable/41447359.

Schnute, J. and Fournier, D. (1980). A new approach to length-frequency analysis: Growth structure. Canadian Journal of Fisheries and Aquatic Science, 37, 1337-1351. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1139/f80-172")}.

The fishgrowth-package help page includes references describing the parameter estimation method.

See Also

gcm, gompertz, gompertzo, richards, richardso, schnute3, vonbert, and vonberto are alternative growth models.

pred_band calculates a prediction band for a fitted growth model.

otoliths_had, otoliths_skj, and tags_skj are example datasets.

fishgrowth-package gives an overview of the package.

Examples

# Model 1: Fit to haddock otoliths

# Explore initial parameter values
plot(len~age, otoliths_had, xlim=c(0,18), ylim=c(0,105), pch=16,
     col="#0080a010")
x <- seq(1, 18, 0.1)
lines(x, vonbert_curve(x, L1=18, L2=67, k=0.1, t1=1, t2=10), lty=3)

# Prepare parameters and data
init <- list(log_L1=log(18), log_L2=log(67), log_k=log(0.1),
             log_sigma_min=log(3), log_sigma_max=log(3))
dat <- list(Aoto=otoliths_had$age, Loto=otoliths_had$len, t1=1, t2=10)
vonbert_objfun(init, dat)

# Fit model
model <- vonbert(init, dat)
fit <- nlminb(model$par, model$fn, model$gr,
              control=list(eval.max=1e4, iter.max=1e4))
report <- model$report()
sdreport <- sdreport(model)

# Plot results
Lhat <- with(report, vonbert_curve(x, L1, L2, k, t1, t2))
lines(x, Lhat, lwd=2, col=2)
legend("bottomright", c("initial curve","model fit"), col=c(1,2), lty=c(3,1),
       lwd=c(1,2), bty="n", inset=0.02, y.intersp=1.25)

# Model summary
report[c("L1", "L2", "k", "sigma_min", "sigma_max")]
fit[-1]
summary(sdreport)

# Plot 95% prediction band
band <- pred_band(x, model)
areaplot::confplot(cbind(lower,upper)~age, band, xlim=c(0,18), ylim=c(0,100),
         ylab="len", col="mistyrose")
points(len~age, otoliths_had, xlim=c(0,18), ylim=c(0,100),
       pch=16, col="#0080a010")
lines(x, Lhat, lwd=2, col=2)
lines(lower~age, band, lty=1, lwd=0.5, col=2)
lines(upper~age, band, lty=1, lwd=0.5, col=2)

#############################################################################

# Model 2: Fit to skipjack otoliths and tags

# Explore initial parameter values
plot(len~age, otoliths_skj, xlim=c(0,4), ylim=c(0,100))
x <- seq(0, 4, 0.1)
points(lenRel~I(lenRel/60), tags_skj, col=4)
points(lenRec~I(lenRel/60+liberty), tags_skj, col=3)
lines(x, vonbert_curve(x, L1=25, L2=75, k=0.8, t1=0, t2=4), lty=2)
legend("bottomright", c("otoliths","tag releases","tac recaptures",
       "initial curve"), col=c(1,4,3,1), pch=c(1,1,1,NA), lty=c(0,0,0,2),
       lwd=c(1.2,1.2,1.2,1), bty="n", inset=0.02, y.intersp=1.25)

# Prepare parameters and data
init <- list(log_L1=log(25), log_L2=log(75), log_k=log(0.8),
             log_sigma_min=log(3), log_sigma_max=log(3),
             log_age=log(tags_skj$lenRel/60))
dat <- list(Aoto=otoliths_skj$age, Loto=otoliths_skj$len,
            Lrel=tags_skj$lenRel, Lrec=tags_skj$lenRec,
            liberty=tags_skj$liberty, t1=0, t2=4)
vonbert_objfun(init, dat)

# Fit model
model <- vonbert(init, dat)
fit <- nlminb(model$par, model$fn, model$gr,
              control=list(eval.max=1e4, iter.max=1e4))
report <- model$report()
sdreport <- sdreport(model)

# Plot results
plot(len~age, otoliths_skj, xlim=c(0,4), ylim=c(0,100))
points(report$age, report$Lrel, col=4)
points(report$age+report$liberty, report$Lrec, col=3)
Lhat <- with(report, vonbert_curve(x, L1, L2, k, t1, t2))
lines(x, Lhat, lwd=2)
legend("bottomright", c("otoliths","tag releases","tac recaptures",
       "model fit"), col=c(1,4,3,1), pch=c(1,1,1,NA), lty=c(0,0,0,1),
       lwd=c(1.2,1.2,1.2,2), bty="n", inset=0.02, y.intersp=1.25)

# Model summary
report[c("L1", "L2", "k", "sigma_min", "sigma_max")]
fit[-1]
head(summary(sdreport), 5)

#############################################################################

# Model 3: Fit to skipjack otoliths only

init <- list(log_L1=log(25), log_L2=log(75), log_k=log(0.8),
             log_sigma_min=log(3), log_sigma_max=log(3))
dat <- list(Aoto=otoliths_skj$age, Loto=otoliths_skj$len, t1=0, t2=4)
model <- vonbert(init, dat)
fit <- nlminb(model$par, model$fn, model$gr,
              control=list(eval.max=1e4, iter.max=1e4))
model$report()[c("L1", "L2", "k", "sigma_min", "sigma_max")]

#############################################################################

# Model 4: Fit to skipjack otoliths only,
# but now estimating constant sigma instead of sigma varying by length

# We do this by omitting log_sigma_max
init <- list(log_L1=log(25), log_L2=log(75), log_k=log(0.8),
             log_sigma_min=log(3))
dat <- list(Aoto=otoliths_skj$age, Loto=otoliths_skj$len, t1=0, t2=4)
model <- vonbert(init, dat)
fit <- nlminb(model$par, model$fn, model$gr,
              control=list(eval.max=1e4, iter.max=1e4))
model$report()[c("L1", "L2", "k", "sigma_min")]

#############################################################################

# Model 5: Fit to skipjack tags only

init <- list(log_L1=log(25), log_L2=log(75), log_k=log(0.8),
             log_sigma_min=log(3), log_sigma_max=log(3),
             log_age=log(tags_skj$lenRel/60))
dat <- list(Lrel=tags_skj$lenRel, Lrec=tags_skj$lenRec,
            liberty=tags_skj$liberty, t1=0, t2=4)
model <- vonbert(init, dat)
fit <- nlminb(model$par, model$fn, model$gr,
              control=list(eval.max=1e4, iter.max=1e4))
model$report()[c("L1", "L2", "k", "sigma_min", "sigma_max")]


fishgrowth documentation built on April 11, 2025, 5:52 p.m.