Testing the stability of the suitability coefficients from an eastern Bering Sea multispecies virtual population analysis

Abstract

Suitability coefficients are important for the estimation of predation mortality in a multispecies virtual population analysis (MSVPA) and subsequent use in the multispecies forecasting model (MSFOR). Testing the assumption of the stability of the suitability coefficients is important in assessing the robustness of the predictions made with MSFOR. We used different statistical methods to partially test this assumption for the eastern Bering Sea MSVPA model with eight species, using stomach content data for the years 1985–1989. Comparison of the estimates from two different sets of stomach content data (set one with all data and set two mainly with data from 1985) suggested that the differences between the two types of estimates were much reduced when the number of predator stomachs sampled increased. In a second approach, we contrasted the residual variances of partial data sets with the results from the fit of the total data set. Results suggested a small increase (∼10.8%) in the variation of the suitability coefficients. Comparison of the means of the suitability coefficients associated with each predator species suggests that only 13 of the 50 possible pairwise contrasts were significantly different (α = 0.05). In general, results suggested that the predator preferences and prey vulnerabilities remained stable over the time period studied. Therefore, MSFOR could be considered as a tool to advise fisheries managers within a multispecies context.

Introduction

Some of the most important parameters estimated in a multispecies virtual population analysis (MSVPA) model are the suitability coefficients. They reflect predator preferences, vulnerability, and availability of prey, which is influenced by the spatial overlap of predators and prey. The suitability coefficients are estimated iteratively in the MSVPA model with the following equation:  

(1)

where U _p,a,i,j is the predator food composition (%) that is assumed equal to the observed stomach content, is the average stock size of the prey p of age a, and W _p,a is the average weight of an individual prey in the stomach of predator i of age j.

Suitability coefficients are an important input for the multispecies forecasting model (MSFOR) because they determine future predation interactions. Therefore, it is important that these coefficients reflect general predation patterns and not the special conditions of any particular year ( Anon., 1989 ). Rice et al . (1991) suggested that in the context of the development of management advice from trophic models, it may be correct (and often necessary) to assume that preferences and vulnerabilities remain constant. However, the prey and predator biomass levels, species age structure, and spatial overlap of predators and prey might have changed between the time when stomach content data were collected and when the suitability estimates were used for providing advice. Therefore, it is important to assess the degree of stability in the suitability coefficients to assure robust results in MSFOR. If parameters are stable, MSVPA models can include trophic interactions in the development of management advice ( Rice et al ., 1991 ). In contrast, if parameters vary from year to year, trophic MSVPA models may have little value for managers ( Rice et al ., 1991 ).

Changes in suitability coefficients based on different years of stomach content data can be explained by random variability in feeding from year to year ( Rice et al ., 1991 ). If there are systematic large changes in suitability coefficients, then there may be significant processes not taken into account by the MSVPA assumptions, such as the spatial overlap of predator/prey or the functional response type. Systematic changes in suitability coefficients could also arise from non-linearities such as prey-switching in the functional feeding response of predators ( Rice et al ., 1991 ). MSVPA commonly assumes a type II functional response, which does not include prey-switching.

Stability of the suitability coefficients in the North Sea was tested using 4 years of stomach content data with a complete set of data for all predators in 1981 and partial information for 3 additional years ( Anon., 1986 , 1989 ; Rice et al ., 1991 ). Those authors used several statistical methods to test the assumption of their stability. Results suggested that predator preferences and prey vulnerabilities were stable over the study period when the effects of changing prey abundance and age structure are accounted for ( Rice et al ., 1991 ). A more extensive “year of the stomach” survey in the North Sea during 1991 allowed a full test of the stability assumption. Analyses suggested considerable shifts in some suitability coefficient estimates depending on which year's stomach data were used. However, with a few exceptions, those shifts did not result in large changes in M2 ( Anon., 1994 ). Those authors found that significant improvements could be obtained in model fit if year effects were added to the log-linear model, called the kernel model ( Equation (2) ). The potential cause could have been prey-switching, but there was no conclusive evidence for this phenomenon ( Anon., 1994 ). Here, we examine the stability assumption for an MSVPA parameterized for the Bering Sea. This assumption can only be partially tested owing to incomplete predator stomach content information, but sufficient data do exist to perform some tests that provide information about the stability of the parameters.

Methods

The eastern Bering Sea MSVPA included eight species: walleye pollock ( Theragra chalcogramma ), Pacific cod ( Gadus macrocephalus ), Greenland turbot ( Reinhardtius hippoglossoides ), and yellowfin sole ( Pleuronectes asper ) as both predator and prey species, rock sole ( Lepidopsetta bilineata ) and Pacific herring ( Clupea pallasi ) as prey species, and arrowtooth flounder ( Atheresthes stomias ) and northern fur seal ( Callorhinus ursinus ) as external predators. Details on age ranges, residual mortality values, catch-at-age data, etc., can be found in Livingston and Jurado (2000) . Technical aspects of the stomach content data are found in Livingston et al . (1993) . Several data sets are required in the estimation of suitability coefficients. However, the predator stomach content data are the main source of information to test the stability assumption. The sampling of stomach contents from the eastern Bering Sea was not uniform across years and quarters for the different predator species ( Table 1 ).

