# Gompertz Growth Model Fit a Gompertz growth model to otoliths and/or tags, using the Schnute parametrization. ```r gompertz(par, data, silent = TRUE, ...) gompertz_curve(t, L1, L2, k, t1, t2) gompertz_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`. ## Returns The `gompertz` function returns a TMB model object, produced by `MakeADFun`. The `gompertz_curve` function returns a numeric vector of predicted length at age. The `gompertz_objfun` function returns the negative log-likelihood as a single number, describing the goodness of fit of `par` to the `data`. ## Details The main function `gompertz` creates a model object, ready for parameter estimation. The auxiliary functions `gompertz_curve` and `gompertz_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. ## Note The Schnute parametrization used in `gompertz` reduces parameter correlation and improves convergence reliability compared to the traditional parametrization used in `gompertzo`. Therefore, the `gompertz` 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`) is preferred over the Schnute parametrization (`L1`, `L2`, `k`). Gompertz is a special case of the Richards (1959) model, where $b=0$. If the best model fit of a `richards` model to a particular dataset involves a very small estimated value of $b$, then the `gompertz` model offers a preferable parametrization, as it produces the same curve using fewer parameters. The Gompertz (1825) growth model, as parametrized by Schnute (1981, Eq. 16) predicts length at age as: $$ \hat L_t ~=~ L_1\exp\!\left[\,\log(L_2/L_1)\,\frac{1-e^{-k(t-t_1)}}{1-e^{-k(t_2-t_1)}}\,\right]Lhat = L1 * exp(log(L2/L1) * (1-exp(-k*(t-t1))) / (1-exp(-k*(t2-t1)))) $$ The variability of length at age increases linearly with length, $$ \sigma_L ~=~ \alpha \,+\, \beta \hat Lsigma = alpha + beta * Lhat $$ where the intercept is c("$alpha =\n$", "$sigma_min - beta * L_min$"), the slope is c("$beta =\n$", "$(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}f= sum(0.5*log(2*pi) + log(sigma_i) + ((L_i-Lhat_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(0.5*log(2*pi) + log(sigma_j) + ((L_j-Lhat_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}+ sum(0.5*log(2*pi) + log(sigma_k) + ((L_k-Lhat_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. ## Examples ```r # 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, gompertz_curve(x, L1=15, L2=70, k=0.4, t1=1, t2=10), lty=3) # Prepare parameters and data init <- list(log_L1=log(20), log_L2=log(70), 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) gompertz_objfun(init, dat) # Fit model model <- gompertz(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, gompertz_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, 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, gompertz_curve(x, L1=28, L2=74, k=1, t1=0, t2=4), 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_L1=log(28), log_L2=log(74), log_k=log(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) gompertz_objfun(init, dat) # Fit model model <- gompertz(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, gompertz_curve(x, L1, L2, k, t1, t2)) 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("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(28), log_L2=log(74), log_k=log(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 <- gompertz(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(28), log_L2=log(74), log_k=log(1), log_sigma_min=log(3)) dat <- list(Aoto=otoliths_skj$age, Loto=otoliths_skj$len, t1=0, t2=4) model <- gompertz(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(28), log_L2=log(74), log_k=log(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 <- gompertz(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")] ``` ## References 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. Schnute, J. (1981). A versatile growth model with statistically stable parameters. **Canadian Journal of Fisheries and Aquatic Science**, 38 , 1128-1140. tools:::Rd_expr_doi("10.1139/f81-153") . The `fishgrowth-package` help page includes references describing the parameter estimation method. ## See Also `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.

gompertz function