A comparison between traditional and measurement-error growth models for weakfish Cynoscion regalis

Data

Age-length data for weakfish Cynoscion regalis were obtained from Wenner & Gregory (2000), with age for the same individual being estimated from sagittal otolith and scale readings. The otolith-scale age comparison database comprised 2,318 weakfish caught intermittently from five states (i.e., New York, Delaware, Maryland, Virginia, and North Carolina) for years 1989, 1992, 1995, and 1996 (Table 1). Individuals were pooled across states and years to fit von Bertalanffy and measurement-error von Bertalanffy growth curves using otolith- or scale-estimated ages. An age-bias plot (Fig. 1A) indicated ageing uncertainty for weakfish, with scale readings tending toward younger age estimates compared to otolith-estimated ages. Also, percent agreement between ageing structures declined with age, suggesting error in the ability of readers to consistently discern age for older fish (i.e., multiplicative ageing uncertainty) (Fig. 1B).

Nonlinear growth models

The von Bertalanffy growth function has a long history in fisheries science and has been used extensively to describe fish growth (i.e., length and weight) as a function of age (Haddon, 2001). Despite criticisms (Roff, 1980), the von Bertalanffy growth curve has been advocated as an appropriate growth model because of its ability to capture observed trends between length and age for a variety of fish species (Chen, Jackson & Harvey, 1992). A recent stock assessment modeled weakfish growth using a von Bertalanffy growth function (Northeast Fisheries Science Center, 2009) that can be written as: (1) $L_i=L_{\infty }\left(1-e^{\left(-k\left({t^{\prime }}_i-t_o\right)\right)}\right)\cdot e^{{\epsilon }_i}$ where L_i is the length-at-age for the ith individual, L_∞ is the asymptotic length, k is the Brody growth coefficient, t_o is the hypothetical length at age 0, and ${t^{\prime }}_i$ is the observed age for the ith individual. Error ε_i is assumed to be independent and normally distributed with mean 0 and variance ${\sigma }_L^2$.

Extending the von Bertalanffy growth model to incorporate measurement error is relatively straight forward, and can be written as: (2) $\begin{array}{l}{L}_i=L_{\infty }\left(1-e^{\left(-k\left(t_i-t_o\right)\right)}\right)\cdot e^{{\epsilon }_i}\\ {t^{\prime }}_i=t_i\cdot e^{{\epsilon }_i}\end{array}$ where t_i is the true age for the ith individual. The logarithm of observed age ${\text{log}}_{\text{e}}\left({t^{\prime }}_i\right)$ is assumed to be independent and normally distributed with mean log_e(t_i) and variance ${\sigma }_A^2$. In order to facilitate the use of a log-normal distribution for observed ages, a small constant (i.e., 10E-05) was added to age-0 individuals during model fitting.

Statistical estimator

A Bayesian estimator was used to construct the joint posterior probability distribution for parameters in the von Bertalanffy and measurement-error von Bertalanffy growth curves. The full conditional distribution for the traditional von Bertalanffy growth model is: (3) $p\left(L_{\infty },k,t_o,{\sigma }_L^2|L\right)\propto {\displaystyle \prod _{i=1}^{n}L}\left(L_i|L_{\infty },k,t_o,{\sigma }_L^2\right)\times \pi \left({\sigma }_L^2\right)\pi \left(L_{\infty }\right)\pi \left(k\right)\pi \left(t_o\right)$

While the full conditional distribution for the measurement-error von Bertalanffy growth model is: (4) $p\left(L_{\infty },k,t_o,{\sigma }_L^2,{\sigma }_A^2,t|L\right)\propto {\displaystyle \prod _{i=1}^{n}L}\left(L_i|L_{\infty },k,t_o,{\sigma }_L^2,t_i\right)L\left({t^{\prime }}_i|t_i,{\sigma }_A^2\right)\times \pi \left({\sigma }_L^2\right)\pi \left({\sigma }_A^2\right)\pi \left(L_{\infty }\right)\pi \left(k\right)\pi \left(t_o\right)\pi \left(t_i\right)$ where p(·) denotes the posterior probability, $\mathcal{L}(\cdot )$ denotes the likelihood function, and $\pi (\cdot )$ denotes the prior distribution.

