richardso: Richards Growth Model (Old Style)
In fishgrowth: Fit Growth Curves to Fish Data
richardso R Documentation
Richards Growth Model (Old Style)
Description
Fit a Richards growth model to otoliths and/or tags, using a traditional
parametrization.
Usage
richardso(par, data, silent = TRUE, ...)
richardso_curve(t, Linf, k, tau, b)
richardso_objfun(par, data)
Arguments
par
a parameter list.
data
a data list.
silent
passed to MakeADFun.
...
passed to MakeADFun.
t
age (vector).
Linf
asymptotic maximum length.
k
growth coefficient.
tau
location parameter.
b
shape parameter.
Details
The main function richardso creates a model object, ready for
parameter estimation. The auxiliary functions richardso_curve and
richardso_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
-
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)
*: 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 richardso function returns a TMB model object, produced by
MakeADFun.
The richardso_curve function returns a numeric vector of predicted
length at age.
The richardso_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 Tjørve and Tjørve (2010,
Eq. 4), predicts length at age as:
\hat L_t ~=~ L_\infty\left(1\,-\,\frac{1}{b}\,
e^{-k(t-\tau)}\right)^{\!b}
The variability of length at age increases linearly with length,
\sigma_L ~=~ \alpha \,+\, \beta \hat L
where the intercept is \alpha=\sigma_{\min} - \beta L_{\min}, the slope is
\beta=(\sigma_{\max}-\sigma_{\min}) / (L_{\max}-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 objective function integrates (sums) the
likelihood components from the otoliths and tags,
f ~=~ \sum_{i=1}^{N_\mathrm{oto}}\, 0.5\log(2\pi) \,+\,
\log\sigma_i \,+\, \frac{(L_i-\hat L_i)^2}{2\sigma_i^2}
~~~~ ~+\,~ \sum_{j=1}^{N_\mathrm{tag}}\, 0.5\log(2\pi) \,+\,
\log\sigma_j \,+\, \frac{(L_j-\hat L_j)^2}{2\sigma_j^2}
~~~~ ~+\,~ \sum_{k=1}^{N_\mathrm{tag}}\, 0.5\log(2\pi) \,+\,
\log\sigma_k \,+\, \frac{(L_k-\hat L_k)^2}{2\sigma_k^2}
where L_i are length measurements from the otolith data, L_j are
length measurements at tag release, and L_k are length measurements at
tag recapture. N_\mathrm{oto} is the number of fish in
the otolith data and N_\mathrm{tag} is the number of fish in the
tagging data.
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.
Tjørve, E. and Tjørve, K.M.C. (2010).
A unified approach to the Richards-model family for use in growth analyses:
Why we need only two model forms.
Journal of Theoretical Biology, 267, 417-425.
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, richardso_curve(x, Linf=100, k=0.1, tau=-1, b=1), lty=3)
# Prepare parameters and data
init <- list(log_Linf=log(100), log_k=log(0.1), tau=-1, b=1,
log_sigma_min=log(3), log_sigma_max=log(3))
dat <- list(Aoto=otoliths_had$age, Loto=otoliths_had$len)
richardso_objfun(init, dat)
# Fit model
model <- richardso(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, richardso_curve(x, Linf, k, tau, b))
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", "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, pch=16, col="#0080a010")
lines(Lhat~age, band, 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, richardso_curve(x, Linf=80, k=0.8, tau=-0.5, b=1), lty=2)
legend("bottomright", c("otoliths","tag releases","tag 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(80), log_k=log(0.8), tau=-0.5, 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)
richardso_objfun(init, dat)
# Fit model
model <- richardso(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
band <- pred_band(x, model)
areaplot::confplot(cbind(lower,upper)~age, band, xlim=c(0,4), ylim=c(0,100),
ylab="len", col="gray90")
points(len~age, otoliths_skj)
points(report$age, report$Lrel, col=4)
points(report$age+report$liberty, report$Lrec, col=3)
Lhat <- with(report, richardso_curve(x, Linf, k, tau, b))
lines(x, Lhat, lwd=2)
legend("bottomright", c("otoliths","tag releases","tag 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", "b", "sigma_min", "sigma_max")]
fit[-1]
head(summary(sdreport), 6)
#############################################################################
# Model 3: Fit to skipjack otoliths only
init <- list(log_Linf=log(80), log_k=log(0.8), tau=-0.5, b=1,
log_sigma_min=log(3), log_sigma_max=log(3))
dat <- list(Aoto=otoliths_skj$age, Loto=otoliths_skj$len)
model <- richardso(init, dat)
fit <- nlminb(model$par, model$fn, model$gr,
control=list(eval.max=1e4, iter.max=1e4))
model$report()[c("Linf", "k", "tau", "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_Linf=log(80), log_k=log(0.8), tau=-0.5, b=1,
log_sigma_min=log(3))
dat <- list(Aoto=otoliths_skj$age, Loto=otoliths_skj$len)
model <- richardso(init, dat)
fit <- nlminb(model$par, model$fn, model$gr,
control=list(eval.max=1e4, iter.max=1e4))
model$report()[c("Linf", "k", "tau", "b", "sigma_min")]
#############################################################################
# Model 5: Fit to skipjack tags only
init <- list(log_Linf=log(80), log_k=log(0.8), tau=-0.5, 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)
model <- richardso(init, dat) # using 1e3 to keep CRAN checks fast,
fit <- nlminb(model$par, model$fn, model$gr, # but try 1e4 to get
control=list(eval.max=1e3, iter.max=1e3)) # better convergence
model$report()[c("Linf", "k", "tau", "b", "sigma_min", "sigma_max")]
fishgrowth documentation built on Nov. 7, 2025, 9:07 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
| richardso | R Documentation |
Richards Growth Model (Old Style)
Description
Fit a Richards growth model to otoliths and/or tags, using a traditional parametrization.
Usage
richardso(par, data, silent = TRUE, ...)
richardso_curve(t, Linf, k, tau, b)
richardso_objfun(par, data)
Arguments
par |
a parameter list. |
data |
a data list. |
silent |
passed to |
... |
passed to |
t |
age (vector). |
Linf |
asymptotic maximum length. |
k |
growth coefficient. |
tau |
location parameter. |
b |
shape parameter. |
Details
The main function richardso creates a model object, ready for
parameter estimation. The auxiliary functions richardso_curve and
richardso_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 -
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)
*: 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 richardso function returns a TMB model object, produced by
MakeADFun.
The richardso_curve function returns a numeric vector of predicted
length at age.
The richardso_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 Tjørve and Tjørve (2010, Eq. 4), predicts length at age as:
\hat L_t ~=~ L_\infty\left(1\,-\,\frac{1}{b}\,
e^{-k(t-\tau)}\right)^{\!b}
The variability of length at age increases linearly with length,
\sigma_L ~=~ \alpha \,+\, \beta \hat L
where the intercept is \alpha=\sigma_{\min} - \beta L_{\min}, the slope is
\beta=(\sigma_{\max}-\sigma_{\min}) / (L_{\max}-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 objective function integrates (sums) the likelihood components from the otoliths and tags,
f ~=~ \sum_{i=1}^{N_\mathrm{oto}}\, 0.5\log(2\pi) \,+\,
\log\sigma_i \,+\, \frac{(L_i-\hat L_i)^2}{2\sigma_i^2}
~~~~ ~+\,~ \sum_{j=1}^{N_\mathrm{tag}}\, 0.5\log(2\pi) \,+\,
\log\sigma_j \,+\, \frac{(L_j-\hat L_j)^2}{2\sigma_j^2}
~~~~ ~+\,~ \sum_{k=1}^{N_\mathrm{tag}}\, 0.5\log(2\pi) \,+\,
\log\sigma_k \,+\, \frac{(L_k-\hat L_k)^2}{2\sigma_k^2}
where L_i are length measurements from the otolith data, L_j are
length measurements at tag release, and L_k are length measurements at
tag recapture. N_\mathrm{oto} is the number of fish in
the otolith data and N_\mathrm{tag} is the number of fish in the
tagging data.
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.
Tjørve, E. and Tjørve, K.M.C. (2010). A unified approach to the Richards-model family for use in growth analyses: Why we need only two model forms. Journal of Theoretical Biology, 267, 417-425.
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, richardso_curve(x, Linf=100, k=0.1, tau=-1, b=1), lty=3)
# Prepare parameters and data
init <- list(log_Linf=log(100), log_k=log(0.1), tau=-1, b=1,
log_sigma_min=log(3), log_sigma_max=log(3))
dat <- list(Aoto=otoliths_had$age, Loto=otoliths_had$len)
richardso_objfun(init, dat)
# Fit model
model <- richardso(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, richardso_curve(x, Linf, k, tau, b))
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", "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, pch=16, col="#0080a010")
lines(Lhat~age, band, 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, richardso_curve(x, Linf=80, k=0.8, tau=-0.5, b=1), lty=2)
legend("bottomright", c("otoliths","tag releases","tag 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(80), log_k=log(0.8), tau=-0.5, 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)
richardso_objfun(init, dat)
# Fit model
model <- richardso(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
band <- pred_band(x, model)
areaplot::confplot(cbind(lower,upper)~age, band, xlim=c(0,4), ylim=c(0,100),
ylab="len", col="gray90")
points(len~age, otoliths_skj)
points(report$age, report$Lrel, col=4)
points(report$age+report$liberty, report$Lrec, col=3)
Lhat <- with(report, richardso_curve(x, Linf, k, tau, b))
lines(x, Lhat, lwd=2)
legend("bottomright", c("otoliths","tag releases","tag 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", "b", "sigma_min", "sigma_max")]
fit[-1]
head(summary(sdreport), 6)
#############################################################################
# Model 3: Fit to skipjack otoliths only
init <- list(log_Linf=log(80), log_k=log(0.8), tau=-0.5, b=1,
log_sigma_min=log(3), log_sigma_max=log(3))
dat <- list(Aoto=otoliths_skj$age, Loto=otoliths_skj$len)
model <- richardso(init, dat)
fit <- nlminb(model$par, model$fn, model$gr,
control=list(eval.max=1e4, iter.max=1e4))
model$report()[c("Linf", "k", "tau", "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_Linf=log(80), log_k=log(0.8), tau=-0.5, b=1,
log_sigma_min=log(3))
dat <- list(Aoto=otoliths_skj$age, Loto=otoliths_skj$len)
model <- richardso(init, dat)
fit <- nlminb(model$par, model$fn, model$gr,
control=list(eval.max=1e4, iter.max=1e4))
model$report()[c("Linf", "k", "tau", "b", "sigma_min")]
#############################################################################
# Model 5: Fit to skipjack tags only
init <- list(log_Linf=log(80), log_k=log(0.8), tau=-0.5, 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)
model <- richardso(init, dat) # using 1e3 to keep CRAN checks fast,
fit <- nlminb(model$par, model$fn, model$gr, # but try 1e4 to get
control=list(eval.max=1e3, iter.max=1e3)) # better convergence
model$report()[c("Linf", "k", "tau", "b", "sigma_min", "sigma_max")]
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.