Table 1

Stomach content data used as input to the MSVPA model for the eastern Bering Sea (PLK – walleye pollock, COD – Pacific cod, GTB – Greenland turbot, YFS – yellowfin sole, ATF – arrowtooth flounder, NFS – northern fur seal, ALL – stomach data available for all predator species).

Year	Quarter I	Quarter II	Quarter III	Quarter IV
1981	–	COD	–	–
1982	NFS	NFS, PLK	NFS	NFS
1983	–	–	–	–
1984	–	–	ATF, COD, GTB, YFS
1985	ALL	COD, GTB, NFS, PLK, YFS	ALL	ALL
1986	–	ATF	ATF, COD, GTB, YFS, PLK	ATF, PLK
1987	PLK, YFS	COD, PLK	ATF, COD, GTB, YFS, PLK	COD, PLK
1988	–	–	ATF, COD, GTB, YFS, PLK	PLK
1989	COD	–	ATF, COD, GTB, YFS, PLK	PLK, YFS
1990	–	ATF, PLK	–	–
1991	–	COD	–	–
1992	–	–	–	–
1993	PLK	–	–	–

Different sets of suitability coefficients were estimated from several MSVPA runs using different stomach content data sets. Two approaches were taken. First, we compared two sets of suitability coefficients estimates from two sets of stomach content data, BSALL (included all stomach contents information) and the data subset BS85 (information from 1985; except for the second quarter for arrowtooth flounder, for which data from the second quarter of 1986 were used). Data from 1985, the only complete year of stomach data, were included in both data sets to make it possible to run the MSVPA. A direct comparison of the two sets of suitability estimates and predation mortalities corresponding to 1985 was made using weighted regression, with the number of stomachs sampled as weight and the estimates from BS85 as response variable. We also compared partial predation mortalities and stock size estimates obtained from both data sets.

In a second approach, we applied the methods used in the North Sea ( Anon., 1989 ; Rice et al ., 1991 ). We built six partial subsets of stomach content data with fixed diet information for quarters I, II, and IV from 1985, and replaced the information for the third quarter for all predators and for the fourth quarter for walleye pollock with information from 1985 to 1989, defining the data sets: BS85, BS86, BS87, BS88, and BS89. Subset BS8589 contained all information from the partial data sets. Although data sets BS85–BS89 differ slightly, they vary in the information corresponding to the third quarter, the main period of predator feeding in the Bering Sea. We ran MSVPA with each data set to obtain a set of suitability coefficients. Most analyses were done using a log-linear model, called the kernel model ( Rice et al ., 1991 ):  

(2)

where S _p,a,i,j is the suitability coefficient. Most of the predictor variables were defined as factors with i representing the predator, p the prey, and q the quarter of the year. The interaction between factors was represented by an “∗” and “ratio” is the ratio of the average weight-at-age of the predator to the average weight-at-age of the prey. The logarithm of the weight ratio and its squared term were defined as covariates. The logarithm of the weight ratio is nested under each predator and prey. Finally, ɛ is a normally distributed error term. Several versions of the kernel model have been used ( Sparre, 1991 ), but the variants never resulted in large changes in the amount of the variability explained. Thus, in this analysis, we only use the model ( Equation (2) ) for the analyses.

Results from the analysis of deviance from the kernel model fit to sets of suitability coefficients (BS8589 and partial data sets) were used to analyse changes in residual variance attributable to the use of different stomach data sets. In addition, to explore the magnitude of change in the suitability coefficients among years, we compared the estimated model coefficients for each predictor variable in the fits of model ( Equation (2) ) to the partial sets ( Anon., 1988 ).

To analyse the impact of the variation in stomach contents in the third quarter (fourth quarter for walleye pollock) among years on estimated suitability coefficients, we adopted the same method as Rice et al . (1991) . The method consists of fitting the kernel model to the data set BS8589 that was restructured (BS8589REST), ensuring that stomach data from each year were treated as coming from a different level of the factor predator. We analyse the changes in variability explained by the model and its predictor variables, as well as the temporal variation of the estimated coefficients associated with each predator factor. We also compared the fitted means of the predator levels for the kernel model fit to the BS8589REST data set. Means of the same predator in different years were contrasted using Tukey's test for multiple comparisons (α = 0.05).

To explore potential prey-switching related to changes in spatial overlap of predators and prey, we regressed the residuals from the kernel model fit to the BS8589 data set against the log-transformed predator and prey biomasses estimated in the MSVPA ( Rice et al ., 1991 ):  