As shown in Eq. (4), observed lengths (L_i) are conditionally independent of observed ages (${t^{\prime }}_i$), with the majority of information about true age (t_i) coming from the prior. Observed length will also inform the estimation of true age through feedback of the likelihoods on the joint posterior distribution, with length often assumed to be a loose proxy for age (e.g., age-length keys). Essentially, measurement-error growth models work to pull observations closer to the median length-at-age, suggesting the need for an informative prior on true age when multiple age determinations are unavailable. If an informative prior on true age is unjustifiable, then multiple age determinations will be necessary to estimate the ageing-error variance(s) or a reference collection will be required, in which true age for a set of individuals is known, so that validation data can help calibrate the model during estimation.

Prior distributions were constructed around historic estimates of weakfish growth, thereby encompassing biological relevancy (Seagraves, 1981; Shepherd & Grimes, 1983; Hawkins, 1988; Villoso, 1989; Lowerre-Barbieri, Chittenden & Barbieri, 1995). Age validation data were unavailable and consequently the latent variable of true age was assumed to follow a uniform distribution, with the lower and upper bound being defined by ±4 of the observed age for each individual, as the largest difference between otolith- and scale-estimated ages was three years (Fig. 1A). The weakfish age-length dataset was characterized by a lack of older, larger-sized individuals compared to the most recent investigation of age and growth (Lowerre-Barbieri, Chittenden & Barbieri, 1995). Changes in weakfish age- and size-structure are most likely a culmination of several factors, including: residual effects of excessive fishing mortality (Northeast Fisheries Science Center, 2009), gear selectivity, and seasonal variation in spatial distribution as a result of differential migration by size (Lowerre-Barbieri, Chittenden & Barbieri, 1995). In order to avoid inflated estimates of asymptotic length (L_∞) and consequent underestimation of the Brody growth coefficient (k), posterior values for L_∞ were bounded by the specified prior distribution. A summary of prior distributions can be found in Table 2.

All models were run with three Markov chains for 100,000 simulations per chain using the software packages WinBUGS version 1.4.3 and R version 2.13.1. Convergence of the Markov chains to the stationary distribution was determined by monitoring trace plots and computing Gelman and Rubin diagnostics. The first 50,000 iterations from each chain were discarded to allow for adequate burn-in and a thinning interval of five was used to reduce autocorrelation among iterative samples and improve computational efficiency. A total number of 30,000 iterations were used to summarize the posterior distribution for each model.

Model selection criteria

Growth is a vital component in discerning the population dynamics of fishes and modeling age-length relationships requires the ability to effectively compare and discriminate among alternative hypotheses that represent biological realism. In this study, model checking and discrimination were conducted using posterior predictive p-values and deviance information criterion (DIC), respectively. While DIC has the potential to identify correct model structure for catch-at-age analyses (Wilberg & Bence, 2008), its ability to select preferred models in an EIV context is less clear (Spiegelhalter et al., 2002; Celeux et al., 2006). To circumvent this issue, posterior predictive model checks and model discrimination statistics were used in an effort to corroborate anecdotal beliefs regarding the applicability of measurement-error models during nonlinear growth curve analyses.

Simulation study

