Application of a generalized additive model (GAM) to reveal relationships between environmental factors and distributions of pelagic fish and krill: a case study in Sendai Bay, Japan

Abstract

Murase, H., Nagashima, H., Yonezaki, S., Matsukura, R., and Kitakado, T. 2009. Application of a generalized additive model (GAM) to reveal relationships between environmental factors and distributions of pelagic fish and krill: a case study in Sendai Bay, Japan. – ICES Journal of Marine Science, 66: 1417–1424.

A generalized additive model (GAM) was applied to fishery-survey data to reveal the influences of environmental factors on the distribution patterns of Japanese anchovy (Engraulis japonicus), sand lance (Ammodytes personatus), and krill (Euphausia pacifica). Echosounder and physical-oceanographic data were collected in Sendai Bay, Japan, in spring 2005. A hierarchical model was used with two spatial strata: (i) presence and absence of each species; and (ii) biomass density of each species, given its presence; and six environmental covariates (surface water temperature, salinity, and chlorophyll, and near-seabed water temperature, salinity, and depth). The results indicate non-linear responses of the two indices to the environmental covariates. In addition, the biomasses estimated by the GAMs were comparable with estimates based on conventional, stratified-random sampling for each species. GAMs are very useful for (i) investigating the effects of environmental factors on the distributions of biological organisms, and (ii) predicting the distributions of animal densities in unsurveyed areas.

abundance estimation, distribution model, echosounder, ecosystem, fish, GAM, habitat model, marine ecology

Introduction

Knowledge of the abundances and dynamics of target fish populations is basic to effective fishery management. It is also important to know how environmental factors affect the distributions of such populations. In 1991, ICES held a workshop on the applications of spatial techniques to acoustic-survey data, where spatial-statistical methods were discussed (ICES, 1993). One of the methods considered was a generalized additive model (GAM). A GAM is a non-parametric, regression technique not restricted by linear relationships, and it is flexible regarding the statistical distribution of the data (Swartzman et al., 1995). GAMs have been applied to acoustic datasets to investigate relationships between environmental factors and horizontal distributions of herring in the Northeastern Atlantic (Maravelias, 1997; Maravelias and Reid, 1997; Bailey et al., 1998; Maravelias et al., 2000a, 2000b), and walleye pollock in the Bering Sea (Swartzman et al., 1994, 1995). GAMs have also been used to depict the effects of environmental factors on the vertical distribution of various fish (Swartzman, 1997; Swartzman et al., 1999; Taylor and Rand, 2003; Winter et al., 2007). However, there have been few applications of GAMs to create maps of the horizontal distribution of fish (Beare et al., 2002). Moreover, GAMs have not been used commonly on echosounder data to estimate the biomasses of commercially fished species and the measurement of the associated uncertainty of the estimate. A stratified, random-sampling method (Jolly and Hampton, 1990) is often used to estimate biomass and its associated uncertainty, but it does not account for the effects of environmental variables. GAMs, however, are flexible enough to model relationships between biomasses and variables describing their environments.

Sendai Bay is located in the northeastern part of Japan, where there are strong seasonal changes in the location of the Subarctic, western-boundary current. In this area, the cold, low-salinity water of the Oyashio Current converges with the warm, high-salinity water of the subtropical, western-boundary Kuroshio Current. Distributions of animals in Sendai Bay change markedly with the seasonal changes in oceanographic conditions. In spring, when the influence of the Oyashio is strong, krill (Euphausia pacifica) and sand lance (Ammodytes personatus) are the dominant pelagic species in the bay. Japanese anchovy (Engraulis japonicus) start migrating from south to north at this time, and some early migrant anchovy are observed in the bay. These three species are important for the commercial fishery there (Nagashima, 2000, 2007; Taki, 2002). For these reasons, Sendai Bay is an ideal region to test the applicability of GAM-based, spatial modelling to reveal relationships between environmental factors and the distributions of these species.

In this paper, GAMs are applied to an acoustic and physical-oceanographic dataset to: (i) examine the effects of environmental variables (e.g. oceanographic conditions) on the distributions of Japanese anchovy, sand lance, and krill; (ii) predict their spatial distributions; and (iii) estimate their biomasses.

Methods