(3)

where B _i,j represents the total biomass of the predator i (age j), and B _p,a the biomass of the prey p (age a). For each predator–prey combination, we fitted separate regression models ( Equation (3) ). We repeated the same process for the residuals from the model fits to the partial data sets (BS85 to BS89).

In order to see whether differences in suitability coefficients obtained from different stomach data subsets (BS85 to BS89) could be due to differences in underlying biomass estimates in the different years, we regressed the difference between sets of the estimated suitability coefficients against the corresponding differences in MSVPA biomasses estimates of predator i and prey p ( Equation (4) ).  

(4)

Finally, we explored the relationships between the suitability coefficient estimates of adult and age-1 walleye pollock, a separation index ( Wespestad et al ., 2000 ), and the eastern Bering Sea bottom trawl index of the age-1 walleye pollock survey ( Wespestad, 1994 ). The separation index is used to compare the distributions of age-1 walleye pollock (length <20 cm) and adult walleye pollock (length >40 cm), and it is defined ( Wespestad et al ., 2000 ) as:  

(5)

where A _i and J _i are the abundance (kg ha ⁻¹ ) of adult and juvenile walleye pollock, respectively. The separation index is zero if the relative distributions of juvenile (age-1) and adult pollock are exactly the same; the index increases when the difference between their distributions increases. More details are found in Wespestad et al . (2000) . We separated the suitability estimates into two groups according to predator age. The first group included adult walleye pollock predators from ages 5 to 7 and the second group from ages 8 to 12+.

Results

Direct comparison of suitability coefficients, predation mortality, and stock estimates

The scatterplot of estimates of suitability coefficients in 1985 from both data sets (BSALL and BS85) shows several outliers ( Figure 1 ) associated with small sample sizes. The weighted regression (estimates from BS85 as response variable) for both sets of suitability coefficients estimates in 1985 was significant (p-value ∼ 0.0) and explained 53.63% of the variability observed. The slope was 0.95 ± 0.02, and the intercept was 0.01 ± 0.002. However, the slope was significantly different from one (H ₀ : β = 1, p-value = 0.01). The fitted residuals from the regression decreased abruptly when sample size increased ( Figure 2 ). Of 1969 fitted residuals, 1597 were >0.01 in absolute value, but only 353 of these had a sample size >200 stomachs.

Figure 1

Comparison of estimates of 1985 suitability coefficients for all predators included derived from two subsets of stomach data (BS85 and BSALL); suit85b – suitability coefficients estimated with fewer than 200 stomachs, suit85 – suitability coefficients estimated with more than 200 stomachs.

Figure 2

Fitted residuals (BSALL vs. BS85) against number of stomachs available for each suitability value (estimates from BSALL data set).

A more detailed analysis with regression models, contrasting estimates from particular predator–prey combinations at a time, showed that, for half the predators, the suitability coefficients of walleye pollock as prey were similar (H ₀ : β = 1; walleye pollock, p-value = 0.64; yellowfin sole, p-value = 0.21; northern fur seal, p-value = 0.94).

We also fitted weighted regression models by predator species. Although most of the models were significant, only three regression models had a slope not significantly different from one (walleye pollock, p-value = 0.40; yellowfin sole, p-value = 0.21; northern fur seal, p-value = 0.98). Only two regression models (Pacific cod and arrowtooth flounder) explained more than 50% of the variability observed.

The weighted regression of estimates of partial M2 (p-value ≈ 0.0), with estimates from the BS8589 data set as response variable, explained 52.72% of the observed variability with a slope significantly bigger than one (H ₀ : β ≤ 1, p-value = 0.0), suggesting that the estimates from the data set BS85 were larger than those from the BSALL data set. Regressions were fitted to predation mortality estimates by predator. Most were significant (p-value ≈ 0.0), except for Greenland turbot, for which there was no relationship between BSALL and BS85 estimates. Half the estimated slopes were >1.0 (walleye pollock, p-value = 0,0; Pacific cod, p-value = 0.0; yellowfin sole, p-value = 0.0). For arrowtooth flounder, the slope was smaller than one (H ₀ : β ≥ 1, p-value = 0.0). The slope corresponding to northern fur seal was one (H ₀ : β = 1, p-value = 0.27).

Finally, in general, the differences between the estimates of stock size from the BS85 and BSALL data sets were <10%, where most of the estimates of stock size from the BS85 data set were larger (not shown). Only walleye pollock and Greenland turbot had the opposite effect. Pacific herring was the species with the greatest deviations between the two estimates.

Analyses with the kernel model

In the second approach, the kernel model was fitted to the BS8589 data set. It explained 84% of the observed variability. The factor for quarter of the year explained the majority of the variability (32.3%). The predator factor and the quarter–predator interaction explained 18.7% and 19.5%, respectively, and the quarter–prey interaction, 8.4%.