To explore the advantages of accounting for ageing error during nonlinear growth curve analyses, we ran simulations that emulated weakfish growth fitting both VBGF and MEVB models to simulated datasets using 1 or 2 age reads. In order to avoid modeling artifacts, we subsetted the original data to only include records collected in 1995 by the state of Delaware, when size truncation was less pronounced. The MEVB model was then fitted to the subsetted dataset using otolith-estimated ages, with the posterior averages of the estimated parameters (i.e., t, L_∞, k, t₀, σ_A, and σ_L) then serving as the true, known values for the simulation. A workflow for the simulation is presented in Fig. 2. Accuracy of the two models was then investigated using the mean absolute error (MAE) and root mean square error (RMSE) statistics, (8) $\text{MAE}=\frac{1}{n}{\displaystyle \sum }_i\left|\widehat{{\theta }_i}-\theta \right|$ (9) $RMSE=\sqrt{\frac{1}{n}{\displaystyle \sum }_i{\left(\widehat{{\theta }_i}-\theta \right)}^2}$ where n is the number of simulations (n = 500), $\widehat{{\theta }_i}$ is the estimated parameter at the ith simulation, and θ is the true, known value used to generate the simulated datasets.

Model discrimination

According to the DIC statistic, traditional von Bertalanffy growth curves outperformed measurement-error growth models for both otolith- and scale-estimated ages (Table 3). Alternatively, posterior predictive p-values for measurement-error growth curves were substantially closer to 0.50 (Table 3; Fig. 3), suggesting improved adequacy of EIV models to reflect observed trends in the age-length relationship for weakfish. All growth curves considered in this study, however, had posterior predictive p-values < 0.50, possibly suggesting underparameterization in the ability of formulated models to partition the overall variance to its respective sources (i.e., variability in age or length). It is also apparent that inclusion of measurement error resulted in correlation between the observed and predicted lengths (Fig. 3), as estimation of true age is being informed, in part, by observed length. Nonetheless, predictive approaches to model comparison may be beneficial for EIV regression, as the utility of information-theoretic-based methods for measurement-error model selection are still circumstantial (Jiao, Reid & Smith, 2009).

von Bertalanffy growth curve parameters

Growth models considering ageing error resulted in higher posterior mean values for L_∞ and t_o (Table 2; Fig. 4), while producing lower posterior mean values for k and σ_L (Table 2; Fig. 4). The posterior mean value for σ_A was higher for scale-estimated ages, although there was substantial overlap between marginal posterior distributions (Fig. 4). Truncation of the joint posterior distribution for L_∞ and k was expected, as specified priors were used to constrain posterior draws to biologically reasonable values. The age-length data for weakfish fail to accurately capture the asymptotic length, leading to unrealistic estimates that are based on extrapolation of the age-length trend (Knight, 1968). As a consequence, measurement-error models demonstrated growth patterns where weakfish grew to reach larger sizes but at slower rates, with traditional von Bertalanffy growth curves overestimating median length-at-age for the observed age range (Fig. 5).

In addition, 95% prediction intervals were wider for traditional von Bertalanffy growth models, compared to their measurement-error analogs (Fig. 5). This is not surprising given that measurement-error models produced lower estimates for variability in predicted lengths (i.e., σ_L), as some of the total variance gets partitioned out for ageing error. Generally, measurement-error growth models produced slightly narrower credible intervals for parameters of the von Bertalanffy growth function, with less difference between posterior mean values for biologically relevant parameter estimates using the different ageing structures (Table 2; Fig. 4).

Simulations

The MAE and RMSE statistics were similar, with only the RMSE being reported for brevity (see Table 4). Results from the simulation study found that the estimated posterior means of the parameters from the MEVB growth models were closer to the “true” values, on average, when ageing error was present (Table 4). This suggests DIC may be inadequate at recommending growth models within EIV contexts, and researchers should pursue other model discrimination statistics when trying to confirm the utility of incorporating ageing uncertainty into growth curve analyses. The simulations also confirmed the intuition gained by using posterior predictive p-values, in that measurement-error models are more adequate at describing growth patterns for weakfish by lessening the degrading effects of ageing error. It is also apparent that multiple age reads can improve estimation of the ageing uncertainty, with simulations also showing that estimated variability in predicted lengths was more biased when ageing error was not considered (Table 4).