Lopez Quintero et al.: Bayesian analysis of the von Bertalanffy growth function 
15 
r 
v + 1 
2 
/ 
r(v/2)(7rv)2 
v+l 
V 
(8) 
z e M is the symmetric Student-^ density with v de- 
grees of freedom, and T{z;v) represents the respec- 
tive cumulative distribution function. In other words, 
we assumed that the original response yi followed a 
log-skew-^ distribution (Marchenko and Genton, 2010), 
which is denoted by 
~LSr(/ii-LL;,af,A,v),f = l,..., n. (9) 
We assumed that v > 1 and considered the first mo- 
ment of the skew-^ distribution (Branco and Dey, 2001); 
therefore, the extra parameter had tobe chosen as 
= 
/7 n(v-i)/2] Agj 
Vtt r(v/2) ^ 1+^2 ’ 
( 10 ) 
so that the transformed errors e- have a zero mean. 
This condition ensured that E{y[) = L[ and allowed us 
to identify the constant in the additive version of the 
regression model. 
Heteroscedasticity was introduced by means of the 
dispersion parameters erf and modeled by using a non- 
negative function m(p; Xj) depending on age x; and a 
heteroscedastic parameter p e M as af = (j^m(p; Xj), 
where > 0. When p = 0, homoscedasticity is recov- 
ered as (jf = cj^m(0; X;) = cr^. In our study, we considered 
2 specific functions for modeling heteroscedasticity: the 
exponential function m(p; X;) = and the power func- 
tion m(p;Xi) = xf^. In both functions, if m(0;xi) = 1, it 
corresponds to the homoscedastic case. 
Asymmetry and heavy tails produced by extreme 
values of length-at-age data were controlled by the pa- 
rameters of shape (A) and degrees of freedom (v). 
Extension to a Bayesian framework 
We advanced a Bayesian analysis for the log-skew- 
t VBGF described in the previous section. Therefore, 
we first noted from the independence assumption and 
Equation 7 that the likelihood function of the unknown 
parameter vector 9 = {(3 ^ is 
fiy' I x,«) = !• + 1), (11) 
where P = {L^,K,Tq)^ are VBGF parameters, y' = {y'i, 
1 x = (xi,...,x„)^, and Zj is as it was previously 
defined. To complete our Bayesian model specification, 
we needed to elicit a prior distribution for the unknown 
parameter vector, say niO). Therefore, the Bayesian in- 
ference on 0 (or function of 9) was based on the posteri- 
or distribution jd^d |x, y') oc/(y'|x, 9)7d9). This posterior 
distribution does not have a closed form (Cancho et al., 
2011), but an estimation could still be calculated by 
using a Markov chain Monte Carlo (MCMC) algorithm 
(Chib and Greenberg, 1995; Cowles and Carlin, 1996; 
Robert and Casella, 2004). 
Given the available methods, we chose to implement 
a hand-tailored component-wise Metropolis-Hasting 
transition scheme; in other words, we selected an ap- 
proach in which 6 is divided into individual pieces that 
are easily updated sequentially with a random walk 
algorithm. Our selection was based on simplicity and 
the need to control all steps in the sampling. Other 
options included the use of variants of BUGS language 
(Lunn et al., 2012), the AD ModelBuilder (Fournier et 
al., 2012), or Stan software (Gelman et al., 2015). Note 
that our selected approach is different from that fol- 
lowed by Siegfried and Sanso (2006), who employed an 
algorithm that included both Gibbs, as well as Metrop- 
olis-Hasting steps. An advantage of using this MCMC 
scheme is that the procedure is subdivided into several 
univariate steps. All proposal distributions were tuned 
to achieve acceptance rates of 25-45% (Robert and Ca- 
sella, 2004). Specifically, we considered the different 
components of as independent (Siegfried and Sanso, 
2006); in other words, idB) becomes the product of the 
marginal prior distributions of ,G^,p,\,v). It follows 
that only these marginal prior distributions must be 
elicited to complete our Bayesian model. 
Prior distributions for the VBGF parameters P were 
chosen as follows. Given that is strictly positive, we 
assumed a left truncated normal distribution with large 
variance (e.g., 100) as the prior distribution for this pa- 
rameter. Other possible and natural choices were log- 
normal, gamma, or even distributions with support in 
a reasonable and restricted interval. Because param- 
eters K and -to are both positive, we used gamma as 
prior distributions for them (Xiao, 1994; Siegfried and 
Sanso, 2006). Contreras-Reyes et al. (2014) reported 
estimates of around 0.16 for southern blue whiting, a 
value that incidentally conforms to the value obtained 
by Siegfried and Sanso (2006) for blue shark {Prionace 
glauca). We used this information to specify that the 
mean of the gamma prior distribution of was around 
0. 15. For -^ 0 , Contreras-Reyes et al. (2014) obtained 
-^0 = 2.5,, which indicates a prior distribution for this 
parameter. 
For the scale parameter cP, we considered the clas- 
sical inverse gamma prior distribution suitable for 
this type of parameter (Zhang et al., 2009). The het- 
eroscedastic parameter p usually takes positive or 
negative values. To give the full power of estimation 
to the data, we chose a noninformative prior, nip) 
1. For the shape parameter A, we set a normal prior 
distribution with a zero mean and large variance. The 
parameter defining the degrees of freedom, v, should 
be strictly larger than 2 to ensure the existence of vari- 
ance in the log-skew-^ model; therefore, we considered 
an exponential prior distribution with mean equaling 
2 (Geweke, 1993; Cancho et al., 2011) and truncated at 
the interval (2, 00 ). These prior specifications are sum- 
marized in Table 1. 
Comparisons and selection of models 
For sake of comparison, we considered 2 additional 
models with constant variance function derived from 
the log-normal distribution. The first one (hereafter, re- 