The model fits to partial data sets (BS85–BS89) consistently explained around 78% of the variability observed, except for the fit to the BS88 data set ( Table 2 ), for which the increase in variability was due to the amount of variability explained by the quarter-of-the-year factor. Once again predator, quarter of the year, and their interaction explained the majority of the variability in the fits.

Table 2

Comparison of the analysis of variance results of the log-linear models for suitability coefficients fitted separately to the partial stomach data sets. SSQ – sum of squares, d.f. – degrees of freedom, r ² – percentage of variability explained by the model, lnratio[pred] – logarithm of the body weight ratio nested under predator, lnratio[prey] – logarithm of the weight ratio nested under prey, lnratio2 – logarithm of the body weight ratio squared.

	BS85		BS86		BS87		BS88		BS89
	SSQ	d.f.	SSQ	d.f.	SSQ	d.f.	SSQ	d.f.	SSQ	d.f.
Quarter	12 556	3	1 330	3	1 023	3	1 704	3	1 416	3
Pred	1 422	5	1 385	5	1 558	5	1 245	5	1 420	5
Prey	289	4	246	4	354	4	322	4	274	4
lnratio2	143	1	69	1	134	1	43	1	46	1
quarter:pred	1 318	10	1 538	10	1 885	10	1 436	10	1 399	10
quarter:prey	620	10	692	10	613	10	703	10	855	10
pred:prey	52	3	56	3	47	3	53	3	56	3
lnratio[pred]	16	5	20	5	29	5	26	5	50	5
lnratio[prey]	63	5	80	5	57	5	53	5	93	5
Residual	1 387	1 590	1 577	1 590	1 684	1 590	1 331	1 590	1 470	1 585
Total	6 566	1 636	6 993	1 636	7 383	1 636	6 916	1 636	7 078	1 631
r 2	0.79	–	0.77	–	0.77	–	0.81	–	0.79	–

The residual sum of squares of individual data sets (BS85–BS89) contains the unexplained within-year variability. If there was no between-year variability in suitability estimates, the residual sum of squares from the BS8589 data set should be approximately equal to the sum of the residual sums of squares of the partial data sets. Any additional variability could then be attributed to changes in suitability coefficients. The residual sum of squares of the partial data sets was 7448.02, whereas the residual sum of squares of the fit to the combined data set (BS8589) was 8352.19. Therefore, the variance part attributable to between-year variability was 10.8% ( Table 2 ). The fit of the kernel model to the restructured data (BS8589REST) explains a little more variability (3%) than the model fitted to the BS8589 data set ( Table 3 ). The percentage of variation explained by the factor predator increased from 19.6% in the fit to the BS8589 set to 27.8% in the model with separate factor levels for each predator and year combination. The addition of extra levels (5–25 degrees of freedom) in the predator factor decreased the amount of variability explained by the interaction quarter–predator ( Table 3 ).

Table 3

Analysis of deviance of the fit to the BS8589built restructured data set; d.f. – degrees of freedom, r ² – percentage of variability explained by the model, % – percentage of variability explained by the predictor variable, lnratio[pred] – logarithm of the body weight ratio nested under predator, lnratio[prey] – logarithm of the weight ratio nested under prey, lnratio2 – logarithm of the body weight ratio squared.

	d.f.	Deviance	Residual d.f.	Residual deviance	Pr(F)	%
NULL	8 179	34 972	–	–	–	–
Quarter	3	6 597	8 176	28 375	0.00	18.86
Predator	25	9 706	8 151	18 669	0.00	27.75
Prey	4	1 808	8 147	16 861	0.00	5.17
lnratio2	1	337	8 146	16 524	0.00	0.96
quarter:pred	14	4 776	8 132	11 748	0.00	13.66
quarter:prey	10	3 275	8 122	8 473	0.00	9.36
pred:prey	11	490	8 111	7 983	0.00	1.40
lnratio[prey]	5	122	8 106	7 861	0.00	0.35
lnratio[pred]	25	682	8 081	7 179	0.00	1.95
r 2	0.79	–	–	–	–	–

The coefficients of the levels of the factor predator of the model fitted to the BS8589REST data set are shown in Table 4 . The coefficients of walleye pollock from 1985 to 1989 were different, but there is consistency in that most were negative and the confidence intervals overlapped. For Pacific cod as predator, coefficients were negative, with less variation than those for walleye pollock. The rest of the species displayed even less variation. The year-specific coefficients for Greenland turbot were similar and negative ( Table 4 ). The coefficients of yellowfin sole as predator were generally negative and differed between years, but their confidence intervals overlapped. The coefficients for arrowtooth flounder were the most consistent. In summary, the majority of year-specific predator coefficients remained fairly consistent with overlapping in the confidence intervals.

Table 4