The survey was conducted in Sendai Bay from 11 to 25 April 2005, as a part of the coastal component of the Japanese Whale Research Programme in the western North Pacific—Phase II (JARPN II). The survey area was stratified into seven blocks (A–G block) based on seabed depth (Figure 1) and a zigzag track line covered each one. In each block, the surveyed distance divided by the square root of the survey area was larger than 6, the number needed for precise estimates of biomass according to Aglen (1989). Acoustic and oceanographic surveys were conducted with the trawler RV “Takuyo maru” (Miyagi prefecture, 120 GT). Track lines were run at 9 knots. A 38- and 120-kHz echosounder system (Simrad EK500) was used for the acoustic survey. The data were recorded with the aid of Echoview (Sonar Data Co., Ltd, Australia). The EK500 was calibrated before the survey using a copper sphere and following the methods of Foote (1982). A trawlnet with a mouth opening 7 m wide by 3.5 m high, and with a 3-mm liner codend, was used to validate the species of target echoes in the survey area. Trawls were conducted in the blocks where scatterer densities were high. Representative echoes were selected, based on knowledge gathered during surveys of this area since 1999. Oceanographic conditions from the sea surface to 1 m above the seabed were recorded by CTD casts at 42 stations. Sea-surface temperature (SST), salinity, and chlorophyll a concentration were recorded every minute along the track line.

Figure 1.

(a) Survey area in Sendai Bay, Japan, with a contour map of surface-water temperature. The area was divided into seven blocks (A–G) to estimate the biomasses of krill, Japanese anchovy, and sand lance using a stratified, random-sampling method. Depths of isobaths are indicated. (b) Contour map of near-seabed water temperature.

Identification of the species contributing echoes was based on the differences of mean volume-backscattering strength (⁠⁠) measured at 120 and 38 kHz (ΔS_v120–38) and the at 120 kHz. Ranges of these values were calculated, based on the results of targeted trawls. Respective ranges of and ΔS_v120–38, which indicate these three species, were: krill, −70 to −45 dB and 10 to 20 dB; adult Japanese anchovy, −50 to −35 dB and −10 to 5 dB; and adult sand lance, −70 to −50 dB and −10 to 10 dB. Nautical-area-backscattering coefficients (s_A) were calculated with the aid of Echoview for each species, for every 0.1 nautical mile transect distance, from the sea surface to the seabed. Theoretical target strength (TS) vs. log-length relationships were derived for krill and sand lance at 120 kHz, based on the distorted-wave Born approximation model (Stanton and Chu, 2000): TS = 50.7 log(L_T) − 150.4 and TS = 49.3 log(L_S) − 160.2, respectively (R. Matsukura, unpublished data), where total length (L_T) and standard length (L_S) are in millimetres. For Japanese anchovy, TS = 20 log(T_L) − 72.5, where T_L is in centimetres (Iversen et al., 1993). Biomass densities (⁠⁠, t nautical mile⁻²) of the three species were estimated as ⁠, where σ_bs is the backscattering cross-sectional area (=10^(0.1×TS)), f_i the frequency, and W_i the weight of the ith length class. As a reference value to evaluate the performance of the GAMs, biomasses of the three species and their variances were estimated based on the stratified, random-sampling method of Jolly and Hampton (1990).

Interactions between the biomasses of the three species and environmental factors were investigated at a longitudinal–latitudinal scale of 30″ × 30″, i.e. the grid-cell size. Six environmental factors [sea surface water temperature (°C), salinity (psu), chlorophyll a concentration (mg m⁻³), and near-seabed (1 m above seabed) water temperature, salinity, and depth (m)] were used as covariates in the spatial modelling. In the survey area, Japanese anchovy was typically distributed in a near-surface layer, whereas krill and sand lance were distributed in both near-surface and near-seabed layers. The value of each environmental factor in each grid cell was estimated using kriging methods with the aid of a GIS, Marine Explorer (Environmental Simulation Laboratory, Japan).

A hierarchical-spatial structure with two strata (first stratum: presence and absence of a biological organism; second stratum: biomass density given its presence) was developed. The two-strata GAMs were used to estimate biomasses and create distribution maps of fish and krill and to examine the effects of environmental factors on their distributions. A GAM with a binomial error distribution that had a logistic-link function was assumed for the first stratum. A GAM with either a Gamma error distribution with the log-link function or a Gaussian error distribution with identical link function was used for the second stratum modelling. In the case of the GAM with a Gaussian error distribution and the identical link function, natural logarithms of were used as dependent variables. All environmental covariates were considered initially for both strata. When the number of observations was few, data for that covariate were pooled into groups of either positive or negative values. Smoothness parameters were estimated with generalized cross-validation (GCV).

