Connecting single-stock assessment models through correlated survival

Connecting single-stock assessment models through correlated survival Abstract Fisheries management is mainly conducted via single-stock assessment models assuming that fish stocks do not interact, except through assumed natural mortalities. Currently, the main alternative is complex ecosystem models which require extensive data, are difficult to calibrate, and have long run times. We propose a simple alternative. In three case studies each with two stocks, we improve the single-stock models, as measured by Akaike information criterion, by adding correlation in the cohort survival. To limit the number of parameters, the correlations are parameterized through the corresponding partial correlations. We consider six models where the partial correlation matrix between stocks follows a band structure ranging from independent assessments to complex correlation structures. Further, a simulation study illustrates the importance of handling correlated data sufficiently by investigating the coverage of confidence intervals for estimated fishing mortality. The results presented will allow managers to evaluate stock statuses based on a more accurate evaluation of model output uncertainty. The methods are directly implementable for stocks with an analytical assessment and do not require any new data sources. Introduction Fisheries management is often based on single-species models for single management areas where the stocks are assumed to have no interaction with surrounding stocks or other species, except via the predetermined natural mortality. Such a simplification is often unrealistic (Eero et al., 2014), and ignoring stock interactions can affect model output (Holsman et al., 2016). These assumptions in single-stock models are, however, convenient as they reduce the data needed on interactions such as migration and predator-prey relations. Collecting fisheries data are costly and time-consuming: The current surveys require expensive scientific cruises and lab processing; tagging studies require catching and manually tagging many individuals with either low probability of recapture, or expensive electronic tags for gathering telemetry data; food-web data such as stomach samples require manual sampling of stomach content and possibly subsequent DNA analysis; spatially explicit catch data require extensive monitoring systems. Further, incorporating such new data sources in stock assessment models requires historic data that covers a long period for reliable statistical results. Several methods ranging in complexity have been developed to improve single-stock assessments by accounting for interactions between stocks (e.g. Gislason and Helgason, 1985; Christensen and Pauly, 1992; Begley and Howel, 2004; Fulton et al., 2004; Lewy and Vinther, 2004; Senina et al., 2008; Methot and Wetzel, 2013). Migration between areas can be explicitly modelled by a transportation matrix that moves, e.g. a proportion of the population to another area (Methot and Wetzel, 2013). These proportions can be informed by tagging studies such as mark-recapture (Quinn et al., 1990; Hannesson et al., 2008) or, for simple models, be estimated from catch and survey data; however, when estimated from catch and survey data, the results may be influenced by other factors such as inaccuracies in the assumed natural mortality or model misspecification. Likewise, predation can be modelled by a rate-of-consumption matrix (Hollowed, 2000b) consisting of factors which describe the predator induced mortality. Estimating predator-induced mortality requires knowledge of complex food webs (Yodzis, 1998; Tarnecki et al., 2016) along with consumption data from stomachs (e.g. Gislason and Helgason, 1985; Daan, 1987; Hollowed, 2000b; Frøysa et al., 2002; Lewy and Vinther, 2004; Richards and Jacobson, 2016). Complex ecosystem or multi-species models often depend on extensive data from different sources such as tagging studies, DNA analysis, environment data or stomach content data (e.g. Christensen and Walters, 2004). Such data sources are not readily available from single-species assessments, are expensive to produce, and may take several years to collect. Data requirements can also include human activity, biochemistry, physiology, geology and oceanography (Link et al., 2010) raising the complexity and data collection cost. Further, ecosystem models often have a vast amount of settings and can be difficult to calibrate and run (Hollowed, 2000a; Chen et al., 2009; Ainsworth and Walters, 2015). Currently, there is a gap between single-stock assessment models and the complex ecosystem models (Collie et al., 2016), with a lack of operational models that combine multiple stocks while only utilizing the data readily available for single-stock assessments, making it directly implementable for fisheries management. The objective of this study is to introduce a method for combining single-stock state-space age-based assessments models to provide a simple alternative to complex ecosystem models. The method is illustrated in three case studies including a total of six European stocks combined two by two. In each case study, we investigate whether combining single-stock models improves the model fit to data compared to individual assessments. We extend the single-stock age-based state-space stock assessment model to allow simultaneous estimation for several stocks through correlated survival processes. Correlations are included in the model through the corresponding partial correlations. Assuming Gaussianity, two random variables from a collection of variables are partially uncorrelated correcting for the remaining variables in the collection if and only if they are conditionally independent given the remaining variables in the collection. Hence, modelling partial correlations allows a clear interpretation of reducing model parameters (i.e. by fixing partial correlations at zero), while still allowing a flexible correlation structure. The multi-stock model introduced does not require further data than the single-stock model. Methods Stock assessment model We implemented age-based state-space stock assessment models for six European stocks (Table 1). The models were implemented with the R-package Template Model Builder (Kristensen et al., 2016). Table 1 Overview of the data sources used in the case studies. Fleet  First year  Last year  First age  Last age  Years with missing  Baltic Cod Case Study            Eastern Baltic Cod (Sub divisions: 25–32)            Commercial  1966  2013  2  8    Survey Q1  1991  2000  2  6    Survey Q1  2001  2013  2  6    Survey Q4  2002  2013  2  5    CPUE tuning index  1997  2013  3  7    Western Baltic Cod (Sub divisions: 23–24)            Commercial  1970  2013  1  7    CPUE tuning index  1997  2013  3  6    Survey Q1  1992  2013  1  6    Survey Q4  1992  2013  0  5    Baltic Herring Case Study            Bothnian Sea Herring (Sub division: 30)            Commercial  1973  2014  1  9    Survey Q4  2007  2014  2  8    Survey Q2  1990  2014  3  9    Central Baltic Herring (Sub divisions: 25–27, 28.2, 29, 32)            Commercial  1974  2014  1  8    Survey Q4  1991  2014  1  8  1993, 1995, 1997  Baltic Mix Case Study            Baltic Sprat (Sub divisions: 22–32)            Commercial  1991  2014  1  8    Survey Q4  1991  2014  1  8  1993, 1995, 1997  Survey Q2  2001  2014  1  8    Survey Q1  1992  2014  1  1  1994, 1996, 1998  Western Baltic Herring (Sub divisions: 22–24)            Commercial  1991  2014  0  8    Survey Q3  1991  2014  1  8  1991, 1999  Survey Q4  1991  2014  0  8  1991–1993, 2001  Survey Q2  1992  2014  0  0    Survey Q1  1991  2014  1  4    Survey Q3  1991  2014  1  4    Fleet  First year  Last year  First age  Last age  Years with missing  Baltic Cod Case Study            Eastern Baltic Cod (Sub divisions: 25–32)            Commercial  1966  2013  2  8    Survey Q1  1991  2000  2  6    Survey Q1  2001  2013  2  6    Survey Q4  2002  2013  2  5    CPUE tuning index  1997  2013  3  7    Western Baltic Cod (Sub divisions: 23–24)            Commercial  1970  2013  1  7    CPUE tuning index  1997  2013  3  6    Survey Q1  1992  2013  1  6    Survey Q4  1992  2013  0  5    Baltic Herring Case Study            Bothnian Sea Herring (Sub division: 30)            Commercial  1973  2014  1  9    Survey Q4  2007  2014  2  8    Survey Q2  1990  2014  3  9    Central Baltic Herring (Sub divisions: 25–27, 28.2, 29, 32)            Commercial  1974  2014  1  8    Survey Q4  1991  2014  1  8  1993, 1995, 1997  Baltic Mix Case Study            Baltic Sprat (Sub divisions: 22–32)            Commercial  1991  2014  1  8    Survey Q4  1991  2014  1  8  1993, 1995, 1997  Survey Q2  2001  2014  1  8    Survey Q1  1992  2014  1  1  1994, 1996, 1998  Western Baltic Herring (Sub divisions: 22–24)            Commercial  1991  2014  0  8    Survey Q3  1991  2014  1  8  1991, 1999  Survey Q4  1991  2014  0  8  1991–1993, 2001  Survey Q2  1992  2014  0  0    Survey Q1  1991  2014  1  4    Survey Q3  1991  2014  1  4    Q1–Q4 indicates at which quarter of the year the survey is conducted. Process model The process model for an individual stock, s, was identical to Nielsen and Berg (2014) model D. The population process followed an exponential decay model,   log ⁡Ns,a,y= log ⁡Ns,a−1,y−1−Fs,a−1,y−1−Ms,a−1,y−1+ηs,a,y,a∈{2,3,…,As+−2,As+−1} with special conditions for the last age group As+ (See Nielsen and Berg, 2014; Albertsen et al., 2017, for details). In the model, Fs,a,y is the fishing mortality, Ms,a,y is the natural mortality, and ηs,a,y is a random Gaussian zero mean increment, which allow inaccuracies in the deterministic part of the model. These inaccuracies can be thought of as random deviations in the natural mortality, which is assumed to be known. While Nielsen and Berg (2014) and Albertsen et al. (2017) model ηs,a,y as independent, we will allow correlations across stocks and ages (see “Connecting stocks” Section), while maintaining independence across years. Recruitment, log ⁡Ns,1,y, was modelled by a random walk. Configuration of model parameters with respect to coupling across ages was done in the same way as for the official advice. Observational model Following Albertsen et al. (2017), the observed log-catch conditional on the process model were modelled by a multivariate normal where the covariance Σs(O) had an AR(1) covariance structure, (Σs(O))i,j=ρo,s|i−j|τo,s,iτo,s,j with ρo,s∈(−1,1) and 0<τo,s,i for all indices i. The catch vector, Cs,y=(Cs,1,y,…,Cs,As+,y), for stock s in year y followed the Baranov catch equation   log ⁡Cs,y= log ⁡Fs,y− log ⁡(Fs,y+Ms,y)+ log ⁡Ns,y+ log ⁡(1− exp ⁡(−Fs,y−Ms,y))+ϵs,y(O) were   ϵs,y(O)∼N(0,Σs(O)),   Ns,y=(Ns,1,y,…,Ns,As+,y),  Fs,y=(Fs,1,y,…,Fs,As+,y), Ms,y=(Ms,1,y,…,Ms,As+,y), and all operations are element-wise. Likewise, survey indices were modelled by   log ⁡Is,y(i)= log ⁡qs(i)+ log ⁡Ns,y−(Fs,y+Ms,y)ds(i)/365+ϵs,y(I,i) where ϵs,y(I,i)∼N(0,Σs(I,i)), Is,y(i) is a vector of age observations for the ith survey in stock s year y, q(i)>0, dS(i) is the number of days into the year the survey is conducted, and Σs(I,i) has an AR(1) structure. For simplicity, years with missing data for single ages were treated as if they were missing all ages for that fleet. Alternatively, missing data could be handled by, for instance, imputing random effects. All observational parameters were estimated; however, they were configured with respect to coupling across ages in the same way as for the official advice. Connecting stocks For many stocks, assuming independence between the error terms, ηs,a,y, is unrealistic. Inaccuracies in the assumed natural mortality will be correlated by connections not included in the model such as environmental changes, migration, or predation. All other things held constant, environmental effects will be seen as positive correlation; migration will be seen as negative correlation; while predation will be seen as negative correlation in a predator-prey relation or as positive correlation if two stocks are influenced by the same predator. To connect the S individual stocks, we temporarily extended the stock log-population processes to have the same number of ages, A, for convenience. As an example, consider a case with two stocks (such as the Baltic cod case study presented in “Case studies” Section) where one stock includes ages 2–8, while the other stock includes ages 0–7. To construct the correlation matrix, both log-population processes have to be temporarily extended to include ages 0–8 to combine the stocks. The individual processes were stacked to one process,   Ny=(N1,1,y,…,N1,A,y,…,NS,1,y,…,NS,A,y)T, for which we modelled the error terms, ηy=(η1,1,y,…,η1,A,y,…,ηS,1,y,…,ηS,A,y)T, by a multivariate normal with covariance ΣN. The covariance was parameterized through the Cholesky decomposition of the precision matrix PN. Defining L to be a lower triangular matrix of parameters with diagonal 1, a positive definite precision matrix was obtained by   PN=LLT If the precision matrix was scaled to have diagonal 1, the off-diagonal elements were the partial correlations with opposite sign. With this construction, the diagonal elements of P could not be estimated freely. To obtain the covariance matrix as the inverse precision, we scaled P to have diagonal 1 and added freely estimated variance parameters. Including variance parameters at this stage maintained the interpretation of the estimated parameters from the single-stock model, where the parameter is directly related to the variance of the specific age. The scaling was done by defining a diagonal scale matrix, S, by   Sij={σi(P−1)iii=j0i≠j where σi>0 were variance parameters. Then the covariance matrix for the population process was obtained by   ΣN=SP−1ST Hence, any covariance matrix, ΣN, can be constructed from a suitable choice of L and σi, i=1,…,S·A. The marginal covariance matrix of the ages included in the assessments was found by using the corresponding indices of ΣN. In the example above, ΣN would be an 18 × 18 matrix where the indices three to nine corresponds to ages 2–8 for the first stock, while indices 10–17 corresponds to ages 0–7 for the second stock. We considered the models Mα, α=0,…,5, defined such that the off-diagonal partial correlation in the error terms of Ns,a,t and Ns′,a′,t was zero if s=s′ or |a−a′|≥α (See Appendix 1 for details). Hence, M0 corresponded to independent assessments while the between-stock correlation has a band structure for M1 to M5. The width of the band increases with the model subscript. When the process is Gaussian as here, a zero partial correlation implies conditional independence; hence, by setting all within-stock partial correlations to zero, any correlation within each stock is a result of partial correlations between the stocks. Therefore, the individual stock models are only changed as an effect of the between-stock connections introduced. Case studies We considered three case studies with two stocks each (Table 1). Two cases represent stocks that are known to migrate between management areas, while the last case represents stocks with intersecting management areas. The Baltic cod (Gadus morhua) case consists of the Eastern and Western Baltic cod stocks; the data was the basis of the 2014 ICES advice (ICES, 2014) for sub divisions 25–32 and 23–24, respectively (Figure 1). The Baltic herring (Clupea harengus membras) case comprised data for the Bothnian Sea and Central Baltic herring stocks used in the 2015 ICES advice (ICES, 2015a) for sub divisions 30 and 25–27, 28.2, 29, 32, respectively (Figure 1). Data for the Baltic mix case were obtained from the 2015 ICES advice for the Baltic sprat (Sprattus sprattus; sub divisions 22–32; ICES, 2015a) and Western Baltic herring (Clupea harengus membras; sub divisions 22–24; ICES, 2015b) stocks. Three of the stocks considered had missing data (Table 1). The Central Baltic Herring and Baltic Sprat data had missing data for all ages in surveys; the Western Baltic Herring was missing data for one age in 1991 for the Q3 survey, while the remaining years with missing data were missing all ages. Process model parameters were coupled between ages in the same way as they were for the advice. Figure 1. View largeDownload slide ICES sub divisions in the North Sea and Baltic Sea (http://geo.ices.dk; accessed 8 February 2017). Figure 1. View largeDownload slide ICES sub divisions in the North Sea and Baltic Sea (http://geo.ices.dk; accessed 8 February 2017). Simulation study We simulated datasets with 50 years from the stock assessment model above with Gaussian observations (See Appendix 2 for details). The datasets were simulated with two stocks combined by partial correlation models M0 to M5. For each model we simulated 150 datasets. To investigate whether the true model could be identified by Akaike information criterion (AIC; Akaike, 1974), and in particular whether AIC would favour overfitting, the models M0 to M5 were fitted to the data. Further, the coverage of CIs was compared between the models. Results Case studies Including correlations in the survival process between stocks improved the model fit, as measured by AIC, in all three case studies (Table 2). For Baltic Cod, M2 had the lowest AIC followed by M3 with an AIC that was 3.04 higher. Model M3 achieved the lowest AIC for Baltic Mix. Here, M2 followed with an AIC difference of 6.59. For Baltic Herring, M1 had the lowest AIC, followed by M3 with a difference of 7.48. Table 2 AIC difference from the best model fit for each case study, with the number of free parameters needed to obtain the correlation structure in parentheses. Model  Baltic cod  Baltic herring  Baltic mix  M0  35.06 (0)  32.99 (0)  14.26 (0)  M1  37.52 (6)  0.00 (8)  24.90 (8)  M2  0.00 (18)  10.96 (23)  6.59 (23)  M3  3.04 (29)  7.48 (36)  0.00 (36)  M4  15.91 (38)  15.24 (47)  16.76 (47)  M5  24.94 (45)  38.27 (56)  13.07 (56)  Model  Baltic cod  Baltic herring  Baltic mix  M0  35.06 (0)  32.99 (0)  14.26 (0)  M1  37.52 (6)  0.00 (8)  24.90 (8)  M2  0.00 (18)  10.96 (23)  6.59 (23)  M3  3.04 (29)  7.48 (36)  0.00 (36)  M4  15.91 (38)  15.24 (47)  16.76 (47)  M5  24.94 (45)  38.27 (56)  13.07 (56)  For the Baltic Cod stocks (Figure 2), the age difference in recruitment was 2 years. With the recruitments 2 years apart, the partial correlation was zero by construction in M3. Further, the recruitment to Eastern Baltic was estimated to be almost partially uncorrelated with Western Baltic age one (Estimate: 0.02; lower: −0.03; upper: 0.05). The resulting correlation in recruitment was also low (Estimate: 0.15; lower: −0.51; upper: 0.50). However, both the partial correlation (Estimate: 0.78; lower: 0.12; upper: 0.97) and correlation (Estimate: 0.93; lower: 0.15; upper: 1.00) between ages two were estimated to be high for the two stocks. For Eastern Baltic Cod, close age groups were estimated to be highly correlated as a result of the partial correlations with Western Baltic Cod. These within-stock correlations between age groups no more than 2 years apart were estimated to be above 0.6, except for ages 6 and 8 where the correlation was 0.48 (0.17; 0.61) and ages 2 and 4 where the correlation was 0.48 (−0.38; 0.92). A majority of the correlations between the stocks, 89%, were estimated to be positive. Of these, 81% were above 0.25, and 53% were above 0.6. With respectively 7 and 8 age groups in the two Baltic Cod stocks, 18 partial correlations were estimated for M3 based on the same number of free parameters in the Choleskey decomposition. Of these estimates, 22% were estimated to be negative, 17% had an absolute value below 0.1, and 28% had an absolute value above 0.7. Figure 2. View largeDownload slide Estimated partial correlation (upper triangle) and corresponding correlation (lower triangle) in non-fishery survival, Σ^N, for model M2 in the Baltic Cod case study between all stocks and ages. By construction, both correlation and partial correlation is 1 in the diagonal. The matrix is extended to have the same ages in both stocks. Ages present in the data are represented by a black rectangle in the diagonal. 95% CIs are included in parentheses (see Appendix 3 for details). A black triangle in the upper right corner indicates that the (partial) correlation is significant. Axis labels indicate the age group, e.g. A3 is age three in the first stock, while B2 is age two in the second stock. Figure 2. View largeDownload slide Estimated partial correlation (upper triangle) and corresponding correlation (lower triangle) in non-fishery survival, Σ^N, for model M2 in the Baltic Cod case study between all stocks and ages. By construction, both correlation and partial correlation is 1 in the diagonal. The matrix is extended to have the same ages in both stocks. Ages present in the data are represented by a black rectangle in the diagonal. 95% CIs are included in parentheses (see Appendix 3 for details). A black triangle in the upper right corner indicates that the (partial) correlation is significant. Axis labels indicate the age group, e.g. A3 is age three in the first stock, while B2 is age two in the second stock. For Baltic Herring, the recruitments were estimated to be highly partially correlated in M1 (Figure 3) with an estimated coefficient of 0.94 (0.58; 0.99). For model M1, the partial correlations are equal to the correlations by construction. As with the recruitment, ages two had a high correlation. This correlation was estimated to be highly negative (Estimate: −0.69) albeit with an uncertain estimate (95% CI: −0.99; 0.96). Ages 5 and 8 were both estimated to be positively correlated with coefficients of 0.74 and 0.55, respectively. However, both 95% CIs included zero (Figure 3). Likewise, the remaining estimated correlations had CIs which included zero. With respectively 8 and 9 ages in the two Baltic Herring stocks, eight partial correlations were estimated for M1 using the same number of free parameters in the Choleskey decomposition. Figure 3. View largeDownload slide Estimated partial correlation (upper triangle) and corresponding correlation (lower triangle) in non-fishery survival, Σ^N, for model M1 in the Baltic Herring case study between all stocks and ages. By construction, both correlation and partial correlation is 1 in the diagonal. The matrix is extended to have the same ages in both stocks. Ages present in the data are represented by a black rectangle in the diagonal. 95% CIs are included in parentheses (see Appendix 3 for details). A black triangle in the upper right corner indicates that the (partial) correlation is significant. Axis labels indicate the age group, e.g. A3 is age three in the first stock, while B2 is age two in the second stock. Figure 3. View largeDownload slide Estimated partial correlation (upper triangle) and corresponding correlation (lower triangle) in non-fishery survival, Σ^N, for model M1 in the Baltic Herring case study between all stocks and ages. By construction, both correlation and partial correlation is 1 in the diagonal. The matrix is extended to have the same ages in both stocks. Ages present in the data are represented by a black rectangle in the diagonal. 95% CIs are included in parentheses (see Appendix 3 for details). A black triangle in the upper right corner indicates that the (partial) correlation is significant. Axis labels indicate the age group, e.g. A3 is age three in the first stock, while B2 is age two in the second stock. In the Baltic Mix case (Figure 4), recruitments were estimated (95% CI) to have a negative partial correlation in M3 of −0.53 (−0.77; −0.07). The resulting correlation in recruitment was also highly negative with an estimated coefficient of −0.65 (−0.94; −0.03). Several age classes were estimated to be highly correlated. Baltic Spart age 3 and Western Baltic Herring age 2 had an estimated correlation of −1.00 (−1.00; −0.49); Sprat age 5 had an estimated correlation with Herring ages 3–5 of 0.91 (−0.85; 0.99), 0.99 (0.59; 1.00), and 0.93 (0.72; 1.00) respectively; and Sprat age 7 was estimated to be correlated with Herring ages 5 and 6 with coefficients of 0.93 (0.55; 0.99) and 0.99 (0.29; 1.00). With 8 and 9 ages respectively in the two Baltic Mix stocks, 36 partial correlations were estimated for M3 using the same number of free parameters in the Choleskey decomposition. Of these partial correlation estimates, 39% were estimated to be negative, 47% had an absolute value below 0.1, and 14% had an absolute value above 0.7. Figure 4. View largeDownload slide Estimated partial correlation (upper triangle) and corresponding correlation (lower triangle) in non-fishery survival, Σ^N, for model M3 in the Baltic Mix case study between all stocks and ages. By construction, both correlation and partial correlation is 1 in the diagonal. The matrix is extended to have the same ages in both stocks. Ages present in the data are represented by a black rectangle in the diagonal. 95% CIs are included in parentheses (see Appendix 3 for details). A black triangle in the upper right corner indicates that the (partial) correlation is significant. Axis labels indicate the age group, e.g. A3 is age three in the first stock, while B2 is age two in the second stock. Figure 4. View largeDownload slide Estimated partial correlation (upper triangle) and corresponding correlation (lower triangle) in non-fishery survival, Σ^N, for model M3 in the Baltic Mix case study between all stocks and ages. By construction, both correlation and partial correlation is 1 in the diagonal. The matrix is extended to have the same ages in both stocks. Ages present in the data are represented by a black rectangle in the diagonal. 95% CIs are included in parentheses (see Appendix 3 for details). A black triangle in the upper right corner indicates that the (partial) correlation is significant. Axis labels indicate the age group, e.g. A3 is age three in the first stock, while B2 is age two in the second stock. Simulation study The simulation study consisted of a total of 900 simulations from the six models M0 to M5. Of these simulations, 81% of the re-estimations identified the correct model from the candidates M0 to M5 (Table 3). Simulations from M3 were the most difficult to identify. Of these, 70% were correctly identified as generated by M3 whereas 21% were identified as M4, 5% as M5, and 3% were identified as M2. Only 7.1% of all simulations obtained the lowest AIC with a smaller model than the true simulation model. Table 3 Number of simulation studies that resulted in the lowest AIC for models M0 to M5.   Best model   True Model  M0  M1  M2  M3  M4  M5  M0  136  12  1  1  0  0  M1  0  140  6  1  2  1  M2  0  0  120  14  8  8  M3  0  0  5  105  32  8  M4  0  0  5  21  107  17  M5  0  0  0  0  33  117    Best model   True Model  M0  M1  M2  M3  M4  M5  M0  136  12  1  1  0  0  M1  0  140  6  1  2  1  M2  0  0  120  14  8  8  M3  0  0  5  105  32  8  M4  0  0  5  21  107  17  M5  0  0  0  0  33  117  Although the AIC was able to identify the true model in most cases, spurious correlations were found by the complex models when the true model did not include correlations. When simulating from M0, all the models M1 to M5 had high estimated correlations. For M5, 67% of the simulations had one or more estimated correlations with an absolute value above 0.5 while 10% had an absolute value above 0.75. The number was reduced with model complexity. The models M1 to M4, had estimated correlations with absolute values above 0.5 for 44, 57, 65, and 67% of the simulations, respectively. Likewise, respectively 3, 9, 11, and 10% of the simulations had estimated correlations with absolute values above 0.75. Besides identifying the true model, the simulations were used to calculate the coverage of Hessian based 95% CIs of the logarithm of the average fishing mortality (Table 4). When simulating from the models M0 and M1, all four estimation models had an average coverage close to 95% based on 150 simulations. This was not the case in the simulations where substantial correlations were included in the true model. When simulating from M2 to M5, estimation model M0 only had a coverage of 91.6, 89.7, 88.5, and 91.9%, respectively. Likewise, when simulating from M5, estimation models M1 to M3 had coverages of 91.7, 89.7, 90.2%, respectively, while models M4 and M5 had higher coverages: 93.4 and 93.2%, respectively. Table 4 Average (over time) coverage of Hessian based 95% CIs of the logarithm of the predicted average (over ages) fishing mortality in the simulation study.   Estimation model   Simulation Model  M0  M1  M2  M3  M4  M5  M0  0.945  0.943  0.941  0.940  0.940  0.940  M1  0.942  0.940  0.939  0.936  0.935  0.934  M2  0.916  0.933  0.939  0.936  0.937  0.936  M3  0.897  0.930  0.939  0.934  0.933  0.931  M4  0.885  0.921  0.929  0.928  0.929  0.931  M5  0.919  0.917  0.897  0.902  0.934  0.932    Estimation model   Simulation Model  M0  M1  M2  M3  M4  M5  M0  0.945  0.943  0.941  0.940  0.940  0.940  M1  0.942  0.940  0.939  0.936  0.935  0.934  M2  0.916  0.933  0.939  0.936  0.937  0.936  M3  0.897  0.930  0.939  0.934  0.933  0.931  M4  0.885  0.921  0.929  0.928  0.929  0.931  M5  0.919  0.917  0.897  0.902  0.934  0.932  Discussion European fisheries management mainly relies on single-stock assessment models, which assumes that stocks are not interconnected. For many stocks, this is not a reasonable assumption, as they are connected by several interactions, such as predation, migration and influences from the environment. Currently, the main alternative to single-stock models is complex ecosystem models, which are data intensive and require extensive tuning and long run times. We introduced a simple and fast alternative: connecting stocks by correlated survival variability. In our case studies, the models were fitted in 1.1–4.5 min. The models only require the data already used in single-stock models, which makes it readily usable for a wide range of stocks. Further, the model can be reduced to independent single-stock models, which allows objective testing of whether connecting stocks is an improvement, through likelihood ratio tests, AIC, or the like. We modelled complicated correlation patterns by a simple structure through the partial correlations; however, other correlation structures could be used instead. Unlike multi-species models and ecosystem models, which try to directly model interactions or common influences of stocks, the present approach estimates whether a good year for one cohort is also a good year for another cohort it may be linked to. Although simplistic and less informative about the ecosystem, this approach is highly operational and easily implemented for any group of stocks with an analytical assessment, since the model does not require any new data sources to be introduced. Further, the approach can be used to model complex correlation structures in other parts of the model such as the fishing mortality process. The case studies presented were chosen to illustrate stocks that are believed to migrate and stocks that share environment. In general, connecting stocks should be considered when: biological knowledge suggests interactions, the management areas intersect, or patterns in model validation tools such as residuals (e.g. Thygesen et al., 2017) or retrospective analysis (e.g. Hurtado-Ferro et al., 2015) are present. If biological knowledge indicates migration or predator-prey relations, it should be tested whether assessments can be improved by correlating the abundance processes. Likewise, when the management areas intersect, the stocks likely share environmental conditions or compete for resources which can correlate the abundance processes. Finally, common or opposite patterns in the residuals or retrospective analysis can indicate missing biological knowledge which can lead to correlated abundance processes. Although modelling partial correlations provides a simple interpretation of parameters, understanding why they arise is not as simple in the present framework. Correlations among stocks or ages can be the result of several processes such as migration, predation or environmental impact. Although further research is needed to evaluate the possibility of extending the framework to directly include any of these processes, it is straight forward to evaluate whether stocks are similar in the sense that they could share parameters in the assessment model. Testing whether parameters are shared among stocks can give a basis for evaluating if stocks should be managed together or separately (Montenegro et al., 2009). For all three case studies, the proposed model improved on the single-stock models as measured by AIC. In general, model M3 performed well for all stocks, although smaller models were preferred for two of the cases. The model allows a flexible correlation structure to be estimated with relatively few parameters. By modelling partial correlations with a band structure between stocks, we obtain very general correlation structures, even within the individual stocks. With this model, correlations within each stock are interpreted as a result of a connection between the stocks. Hence, the model only changes the single-stock model as an effect of the connection between stocks. In the simulation study we showed that the true model can be identified by AIC, although with a risk of choosing a model including too many partial correlations. It is well-known that AIC can favour overfitting, when the true model has a finite parameter space, and is among the candidate models (Burnham and Anderson, 2003). However, for stock assessment models the true data generating system is too complex to be described by a simple model; observed data will depend on, amongst other things, the hydrological systems, individual fish behaviour, individual fisherman behaviour, and a wide array of data accumulation procedures. Any practically applicable candidate model will require strong simplifications; therefore, the AIC does not necessarily select the true model in the case studies, but simply the model that has the best fit to the data, while penalizing for the additional parameters used. Although the selected model provides the best fit to data, interpretation of individual correlation parameters must be done with caution. Estimated correlations in the abundance process will be a result of several effects such as environmental conditions, predation and migration, which cannot be separated in the presented model. Further, if the correlation between e.g. observations are not adequately modelled, the state-space model can somewhat compensate to increase the model fit using the correlations in the abundance process (Thygesen et al., 2017). Finally, the typically short time series available for stock assessment models can result in imprecise individual correlation estimates, which clutters the overall pattern seen. Ignoring correlations in data result in unreliable CIs for model output. In a simulation study, we investigated the CIs of the estimated fishing mortality; both when correlations were present and absent. We found that when severe correlations are present in the data, (i.e. simulating from models M3 to M5), ignoring the correlations results in CIs that are too narrow. When no correlations were estimated, Hessian based 95% CIs only had a coverage of 89.7, 88.5, and 91.9%, respectively, when simulating from models M3 to M5, while the more flexible correlation structures provided CIs close to the correct coverage. Accurate assessments of model output uncertainty is important for sustainable fisheries management (Chaput, 2004; Cadrin and Dickey-Collas, 2014). Underestimating the uncertainty of model output used for management will increase the probability of overfishing (Punt et al., 2012). Hence, if correlations in data are ignored, a fishery perceived as sustainable may in fact have a non-sustainable risk of depletion. These results will allow managers and policymakers to make informed decisions based on a more accurate assessment of model output uncertainty. As a simple extension of single-species assessments, the model output can be evaluated per stock, while giving an accurate evaluation of uncertainty. Likewise, connecting stocks may increase the accuracy of forecasts. Together with the estimated correlations, these improvements can provide guidance on managing single stocks and the potential influence on surrounding stocks. The procedure outlined is directly applicable without additional data sources. Acknowledgements The authors thank handling editor Dr Shijie Zhou and two anonymous reviewers for their valuable comments to improve this article. Funding Funding for this study was provided by a grant from the European Fisheries and Maritime Fund (33113-B-15-003, Ministry of Environment and Food in Denmark). References Ainsworth C., Walters C. 2015. Ten common mistakes made in ecopath with ecosim modelling. Ecological Modelling , 308: 14– 17. Google Scholar CrossRef Search ADS   Akaike H. 1974. A new look at the statistical model identification. IEEE Transactions on Automatic Control , 19: 716– 723. Google Scholar CrossRef Search ADS   Albertsen C. M., Nielsen A., Thygesen U. H. 2017. Choosing the observational likelihood in state-space stock assessment models. Canadian Journal of Fisheries and Aquatic Sciences , 74: 779– 789. Google Scholar CrossRef Search ADS   Begley J., Howel D. 2004. An overview of gadget, the globally applicable area-disaggregated general ecosystem toolbox. Technical Report, ICES document CM/FF:13. Burnham K., Anderson D. 2003. Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach . Springer, New York. Cadrin S. X., Dickey-Collas M. 2014. Stock assessment methods for sustainable fisheries. ICES Journal of Marine Science , 72: 1– 6. Google Scholar CrossRef Search ADS   Chaput G. 2004. Considerations for using spawner reference levels for managing single- and mixed-stock fisheries of atlantic salmon. ICES Journal of Marine Science , 61: 1379– 1388. Google Scholar CrossRef Search ADS   Chen Z., Xu S., Qiu Y., Lin Z., Jia X. 2009. Modeling the effects of fishery management and marine protected areas on the beibu gulf using spatial ecosystem simulation. Fisheries Research , 100: 222– 229. Google Scholar CrossRef Search ADS   Christensen V., Pauly D. 1992. ECOPATH II — a software for balancing steady-state ecosystem models and calculating network characteristics. Ecological Modelling , 61: 169– 185. Google Scholar CrossRef Search ADS   Christensen V., Walters C. J. 2004. Ecopath with ecosim: methods, capabilities and limitations. Ecological Modelling , 172: 109– 139. Google Scholar CrossRef Search ADS   Collie J. S., Botsford L. W., Hastings A., Kaplan I. C., Largier J. L., Livingston P. A., Plagányi É., et al.   2016. Ecosystem models for fisheries management: finding the sweet spot. Fish and Fisheries , 17: 101– 125. Google Scholar CrossRef Search ADS   Daan N. 1987. Multispecies versus single-species assessment of north sea fish stocks. Canadian Journal of Fisheries and Aquatic Sciences , 44: s360– s370. doi:10.1139/f87-337. Google Scholar CrossRef Search ADS   Eero M., Hemmer-Hansen J., Hussy K. 2014. Implications of stock recovery for a neighbouring management unit: experience from the baltic cod. ICES Journal of Marine Science , 71: 1458– 1466. Google Scholar CrossRef Search ADS   Frøysa K. G., Bogstad B., Skagen D. W. 2002. Fleksibest—an age-length structured fish stock assessment model. Fisheries Research , 55: 87– 101. Google Scholar CrossRef Search ADS   Fulton E. A., Smith A. D. M., Punt A. E. 2004. Ecological indicators of the ecosystem effects of fishing: Final report. Report no. r99/1546. Technical Report, Australian Fisheries Management Authority, Canberra. Gislason H., Helgason T. 1985. Species interaction in assessment of fish stocks with special application to the north sea. Dana  , 5: 1– 44. Hannesson S., Jakobsdottir A., Begley J., Taylor L., Stefansson G. 2008. On the use of tagging data in statistical multispecies multi-area models of marine populations. ICES Journal of Marine Science , 65: 1762– 1772. Google Scholar CrossRef Search ADS   Hollowed A. 2000a. Are multispecies models an improvement on single-species models for measuring fishing impacts on marine ecosystems?. ICES Journal of Marine Science , 57: 707– 719. Google Scholar CrossRef Search ADS   Hollowed A. 2000b. Including predation mortality in stock assessments: a case study for gulf of alaska walleye pollock. ICES Journal of Marine Science , 57: 279– 293. Google Scholar CrossRef Search ADS   Holsman K. K., Ianelli J., Aydin K., Punt A. E., Moffitt E. A. 2016. A comparison of fisheries biological reference points estimated from temperature-specific multi-species and single-species climate-enhanced stock assessment models. Deep Sea Research Part II: Topical Studies in Oceanography , 134: 360– 378. Google Scholar CrossRef Search ADS   Hurtado-Ferro F., Szuwalski C. S., Valero J. L., Anderson S. C., Cunningham C. J., Johnson K. F., Licandeo R., et al.   2015. Looking in the rear-view mirror: bias and retrospective patterns in integrated, age-structured stock assessment models. ICES Journal of Marine Science , 72: 99– 110. Google Scholar CrossRef Search ADS   ICES 2014. Report of the baltic fisheries assessment working group (WGBFAS), 3-10 April 2014, ICES HQ, Copenhagen, Denmark. ICES CM 2014/ACOM:10, ICES. ICES 2015a. Report of the Baltic Fisheries Assessment Working Group (WGBFAS), 14-21 April 2014, ICES HQ, Copenhagen, Denmark. ICES CM 2015/ACOM:10, ICES. ICES 2015b. Report of the Herring Assessment Working Group for the Area South of 62°N (HAWG), 10-19 March 2015, ICES HQ, Copenhagen, Denmark. ICES CM 2015/ACOM:06, ICES. Kristensen K., Nielsen A., Berg C. W., Skaug H., Bell B. 2016. Tmb: automatic differentiation and laplace approximation. Journal of Statistical Software , 70: 1– 21. Google Scholar CrossRef Search ADS   Lewy P., Vinther M. 2004. A stochastic age-length-structured multispecies model applied to north sea stocks. Technical Report, ICES document CM/FF:20. Link J. S., Fulton E. A., Gamble R. J. 2010. The northeast US application of ATLANTIS: A full system model exploring marine ecosystem dynamics in a living marine resource management context. Progress in Oceanography , 87: 214– 234. Google Scholar CrossRef Search ADS   Methot R. D., Wetzel C. R. 2013. Stock synthesis: A biological and statistical framework for fish stock assessment and fishery management. Fisheries Research , 142: 86– 99. Google Scholar CrossRef Search ADS   Montenegro C., Maunder M. N., Zilleruelo M. 2009. Improving management advice through spatially explicit models and sharing information. Fisheries Research , 100: 191– 199. Google Scholar CrossRef Search ADS   Nielsen A., Berg C. 2014. Estimation of time-varying selectivity in stock assessments using state-space models. Fisheries Research , 158: 96– 101. Google Scholar CrossRef Search ADS   Punt A. E., Siddeek M. S. M., Garber-Yonts B., Dalton M., Rugolo L., Stram D., Turnock B. J., et al.   2012. Evaluating the impact of buffers to account for scientific uncertainty when setting TACs: application to red king crab in bristol bay, alaska. ICES Journal of Marine Science , 69: 624– 634. Google Scholar CrossRef Search ADS   Quinn T. J., Deriso R. B., Neal P. R. 1990. Migratory catch-age analysis. Canadian Journal of Fisheries and Aquatic Sciences , 47: 2315– 2327. Google Scholar CrossRef Search ADS   Richards R. A., Jacobson L. D. 2016. A simple predation pressure index for modeling changes in natural mortality: application to gulf of maine northern shrimp stock assessment. Fisheries Research , 179: 224– 236. Google Scholar CrossRef Search ADS   Senina I., Sibert J., Lehodey P. 2008. Parameter estimation for basin-scale ecosystem-linked population models of large pelagic predators: application to skipjack tuna. Progress in Oceanography , 78: 319– 335. Google Scholar CrossRef Search ADS   Tarnecki J. H., Wallace A. A., Simons J. D., Ainsworth C. H. 2016. Progression of a Gulf of Mexico food web supporting Atlantis ecosystem model development. Fisheries Research , 179: 237– 250. Google Scholar CrossRef Search ADS   Thygesen U. H., Albertsen C. M., Berg C. W., Kristensen K., Nielsen A. 2017. Validation of ecological state space models using the laplace approximation. Environmental and Ecological Statistics , 24: 317–339. Yodzis P. 1998. Local trophodynamics and the interaction of marine mammals and fisheries in the benguela ecosystem. Journal of Animal Ecology , 67: 635– 658. Google Scholar CrossRef Search ADS   Appendix 1: Unstructured covariance matrix Let L be a lower triangular matrix of parameters with diagonal 1, and normalized such that ⟨Li·,Li·⟩=1 for all i. Then Σ=LLT is an arbitrarily scaled covariance matrix. Off-diagonal elements of the covariance matrix can be forced to be zero (e.g. https://github.com/admb-project/examples/tree/master/admb-tricks/parameterization/covariance-matrices; accessed: 06 March 2017) by noting that   Σjk=⟨Lj·,Lk·⟩. Hence, by defining L˜ such that L˜j·=Lj·−⟨Lj·,Lk·⟩Lk· and L˜i·=Li· for i≠j, then for the covariance matrix Σ˜=L˜L˜T,   Σ˜jk=⟨L˜j·,L˜k·⟩=⟨Lj·−⟨Lj·,Lk·⟩Lk·,Lk·⟩=⟨Lj·,Lk·⟩−⟨Lj·,Lk·⟩⟨Lk·,Lk·⟩=⟨Lj·,Lk·⟩−⟨Lj·,Lk·⟩=0 Although this method allows us to impose constraints on the covariance matrix, it is, in practice, difficult to pre-determine which parameters will be redundant. Instead, the gradient and hessian of the resulting likelihood function can be examined to determine parameters that do not affect the likelihood value. These parameters must be fixed in optimization. Appendix 2: Simulation study details Stock structure We simulated 100 years from the model and discarded the first half of the time series. Two stocks were simulated with 5 age groups each. The last age group was considered a plus group. Natural mortality was fixed at 0.2. The initial log-population size was set to 16 for the first age group and was 0.5 less for each subsequent age group. Fishing mortality log ⁡-fishing mortality for stock s, age a, and year y was simulated by   log ⁡Fs,a,y= cos ⁡(y10+10·1{2}(s))−1+ϵs,a,y where   1{2}(s)={1 if s=20 if s≠2, and ϵs,a,y∼N(0,0.012). Log-population process The log ⁡-population given fishing mortality was simulated from the models M0 to M5 above. The off-diagonal non-zero entries, Pi,j, of the precision, P, matrix for Mα was constructed as −(α+1)/(α−|ai−aj|+1)2, where ai is the age corresponding to the ith index of the process. This construction ensures non-negligible partial correlations in the elements not included in the smaller models. Diagonal elements, Pi,i, were set to ∑j≠iPi,j+ exp ⁡(−α) to ensure the precision matrix was positive-definite. The corresponding covariance matrix was scaled to have diagonal elements 0.12. Observations log ⁡-observations given fishing mortality and log ⁡-population were simulated from a multivariate normal with AR(1) covariance structure. The AR(1) parameter was set to 0. For commercial data, the variances for all ages were 1. A first quarter survey was simulated from a multivariate normal with AR(1) covariance structure with parameter 0 and variances exp ⁡(−4). The catchability was set to exp ⁡(−7). Re-estimation Re-estimation was performed assuming an AR(1) structure for the correlations in the log ⁡-F random walk, random walk recruitment, and univariate log ⁡-normal observations. Appendix 3: CI for correlation parameters Simulation based CIs for correlation parameters were constructed through the profile likelihood. Splitting the parameter vector θ into correlation parameters θc and other parameters θo, we simulated proposal parameters θc(i), i=1,…,N, from the hessian based Gaussian approximation of the distribution of the estimator. The proposal was accepted if the hypothesis θc=θc(i) could not be rejected by a likelihood ratio test with significance level α. Each accepted simulated correlation parameter vector was transformed to a partial correlation and a correlation matrix. For each matrix element, minimum and maximum values were used as (1−α)·100 % CI. We simulated N=1000 proposal parameter vectors for each case study, and obtained acceptance rates of 0.80, 0.41, and 0.63 for Baltic Herring, Mix and Cod, respectively. For Baltic Mix and Cod, the variance of the proposal distribution was reduced to increase acceptance. © International Council for the Exploration of the Sea 2017. All rights reserved. For Permissions, please email: journals.permissions@oup.com http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png ICES Journal of Marine Science Oxford University Press

Connecting single-stock assessment models through correlated survival

Loading next page...
 
/lp/ou_press/connecting-single-stock-assessment-models-through-correlated-survival-NKR9vbFRV0
Publisher
ICES/CIEM
Copyright
© International Council for the Exploration of the Sea 2017. All rights reserved. For Permissions, please email: journals.permissions@oup.com
ISSN
1054-3139
eISSN
1095-9289
D.O.I.
10.1093/icesjms/fsx114
Publisher site
See Article on Publisher Site

Abstract

Abstract Fisheries management is mainly conducted via single-stock assessment models assuming that fish stocks do not interact, except through assumed natural mortalities. Currently, the main alternative is complex ecosystem models which require extensive data, are difficult to calibrate, and have long run times. We propose a simple alternative. In three case studies each with two stocks, we improve the single-stock models, as measured by Akaike information criterion, by adding correlation in the cohort survival. To limit the number of parameters, the correlations are parameterized through the corresponding partial correlations. We consider six models where the partial correlation matrix between stocks follows a band structure ranging from independent assessments to complex correlation structures. Further, a simulation study illustrates the importance of handling correlated data sufficiently by investigating the coverage of confidence intervals for estimated fishing mortality. The results presented will allow managers to evaluate stock statuses based on a more accurate evaluation of model output uncertainty. The methods are directly implementable for stocks with an analytical assessment and do not require any new data sources. Introduction Fisheries management is often based on single-species models for single management areas where the stocks are assumed to have no interaction with surrounding stocks or other species, except via the predetermined natural mortality. Such a simplification is often unrealistic (Eero et al., 2014), and ignoring stock interactions can affect model output (Holsman et al., 2016). These assumptions in single-stock models are, however, convenient as they reduce the data needed on interactions such as migration and predator-prey relations. Collecting fisheries data are costly and time-consuming: The current surveys require expensive scientific cruises and lab processing; tagging studies require catching and manually tagging many individuals with either low probability of recapture, or expensive electronic tags for gathering telemetry data; food-web data such as stomach samples require manual sampling of stomach content and possibly subsequent DNA analysis; spatially explicit catch data require extensive monitoring systems. Further, incorporating such new data sources in stock assessment models requires historic data that covers a long period for reliable statistical results. Several methods ranging in complexity have been developed to improve single-stock assessments by accounting for interactions between stocks (e.g. Gislason and Helgason, 1985; Christensen and Pauly, 1992; Begley and Howel, 2004; Fulton et al., 2004; Lewy and Vinther, 2004; Senina et al., 2008; Methot and Wetzel, 2013). Migration between areas can be explicitly modelled by a transportation matrix that moves, e.g. a proportion of the population to another area (Methot and Wetzel, 2013). These proportions can be informed by tagging studies such as mark-recapture (Quinn et al., 1990; Hannesson et al., 2008) or, for simple models, be estimated from catch and survey data; however, when estimated from catch and survey data, the results may be influenced by other factors such as inaccuracies in the assumed natural mortality or model misspecification. Likewise, predation can be modelled by a rate-of-consumption matrix (Hollowed, 2000b) consisting of factors which describe the predator induced mortality. Estimating predator-induced mortality requires knowledge of complex food webs (Yodzis, 1998; Tarnecki et al., 2016) along with consumption data from stomachs (e.g. Gislason and Helgason, 1985; Daan, 1987; Hollowed, 2000b; Frøysa et al., 2002; Lewy and Vinther, 2004; Richards and Jacobson, 2016). Complex ecosystem or multi-species models often depend on extensive data from different sources such as tagging studies, DNA analysis, environment data or stomach content data (e.g. Christensen and Walters, 2004). Such data sources are not readily available from single-species assessments, are expensive to produce, and may take several years to collect. Data requirements can also include human activity, biochemistry, physiology, geology and oceanography (Link et al., 2010) raising the complexity and data collection cost. Further, ecosystem models often have a vast amount of settings and can be difficult to calibrate and run (Hollowed, 2000a; Chen et al., 2009; Ainsworth and Walters, 2015). Currently, there is a gap between single-stock assessment models and the complex ecosystem models (Collie et al., 2016), with a lack of operational models that combine multiple stocks while only utilizing the data readily available for single-stock assessments, making it directly implementable for fisheries management. The objective of this study is to introduce a method for combining single-stock state-space age-based assessments models to provide a simple alternative to complex ecosystem models. The method is illustrated in three case studies including a total of six European stocks combined two by two. In each case study, we investigate whether combining single-stock models improves the model fit to data compared to individual assessments. We extend the single-stock age-based state-space stock assessment model to allow simultaneous estimation for several stocks through correlated survival processes. Correlations are included in the model through the corresponding partial correlations. Assuming Gaussianity, two random variables from a collection of variables are partially uncorrelated correcting for the remaining variables in the collection if and only if they are conditionally independent given the remaining variables in the collection. Hence, modelling partial correlations allows a clear interpretation of reducing model parameters (i.e. by fixing partial correlations at zero), while still allowing a flexible correlation structure. The multi-stock model introduced does not require further data than the single-stock model. Methods Stock assessment model We implemented age-based state-space stock assessment models for six European stocks (Table 1). The models were implemented with the R-package Template Model Builder (Kristensen et al., 2016). Table 1 Overview of the data sources used in the case studies. Fleet  First year  Last year  First age  Last age  Years with missing  Baltic Cod Case Study            Eastern Baltic Cod (Sub divisions: 25–32)            Commercial  1966  2013  2  8    Survey Q1  1991  2000  2  6    Survey Q1  2001  2013  2  6    Survey Q4  2002  2013  2  5    CPUE tuning index  1997  2013  3  7    Western Baltic Cod (Sub divisions: 23–24)            Commercial  1970  2013  1  7    CPUE tuning index  1997  2013  3  6    Survey Q1  1992  2013  1  6    Survey Q4  1992  2013  0  5    Baltic Herring Case Study            Bothnian Sea Herring (Sub division: 30)            Commercial  1973  2014  1  9    Survey Q4  2007  2014  2  8    Survey Q2  1990  2014  3  9    Central Baltic Herring (Sub divisions: 25–27, 28.2, 29, 32)            Commercial  1974  2014  1  8    Survey Q4  1991  2014  1  8  1993, 1995, 1997  Baltic Mix Case Study            Baltic Sprat (Sub divisions: 22–32)            Commercial  1991  2014  1  8    Survey Q4  1991  2014  1  8  1993, 1995, 1997  Survey Q2  2001  2014  1  8    Survey Q1  1992  2014  1  1  1994, 1996, 1998  Western Baltic Herring (Sub divisions: 22–24)            Commercial  1991  2014  0  8    Survey Q3  1991  2014  1  8  1991, 1999  Survey Q4  1991  2014  0  8  1991–1993, 2001  Survey Q2  1992  2014  0  0    Survey Q1  1991  2014  1  4    Survey Q3  1991  2014  1  4    Fleet  First year  Last year  First age  Last age  Years with missing  Baltic Cod Case Study            Eastern Baltic Cod (Sub divisions: 25–32)            Commercial  1966  2013  2  8    Survey Q1  1991  2000  2  6    Survey Q1  2001  2013  2  6    Survey Q4  2002  2013  2  5    CPUE tuning index  1997  2013  3  7    Western Baltic Cod (Sub divisions: 23–24)            Commercial  1970  2013  1  7    CPUE tuning index  1997  2013  3  6    Survey Q1  1992  2013  1  6    Survey Q4  1992  2013  0  5    Baltic Herring Case Study            Bothnian Sea Herring (Sub division: 30)            Commercial  1973  2014  1  9    Survey Q4  2007  2014  2  8    Survey Q2  1990  2014  3  9    Central Baltic Herring (Sub divisions: 25–27, 28.2, 29, 32)            Commercial  1974  2014  1  8    Survey Q4  1991  2014  1  8  1993, 1995, 1997  Baltic Mix Case Study            Baltic Sprat (Sub divisions: 22–32)            Commercial  1991  2014  1  8    Survey Q4  1991  2014  1  8  1993, 1995, 1997  Survey Q2  2001  2014  1  8    Survey Q1  1992  2014  1  1  1994, 1996, 1998  Western Baltic Herring (Sub divisions: 22–24)            Commercial  1991  2014  0  8    Survey Q3  1991  2014  1  8  1991, 1999  Survey Q4  1991  2014  0  8  1991–1993, 2001  Survey Q2  1992  2014  0  0    Survey Q1  1991  2014  1  4    Survey Q3  1991  2014  1  4    Q1–Q4 indicates at which quarter of the year the survey is conducted. Process model The process model for an individual stock, s, was identical to Nielsen and Berg (2014) model D. The population process followed an exponential decay model,   log ⁡Ns,a,y= log ⁡Ns,a−1,y−1−Fs,a−1,y−1−Ms,a−1,y−1+ηs,a,y,a∈{2,3,…,As+−2,As+−1} with special conditions for the last age group As+ (See Nielsen and Berg, 2014; Albertsen et al., 2017, for details). In the model, Fs,a,y is the fishing mortality, Ms,a,y is the natural mortality, and ηs,a,y is a random Gaussian zero mean increment, which allow inaccuracies in the deterministic part of the model. These inaccuracies can be thought of as random deviations in the natural mortality, which is assumed to be known. While Nielsen and Berg (2014) and Albertsen et al. (2017) model ηs,a,y as independent, we will allow correlations across stocks and ages (see “Connecting stocks” Section), while maintaining independence across years. Recruitment, log ⁡Ns,1,y, was modelled by a random walk. Configuration of model parameters with respect to coupling across ages was done in the same way as for the official advice. Observational model Following Albertsen et al. (2017), the observed log-catch conditional on the process model were modelled by a multivariate normal where the covariance Σs(O) had an AR(1) covariance structure, (Σs(O))i,j=ρo,s|i−j|τo,s,iτo,s,j with ρo,s∈(−1,1) and 0<τo,s,i for all indices i. The catch vector, Cs,y=(Cs,1,y,…,Cs,As+,y), for stock s in year y followed the Baranov catch equation   log ⁡Cs,y= log ⁡Fs,y− log ⁡(Fs,y+Ms,y)+ log ⁡Ns,y+ log ⁡(1− exp ⁡(−Fs,y−Ms,y))+ϵs,y(O) were   ϵs,y(O)∼N(0,Σs(O)),   Ns,y=(Ns,1,y,…,Ns,As+,y),  Fs,y=(Fs,1,y,…,Fs,As+,y), Ms,y=(Ms,1,y,…,Ms,As+,y), and all operations are element-wise. Likewise, survey indices were modelled by   log ⁡Is,y(i)= log ⁡qs(i)+ log ⁡Ns,y−(Fs,y+Ms,y)ds(i)/365+ϵs,y(I,i) where ϵs,y(I,i)∼N(0,Σs(I,i)), Is,y(i) is a vector of age observations for the ith survey in stock s year y, q(i)>0, dS(i) is the number of days into the year the survey is conducted, and Σs(I,i) has an AR(1) structure. For simplicity, years with missing data for single ages were treated as if they were missing all ages for that fleet. Alternatively, missing data could be handled by, for instance, imputing random effects. All observational parameters were estimated; however, they were configured with respect to coupling across ages in the same way as for the official advice. Connecting stocks For many stocks, assuming independence between the error terms, ηs,a,y, is unrealistic. Inaccuracies in the assumed natural mortality will be correlated by connections not included in the model such as environmental changes, migration, or predation. All other things held constant, environmental effects will be seen as positive correlation; migration will be seen as negative correlation; while predation will be seen as negative correlation in a predator-prey relation or as positive correlation if two stocks are influenced by the same predator. To connect the S individual stocks, we temporarily extended the stock log-population processes to have the same number of ages, A, for convenience. As an example, consider a case with two stocks (such as the Baltic cod case study presented in “Case studies” Section) where one stock includes ages 2–8, while the other stock includes ages 0–7. To construct the correlation matrix, both log-population processes have to be temporarily extended to include ages 0–8 to combine the stocks. The individual processes were stacked to one process,   Ny=(N1,1,y,…,N1,A,y,…,NS,1,y,…,NS,A,y)T, for which we modelled the error terms, ηy=(η1,1,y,…,η1,A,y,…,ηS,1,y,…,ηS,A,y)T, by a multivariate normal with covariance ΣN. The covariance was parameterized through the Cholesky decomposition of the precision matrix PN. Defining L to be a lower triangular matrix of parameters with diagonal 1, a positive definite precision matrix was obtained by   PN=LLT If the precision matrix was scaled to have diagonal 1, the off-diagonal elements were the partial correlations with opposite sign. With this construction, the diagonal elements of P could not be estimated freely. To obtain the covariance matrix as the inverse precision, we scaled P to have diagonal 1 and added freely estimated variance parameters. Including variance parameters at this stage maintained the interpretation of the estimated parameters from the single-stock model, where the parameter is directly related to the variance of the specific age. The scaling was done by defining a diagonal scale matrix, S, by   Sij={σi(P−1)iii=j0i≠j where σi>0 were variance parameters. Then the covariance matrix for the population process was obtained by   ΣN=SP−1ST Hence, any covariance matrix, ΣN, can be constructed from a suitable choice of L and σi, i=1,…,S·A. The marginal covariance matrix of the ages included in the assessments was found by using the corresponding indices of ΣN. In the example above, ΣN would be an 18 × 18 matrix where the indices three to nine corresponds to ages 2–8 for the first stock, while indices 10–17 corresponds to ages 0–7 for the second stock. We considered the models Mα, α=0,…,5, defined such that the off-diagonal partial correlation in the error terms of Ns,a,t and Ns′,a′,t was zero if s=s′ or |a−a′|≥α (See Appendix 1 for details). Hence, M0 corresponded to independent assessments while the between-stock correlation has a band structure for M1 to M5. The width of the band increases with the model subscript. When the process is Gaussian as here, a zero partial correlation implies conditional independence; hence, by setting all within-stock partial correlations to zero, any correlation within each stock is a result of partial correlations between the stocks. Therefore, the individual stock models are only changed as an effect of the between-stock connections introduced. Case studies We considered three case studies with two stocks each (Table 1). Two cases represent stocks that are known to migrate between management areas, while the last case represents stocks with intersecting management areas. The Baltic cod (Gadus morhua) case consists of the Eastern and Western Baltic cod stocks; the data was the basis of the 2014 ICES advice (ICES, 2014) for sub divisions 25–32 and 23–24, respectively (Figure 1). The Baltic herring (Clupea harengus membras) case comprised data for the Bothnian Sea and Central Baltic herring stocks used in the 2015 ICES advice (ICES, 2015a) for sub divisions 30 and 25–27, 28.2, 29, 32, respectively (Figure 1). Data for the Baltic mix case were obtained from the 2015 ICES advice for the Baltic sprat (Sprattus sprattus; sub divisions 22–32; ICES, 2015a) and Western Baltic herring (Clupea harengus membras; sub divisions 22–24; ICES, 2015b) stocks. Three of the stocks considered had missing data (Table 1). The Central Baltic Herring and Baltic Sprat data had missing data for all ages in surveys; the Western Baltic Herring was missing data for one age in 1991 for the Q3 survey, while the remaining years with missing data were missing all ages. Process model parameters were coupled between ages in the same way as they were for the advice. Figure 1. View largeDownload slide ICES sub divisions in the North Sea and Baltic Sea (http://geo.ices.dk; accessed 8 February 2017). Figure 1. View largeDownload slide ICES sub divisions in the North Sea and Baltic Sea (http://geo.ices.dk; accessed 8 February 2017). Simulation study We simulated datasets with 50 years from the stock assessment model above with Gaussian observations (See Appendix 2 for details). The datasets were simulated with two stocks combined by partial correlation models M0 to M5. For each model we simulated 150 datasets. To investigate whether the true model could be identified by Akaike information criterion (AIC; Akaike, 1974), and in particular whether AIC would favour overfitting, the models M0 to M5 were fitted to the data. Further, the coverage of CIs was compared between the models. Results Case studies Including correlations in the survival process between stocks improved the model fit, as measured by AIC, in all three case studies (Table 2). For Baltic Cod, M2 had the lowest AIC followed by M3 with an AIC that was 3.04 higher. Model M3 achieved the lowest AIC for Baltic Mix. Here, M2 followed with an AIC difference of 6.59. For Baltic Herring, M1 had the lowest AIC, followed by M3 with a difference of 7.48. Table 2 AIC difference from the best model fit for each case study, with the number of free parameters needed to obtain the correlation structure in parentheses. Model  Baltic cod  Baltic herring  Baltic mix  M0  35.06 (0)  32.99 (0)  14.26 (0)  M1  37.52 (6)  0.00 (8)  24.90 (8)  M2  0.00 (18)  10.96 (23)  6.59 (23)  M3  3.04 (29)  7.48 (36)  0.00 (36)  M4  15.91 (38)  15.24 (47)  16.76 (47)  M5  24.94 (45)  38.27 (56)  13.07 (56)  Model  Baltic cod  Baltic herring  Baltic mix  M0  35.06 (0)  32.99 (0)  14.26 (0)  M1  37.52 (6)  0.00 (8)  24.90 (8)  M2  0.00 (18)  10.96 (23)  6.59 (23)  M3  3.04 (29)  7.48 (36)  0.00 (36)  M4  15.91 (38)  15.24 (47)  16.76 (47)  M5  24.94 (45)  38.27 (56)  13.07 (56)  For the Baltic Cod stocks (Figure 2), the age difference in recruitment was 2 years. With the recruitments 2 years apart, the partial correlation was zero by construction in M3. Further, the recruitment to Eastern Baltic was estimated to be almost partially uncorrelated with Western Baltic age one (Estimate: 0.02; lower: −0.03; upper: 0.05). The resulting correlation in recruitment was also low (Estimate: 0.15; lower: −0.51; upper: 0.50). However, both the partial correlation (Estimate: 0.78; lower: 0.12; upper: 0.97) and correlation (Estimate: 0.93; lower: 0.15; upper: 1.00) between ages two were estimated to be high for the two stocks. For Eastern Baltic Cod, close age groups were estimated to be highly correlated as a result of the partial correlations with Western Baltic Cod. These within-stock correlations between age groups no more than 2 years apart were estimated to be above 0.6, except for ages 6 and 8 where the correlation was 0.48 (0.17; 0.61) and ages 2 and 4 where the correlation was 0.48 (−0.38; 0.92). A majority of the correlations between the stocks, 89%, were estimated to be positive. Of these, 81% were above 0.25, and 53% were above 0.6. With respectively 7 and 8 age groups in the two Baltic Cod stocks, 18 partial correlations were estimated for M3 based on the same number of free parameters in the Choleskey decomposition. Of these estimates, 22% were estimated to be negative, 17% had an absolute value below 0.1, and 28% had an absolute value above 0.7. Figure 2. View largeDownload slide Estimated partial correlation (upper triangle) and corresponding correlation (lower triangle) in non-fishery survival, Σ^N, for model M2 in the Baltic Cod case study between all stocks and ages. By construction, both correlation and partial correlation is 1 in the diagonal. The matrix is extended to have the same ages in both stocks. Ages present in the data are represented by a black rectangle in the diagonal. 95% CIs are included in parentheses (see Appendix 3 for details). A black triangle in the upper right corner indicates that the (partial) correlation is significant. Axis labels indicate the age group, e.g. A3 is age three in the first stock, while B2 is age two in the second stock. Figure 2. View largeDownload slide Estimated partial correlation (upper triangle) and corresponding correlation (lower triangle) in non-fishery survival, Σ^N, for model M2 in the Baltic Cod case study between all stocks and ages. By construction, both correlation and partial correlation is 1 in the diagonal. The matrix is extended to have the same ages in both stocks. Ages present in the data are represented by a black rectangle in the diagonal. 95% CIs are included in parentheses (see Appendix 3 for details). A black triangle in the upper right corner indicates that the (partial) correlation is significant. Axis labels indicate the age group, e.g. A3 is age three in the first stock, while B2 is age two in the second stock. For Baltic Herring, the recruitments were estimated to be highly partially correlated in M1 (Figure 3) with an estimated coefficient of 0.94 (0.58; 0.99). For model M1, the partial correlations are equal to the correlations by construction. As with the recruitment, ages two had a high correlation. This correlation was estimated to be highly negative (Estimate: −0.69) albeit with an uncertain estimate (95% CI: −0.99; 0.96). Ages 5 and 8 were both estimated to be positively correlated with coefficients of 0.74 and 0.55, respectively. However, both 95% CIs included zero (Figure 3). Likewise, the remaining estimated correlations had CIs which included zero. With respectively 8 and 9 ages in the two Baltic Herring stocks, eight partial correlations were estimated for M1 using the same number of free parameters in the Choleskey decomposition. Figure 3. View largeDownload slide Estimated partial correlation (upper triangle) and corresponding correlation (lower triangle) in non-fishery survival, Σ^N, for model M1 in the Baltic Herring case study between all stocks and ages. By construction, both correlation and partial correlation is 1 in the diagonal. The matrix is extended to have the same ages in both stocks. Ages present in the data are represented by a black rectangle in the diagonal. 95% CIs are included in parentheses (see Appendix 3 for details). A black triangle in the upper right corner indicates that the (partial) correlation is significant. Axis labels indicate the age group, e.g. A3 is age three in the first stock, while B2 is age two in the second stock. Figure 3. View largeDownload slide Estimated partial correlation (upper triangle) and corresponding correlation (lower triangle) in non-fishery survival, Σ^N, for model M1 in the Baltic Herring case study between all stocks and ages. By construction, both correlation and partial correlation is 1 in the diagonal. The matrix is extended to have the same ages in both stocks. Ages present in the data are represented by a black rectangle in the diagonal. 95% CIs are included in parentheses (see Appendix 3 for details). A black triangle in the upper right corner indicates that the (partial) correlation is significant. Axis labels indicate the age group, e.g. A3 is age three in the first stock, while B2 is age two in the second stock. In the Baltic Mix case (Figure 4), recruitments were estimated (95% CI) to have a negative partial correlation in M3 of −0.53 (−0.77; −0.07). The resulting correlation in recruitment was also highly negative with an estimated coefficient of −0.65 (−0.94; −0.03). Several age classes were estimated to be highly correlated. Baltic Spart age 3 and Western Baltic Herring age 2 had an estimated correlation of −1.00 (−1.00; −0.49); Sprat age 5 had an estimated correlation with Herring ages 3–5 of 0.91 (−0.85; 0.99), 0.99 (0.59; 1.00), and 0.93 (0.72; 1.00) respectively; and Sprat age 7 was estimated to be correlated with Herring ages 5 and 6 with coefficients of 0.93 (0.55; 0.99) and 0.99 (0.29; 1.00). With 8 and 9 ages respectively in the two Baltic Mix stocks, 36 partial correlations were estimated for M3 using the same number of free parameters in the Choleskey decomposition. Of these partial correlation estimates, 39% were estimated to be negative, 47% had an absolute value below 0.1, and 14% had an absolute value above 0.7. Figure 4. View largeDownload slide Estimated partial correlation (upper triangle) and corresponding correlation (lower triangle) in non-fishery survival, Σ^N, for model M3 in the Baltic Mix case study between all stocks and ages. By construction, both correlation and partial correlation is 1 in the diagonal. The matrix is extended to have the same ages in both stocks. Ages present in the data are represented by a black rectangle in the diagonal. 95% CIs are included in parentheses (see Appendix 3 for details). A black triangle in the upper right corner indicates that the (partial) correlation is significant. Axis labels indicate the age group, e.g. A3 is age three in the first stock, while B2 is age two in the second stock. Figure 4. View largeDownload slide Estimated partial correlation (upper triangle) and corresponding correlation (lower triangle) in non-fishery survival, Σ^N, for model M3 in the Baltic Mix case study between all stocks and ages. By construction, both correlation and partial correlation is 1 in the diagonal. The matrix is extended to have the same ages in both stocks. Ages present in the data are represented by a black rectangle in the diagonal. 95% CIs are included in parentheses (see Appendix 3 for details). A black triangle in the upper right corner indicates that the (partial) correlation is significant. Axis labels indicate the age group, e.g. A3 is age three in the first stock, while B2 is age two in the second stock. Simulation study The simulation study consisted of a total of 900 simulations from the six models M0 to M5. Of these simulations, 81% of the re-estimations identified the correct model from the candidates M0 to M5 (Table 3). Simulations from M3 were the most difficult to identify. Of these, 70% were correctly identified as generated by M3 whereas 21% were identified as M4, 5% as M5, and 3% were identified as M2. Only 7.1% of all simulations obtained the lowest AIC with a smaller model than the true simulation model. Table 3 Number of simulation studies that resulted in the lowest AIC for models M0 to M5.   Best model   True Model  M0  M1  M2  M3  M4  M5  M0  136  12  1  1  0  0  M1  0  140  6  1  2  1  M2  0  0  120  14  8  8  M3  0  0  5  105  32  8  M4  0  0  5  21  107  17  M5  0  0  0  0  33  117    Best model   True Model  M0  M1  M2  M3  M4  M5  M0  136  12  1  1  0  0  M1  0  140  6  1  2  1  M2  0  0  120  14  8  8  M3  0  0  5  105  32  8  M4  0  0  5  21  107  17  M5  0  0  0  0  33  117  Although the AIC was able to identify the true model in most cases, spurious correlations were found by the complex models when the true model did not include correlations. When simulating from M0, all the models M1 to M5 had high estimated correlations. For M5, 67% of the simulations had one or more estimated correlations with an absolute value above 0.5 while 10% had an absolute value above 0.75. The number was reduced with model complexity. The models M1 to M4, had estimated correlations with absolute values above 0.5 for 44, 57, 65, and 67% of the simulations, respectively. Likewise, respectively 3, 9, 11, and 10% of the simulations had estimated correlations with absolute values above 0.75. Besides identifying the true model, the simulations were used to calculate the coverage of Hessian based 95% CIs of the logarithm of the average fishing mortality (Table 4). When simulating from the models M0 and M1, all four estimation models had an average coverage close to 95% based on 150 simulations. This was not the case in the simulations where substantial correlations were included in the true model. When simulating from M2 to M5, estimation model M0 only had a coverage of 91.6, 89.7, 88.5, and 91.9%, respectively. Likewise, when simulating from M5, estimation models M1 to M3 had coverages of 91.7, 89.7, 90.2%, respectively, while models M4 and M5 had higher coverages: 93.4 and 93.2%, respectively. Table 4 Average (over time) coverage of Hessian based 95% CIs of the logarithm of the predicted average (over ages) fishing mortality in the simulation study.   Estimation model   Simulation Model  M0  M1  M2  M3  M4  M5  M0  0.945  0.943  0.941  0.940  0.940  0.940  M1  0.942  0.940  0.939  0.936  0.935  0.934  M2  0.916  0.933  0.939  0.936  0.937  0.936  M3  0.897  0.930  0.939  0.934  0.933  0.931  M4  0.885  0.921  0.929  0.928  0.929  0.931  M5  0.919  0.917  0.897  0.902  0.934  0.932    Estimation model   Simulation Model  M0  M1  M2  M3  M4  M5  M0  0.945  0.943  0.941  0.940  0.940  0.940  M1  0.942  0.940  0.939  0.936  0.935  0.934  M2  0.916  0.933  0.939  0.936  0.937  0.936  M3  0.897  0.930  0.939  0.934  0.933  0.931  M4  0.885  0.921  0.929  0.928  0.929  0.931  M5  0.919  0.917  0.897  0.902  0.934  0.932  Discussion European fisheries management mainly relies on single-stock assessment models, which assumes that stocks are not interconnected. For many stocks, this is not a reasonable assumption, as they are connected by several interactions, such as predation, migration and influences from the environment. Currently, the main alternative to single-stock models is complex ecosystem models, which are data intensive and require extensive tuning and long run times. We introduced a simple and fast alternative: connecting stocks by correlated survival variability. In our case studies, the models were fitted in 1.1–4.5 min. The models only require the data already used in single-stock models, which makes it readily usable for a wide range of stocks. Further, the model can be reduced to independent single-stock models, which allows objective testing of whether connecting stocks is an improvement, through likelihood ratio tests, AIC, or the like. We modelled complicated correlation patterns by a simple structure through the partial correlations; however, other correlation structures could be used instead. Unlike multi-species models and ecosystem models, which try to directly model interactions or common influences of stocks, the present approach estimates whether a good year for one cohort is also a good year for another cohort it may be linked to. Although simplistic and less informative about the ecosystem, this approach is highly operational and easily implemented for any group of stocks with an analytical assessment, since the model does not require any new data sources to be introduced. Further, the approach can be used to model complex correlation structures in other parts of the model such as the fishing mortality process. The case studies presented were chosen to illustrate stocks that are believed to migrate and stocks that share environment. In general, connecting stocks should be considered when: biological knowledge suggests interactions, the management areas intersect, or patterns in model validation tools such as residuals (e.g. Thygesen et al., 2017) or retrospective analysis (e.g. Hurtado-Ferro et al., 2015) are present. If biological knowledge indicates migration or predator-prey relations, it should be tested whether assessments can be improved by correlating the abundance processes. Likewise, when the management areas intersect, the stocks likely share environmental conditions or compete for resources which can correlate the abundance processes. Finally, common or opposite patterns in the residuals or retrospective analysis can indicate missing biological knowledge which can lead to correlated abundance processes. Although modelling partial correlations provides a simple interpretation of parameters, understanding why they arise is not as simple in the present framework. Correlations among stocks or ages can be the result of several processes such as migration, predation or environmental impact. Although further research is needed to evaluate the possibility of extending the framework to directly include any of these processes, it is straight forward to evaluate whether stocks are similar in the sense that they could share parameters in the assessment model. Testing whether parameters are shared among stocks can give a basis for evaluating if stocks should be managed together or separately (Montenegro et al., 2009). For all three case studies, the proposed model improved on the single-stock models as measured by AIC. In general, model M3 performed well for all stocks, although smaller models were preferred for two of the cases. The model allows a flexible correlation structure to be estimated with relatively few parameters. By modelling partial correlations with a band structure between stocks, we obtain very general correlation structures, even within the individual stocks. With this model, correlations within each stock are interpreted as a result of a connection between the stocks. Hence, the model only changes the single-stock model as an effect of the connection between stocks. In the simulation study we showed that the true model can be identified by AIC, although with a risk of choosing a model including too many partial correlations. It is well-known that AIC can favour overfitting, when the true model has a finite parameter space, and is among the candidate models (Burnham and Anderson, 2003). However, for stock assessment models the true data generating system is too complex to be described by a simple model; observed data will depend on, amongst other things, the hydrological systems, individual fish behaviour, individual fisherman behaviour, and a wide array of data accumulation procedures. Any practically applicable candidate model will require strong simplifications; therefore, the AIC does not necessarily select the true model in the case studies, but simply the model that has the best fit to the data, while penalizing for the additional parameters used. Although the selected model provides the best fit to data, interpretation of individual correlation parameters must be done with caution. Estimated correlations in the abundance process will be a result of several effects such as environmental conditions, predation and migration, which cannot be separated in the presented model. Further, if the correlation between e.g. observations are not adequately modelled, the state-space model can somewhat compensate to increase the model fit using the correlations in the abundance process (Thygesen et al., 2017). Finally, the typically short time series available for stock assessment models can result in imprecise individual correlation estimates, which clutters the overall pattern seen. Ignoring correlations in data result in unreliable CIs for model output. In a simulation study, we investigated the CIs of the estimated fishing mortality; both when correlations were present and absent. We found that when severe correlations are present in the data, (i.e. simulating from models M3 to M5), ignoring the correlations results in CIs that are too narrow. When no correlations were estimated, Hessian based 95% CIs only had a coverage of 89.7, 88.5, and 91.9%, respectively, when simulating from models M3 to M5, while the more flexible correlation structures provided CIs close to the correct coverage. Accurate assessments of model output uncertainty is important for sustainable fisheries management (Chaput, 2004; Cadrin and Dickey-Collas, 2014). Underestimating the uncertainty of model output used for management will increase the probability of overfishing (Punt et al., 2012). Hence, if correlations in data are ignored, a fishery perceived as sustainable may in fact have a non-sustainable risk of depletion. These results will allow managers and policymakers to make informed decisions based on a more accurate assessment of model output uncertainty. As a simple extension of single-species assessments, the model output can be evaluated per stock, while giving an accurate evaluation of uncertainty. Likewise, connecting stocks may increase the accuracy of forecasts. Together with the estimated correlations, these improvements can provide guidance on managing single stocks and the potential influence on surrounding stocks. The procedure outlined is directly applicable without additional data sources. Acknowledgements The authors thank handling editor Dr Shijie Zhou and two anonymous reviewers for their valuable comments to improve this article. Funding Funding for this study was provided by a grant from the European Fisheries and Maritime Fund (33113-B-15-003, Ministry of Environment and Food in Denmark). References Ainsworth C., Walters C. 2015. Ten common mistakes made in ecopath with ecosim modelling. Ecological Modelling , 308: 14– 17. Google Scholar CrossRef Search ADS   Akaike H. 1974. A new look at the statistical model identification. IEEE Transactions on Automatic Control , 19: 716– 723. Google Scholar CrossRef Search ADS   Albertsen C. M., Nielsen A., Thygesen U. H. 2017. Choosing the observational likelihood in state-space stock assessment models. Canadian Journal of Fisheries and Aquatic Sciences , 74: 779– 789. Google Scholar CrossRef Search ADS   Begley J., Howel D. 2004. An overview of gadget, the globally applicable area-disaggregated general ecosystem toolbox. Technical Report, ICES document CM/FF:13. Burnham K., Anderson D. 2003. Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach . Springer, New York. Cadrin S. X., Dickey-Collas M. 2014. Stock assessment methods for sustainable fisheries. ICES Journal of Marine Science , 72: 1– 6. Google Scholar CrossRef Search ADS   Chaput G. 2004. Considerations for using spawner reference levels for managing single- and mixed-stock fisheries of atlantic salmon. ICES Journal of Marine Science , 61: 1379– 1388. Google Scholar CrossRef Search ADS   Chen Z., Xu S., Qiu Y., Lin Z., Jia X. 2009. Modeling the effects of fishery management and marine protected areas on the beibu gulf using spatial ecosystem simulation. Fisheries Research , 100: 222– 229. Google Scholar CrossRef Search ADS   Christensen V., Pauly D. 1992. ECOPATH II — a software for balancing steady-state ecosystem models and calculating network characteristics. Ecological Modelling , 61: 169– 185. Google Scholar CrossRef Search ADS   Christensen V., Walters C. J. 2004. Ecopath with ecosim: methods, capabilities and limitations. Ecological Modelling , 172: 109– 139. Google Scholar CrossRef Search ADS   Collie J. S., Botsford L. W., Hastings A., Kaplan I. C., Largier J. L., Livingston P. A., Plagányi É., et al.   2016. Ecosystem models for fisheries management: finding the sweet spot. Fish and Fisheries , 17: 101– 125. Google Scholar CrossRef Search ADS   Daan N. 1987. Multispecies versus single-species assessment of north sea fish stocks. Canadian Journal of Fisheries and Aquatic Sciences , 44: s360– s370. doi:10.1139/f87-337. Google Scholar CrossRef Search ADS   Eero M., Hemmer-Hansen J., Hussy K. 2014. Implications of stock recovery for a neighbouring management unit: experience from the baltic cod. ICES Journal of Marine Science , 71: 1458– 1466. Google Scholar CrossRef Search ADS   Frøysa K. G., Bogstad B., Skagen D. W. 2002. Fleksibest—an age-length structured fish stock assessment model. Fisheries Research , 55: 87– 101. Google Scholar CrossRef Search ADS   Fulton E. A., Smith A. D. M., Punt A. E. 2004. Ecological indicators of the ecosystem effects of fishing: Final report. Report no. r99/1546. Technical Report, Australian Fisheries Management Authority, Canberra. Gislason H., Helgason T. 1985. Species interaction in assessment of fish stocks with special application to the north sea. Dana  , 5: 1– 44. Hannesson S., Jakobsdottir A., Begley J., Taylor L., Stefansson G. 2008. On the use of tagging data in statistical multispecies multi-area models of marine populations. ICES Journal of Marine Science , 65: 1762– 1772. Google Scholar CrossRef Search ADS   Hollowed A. 2000a. Are multispecies models an improvement on single-species models for measuring fishing impacts on marine ecosystems?. ICES Journal of Marine Science , 57: 707– 719. Google Scholar CrossRef Search ADS   Hollowed A. 2000b. Including predation mortality in stock assessments: a case study for gulf of alaska walleye pollock. ICES Journal of Marine Science , 57: 279– 293. Google Scholar CrossRef Search ADS   Holsman K. K., Ianelli J., Aydin K., Punt A. E., Moffitt E. A. 2016. A comparison of fisheries biological reference points estimated from temperature-specific multi-species and single-species climate-enhanced stock assessment models. Deep Sea Research Part II: Topical Studies in Oceanography , 134: 360– 378. Google Scholar CrossRef Search ADS   Hurtado-Ferro F., Szuwalski C. S., Valero J. L., Anderson S. C., Cunningham C. J., Johnson K. F., Licandeo R., et al.   2015. Looking in the rear-view mirror: bias and retrospective patterns in integrated, age-structured stock assessment models. ICES Journal of Marine Science , 72: 99– 110. Google Scholar CrossRef Search ADS   ICES 2014. Report of the baltic fisheries assessment working group (WGBFAS), 3-10 April 2014, ICES HQ, Copenhagen, Denmark. ICES CM 2014/ACOM:10, ICES. ICES 2015a. Report of the Baltic Fisheries Assessment Working Group (WGBFAS), 14-21 April 2014, ICES HQ, Copenhagen, Denmark. ICES CM 2015/ACOM:10, ICES. ICES 2015b. Report of the Herring Assessment Working Group for the Area South of 62°N (HAWG), 10-19 March 2015, ICES HQ, Copenhagen, Denmark. ICES CM 2015/ACOM:06, ICES. Kristensen K., Nielsen A., Berg C. W., Skaug H., Bell B. 2016. Tmb: automatic differentiation and laplace approximation. Journal of Statistical Software , 70: 1– 21. Google Scholar CrossRef Search ADS   Lewy P., Vinther M. 2004. A stochastic age-length-structured multispecies model applied to north sea stocks. Technical Report, ICES document CM/FF:20. Link J. S., Fulton E. A., Gamble R. J. 2010. The northeast US application of ATLANTIS: A full system model exploring marine ecosystem dynamics in a living marine resource management context. Progress in Oceanography , 87: 214– 234. Google Scholar CrossRef Search ADS   Methot R. D., Wetzel C. R. 2013. Stock synthesis: A biological and statistical framework for fish stock assessment and fishery management. Fisheries Research , 142: 86– 99. Google Scholar CrossRef Search ADS   Montenegro C., Maunder M. N., Zilleruelo M. 2009. Improving management advice through spatially explicit models and sharing information. Fisheries Research , 100: 191– 199. Google Scholar CrossRef Search ADS   Nielsen A., Berg C. 2014. Estimation of time-varying selectivity in stock assessments using state-space models. Fisheries Research , 158: 96– 101. Google Scholar CrossRef Search ADS   Punt A. E., Siddeek M. S. M., Garber-Yonts B., Dalton M., Rugolo L., Stram D., Turnock B. J., et al.   2012. Evaluating the impact of buffers to account for scientific uncertainty when setting TACs: application to red king crab in bristol bay, alaska. ICES Journal of Marine Science , 69: 624– 634. Google Scholar CrossRef Search ADS   Quinn T. J., Deriso R. B., Neal P. R. 1990. Migratory catch-age analysis. Canadian Journal of Fisheries and Aquatic Sciences , 47: 2315– 2327. Google Scholar CrossRef Search ADS   Richards R. A., Jacobson L. D. 2016. A simple predation pressure index for modeling changes in natural mortality: application to gulf of maine northern shrimp stock assessment. Fisheries Research , 179: 224– 236. Google Scholar CrossRef Search ADS   Senina I., Sibert J., Lehodey P. 2008. Parameter estimation for basin-scale ecosystem-linked population models of large pelagic predators: application to skipjack tuna. Progress in Oceanography , 78: 319– 335. Google Scholar CrossRef Search ADS   Tarnecki J. H., Wallace A. A., Simons J. D., Ainsworth C. H. 2016. Progression of a Gulf of Mexico food web supporting Atlantis ecosystem model development. Fisheries Research , 179: 237– 250. Google Scholar CrossRef Search ADS   Thygesen U. H., Albertsen C. M., Berg C. W., Kristensen K., Nielsen A. 2017. Validation of ecological state space models using the laplace approximation. Environmental and Ecological Statistics , 24: 317–339. Yodzis P. 1998. Local trophodynamics and the interaction of marine mammals and fisheries in the benguela ecosystem. Journal of Animal Ecology , 67: 635– 658. Google Scholar CrossRef Search ADS   Appendix 1: Unstructured covariance matrix Let L be a lower triangular matrix of parameters with diagonal 1, and normalized such that ⟨Li·,Li·⟩=1 for all i. Then Σ=LLT is an arbitrarily scaled covariance matrix. Off-diagonal elements of the covariance matrix can be forced to be zero (e.g. https://github.com/admb-project/examples/tree/master/admb-tricks/parameterization/covariance-matrices; accessed: 06 March 2017) by noting that   Σjk=⟨Lj·,Lk·⟩. Hence, by defining L˜ such that L˜j·=Lj·−⟨Lj·,Lk·⟩Lk· and L˜i·=Li· for i≠j, then for the covariance matrix Σ˜=L˜L˜T,   Σ˜jk=⟨L˜j·,L˜k·⟩=⟨Lj·−⟨Lj·,Lk·⟩Lk·,Lk·⟩=⟨Lj·,Lk·⟩−⟨Lj·,Lk·⟩⟨Lk·,Lk·⟩=⟨Lj·,Lk·⟩−⟨Lj·,Lk·⟩=0 Although this method allows us to impose constraints on the covariance matrix, it is, in practice, difficult to pre-determine which parameters will be redundant. Instead, the gradient and hessian of the resulting likelihood function can be examined to determine parameters that do not affect the likelihood value. These parameters must be fixed in optimization. Appendix 2: Simulation study details Stock structure We simulated 100 years from the model and discarded the first half of the time series. Two stocks were simulated with 5 age groups each. The last age group was considered a plus group. Natural mortality was fixed at 0.2. The initial log-population size was set to 16 for the first age group and was 0.5 less for each subsequent age group. Fishing mortality log ⁡-fishing mortality for stock s, age a, and year y was simulated by   log ⁡Fs,a,y= cos ⁡(y10+10·1{2}(s))−1+ϵs,a,y where   1{2}(s)={1 if s=20 if s≠2, and ϵs,a,y∼N(0,0.012). Log-population process The log ⁡-population given fishing mortality was simulated from the models M0 to M5 above. The off-diagonal non-zero entries, Pi,j, of the precision, P, matrix for Mα was constructed as −(α+1)/(α−|ai−aj|+1)2, where ai is the age corresponding to the ith index of the process. This construction ensures non-negligible partial correlations in the elements not included in the smaller models. Diagonal elements, Pi,i, were set to ∑j≠iPi,j+ exp ⁡(−α) to ensure the precision matrix was positive-definite. The corresponding covariance matrix was scaled to have diagonal elements 0.12. Observations log ⁡-observations given fishing mortality and log ⁡-population were simulated from a multivariate normal with AR(1) covariance structure. The AR(1) parameter was set to 0. For commercial data, the variances for all ages were 1. A first quarter survey was simulated from a multivariate normal with AR(1) covariance structure with parameter 0 and variances exp ⁡(−4). The catchability was set to exp ⁡(−7). Re-estimation Re-estimation was performed assuming an AR(1) structure for the correlations in the log ⁡-F random walk, random walk recruitment, and univariate log ⁡-normal observations. Appendix 3: CI for correlation parameters Simulation based CIs for correlation parameters were constructed through the profile likelihood. Splitting the parameter vector θ into correlation parameters θc and other parameters θo, we simulated proposal parameters θc(i), i=1,…,N, from the hessian based Gaussian approximation of the distribution of the estimator. The proposal was accepted if the hypothesis θc=θc(i) could not be rejected by a likelihood ratio test with significance level α. Each accepted simulated correlation parameter vector was transformed to a partial correlation and a correlation matrix. For each matrix element, minimum and maximum values were used as (1−α)·100 % CI. We simulated N=1000 proposal parameter vectors for each case study, and obtained acceptance rates of 0.80, 0.41, and 0.63 for Baltic Herring, Mix and Cod, respectively. For Baltic Mix and Cod, the variance of the proposal distribution was reduced to increase acceptance. © International Council for the Exploration of the Sea 2017. All rights reserved. For Permissions, please email: journals.permissions@oup.com

Journal

ICES Journal of Marine ScienceOxford University Press

Published: Jan 1, 2018

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off