Coefficients for the levels of the factor predator as a main effect obtained from the fit of the log-linear model to the BS8589REST data set. s.e. – standard error.

Level of predictor variable	Parameter estimate	s.e.
Walleye pollock 1985	−1.10	0.80
Walleye pollock 1986	−0.41	0.81
Walleye pollock 1987	0.08	0.87
Walleye pollock 1988	−1.24	0.82
Walleye pollock 1989	−2.54	0.75
Pacific cod 1985	−2.69	0.36
Pacific cod 1986	−2.30	0.36
Pacific cod 1987	−3.23	0.36
Pacific cod 1988	−2.23	0.36
Pacific cod 1989	−2.00	0.35
Greenland turbot 1985	−4.27	0.40
Greenland turbot 1986	−4.11	0.41
Greenland turbot 1987	−4.22	0.41
Greenland turbot 1988	−3.22	0.41
Greenland turbot 1989	−3.34	0.41
Yellowfin sole 1985	−0.08	3.31
Yellowfin sole 1986	−4.51	1.74
Yellowfin sole 1987	−1.49	2.10
Yellowfin sole 1988	0.61	3.12
Yellowfin sole 1989	−0.95	2.83
Arrowtooth flounder 1985	0.38	0.40
Arrowtooth flounder 1986	0.57	0.40
Arrowtooth flounder 1987	0.57	0.40
Arrowtooth flounder 1988	0.49	0.39
Arrowtooth flounder 1989	0.51	0.37

Owing to the unbalanced design in the sampling of the stomach data, only the coefficients of the combination walleye pollock as prey and Pacific cod as predator ( Table 5 ) were estimated. For walleye pollock, the coefficients were consistent for the majority of predators and years, except for Greenland turbot in 1988 and 1989 ( Table 5 ). For Pacific cod, the estimates of the interactions Pacific cod–Pacific cod and Pacific cod–yellowfin sole were similar, i.e. positive with overlapping confidence intervals. The estimates of Pacific cod–rock sole were also consistent, i.e. generally negative with absolute values smaller than 0.2 and confidence intervals overlapping. The slopes of the tendencies of size preferences (natural logarithm of the weight ratio nested under predator) were negative with absolute value smaller than one. The majority of the confidence intervals overlapped.

Table 5

Coefficient estimates for the predator–prey interactions from the fit of the log-linear model to the BS8589REST data set. * denotes the species as prey.

	Species-year	Slope estimate	s.e.
Walleye pollock as prey	plk 85	−2.01	0.50
	plk 86	−2.51	0.51
	plk 87	−2.57	0.54
	plk 88	−1.97	0.52
	plk 89	−1.60	0.48
Walleye pollock as prey	cod 85	−0.36	0.12
	cod 86	−0.31	0.12
	cod 87	−0.37	0.12
	cod 88	−0.28	0.12
	cod 89	−0.21	0.11
Walleye pollock as prey	gtb 85	0.08	0.13
	gtb 86	0.00	0.13
	gtb 87	−0.02	0.13
	gtb 88	−0.33	0.13
	gtb 89	−0.28	0.13
Walleye pollock as prey	atf 85	0.10	0.12
	atf 86	0.08	0.12
	atf 87	0.07	0.12
	atf 88	0.03	0.12
Pacific cod as predator	cod 85*	0.06	0.16
	cod 86*	0.09	0.16
	cod 87*	0.05	0.16
	cod 88*	0.01	0.16
Pacific cod as predator	yfs 85*	0.20	0.08
	yfs 86*	0.16	0.08
	yfs 87*	0.09	0.08
	yfs 88*	0.06	0.08
Pacific cod as predator	sole 85*	−0.02	0.11
	sole 86*	−0.14	0.11
	sole 87*	0.01	0.11
	sole 88*	−0.02	0.11

Also, as a consequence of the unbalanced design of the stomach sampling, only the fitted means of the log-transformed suitability coefficients of the factor predator were estimable. The multiple comparisons Tukey's test suggested that 37 of the 50 pairwise contrasts were not significantly different (α = 0.05). In particular, the walleye pollock means for the BS85–BS88 data sets were not significantly different. Also, the mean from the BS89 data set was not significantly different from the mean from the BS86 data set. For Pacific cod, only the mean from the BS87 data set differed significantly from the rest of the means. In the case of Greenland turbot the means did not differ from each other (α = 0.05). The yellowfin sole means from the 1985, 1987, and 1989 data sets did not differ among themselves. No significant differences were found (α = 0.05) for arrowtooth flounder. In summary, 74% of the pairwise comparisons from different years of stomach data were not significant, implying that the majority of the estimates of suitability coefficients from these individual years of stomach data were not significantly different.

