gcm: Growth Cessation Model
In fishgrowth: Fit Growth Curves to Fish Data

gcm	R Documentation

Growth Cessation Model

Description

Fit a growth cessation model (GCM) to otoliths and/or tags.

Usage

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

gcm_curve(t, L0, rmax, k, t50)

gcm_objfun(par, data)

Arguments

`par`	a parameter list.
`data`	a data list.
`silent`	passed to `MakeADFun`.
`...`	passed to `MakeADFun`.
`t`	age (vector).
`L0`	predicted length at age 0.
`rmax`	shape parameter that determines the initial slope.
`k`	shape parameter that determines how quickly the growth curve reaches the asymptotic maximum.
`t50`	shape parameter that determines the logistic function midpoint.

Details

The main function gcm creates a model object, ready for parameter estimation. The auxiliary functions gcm_curve and gcm_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:

L0, predicted length at age 0
log_rmax, shape parameter that determines the initial slope
log_k, shape parameter that determines how quickly the growth curve reaches the asymptotic maximum
t50, shape parameter that determines the logistic function midpoint
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)

*: 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 gcm function returns a TMB model object, produced by MakeADFun.

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

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

Note

The growth cessation model (Maunder et al. 2018) predicts length at age as:

\hat L_t ~=~ L_0 ~+~ r_{\max}\!\left[\,\frac{\log\left(1+e^{-kt_{50}} \right) \;-\;\log\left(1+e^{k(t-t_{50})}\right)}{k}\;+\;t\:\right]

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

Maunder, M.N., Deriso, R.B., Schaefer, K.M., Fuller, D.W., Aires-da-Silva, A.M., Minte-Vera, C.V., and Campana, S.E. (2018). The growth cessation model: a growth model for species showing a near cessation in growth with application to bigeye tuna (Thunnus obesus). Marine Biology, 165, 76. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1007/s00227-018-3336-9")}.

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

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, gcm_curve(x, L0=5, rmax=20, k=0.15, t50=0), lty=3)

# Prepare parameters and data
init <- list(L0=5, log_rmax=log(20), log_k=log(0.15), t50=-1,
             log_sigma_min=log(3), log_sigma_max=log(3))
dat <- list(Aoto=otoliths_had$age, Loto=otoliths_had$len)
gcm_objfun(init, dat)

# Fit model
model <- gcm(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, gcm_curve(x, L0, rmax, k, t50))
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("L0", "rmax", "k", "t50", "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, gcm_curve(x, L0=20, rmax=120, k=2, t50=0), 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(L0=20, log_rmax=log(120), log_k=log(4), t50=0,
             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)
gcm_objfun(init, dat)

# Fit model
model <- gcm(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, gcm_curve(x, L0, rmax, k, t50))
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("L0", "rmax", "k", "t50", "sigma_min", "sigma_max")]
fit[-1]
head(summary(sdreport), 6)

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

# Model 3: Stepwise estimation procedure, described by Maunder et al. (2018)
# - estimate L0 and rmax using linear regression on younger fish
# - estimate k and t50 using GCM and all data, keeping L0 and rmax fixed

# Estimate L0 and rmax
plot(otoliths_skj, xlim=c(0,4), ylim=c(0,100))
fm <- lm(len~age, otoliths_skj)
abline(fm)
L0 <- coef(fm)[[1]]
rmax <- coef(fm)[[2]]

# Explore initial parameter values (k, t50, age)
t <- seq(0, 4, by=0.1)
points(t, gcm_curve(t, L0, rmax, k=3, t50=2), col="gray")
points(lenRel~I(lenRel/50), tags_skj, col=4)
points(lenRec~I(lenRel/50+liberty), tags_skj, col=3)
legend("bottomright", c("otoliths","tag releases","tac recaptures",
       "linear regression (otoliths)"), 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)

# Prepare parameters
init <- list(L0=L0, log_rmax=log(rmax), log_k=log(3), t50=2,
             log_sigma_min=log(3), log_sigma_max=log(3),
             log_age=log(tags_skj$lenRel/50))

# Fit model
map <- list(L0=factor(NA), log_rmax=factor(NA))  # fix L0 and rmax
model <- gcm(init, dat, map=map)
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, gcm_curve(x, L0, rmax, k, t50))
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("L0", "rmax", "k", "t50", "sigma_min", "sigma_max")]
fit[-1]
head(summary(sdreport), 6)

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