Effective covariates were selected following Wood (2001). Covariates were deleted from the models if the following three criteria were met: (i) the estimated degrees of freedom was close to 1; (ii) the confidence interval was zero everywhere; and (iii) the GCV score decreased when the term was dropped. Models with the lowest GCV scores were selected. However, if the use of a selected model resulted in: (i) unrealistically large biomasses relative to those predicted by the stratified, random-sampling method, or (ii) convergence errors in the estimation of uncertainty of biomass, the covariates were deleted until approximate significance levels of all smoother terms were <0.001. Wood (2001) indicated that the deletion of terms is sometimes subjective, and that if the deletion of a term results in a small change in the GCV score, the term should be deleted.

This GAM-based analysis used the “mgcv” package (version 1.30-27) of the R program (R Development Core Team, 2007). “Deviance explained” (analogous to variance in a linear regression), adjusted r², and GCV scores were calculated. The shapes of the functional forms for the selected covariates were plotted. When the slopes of the functional forms are positive, the covariates are related positively to the dependent variables, or vice versa. Selected models allowed prediction of in unsurveyed cells. Biomass density for each species was estimated as the product of the results from the two strata. The coefficient of variation, CV of ⁠, was estimated using a parametric bootstrap with 1000 iterations, assuming lognormally distributed errors. GAM-estimated was compared with the echosounder-estimated ⁠.

Results

Contour maps of the surface and near-seabed water temperatures are presented in Figure 1. Generally, both surface water and near-seabed water temperatures were high in coastal areas and low offshore. However, the distribution of surface water temperatures was heterogeneous compared with near-seabed temperatures.

Target-identification trawls were conducted at six stations (Table 1). Echoes from krill and sand lance were each identified at three stations. Because there was no opportunity to identify echoes of Japanese anchovy during the survey, their identification was based on previous results (Nagashima, 2006).

Table 1.

Summary of target trawls.

Block	Station #	Positions	Total catch (kg)	Species (% in total catch)
B	1	38°05′N 141°08′E	0.87	Sand lance (100)
C	1	38°03′N 141°22′E	2	Krill (100)
	2	38°13′N 141°26′E	0.65	Sand lance (100)
D	1	37°58′N 141°32′E	2.4	Krill (≈100)
	2	38°00′N 141°45′E	3.3	Krill (≈100)
F	1	37–51°N 141°16′E	1.5	Sand lance (94.4)

Krill were found in all the blocks, whereas Japanese anchovy and sand lance were found only in Blocks B, C, E, and F, and in B–F, respectively (Figure 2). Initially, GAMs were calculated using data from all the blocks. However, because a convergence error was encountered when estimating uncertainty in the estimate of biomass, GAMs were recalculated using only blocks with a species present. Selected GAMs for each species are summarized in Table 2. Selected covariates differed among the various species.

Figure 2.

Distributions of krill (top), Japanese anchovy (middle), and sand lance (near-seabed) observed along the track line by the echosounder (left) and predicted by GAMs (right). Circles in the plots on the left represent observed densities (t nautical mile⁻²). Contours of biomasses (t) estimated by GAMs are drawn in the plots on the right.

Table 2.

Selected GAM-based, spatial models for krill, Japanese anchovy, and sand lance.

	Krill				Japanese anchovy				Sand lance
Parameter	First stratum		Second stratum		First stratum		Second stratum		First stratum		Second stratum
Family	Binomial		Gamma		Binomial		Gaussian		Binomial		Gamma
Link function	Logit		Log		Logit		Identity		Logit		Log
Adjusted r2	0.86		0.15		0.06		0.29		0.23		0.07
Deviance explained (%)	81.0		22.2		13.9		35.0		19.1		27.6
GCV score	0.27		3.39		0.49		2.91		1.13		2.62
	d.f.	p-value	d.f.	p-value	d.f.	p-value	d.f.	p-value	d.f.	p-value	d.f.	p-value
Covariates
Near-seabed temperature	8.34	<0.001	6.43	<0.001	–	–	–	–	–	–	7.30	<0.001
Near-seabed salinity	7.37	<0.001	2.57	<0.001	1.81	<0.01	–	–	–	–	–	–
SST	8.65	<0.001	3.43	<0.001	3.08	<0.1	–	–	–	–	–	–
Sea surface salinity	8.84	<0.001	–	–	8.55	<0.001	–	–	–	–	–	–
Sea surface chlorophyll a concentrations	8.96	<0.001	4.41	<0.001	8.64	<0.001	–	–	6.96	<0.001	5.36	<0.001
Depth	8.69	<0.001	7.87	<0.001	1.07	<0.001	6.86	<0.001	5.22	<0.001	8.66	<0.001