We found that 50% of the regressions of the residuals of the kernel model (BS8589 data set) to predator and prey biomass were statistically significant ( Equation (3) ). However, predator and prey biomass explained little residual variability, with the highest r ² = 0.199 ( Table 6 ). This suggests that neither the biomass of predator nor the biomass of the prey can explain the differences between the estimates from the MSVPA and the kernel model.

Table 6

Results for the regression of the log-linear model residuals (BS8589 data set) against the logged predator and prey biomass.

Species	p-value	r 2	Predator biomass	p-value	Prey biomass	p-value
Predator
Walleye pollock	0.246	0.006	0.086	0.141	−0.021	0.420
Pacific cod	0.000	0.008	0.065	0.000	−0.005	0.661
Greenland turbot	0.000	0.021	−0.174	0.000	0.026	0.210
Yellowfin sole	0.071	0.029	−0.006	0.939	−0.027	0.022
Prey
Walleye pollock	0.102	0.002	−0.028	0.051	0.027	0.222
Pacific cod	0.020	0.050	0.201	0.007	−0.001	0.989
Greenland turbot	0.616	0.045	−0.205	0.536	−0.253	0.449
Yellowfin sole	0.000	0.023	0.064	0.005	0.106	0.001
Rock sole	0.047	0.011	0.098	0.029	−0.026	0.282
Pacific herring	0.000	0.199	0.116	0.001	−0.718	0.000

In many cases, differences in predator and prey biomass explained observed differences in suitability coefficients obtained from the partial data sets ( Table 7 ). For walleye pollock as predator, the majority of the slopes of the change in prey biomass were positive. When selected as prey, the majority of the slope estimates were negative. Thus, for walleye pollock as predator, prey suitability estimates were higher when prey biomass increased.

Table 7

Results of the regression of the pairwise difference in suitability coefficients of walleye pollock against differences in predator and prey biomass.

Years	p-value	r 2	Predator biomass difference slope	p-value	Prey biomass difference slope	p-value
Walleye pollock as predator
85–86	0.025	0.22	−5.82E−08	0.084	2.85E−08	0.016
85–87	6.69E−07	0.71	3.21E−08	0.073	2.45E−08	0.000
85–88	0.165	0.21	−4.08E−08	0.453	−1.87E−08	0.227
85–89	0.000	0.35	−9.18E−09	0.592	3.66E−08	0.000
86–87	0.738	0.02	−1.98E−08	0.600	−2.00E−08	0.567
86–88	4.95E−07	0.86	2.01E−09	0.890	−1.15E−07	0.000
86–89	0.111	0.09	−3.33E−08	0.174	2.63E−08	0.138
87–88	5.84E−10	0.99	−2.46E−09	0.812	7.45E−07	0.000
87–89	0.005	0.35	−5.14E−08	0.203	3.14E−08	0.005
88–89	1.07E−06	0.84	−2.36E−10	0.989	−6.01E−08	0.000
Walleye pollock as prey
85–86	0.059	0.02	−3.57E−07	0.313	1.21E−07	0.03
85–87	0.060	0.02	5.49E−07	0.429	−1.26E−07	0.02
85–88	7.40E−11	0.16	3.95E−07	0.326	−2.34E−07	0.00
85–89	1.08E−06	0.09	−2.79E−07	0.258	−3.68E−07	0.00
86–87	0.040	0.03	−2.60E−08	0.956	3.14E−07	0.01
86–88	0.642	0.00	5.18E−08	0.904	8.44E−08	0.35
86–89	0.001	0.06	2.98E−07	0.235	−3.40E−07	0.00
87–88	0.084	0.02	−9.14E−08	0.877	−3.20E−07	0.03
87–89	0.147	0.02	3.74E−07	0.296	−1.30E−07	0.08
88–89	0.571	0.00	2.26E−07	0.304	5.22E−09	0.92

There was no relationship between the suitability coefficients with the separation index for the first group (ages 5–7). Results were similar for the comparison of the age-1 survey and the suitability coefficients. For the second group (ages 8–12+), the relationship between the suitability coefficients and the separation index was not significant, but their correlation was 0.56 ( Figure 3 ). No relationship was found between the age-1 walleye pollock survey and the corresponding suitability coefficients. Therefore, neither the separation index nor the walleye pollock survey explained the observed variation in the suitability coefficients.

Figure 3

Comparison of suitability coefficients, ages 5–7 walleye pollock as predator and age-1 pollock as prey, against the separation index.

Discussion

The recognition of the importance of species interactions in the dynamics of fish populations is a key step towards a holistic approach to fisheries management. MSVPA and MSFOR incorporate predation interactions, and assume that the Andersen and Ursin (1977) formulation captures components such as prey preference, season, and size preference ( Rice et al ., 1991 ). MSFOR assumes stability of the suitability coefficients. If this holds, MSFOR using a reliable future recruitment assumption (stock-recruitment model for each species) and potential correlations between species recruitment could be considered as a consistent model to explore the indirect effects of fishing on the dynamics of species attributable to trophic links. Therefore, MSVPA and MSFOR may be considered a first step in the development of tools that provide advice to fisheries managers in a multispecies context.

