Last chance! 50% off unlimited learning
Sale ends in
Fit a Gompertz growth model to otoliths and/or tags, using a traditional parametrization.
gompertzo(par, data, silent = TRUE, ...)gompertzo_curve(t, Linf, k, tau)
gompertzo_objfun(par, data)
The gompertzo function returns a TMB model object, produced by
MakeADFun.
The gompertzo_curve function returns a numeric vector of predicted
length at age.
The gompertzo_objfun function returns the negative log-likelihood as a
single number, describing the goodness of fit of par to the
data.
a parameter list.
a data list.
passed to MakeADFun.
passed to MakeADFun.
age (vector).
asymptotic maximum length.
growth coefficient.
location parameter.
The main function gompertzo creates a model object, ready for
parameter estimation. The auxiliary functions gompertzo_curve and
gompertzo_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_Linf, asymptotic maximum length
log_k, growth coefficient
tau, location parameter
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.
Gompertz, B. (1825). On the nature of the function expressive of the law of human mortality, and on a new mode of determining the value of life contingencies. Philosophical Transactions of the Royal Society, 115, 513-583.
Ricker, W.E. (1979). Growth rates and models. In: W.S. Hoar et al. (eds.) Fish physiology 8: Bioenergetics and growth. New York: Academic Press, pp. 677-743.
The fishgrowth-package help page includes references describing
the parameter estimation method.
gcm, gompertz, gompertzo,
richards, richards, 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.
# 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, gompertzo_curve(x, Linf=73, k=0.4, tau=2), lty=3)
# Prepare parameters and data
init <- list(log_Linf=log(73), log_k=log(0.4), tau=2,
log_sigma_min=log(3), log_sigma_max=log(3))
dat <- list(Aoto=otoliths_had$age, Loto=otoliths_had$len)
gompertzo_objfun(init, dat)
# Fit model
model <- gompertzo(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, gompertzo_curve(x, Linf, k, tau))
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("Linf", "k", "tau", "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, gompertzo_curve(x, Linf=75, k=1, tau=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(log_Linf=log(75), log_k=log(1), tau=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)
gompertzo_objfun(init, dat)
# Fit model
model <- gompertzo(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, gompertzo_curve(x, Linf, k, tau))
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("Linf", "k", "tau", "sigma_min", "sigma_max")]
fit[-1]
head(summary(sdreport), 5)
#############################################################################
# Model 3: Fit to skipjack otoliths only
init <- list(log_Linf=log(75), log_k=log(1), tau=0,
log_sigma_min=log(3), log_sigma_max=log(3))
dat <- list(Aoto=otoliths_skj$age, Loto=otoliths_skj$len)
model <- gompertzo(init, dat)
fit <- nlminb(model$par, model$fn, model$gr,
control=list(eval.max=1e4, iter.max=1e4))
model$report()[c("Linf", "k", "tau", "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_Linf=log(75), log_k=log(1), tau=0,
log_sigma_min=log(3))
dat <- list(Aoto=otoliths_skj$age, Loto=otoliths_skj$len)
model <- gompertzo(init, dat)
fit <- nlminb(model$par, model$fn, model$gr,
control=list(eval.max=1e4, iter.max=1e4))
model$report()[c("Linf", "k", "tau", "sigma_min")]
#############################################################################
# Model 5: Fit to skipjack tags only
init <- list(log_Linf=log(75), log_k=log(1), tau=0,
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)
model <- gompertzo(init, dat)
fit <- nlminb(model$par, model$fn, model$gr,
control=list(eval.max=1e4, iter.max=1e4))
model$report()[c("Linf", "k", "tau", "sigma_min", "sigma_max")]
Run the code above in your browser using DataLab