Approximate significance levels (p-value) and degrees of freedom (d.f.) are displayed for each of the covariates.

The shapes of the functional forms for selected covariates of the first and second strata are illustrated in Figures 3–5. These indicate that the three species displayed non-linear responses to the covariates. For instance, krill presence is related non-linearly to near-seabed temperature, whereas its biomass density exhibits a monotonic, negative response. Functional forms of near-seabed salinity were similar for both the first and the second strata for krill. All environmental covariates were selected for both the first and the second strata for krill, except for sea-surface salinity in the second stratum. “Deviance explained” for the first and the second strata was 81.0 and 22.2%, respectively.

Figure 3.

Smoothed fits of covariates modelling (a) the presence–absence and (b) the biomass density of krill. Tick marks on the x-axis are observed datapoints. The y-axis represents the spline function. Dashed lines indicate 95% confidence bounds.

Figure 4.

Smoothed fits of covariates modelling (a) the presence–absence and (b) the biomass density of Japanese anchovy. Tick marks on the x-axis are observed datapoints. The y-axis represents the spline function. Dashed lines indicate 95% confidence bounds.

Figure 5.

Smoothed fits of covariates modelling (a) the presence–absence and (b) the biomass density of sand lance. Tick marks on the x-axis are observed datapoints. The y-axis represents the spline function. Dashed lines indicate 95% confidence bounds.

For Japanese anchovy, different covariates were selected for models of the first and the second strata. For the first stratum, they included SST, salinity, and chlorophyll a, and near-seabed salinity and depth. For the second stratum, only depth was selected. “Deviance explained” for the first and the second strata was 13.9 and 35.0%, respectively.

For sand lance, the covariates for the first stratum included sea-surface chlorophyll a and depth. The covariates for the second stratum included sea-surface chlorophyll a, near-seabed temperature, and depth. The shape of the functional forms indicated that biomass density of sand lance was positively related to near-seabed temperature. The shapes of the functional forms for sea-surface chlorophyll a were similar for the first and the second strata. Although depth was selected for the first and the second strata, the responses differed. The probability of presence peaked at ∼50 m, decreased down to ∼150 m, then increased towards 200 m. Biomass density peaked at 100 m and declined at both shallow and deep seabed depths. “Deviance explained” for the first and the second strata was, respectively, 19.1 and 27.6%.

Biomass estimates of krill, Japanese anchovy, and sand lance based on GAMs were, respectively, 1.49 × 10⁶ t (CV = 0.07), 0.71 × 10³ t (CV = 0.48), and 7.54 × 10³ t (CV = 0.02). Of the bootstrap seeds, 3% were discarded as outliers during the calculation of uncertainty in the krill-biomass estimate. For comparison, biomass estimates of krill, Japanese anchovy, and sand lance based on the stratified, random-sampling method were, respectively, 1.56 × 10⁶ t (CV = 0.14), 1.32 × 10³ t (CV = 0.33), and 7.68 × 10³ t (CV = 0.19). The estimated values of were similar to the observed ⁠, indicating that the models fitted the data well.

Predicted spatial distributions are illustrated in Figure 2. There are large discrepancies in the cells with large observed ⁠. In fact, estimated by GAMs tended to be underestimated in cells with high observed ⁠.

Discussion

The study demonstrated that GAM-based, spatial modelling could be used to create plausible fish- and zooplankton-distribution maps from acoustic-survey data. The resulting distribution maps of krill, Japanese anchovy, and sand lance match the animal distributions observed by the echosounder along the track lines. Although geostatistical methods, especially kriging, have often been used to create animal-distribution maps using acoustic-survey data (Rivoirard et al., 2000, for reference), they do not account for the effects of environmental factors on the distributions of such organisms. In contrast, GAMs can take account of such factors by treating them as covariates.

