richards {fishgrowth} | R Documentation |
Richards Growth Model
Description
Fit a Richards growth model to otoliths and/or tags, using the Schnute parametrization.
Usage
richards(par, data, silent = TRUE, ...)
richards_curve(t, L1, L2, k, b, t1, t2)
richards_objfun(par, data)
Arguments
par |
a parameter list. |
data |
a data list. |
silent |
passed to |
... |
passed to |
t |
age (vector). |
L1 |
predicted length at age |
L2 |
predicted length at age |
k |
growth coefficient. |
b |
shape parameter. |
t1 |
age where predicted length is |
t2 |
age where predicted length is |
Details
The main function richards
creates a model object, ready for parameter
estimation. The auxiliary functions richards_curve
and
richards_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 aget1
-
log_L2
, predicted length at aget2
-
log_k
, growth coefficient -
b
, shape 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) -
t1
, age where predicted length isL1
-
t2
, age where predicted length isL2
*: 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 richards
function returns a TMB model object, produced by
MakeADFun
.
The richards_curve
function returns a numeric vector of predicted
length at age.
The richards_objfun
function returns the negative log-likelihood as a
single number, describing the goodness of fit of par
to the
data
.
Note
The Schnute parametrization used in richards
reduces parameter
correlation and improves convergence reliability compared to the traditional
parametrization used in richardso
. Therefore, the
richards
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
,
tau
, b
) is preferred over the Schnute parametrization
(L1
, L2
, k
, b
).
The Richards (1959) growth model, as parametrized by Schnute (1981, Eq. 15), predicts length at age as:
\hat L_t ~=~ \left[\:L_1^b\;+\;(L_2^b-L_1^b)\,
\frac{1-e^{-k(t-t_1)}}{1-e^{-k(t_2-t_1)}}\,\right]^{1/b}
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
Richards, F.J. (1959). A flexible growth function for empirical use. Journal of Experimental Botany, 10, 290-300. https://www.jstor.org/stable/23686557.
Schnute, J. (1981). A versatile growth model with statistically stable parameters. Canadian Journal of Fisheries and Aquatic Science, 38, 1128-1140. doi:10.1139/f81-153.
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, richards_curve(x, L1=18, L2=67, k=0.1, b=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), b=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)
richards_objfun(init, dat)
# Fit model
model <- richards(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, richards_curve(x, L1, L2, k, b, 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", "b", "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, richards_curve(x, L1=25, L2=75, k=0.8, b=1, 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), b=1,
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)
richards_objfun(init, dat)
# Fit model
model <- richards(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, richards_curve(x, L1, L2, k, b, 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", "b", "sigma_min", "sigma_max")]
fit[-1]
head(summary(sdreport), 6)
#############################################################################
# Model 3: Fit to skipjack otoliths only
init <- list(log_L1=log(25), log_L2=log(75), log_k=log(0.8), b=1,
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 <- richards(init, dat)
fit <- nlminb(model$par, model$fn, model$gr,
control=list(eval.max=1e4, iter.max=1e4))
model$report()[c("L1", "L2", "k", "b", "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), b=1,
log_sigma_min=log(3))
dat <- list(Aoto=otoliths_skj$age, Loto=otoliths_skj$len, t1=0, t2=4)
model <- richards(init, dat)
fit <- nlminb(model$par, model$fn, model$gr,
control=list(eval.max=1e4, iter.max=1e4))
model$report()[c("L1", "L2", "k", "b", "sigma_min")]
#############################################################################
# Model 5: Fit to skipjack tags only
init <- list(log_L1=log(25), log_L2=log(75), log_k=log(0.8), b=1,
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 <- richards(init, dat)
fit <- nlminb(model$par, model$fn, model$gr,
control=list(eval.max=1e4, iter.max=1e4))
model$report()[c("L1", "L2", "k", "b", "sigma_min", "sigma_max")]