42 
Fishery Bulletin 116(1) 
Model implementation and posterior analysis 
The model was run in the Centre of Supercomputing of 
Galicia (CESGA) by using the Galician virtual super¬ 
computer (SVG) cluster with an Ivy Bridge microarchi¬ 
tecture (Intel Corp., Santa Clara, CA). The software 
used was R, vers. 3.2.0 (R Core Team, 2015) and JAGS 
4.2.0 (Just Another Gibbs Sampler) (Plummer, 2012a). 
The R package rjags links the R code with the MCMC 
sampler (Plummer, 2012b; see the supplementary text 
(online only) for the JAGS code used in this study). The 
cluster required 311.28 hours for 170,000 iterations of 
two chains running in different cores. The first 10,000 
iterations were for adaptation without thinning, and 
then a thinning of 100 iterations was applied. After 
a burn-in period of 80,000 iterations, potential scale- 
reduction factors (Gelman et al., 2013) and the Heidel¬ 
berg and Welch convergence diagnostic (Heidelberger 
and Welch, 1981, 1983) were calculated for all monitored 
parameters by using the coda package (Plummer et al., 
2006). The potential-scale reduction factor remained be¬ 
low 1.05 for all of them. Heidelberg and Welch tests were 
passed for both chains and all parameters. 
We used the Kullback-Leibler (KL) divergence (Kull- 
back and Leibler, 1951) as a measure of statistical 
distance between the obtained distributions. The KL 
divergence measures the loss of information when the 
predicted distribution is used instead of the real one. It 
is also considered as a basis for valid inference in eco¬ 
logical Bayesian approaches (Burnham and Anderson, 
2001). It was calculated by using the KLdiv function 
from the R package flexmix (Griin and Leisch, 2008). 
In order to compare the estimated time series with 
available data, the expected value of the quadratic loss 
and its square root, i.e., mean squared error and root 
mean-squared error (MSE and RMSE) were calculated. 
Standardization was necessary to compare yearly mod¬ 
el outputs and in situ data for juveniles because they 
express juvenile abundance at different scales. 
Results 
Posterior distributions and the joint posterior of the 
parameters are displayed in Figure 2 and Supplemen¬ 
tary Figure 1 (online only), respectively. The posterior 
for initial population N* is coherent with the initial 
modeled population in Ruiz et al. (2009), where an un¬ 
informative prior was also implemented. The annual 
fishing-induced mortality mode is higher than the an¬ 
nual natural mortality mode; that finding is consistent 
with the restrictions placed on the priors, but natural 
mortality values are very close to the boundary. The 
mode for parameters L„, and growth rate g for the von 
Berttalanfy growth model are close to the middle point 
of the limits established by the priors according to Bel- 
lido et al. (2000). The easterlies parameter (A) is con¬ 
centrated in values between one and two, and the dis¬ 
charge parameter (p) remains in low values with 0.005 
as the mode of the distribution. 
In Figure 3A, a temporary decay in the modeled re¬ 
cruitment time series from 1992 to 1995 is observed to¬ 
gether with a sudden recovery in the next year. Figure 
3B shows a consistent evolution of catches compared 
with the modeled stock population size, a result that 
is coherent with the assumption of catches proportional 
to abundance. This proportionality is derived by con¬ 
sidering a constant mortality in Baranov’s catch equa¬ 
tion which is equivalent to the expected value of c t in 
Equation 9. 
A detailed representation of the juvenile population 
is presented in Figure 4 to assess the performance of 
the higher temporal resolution for the first life stages. 
The dual-time resolution results in a consistent rela¬ 
tion between the model output for the juvenile popu¬ 
lation (R t (3)+R t (4), t= 1,..., T) and in situ data (Drake 
et al., 2007) and offers a smoother version of the for¬ 
mer, although the way these data are recorded is prone 
to sampling noise. There are also some differences be¬ 
tween both time series—the most remarkable in year 
2000, which shows a smaller estimated number of ju¬ 
veniles than the number reported from in situ data, as 
detailed in Table 3. 
The model seems to offer a smoother version of the 
in situ data although the way these data were recorded 
(Drake et al., 2007) is prone to sampling noise. 
An additional contrast between the output and the 
data not used in the model formulation or parametriza- 
tion comes from the length structure of the population. 
The growth module, included through the matrix G, 
allows the simulation of the length frequency in the 
stock for each month (l t ). The KL divergence between 
model outputs and ICES data was smaller than 2 in 
80.63 of the months considered (Fig. 5), showing that 
a von Berttalanfy-based matrix could model the length 
frequency of the population without a significant loss of 
information in many cases. A lower KL divergence was 
obtained in the later years, where the sensitivity to 
initial values was reduced. Main differences are in the 
first two length classes, as can be seen in the compari¬ 
son between modeled and observed length distributions 
for each month in Supplementary Figure 2 (online only). 
The lengths reported by ICES are in the rank of the es¬ 
timated length (between the 5th and 95th percentiles) 
for the first and second length classes in the 62% and 
59% of all the 222 months, respectively, compared with 
63%, 63%, 81%, and 97% for the other length classes. 
Discussion 
Bayesian models have proven their value for assessing 
the underlying stock-recruitment relationship of ex¬ 
ploited species, along with the associated uncertainty 
(Dorn, 2002; Mantyniemi and Romakkaniemi, 2002; 
Michielsens and McAllister, 2004). In many stocks 
and particularly in small pelagic fish, the sensitivity 
of recruitment to the environment impedes the use of 
traditional stock-recruitment relationships (e.g., Rick¬ 
er or Beverton-Holt curves) as the unique driver of 