Maravelias (1997) used two-strata GAMs (first stratum: presence and absence; second stratum: number of individuals given presence) to study trends in abundance and geographic distribution of North Sea herring. In the current study, two-strata GAMs were used to estimate biomasses, create distribution maps, and explore the effects of environmental factors. Unfortunately, the low values obtained for “Deviance explained”, especially for Japanese anchovy and sand lance, indicate that these GAMs need improvement before they can be applied successfully to predict biomasses of these species in Sendai Bay in spring. The inclusion of different environmental covariates might improve the results. For instance, because sand lance have strong preferences for certain seabed types (Kobayashi et al., 1995), the inclusion of seabed metrics as covariates might improve these GAM results. It is therefore important to collect data on potentially influential covariates.

Taki et al. (1996) and Taki and Ogishima (1997) used net samples to demonstrate that distributions of E. pacifica are related to water temperature ranging from 7 to 14°C and salinity <34 psu, and that their biomass was high at depths of 100–150 m in April. Finer-scale interactions between animal distributions and environmental factors have been revealed by GAMs (Swartzman, 1997; Swartzman et al., 1999; Winter et al., 2007). In our results, the GAMs confirmed that krill presence is related positively to near-seabed temperature, but also revealed that its biomass density is related negatively to temperature. Differences in the functional forms between the first and the second strata suggest that the densities of krill schools may be high in areas with high water temperature, and low school densities may be found in areas with low water temperature.

Previous qualitative studies have indicated that oceanographic conditions have varying effects on the distributions of Japanese anchovy, depending on the season (Nagashima, 2006). Because April is the time of the early migration of Japanese anchovy from the south to Sendai Bay, data from one year are insufficient to assess any environmental linkages to their biomass or distribution. However, a GAM analysis of a multiyear dataset might reveal the environmental parameters that control the timing of the northerly migration of Japanese anchovy.

Kobayashi et al. (1995) reported that the biomass of adult sand lance was high in the central region of Sendai Bay, where the water depth was 40–70 m and the seabed was medium to large gravel. Small sand lance (10–13 cm) dominated where water depth was 40–50 m, whereas large sand lance (>13 cm) dominated where water depth was 50–70 m. The amount of commercially caught sand lance was related positively to the mean SST in March. The shape of the functional forms of the near-seabed temperature for the second stratum in our study also indicated a positive relationship between the biomass of sand lance and temperature. Comparing the shapes of the functional forms between the first and the second strata might differentiate the distribution patterns of sand lance according to their length.

Perhaps different combinations of selected covariates and shapes of functional forms for the three species could identify different ecological niches and habitat requirements. In future, ecological theory should be incorporated explicitly into the GAM analysis, using appropriate statistical methods to select the environmental covariates to be tested (Austin, 2007).

Relatively low CVs for krill and sand lance might have resulted from accounting for covariates in the GAMs. The relatively high CV for Japanese anchovy could have resulted from poor GAM performance when applied to data from patchily distributed animals. The residuals indicate that it could be difficult to get GAMs to predict cells with large biomasses.

The distribution maps produced with the GAMs are in good agreement with the observed distributions of krill, Japanese anchovy, and sand lance along the survey tracks. An objective criterion for evaluating such GAM maps is being sought. In addition, for future work, many other predictive habitat-distribution models could be considered (see Guisan and Zimmermann, 2000; Austin, 2002, 2007; Leathwick et al., 2006, for reviews).

The accuracies of the species-identification algorithms can be improved (Nagashima et al., 2008). Other sources of measurement uncertainty, such as those stemming from the calibration, target-strength model, length-to-weight model, bubble attenuation, signal thresholding, and survey-area definition, could also affect the accuracy of the biomass estimation. The effect of these uncertainties on biomass estimation should also be investigated (Demer, 2004).

This study demonstrates that GAM-based, spatial modelling is useful to create plausible distribution maps and to estimate biomass, while considering environmental factors. Such GAM-based methods could even be extended to studies of predator–prey relationships (Smout and Lindstrøm, 2007).

Acknowledgements

We thank the captain of RV “Takuyo-Maru”, Hiroaki Kimura, and his crew who assisted us in collecting this valuable dataset. The study was conducted as a part of the coastal component of JARPN II; we thank coordinators Hidehiro Kato (Tokyo University of Marine Science and Technology) and Yoshihiro Fujise (the Institute of Cetacean Research) for their valuable comments on the survey. Kazushi Miyashita (Hokkaido University) kindly provided critical comments on the acoustic-data analysis methods. The survey was supported by the Fisheries Agency of Japan, the Fisheries Research Agency of Japan, the Miyagi Prefecture, and the Institute of Cetacean Research. We thank these institutions for their support.