An analysis of the stability of the suitability coefficients requires at least two sets of independent data ( Sparre, 1991 ). However, two complete sets of stomach content data for the Bering Sea do not exist. Therefore, the hypothesis of the stability of the suitability coefficients was only partially tested. A complete test of this assumption would require comparison of the mean suitability coefficient for each combination of predator–prey–age–quarter. This task is difficult because we cannot individually analyse all combinations produced by MSVPA. Instead, we used indirect methods such as the direct comparison of estimates ( Anon., 1994 ), and the methods based on the kernel model ( Sparre, 1991 ).

The ICES multispecies working group had two complete sets of stomach content data available for the North Sea (1981 and 1991). They compared the two corresponding sets of suitability coefficients, and also found that although there were considerable changes in the suitability coefficients of prey species, the changes did not produce substantial changes in the partial predation mortality estimated with the two runs ( Anon., 1994 ).

Although the two sets (BSALL and BS85) of walleye pollock suitability coefficients for 1985 were comparable, the corresponding predation mortality had small changes producing differences in the estimates of the stock size during the earlier years of the stock assessment. These differences could arise as a result of differences between suitability coefficients for a particular predator–prey combination in a given quarter. In particular, feeding takes place mainly in the third and fourth quarters, so small changes in suitability coefficients could influence the estimation of predation mortality and consequently the stock size estimates. The majority of the species had small changes in the suitability coefficients, with no visible effect on estimates of predation mortality and stock size, with the exception of Pacific herring. Changes in Pacific herring suitability coefficients as prey suggest that this species was more suitable as a prey for Pacific cod, yellowfin sole, and arrowtooth flounder when the BS85 data set was used. Therefore, despite some discrepancies in the estimates of the stock size of some species, overall the majority of stock size estimates were similar.

The results suggest that the variability of the deviations of suitability coefficients was reduced when the number of stomachs sampled increased, a result also found for the North Sea (A. Kempf, Institut für Hydrobiologie und Fischereiwissenschaften, Hamburg, pers. comm.). Therefore, it is recommended that stomach data from small sample sizes be excluded in future MSVPA runs.

Results from the kernel model showed a small interannual overall variability of the suitability coefficients (10.8%). This variability included more than systematic changes in preferences and vulnerabilities, but also sampling error over years ( Sparre, 1991 ). However, this is a conservative estimate, because all partial data sets contained the same stomach data from 1985 for three, but one quarter. Had the partial data sets contained no common data, which de facto introduces between-year variance estimates, the residual sum of squares of individual partial data sets might have been smaller and hence the between-year variance part >10.8%. This result agrees with previous analyses for the North Sea, in which a comparable interannual variation was found ( Anon., 1989 , 1994 ; Sparre, 1991 ).

Results from analyses with the BS8589REST data and the kernel model showed a small increase in the variability explained, suggesting a small temporal variability. A comparison of the estimated coefficients of the factor predator and the predator–prey interactions showed remarkable consistency. The parameters that displayed the greatest variability were associated with walleye pollock. However, despite its complex dynamics as a result of its important role as both predator and prey, its parameters were still similar with regard to the sign of the coefficients.

The pairwise comparison of the fitted means of logged suitability coefficients showed more similarities than differences, suggesting that the suitability coefficients remain consistent over time. However, we assumed that the suitability coefficients of each partial data set are independent. It is difficult to assess the validity of this hypothesis, because stomach information from 1985 (quarters I, II, and IV) was used to estimate the suitability coefficients of other years. Therefore, although the variance of various factors may have been underestimated, the degrees of freedom are overestimated. It is difficult to evaluate whether these effects cancel each other out ( Anon., 1988 ).

The overall results from the kernel model revealed some interannual changes in the suitability coefficients. However, there are fewer differences than similarities, suggesting that the suitability coefficients are stable over time, perhaps reflecting general aspects of feeding patterns of predators and not the special conditions of any particular year ( Anon., 1989 ). Similar results have also been found in at least two separate analyses for the North Sea ( Anon., 1988 , 1994 ; Sparre, 1991 ). Results suggest that MSVPA correctly represents the predation interactions. With adequate data to parameterize the model, the suitability coefficients remain stable long enough to be useful for fisheries management ( Sparre, 1991 ). Recent preliminary simulation results (A. Kempf, pers. comm.) suggest that MSFOR would produce important changes in the long-term estimates of stock biomass when different sets of stomach content data are used. The present work focused on the stock assessment stage, so no simulations were done. Therefore, we did not explore the implications of different sets of suitability coefficients in the long-term estimates of stock size. Such sets of simulations need to be carried out for the Bering Sea.

