8 Growth
Collins 1889.
Growth is a basic aspect of a fish’s life history, but it is also an important process for managing fish populations. For example, a minimum size for harvest might be used to increase the proportion of fish surviving to the size (and age) of sexual maturity. Another consideration in developing a harvest strategy is finding a balance between the growth rate (increasing size of each individual in a cohort) and the loss rate due to natural mortality (reducing the number of individuals in a cohort).
There are typically two aspects to growth modeling: age versus length and length versus weight. With those two equations, it is possible to predict the length or weight at any age, and to examine harvest strategies using minimum size or slot regulations. The focus in commercial marine fisheries is typically on weight because income from commercial fishing is a function of weight of the catch. Length gets more emphasis in freshwater fisheries management because anglers are not selling their catches (and trophy designations are generally length-based).
8.1 Growth in length
The relationship between age and length for adult fish usually is a curve of decreasing slope; that is, growth gradually slows with age as individuals approach some average maximum size. The equation most often used for age and length is the von Bertalanffy growth function: \(L_t=L_\infty*(1-e^{-k*(t-t_0)})\). This equation for size at time t has three parameters: the asymptotic size or average maximum length (\(L_\infty\)), growth rate (\(k\)), and time (age) at length of 0 (\(t_0\)). The last parameter is a theoretical one rather than something that is observed. It can take on positive or negative values, and can be thought of as an extrapolation from sizes at older ages.
There are two usual sources of data for fitting a von Bertalanffy growth curve: paired age-length observations from a laboratory examination of otoliths or other hard parts that form annual marks (Campana 2001); or a tagging study. We consider both approaches in the following sections.
8.1.1 Age-length
Age:length observations may come from a survey or from sampling a fishery. Some mis-ageing invariably occurs when interpreting otoliths or other hard parts (Campana 2001), but for simplicity we assume that the ages are determined without error. We further assume that fish are not sampled randomly but are instead a systematic sample by age. This ensures that the data set includes less common ages (sizes), particularly older fish.
The first step in fitting the model is to generate a simulated age:length data set.
# Fitting von Bertalanffy curve to simulated data
rm(list=ls()) # Clear Environment
# Choose von Bertalanffy parameter values and other simulation settings
L_inf <- 80 # e.g., 80-cm fish
k <- 0.3 # Slope
t_0 <- -0.5 # Age at length 0
MaxAge <- 12
SystSample <- 5 # Number per age group
N.AgeLen <- SystSample * MaxAge
Age <- rep(1:MaxAge, each=SystSample)
Var_L <- 10 # Variance about age-length relationship
Len <- rnorm(n=N.AgeLen,L_inf*(1-exp(-k*(Age-t_0))), sd=sqrt(Var_L))
plot(Age, Len, xlab="Age", ylab="Length (cm)")The parameter values and other simulation settings provide a reasonable sample size and well-defined growth pattern and asymptote. The latter characteristic is critical in obtaining reliable results. For real data sets, it is helpful to look at the shape of the scatter plot of length versus age. Data sets with few young fish or lacking a well-defined asymptote may result in poor estimates. Note also that we are using integer ages. If there was good information about the timing of spawning, ages could be calculated as the integer age plus the fraction of the year since the spawning season. Using fractional ages would be more important for short-lived species. The model fitting process is the same whether using integer or fractional (real valued) ages.
We assume an additive error term; i.e., variability in length at age is constant. Whether an additive or multiplicative error term is appropriate can often be seen by examining a scatter plot of length versus age. A multiplicative error term would be called for when points have greater spread for older fish.
The JAGS code for fitting the model is generally similar to the simulation, except that the dnorm() function for the likelihood requires precision, defined as the inverse of the variance. Vague uniform priors are used for all four parameters. The prior distribution ranges and initial values may need to be adjusted if the range of simulated (or observed) lengths is varied. Note that max(Len) is used as a convenient upper bound for initial values of \(L_\infty\), but a fixed value greater than observed lengths could also be used.
# Load necessary library
library(rjags)
library(R2jags)
# JAGS code
sink("GrowthCurve.txt")
cat("
model{
# Priors
L_inf.est ~ dunif(0, 200)
k.est ~ dunif(0, 2)
t0.est ~ dunif(-5, 5)
Var_L.est ~ dunif(0,100) # Variability in length at age
# Calculated value
tau.est <- 1/Var_L.est
# Likelihood
for (i in 1:N.AgeLen) {
Len_hat[i] <- L_inf.est*(1-exp(-k.est*(Age[i]-t0.est)))
Len[i] ~ dnorm(Len_hat[i], tau.est)
} #i
}
",fill=TRUE)
sink()
# Bundle data
jags.data <- list("N.AgeLen", "Age", "Len")
# Initial values
jags.inits <- function(){ list(L_inf.est=runif(n=1, min=0, max=max(Len)),
k.est=runif(n=1, min=0, max=2),
t0.est=runif(n=1, min=-5, max=5),
Var_L.est=runif(n=1, min=0, max=100))}
model.file <- 'GrowthCurve.txt'
# Parameters monitored
jags.params <- c("L_inf.est", "k.est", "t0.est", "Var_L.est")
# Call JAGS from R
jagsfit <- jags(data=jags.data, jags.params, inits=jags.inits,
n.chains = 3, n.thin=1, n.iter = 10000,
model.file)
print(jagsfit)
plot(jagsfit)Convergence is generally rapid and the estimates from simulated data are reliable. Try varying some of the simulation settings and parameter values to see how the model performs under different conditions. For example, note how the growth curve changes as you try different values for the von Bertalanffy parameters and the variance about the curve. (One caution is to ensure that the assumed additive error term does not produce negative length values. If that occurs, one solution would be to switch to a multiplicative error term, using a lognormal rather than a normal distribution (see Section 8.2.) Another setting to vary is the maximum age. A low value, as might occur if the population was heavily exploited, can make it difficult to estimate \(L_\infty\). One Bayesian approach for dealing with limited data is the use of informative prior distributions (Doll and Jacquemin 2018).
8.1.2 Using tagging data
Growth information can also be obtained from a tagging study (Wang et al. 1995; Zhang et al. 2009; Scherrer et al. 2021). Fish seen on two or more occasions provide information on growth as a function of size at initial capture, time-at-large (interval between tagging and a subsequent recapture) and the growth increment. This approach differs from using age and length data because age at tagging is not known. However, it can be made compatible with a traditional age:length analysis by estimating relative age (true age - \(t_0\)) of all individuals as additional parameters (Wang et al. 1995; Zhang et al. 2009; Scherrer et al. 2021). This allows growth to be modeled as a function of age rather than size. Flexible versions of these models can allow for separate growth parameters (\(L_\infty\), k) for each individual. Here we fit a simpler model (Model 4 of Scherrer et al. (2021)) with all individuals sharing \(L_\infty\) and k. It is straightforward to modify the code to allow for either or both parameters to vary among individuals (Scherrer et al. 2021).
One important practical question is how to obtain length data when fish are encountered after tagging. Using lengths from fishers provides the greatest sample size but runs the risk of substantial measurement error (or outright bias if an angler “estimates” the length). Length data obtained by biologists limits the sample size but ensures data quality, and could be obtained through field survey sampling or fishery monitoring by creel clerks or at-sea observers. In order for our results to apply to the whole population, we also assume that growth is the same for tagged and untagged fish. The model allows for measurement error and other sources of variability (e.g., environmental, genetic) that cause individual fish lengths to vary about the fitted curve. A constant additive error is assumed here but the model could be modified to allow for a multiplicative error term. We begin with simulation code for generating the observations:
rm(list=ls()) # Clear Environment
L_inf <- 80 # e.g., 80-cm fish
k <- 0.3 # Slope
Var_L <- 10 # Variance about age-length relationship
N.Tag <- 120 # Number of fish recaptured at least once
n <- rpois(N.Tag, lambda=0.1)+2 # Number of encounters for each recaptured individual, >2 uncommon
# Generate vector for relative age at initial capture (true age + t0)
Shape <- 5 # Parameters for gamma distribution generating relative ages
Rate <- 2
A <- rgamma(n=N.Tag, shape=Shape, rate=Rate) # Relative age vector
hist(A, xlab="Relative age", main="")
L <- array(data=NA, dim=c(N.Tag, max(n)))
dt <- array(data=NA, dim=c(N.Tag, max(n)-1))
for (i in 1:N.Tag){
for (j in 1:(n[i]-1)){
dt[i,j] <- rgamma(n=1, shape=6, rate=2) # Random dist of times at large
}
}
hist(dt, xlab="Time at large", main="")
for (i in 1:N.Tag){
L[i,1] <- rnorm(n=1, mean=L_inf *(1.0 - exp(-k*A[i])), sqrt(Var_L))
for (j in 2:n[i]){
L[i, j] <- rnorm(n=1, L_inf*(1.0 - exp(-k*(A[i]+dt[i, j-1]))), sqrt(Var_L))
} #j
} #i
plot(A, L[,1], xlab="Relative age", ylab="Length at first capture")
plot(dt[,1], (L[,2]-L[,1]), xlab="Time at large", ylab="Growth increment")
plot(L[,1], (L[,2]-L[,1]), xlab="Length at first capture", ylab="Growth increment")Choices for the growth parameters are arbitrary but should produce a well-defined asymptote when combined with the other settings (time-at-large, relative age). The sample size of fish recaptured at least once is relatively robust but can of course be varied to judge model performance given fewer or more individuals. To allow for random variation in the tagging data set, a Poisson distribution is used to generate the number of encounters for each individual. Each individual is seen at least twice but the lambda parameter can be adjusted to vary the typical number of encounters (use hist(n) to see distribution). As in previous studies (Wang et al. 1995; Zhang et al. 2009; Scherrer et al. 2021), we use a gamma distribution for relative age. The two gamma parameters (shape and rate) can be varied to obtain relative age distributions with desired characteristics (as observed on the hist() plot). Here they are chosen to provide a data set of mostly younger fish. Times-at-large are also generated using a gamma distribution, with shape and rate settings chosen to obtain a realistic pattern (mostly shorter times). Three plots are used to examine the tagging data set. The first (relative age versus length) may not have a well-defined asymptote if mostly younger fish are tagged. The plot can be examined to ensure that the chosen variance in length at age is reasonable (and not generating negative lengths). The second and third plots are noisy and difficult to interpret because growth increment is a function of both time-at-large and size at tagging. They are nevertheless useful for ensuring that the data set encompasses the early faster-growing period and slowed growth near the asymptote.
The JAGS code for fitting the model uses uninformative uniform prior distributions for the two growth parameters, variance in length-at-age, and the two gamma hyperparameters that describe the distribution of relative ages.
# Load necessary library
library(rjags)
library(R2jags)
# JAGS code
# Code modified from Scherrer et al. 2021: model 4
sink("vonBert_GrowthInc.txt")
cat("
model{
# Priors
Var_L.est ~ dunif(0, 100) # Variance about expected length
tau.est <- 1/Var_L.est # Precision
k.est ~ dunif(0, 2)
L_inf.est ~ dunif(0, 200)
Shape.est ~ dunif(0, 100) # Gamma hyperparameters for relative age
Rate.est ~ dunif(0, 100)
# Likelihood
for (i in 1:N.Tag) {
A.est[i] ~ dgamma(Shape.est, Rate.est) # Relative age (true age + t0)
L_Exp[i, 1] <- L_inf.est *(1.0 - exp(-k.est*A.est[i])) # Expected length at capture
L[i, 1] ~ dnorm(L_Exp[i, 1], tau.est) # Length at initial capture
for (j in 2:n[i]) { # Recapture lengths
L_Exp[i, j] <- L_inf.est*(1.0 - exp(-k.est*(A.est[i]+dt[i, j-1])))
L[i, j] ~ dnorm(L_Exp[i, j], tau.est)
} #j
} #i
}
",fill=TRUE)
sink()
# Bundle data
jags.data <- list("N.Tag", "n", "L", "dt")
# Initial values
jags.inits <- function(){ list(L_inf.est=runif(n=1, min=0, max=200),
k.est=runif(n=1, min=0, max=2),
Var_L.est=runif(n=1, min=0, max=100),
Shape.est=runif(n=1, min=0, max=100),
Rate.est=runif(n=1, min=0, max=100)
)}
model.file <- 'vonBert_GrowthInc.txt'
# Parameters monitored
jags.params <- c("L_inf.est", "k.est", "Var_L.est" ,"Shape.est", "Rate.est"
)
# Call JAGS from R
jagsfit <- jags(data=jags.data, jags.params, inits=jags.inits,
n.chains = 3, n.thin=1, n.iter = 20000,
model.file)
print(jagsfit)
plot(jagsfit)The model requires a considerable number of updates to achieve convergence, which is understandable given the large number of parameters (due to estimating relative age for each individual). The estimates tend to be reliable, which can be attributed to the large sample size and good ranges for length-at-tagging and time-at-large. Varying the simulation settings for time-at-large and relative age are helpful for understanding this approach to growth analysis. It should not matter for this model whether observations are from more individuals or more captures per individual. Increasing the number of observations per individual (Poisson lambda) would be relevant if fitting a model with growth parameters that varied among individuals.
For a given sample size, precision is lower for tagging compared to age:length data. This seems reasonable given that relative age must be estimated in a tagging study, whereas absolute ages are (assumed to be) known in an age:length analysis. Nevertheless, there are some situations where tagging data are preferred to age:length, such as for tropical fishes that do not show annual marks on hard parts.
It is straightforward to carry out a combined analysis using age:length and tagging data. The age:length data provide information on all three growth parameters (\(L_\infty\), \(k\), \(t_0\)) as well as the variance term. The tagging data provide information on all parameters except \(t_0\). A combined analysis might be preferable if multiple methods provide better coverage over the entire age/size range. For example, Scherrer et al. (2021) demonstrated the benefits of a combined analysis using age:length, tagging, and length frequency data.
8.2 Growth in weight
The length-weight relationship for fish is typically nonlinear: \(W=a*L^b\). The ‘a’ parameter is a scaling coefficient and varies depending on the units chosen for measuring length and weight (e.g., cm vs mm or kg vs g). The ‘b’ parameter is close to 3 for most fishes, because weight is a three-dimensional value, increasing approximately as a cube of length.
The length-weight relationship is useful for several reasons. Length can be measured more quickly and accurately than weight, especially on the water. Fortunately, there tends to be a very strong relationship between length and weight, so it is often advantageous to measure length and predict weight. In addition, the length-weight relationship provides information about fish condition (shape). Fish in poor condition may be overcrowded (e.g., in an pond or aquaculture setting) or lacking in suitable forage.
We begin with simulation code, using an arbitrary range of 20 to 80 for length (e.g., measured in cm). For our example, the scaling parameter ‘a’ will be around 1E-6, so we instead work with the ln-scale value (-12). The model can be fitted using either ‘a’ or the ln-scale transformed parameter. Any assumed value around 3 is fine for the ‘b’ parameter. A multiplicative error term is assumed, as variability in weight typically increases as a function of length (Hayes and Brodziack (1997)).
# Fitting length:weight relationship to simulated data
rm(list=ls()) # Clear Environment
# Choose parameter values and other simulation settings
ln_L_W_a <- -12 # e.g., L in cm, weight in kg
L_W_b <- 3.2
MinL <- 20
MaxL <- 80
SystSample <- 3 # Number per length group
Len <- rep(seq(from=MinL, to=MaxL, by=2), each=SystSample)
SampSize <- length(Len)
W_Var <- 0.1 # Ln-scale variance in weight at length
Wgt <- rlnorm(n=SampSize, mean=(ln_L_W_a+log(Len^(L_W_b))),
sd=sqrt(W_Var))
par(mfcol=c(1,1)) # Ensure default of 1 panel (changes in trace plots)
par(mar=c(5.1, 4.1, 4.1, 2.1)) # Ensure default margins
plot(Len, Wgt, xlab="Length (cm)", ylab="Weight (kg)")The simulation uses a systematic sample for length, using the seq() function nested within rep(). Use Help for the rep() function and try each function separately in the Console to understand how the two functions work together to generate the length vector. The plot par() function is used to undo formatting changes in the trace plots generated below. As the plot illustrates, variability around the underlying length-weight relationship increases as a function of length.
The JAGS code is relatively short and fitting appears straightforward, although appearances can be deceiving! This simple model is probably the most difficult one to fit in this book, because the two parameters are highly correlated. The high correlation means that offsetting changes in ‘a’ and ‘b’ can provide essentially the same fit, so the updating process is very slow. This model provides a good opportunity to learn about the updating process and using Rhat and trace plots to judge convergence.
The JAGS code for fitting the model begins with broad uniform ranges for the prior distributions. Note that we use same lognormal distribution for the likelihood as in the simulation code (dlnorm() in JAGS). The lognormal distribution assumes a constant, additive error term which is a multiplicative error in the back-transformed arithmetic scale. It is convenient to work with the ln-scale transformed value (\(log(a)+log(L^b)\)) but the model can also be fitted working directly with ‘a’ parameter (\(log(aL^b)\)).
# Load necessary library
library(rjags)
library(R2jags)
# JAGS code
sink("LenWgt.txt")
cat("
model{
# Priors
ln_L_W_a.est ~ dunif(-20, 20)
L_W_b.est ~ dunif(1, 4)
W_Var.est ~ dunif(0,0.5) # Multiplicative error. Variability around ln-scale line
# Calculated value
precision <- 1/W_Var.est
# Likelihood
for (i in 1:SampSize) {
lnWgt_hat[i] <- ln_L_W_a.est+log(Len[i]^L_W_b.est)
Wgt[i] ~ dlnorm(lnWgt_hat[i], precision)
} #i
}
",fill=TRUE)
sink()
# Bundle data
jags.data <- list("SampSize", "Wgt", "Len")
# Initial values
jags.inits <- function(){ list(ln_L_W_a.est=runif(n=1, min=-20, 20),
L_W_b.est=runif(n=1, min=1, max=4),
W_Var.est=runif(n=1, min=0, max=0.5))}
model.file <- 'LenWgt.txt'
# Parameters monitored
jags.params <- c("ln_L_W_a.est", "L_W_b.est", "W_Var.est")
# Call JAGS from R
jagsfit <- jags(data=jags.data, jags.params, inits=jags.inits,
n.chains = 3, n.thin=1, n.iter = 5000,
model.file)
print(jagsfit)
plot(jagsfit)
par(mar=c(5.1, 4.1, 4.1, 2.1)) # Set default margins
plot(jagsfit$BUGSoutput$sims.list$ln_L_W_a.est,
jagsfit$BUGSoutput$sims.list$L_W_b.est)
cor(jagsfit$BUGSoutput$sims.list$ln_L_W_a.est,
jagsfit$BUGSoutput$sims.list$L_W_b.est)
jagsfit.mcmc <- as.mcmc(jagsfit) # Creates an MCMC object for plotting
par(mar=c(1,1,1,1)) # Avoid error for "figure margins too large"
plot(jagsfit.mcmc) # Trace and density plots
autocorr.diag(jagsfit.mcmc)With 5,000 MCMC iterations, the estimates are generally close to the true values but Rhat is inflated. Plotting the sequence of MCMC estimates for ln-scale ‘a’ versus ‘b’ makes clear the extremely high correlation between the two parameters. A numerical estimate of the correlation (~0.99) is obtained in the Console window by using the cor() function.
The trace plots show clearly that many additional updates are needed for the length-weight parameters. The high autocorrelations for ‘a’ and ‘b’ cause the trace plots to look more like sea serpents than grassy lawns. Compare this result to a longer run using 50,000 updates. The trace plots for ‘a’ and ‘b’ now have a much better appearance and we have good low estimates for Rhat. The autocorrelations can be reduced by thinning (e.g., try n.thin=10). This improves the appearance of trace plot, but the reduced sample size slightly reduces precision and is not generally recommended (Link and Eaton 2012).
It can be helpful to add code for displaying the length-weight observations and the fitted curve. The as.numeric() function converts each array (of length 1) to a scalar for use in generating the fitted curve, which is added to the plot using the points function.
par(mar=c(5.1, 4.1, 4.1, 2.1)) # Restore default margins
plot(Len, Wgt, xlab="Length (cm)", ylab="Weight (kg)")
par_a <- exp(as.numeric(jagsfit$BUGSoutput$mean$ln_L_W_a.est))
par_b <- as.numeric(jagsfit$BUGSoutput$mean$L_W_b.est)
Wgt_hat <- par_a*(Len^par_b)
points(Len, Wgt_hat, type="l", col="red")8.3 Exercises
Before carrying out simulation runs, describe the expected outcome in fitting the von Bertalanffy growth curve to age:length data if you changed the growth rate (k) to 0.05. Report the results for asymptotic size estimates from several runs. What other simulation setting would need to change in order to accommodate such a slow growth rate?
Modify the code from Section 8.1.1 to add the fitted curve to the age:length scatter plot. As a hint, you can either redo the scatter plot after the
plot(jagsfit)line or add a # to turn off the jagsfit plot.Before carrying out simulation runs, describe the expected outcome in fitting the von Bertalanffy growth curve to age:length data if you changed the maximum age to 4 (e.g., if the stock was heavily exploited). Report the results for asymptotic size estimates from several runs.
Modify the code for the von Bertalanffy age:length analysis to obtain the Age vector using the rgamma function that provided relative ages for the tagging analysis (shape=5, rate=2). How does that change affect estimates of asymptotic size? Find a pair of values for shape and rate that result in reliable estimates of the growth parameters.