The interannual variability of the suitability coefficients may be related to spatial overlap and changes in the stock biomass of the prey populations ( Sparre, 1991 ; Anon., 1994 ). Results of the analyses for the Bering Sea confirm that the changes in suitability coefficients could be related to changes in the prey biomass, and that this factor may be included in the kernel model as a predictor variable.

Similar results caused the ICES multispecies working group to explore the possibility of prey-switching by modifying the kernel model. However, they did not find strong evidence of prey-switching in the North Sea. They suggested that general trends in prey biomass may be confounded with the year effect, and so act as a proxy rather than as a cause of the year effect ( Anon., 1994 ). Larsen and Gislason (1992) found that including prey-switching did not significantly improve the correlation between MSVPA estimates of year-class strength and independent stock estimates for the North Sea. Their results suggested that the hypothesis of constant suitability coefficients could not be rejected, but that if prey-switching occurs, it would be negative prey-switching, implying that the suitability to a prey species declines when its abundance increases. Other authors ( Rindorf et al ., 1998 ) have even suggested that the MSVPA suitability model is not appropriate to describe a predator's choice of prey. Therefore, more analyses are needed to test the stability assumption; however, it is necessary to make them in a new statistical framework, allowing the use of standard tools such as maximum likelihood estimation and likelihood ratios for model comparison.

Finally, there was no clear relationship between the suitability coefficients and the separation index. Similarly, no relationship was found when compared against the age-1 pollock survey. These comparisons included only five points (years 1985–1989) owing to the limited availability of stomach content data. These data do not allow conclusions to be drawn about the relationships explored. This aspect of the analysis will be improved when additional stomach content data are added to the Bering Sea MSVPA model.

All methods used in these analyses suggest that the suitability coefficients remain fairly constant over time. The methods analyse the change of suitability coefficients in an indirect way from a global perspective. It is difficult to track the temporal changes in the individual suitability coefficients for each predator–prey–age–quarter combination. A quick review of the temporal trend of some of these combinations would reveal temporal changes in suitability coefficients. It seems that, in an overall context, those changes could cancel each other out in such a way that the global changes are minimized.

It is difficult to imagine that in a system as dynamic as the Bering Sea, the suitability coefficients could remain constant over time. Most likely, changes in walleye pollock biomass and in the spatial overlap of adult and juvenile populations ( Wespestad et al ., 2000 ), together with environmental change (e.g. climate regime shifts), and population fluctuations of several species observed in the last few decades ( NRC, 1996 ), could potentially affect the relative constancy of the suitability coefficients. The suitability coefficients calculated in the MSVPA may contain a spatial overlap component that cannot be eliminated within the model ( Anon., 1989 ). Although the overlapping term did not explain a large amount of the variability in the models of the North Sea ( Anon., 1988 ), in the Bering Sea there is a suggestion that the spatial overlap of adult and juvenile populations determines the intensity of walleye pollock cannibalism, an important controlling element of interannual recruitment variability ( Wespestad et al ., 2000 ).

As spatial overlap depends on climate variability in the Bering Sea ( Wespestad et al ., 2000 ), it is necessary to improve our stomach content database to include spatial components for future analysis. The random component of the climate variability influencing the spatial overlap of juvenile and adult populations and the intensity of cannibalism suggests that representing the suitability coefficients with a statistical distribution would be an important improvement. However, with current MSVPA technology, estimation of the suitability coefficients distributions is not possible. Therefore, it is necessary to introduce statistical assumptions into a multispecies framework. This task requires a new approach to assess the dynamics of populations. MSVPA is based on Gulland's method (1965) , which is a backward non-linear sequential method that estimates stock size based mainly on catch-at-age data, diet files, and predator ration. The model is deterministic and does not provide a measure of how well the parameters were estimated ( Megrey, 1989 ). Therefore, we have developed a multispecies statistical model (MSM) that includes the separable fishing mortality assumption ( Doubleday, 1976 ) in a multispecies scenario. In this way, the common procedures of minimizing the sum of squares and/or maximizing the likelihood can be used. This approach would allow estimation of the posterior distribution of suitability coefficients and other model parameters, a task not possible with the current MSVPA model. This approach is potentially an important advance in providing advice to fisheries managers, because it allows consideration of uncertainty in a multispecies modelling context and will help to establish useful scenarios for the management of groundfish resources in the Bering Sea.

We thank Ivonne Ortiz (University of Washington) and Kerim Aydim (AFSC) for reviewing a draft of this manuscript, and two anonymous reviewers whose comments helped to improve the final product. We also thank the AFSC Food Habits laboratory for their work in the analysis of an immense number of stomachs that provided the data for carrying out this research. The publication was partially funded by the Joint Institute for the Study of the Atmosphere and Ocean (JISAO) under National Oceanic and Atmospheric Administration cooperative agreement number NA67RJ01555, contribution 1070.