TY - JOUR
AU - Barnett, Lewis A. K.
AB - Several approaches have been developed over the last decade to simultaneously estimate distribution or density for multiple species (e.g. “joint species distribution” or “multispecies occupancy” models). However, there has been little research comparing estimates of abundance trends or distribution shifts from these multispecies models with similar single-species estimates. We seek to determine whether a model including correlations among species (and particularly species that may affect habitat quality, termed “biogenic habitat”) improves predictive performance or decreases standard errors for estimates of total biomass and distribution shift relative to similar single-species models. To accomplish this objective, we apply a vector-autoregressive spatio-temporal (VAST) model that simultaneously estimates spatio-temporal variation in density for multiple species, and present an application of this model using data for eight US Pacific Coast rockfishes (Sebastes spp.), thornyheads (Sebastolobus spp.), and structure-forming invertebrates (SFIs). We identified three fish groups having similar spatial distribution (northern Sebastes, coastwide Sebastes, and Sebastolobus species), and estimated differences among groups in their association with SFI. The multispecies model was more parsimonious and had better predictive performance than fitting a single-species model to each taxon individually, and estimated fine-scale variation in density even for species with relatively few encounters (which the single-species model was unable to do). However, the single-species models showed similar abundance trends and distribution shifts to those of the multispecies model, with slightly smaller standard errors. Therefore, we conclude that spatial variation in density (and annual variation in these patterns) is correlated among fishes and SFI, with congeneric fishes more correlated than species from different genera. However, explicitly modelling correlations among fishes and biogenic habitat does not seem to improve precision for estimates of abundance trends or distribution shifts for these fishes. Introduction There are several benefits to simultaneously analysing the distribution and density of multiple species within a natural community. Multispecies models of spatial distribution can estimate associations among species (Latimer et al., 2009; Ovaskainen et al., 2016; Thorson et al., 2015a, 2016a), such that the presence or absence of a given species can be used as an indicator of habitat for other species when reliable habitat variables are otherwise lacking (Ovaskainen et al., 2010). Multispecies models fitted to presence/absence data (termed “multispecies occupancy models”) can also be used in some cases to identify the impact of management actions more efficiently than using single-species occupancy models (Zipkin et al., 2010). Furthermore, research shows that estimating the distribution for each species individually and then summarizing community-level properties by stacking results from single-species analyses can result in improper inference about ecological communities (Clark et al., 2014). The predictive performance of species distribution models is often improved when available covariates are included that are informative about habitat quality. Unfortunately, environmental variables associated with habitat quality are difficult to measure for many species, including demersal marine fishes. To overcome this difficulty, new species distribution modelling (SDM) techniques may allow differences in habitat to be inferred from spatial variation in the density of species with similar habitat requirements (Latimer et al., 2009; Ovaskainen et al., 2010). For example, joint species distribution models have previously been used to show strong covariation in population density among US Pacific Coast rockfishes and thornyheads (Sebastes and Sebastolobus spp.), and these correlations imply that the population density of one species is informative about the density of correlated species (Thorson et al., 2015a). Similarly, joint dynamic species distribution models (JDSDMs) can estimate abundance trends for infrequently encountered species and have revealed similarities in spatio-temporal dynamics among related butterfly species (Thorson et al., 2016a,,b). However, JDSDMs have not previously been used to explore associations between fishes and species that are associated with specific habitat features (e.g. structure-forming invertebrates, SFI). Marine fishes are intensively managed in many parts of the developed world, and the management of marine fisheries is strongly linked to estimates of population status and productivity from population models (termed “stock assessment models”) throughout North America and Europe (Methot, 2009; Maunder and Punt, 2013). Although these stock assessment models often integrate many different types of information, time-series that are proportional to population abundance (“abundance indices”) are often among the most critical (Francis, 2011). For this reason, there is considerable research regarding best practices for minimizing error when estimating abundance indices for fishes from survey data (Walters, 2003; Maunder and Punt, 2004; Shelton et al., 2014). Similarly, survey data are increasingly used to estimate shifts in fish distribution over time (e.g. due to climate change), and distribution shifts are often measured by estimating the centroid of the population's distribution and shifts in this centroid over time (Perry et al., 2005; Pinsky et al., 2013). Research suggests that spatio-temporal models are statistically efficient and can improve precision for estimates of abundance indices or distribution shifts relative to nonspatial models, given limited available data (Thorson et al., 2015b, 2016b). Recently, novel methods have been proposed for estimating abundance indices by simultaneously fitting a JDSDM to data for multiple species (Thorson et al., 2016). However, there is little research comparing the single- and multispecies approaches to estimating abundance indices for marine fishes. For three reasons, Pacific rockfishes and their close relatives provide an interesting example when studying associations between fishes and species that affect habitat suitability (“biogenic habitat”) or the potential benefit of these associations when estimating abundance indices or distribution shifts. Most importantly, Pacific rockfishes manifest an astounding diversity of species, with more than 65 species co-occurring in the Northeast Pacific (Hyde and Vetter, 2007) and exhibit a wide range of life history strategies (Love et al., 2002; Mangel et al., 2007). Given this life history diversity, rockfishes likely include species whose spatial distributions are both strongly correlated and relatively uncorrelated with SFI. Second, Pacific rockfishes differ in functional traits related to the feeding type and efficiency [eye and gill raker size, Ingram and Shurin, 2009], so species with similar spatial distribution and feeding types might exhibit correlated changes in productivity over time in response to variable food supply. Therefore, bottom-up drivers of abundance or distribution changes would result in correlated abundance or distribution changes over time for species with similar feeding types. Third, many Pacific rockfishes have low and extremely variable population densities (Thorson et al., 2011), such that single-species estimates of trends in population abundance or population distribution are frequently imprecise (Thorson et al., 2015b, 2016b). Given these characteristics of the rockfish assemblage, the inclusion of information about species associations and biogenic habitat when estimating population abundance may increase precision and thereby improve stock assessments. Given the potential benefit of estimating habitat quality from the density of co-occurring marine species when estimating abundance indices, we seek to simultaneously estimate the density of Pacific rockfishes and structure-forming invertebrates at a coastwide scale. Specifically, we seek to answer three questions: (i) do Pacific rockfishes have an association with structure-forming invertebrates on the US West Coast? (ii) Is this association similar or variable among rockfish species? and (iii) Does the inclusion of information regarding co-occurrence (either among rockfishes or between rockfish and SFI) improve predictions of local rockfish density or increase precision when estimating rockfish abundance trends or population distribution? To address these questions, we develop a vector-autoregressive spatio-temporal (VAST) model for jointly analysing catch-rate data for fish and structure-forming invertebrates and apply the model to data for eight rockfish species and SFI during 2003–2014. Methods Pacific rockfishes Pacific rockfishes (genus Sebastes) and thornyheads (genus Sebastolobus), hereafter collectively called “rockfishes”, are one of the dominant species groups within the assemblage of bottom-associated fishes off the US West Coast. Pacific rockfishes in this region are monitored by the West Coast groundfish bottom trawl survey (WCGBTS) conducted annually by the Northwest Fisheries Science Center since 2003 (Bradburn et al., 2011). The WCGBTS covers areas between the Canada and Mexico borders in 55–1280 m depth, and survey stations for each year are chosen at random within strata defined by depth and latitude (two regions divided at Point Conception, CA). Four commercial vessels (20–28 m length) are chartered each year to sample from mid-May to late October, conducting ca. 15-min tows at a speed of 2 knots using a standard Aberdeen-type trawl with a 3.8-cm mesh codend liner, 25.9-m headrope, and 31.7-m footrope. All fishes and invertebrates are sorted at sea to the lowest possible taxon, and their wet weight is measured. For the purposes of our analysis, we take the midpoint of each haul to represent the location of each biological sample. We analyse these survey data between the years 2003 and 2014, focusing on structure-forming invertebrates and eight species of Pacific rockfish (Table 1) that are frequently captured within the survey and for which there was previous documentation of association with structure-forming invertebrates at fine spatial scales (Love et al., 2002). We aggregate the structure-forming invertebrate taxa into a single grouping to obtain adequate encounter rates for estimating the distribution for structure-forming invertebrates. This SFI group primarily consists of sponges (phylum Porifera), anemones (order Actiniaria), and sea pens (order Pennatulacea), along with fewer observations of true corals (subclass Hexacorallia) and other soft corals (subclass Octocorallia). Although the survey is primarily designed to capture demersal fishes and is not as effective as visual methods for assessing structure-forming invertebrates, it is the primary source of data for estimating spatio-temporal associations between demersal fishes and biogenic habitat at large spatial and temporal scales off the US West Coast. Bottom-trawl samples have been shown to be a good predictor of biogenic habitat distribution in areas such as the eastern Bering Sea based on validation using camera surveys (Rooper et al., 2016). Table 1 List of taxa (common and scientific name), the abbreviations used to indicate taxa in plots, and the total number of encounters during 2003–2014. Common name . Scientific name . Plotting code . Encounters . Structure-forming invertebrates (SFIs) — SFI 6383 Longspine thornyhead Sebastolobus altivelis L. spine 2758 Shortspine thornyhead Sebastolobus alascanus S. spine 3891 Darkblotched rockfish Sebastes crameri Dark 1338 Pacific Ocean perch Sebastes alutus POP 547 Sharpchin rockfish Sebastes zacentrus Sharp 490 Splitnose rockfish Sebastes diploproa Split 1619 Stripetail rockfish Sebastes saxicola Stripe 1630 Greenspotted rockfish Sebastes chlorostictus Green 434 Common name . Scientific name . Plotting code . Encounters . Structure-forming invertebrates (SFIs) — SFI 6383 Longspine thornyhead Sebastolobus altivelis L. spine 2758 Shortspine thornyhead Sebastolobus alascanus S. spine 3891 Darkblotched rockfish Sebastes crameri Dark 1338 Pacific Ocean perch Sebastes alutus POP 547 Sharpchin rockfish Sebastes zacentrus Sharp 490 Splitnose rockfish Sebastes diploproa Split 1619 Stripetail rockfish Sebastes saxicola Stripe 1630 Greenspotted rockfish Sebastes chlorostictus Green 434 Open in new tab Table 1 List of taxa (common and scientific name), the abbreviations used to indicate taxa in plots, and the total number of encounters during 2003–2014. Common name . Scientific name . Plotting code . Encounters . Structure-forming invertebrates (SFIs) — SFI 6383 Longspine thornyhead Sebastolobus altivelis L. spine 2758 Shortspine thornyhead Sebastolobus alascanus S. spine 3891 Darkblotched rockfish Sebastes crameri Dark 1338 Pacific Ocean perch Sebastes alutus POP 547 Sharpchin rockfish Sebastes zacentrus Sharp 490 Splitnose rockfish Sebastes diploproa Split 1619 Stripetail rockfish Sebastes saxicola Stripe 1630 Greenspotted rockfish Sebastes chlorostictus Green 434 Common name . Scientific name . Plotting code . Encounters . Structure-forming invertebrates (SFIs) — SFI 6383 Longspine thornyhead Sebastolobus altivelis L. spine 2758 Shortspine thornyhead Sebastolobus alascanus S. spine 3891 Darkblotched rockfish Sebastes crameri Dark 1338 Pacific Ocean perch Sebastes alutus POP 547 Sharpchin rockfish Sebastes zacentrus Sharp 490 Splitnose rockfish Sebastes diploproa Split 1619 Stripetail rockfish Sebastes saxicola Stripe 1630 Greenspotted rockfish Sebastes chlorostictus Green 434 Open in new tab Vector-autoregressive spatio-temporal (VAST) model We seek to estimate the association among fishes and structure-forming invertebrates and, therefore, model correlations among density d(s,c,t) for each taxon c (indicating fish species or the SFI group) at location s and time t (all symbols are summarized in Table 2). To do so, we build upon recent research regarding JDSDMs. In particular, we propose a VAST model, where the probability distribution for catch data bi is decomposed into two components representing (i) the probability of encounter psi,ci,ti for the location si ⁠, taxon ci ⁠, and year ti of the ith sample, and (ii) the expected catch rate rsi,ci,ti ⁠, given that taxon ci is encountered. Decomposing catch rates into encounter-probability p and positive catch rates r is commonly conducted using delta models (Maunder and Punt, 2004; Martin et al., 2005), although delta models have not previously been used within JDSDMs. Using a delta model allows us to separately identify species with similar distribution (similarities in occupied habitat) vs. similar density (similarities in hotspots within their distribution). Therefore, we specify: Pr⁡(bi=B)={1−p(si,ci,ti) if B=0p(si,ci,ti)×Lognormal{B|log[wi×r(si,ci,ti)],σc2} if B>0(1) where Lognormal(x|μ,σ2) is a lognormal probability distribution function for value x ⁠, given a log-mean of μ and a variance of σ2 ⁠, and wi is the area swept for the ith sample. Table 2 List of symbols representing indices, data, fixed effects, random effects, and derived quantities defined in the main text. Name . Symbol . Dimension . Type . Sample i — Index Taxon c — Index Location s — Index Year t — Index Taxon (species or species-group) c — Index Vessel v — Index Probability distribution for catches B — Random variable Catch data bi ni Data Area-swept for each sample wi ni Data Area associated with each location a(s) ns Data Statistic associated with each location x(s) ns Data Variance in positive catch rates σ2 — Fixed effect Intercept for p γpc,t nc×nt Fixed effect Intercept for r γrc,t nc×nt Fixed effect Decorrelation distance in ɛps,c,t κp — Fixed effect Decorrelation distance in ɛrs,c,t κr — Fixed effect Geometric anisotropy H 2×2 Fixed effect Factor approximation to Vɛp Lɛp nc×nf Fixed effect Factor approximation to Vɛr Lɛr nc×nf Fixed effect Factor approximation to Vδp Lδp nc×nf Fixed effect Factor approximation to Vδr Lδr nc×nf Fixed effect Spatio-temporal variation in p ɛps,c,t ns×nc×nt Random effect Spatio-temporal variation in r ɛrs,c,t ns×nc×nt Random effect Vessel effect for p δp(c,v) nc×nv Random effect Vessel effect for r δr(c,v) nc×nv Random effect Encounter probability p(s,c,t) ns×nc×nt Derived quantity Positive catch rates r(s,c,t) ns×nc×nt Derived quantity Local density d(s,c,t) ns×nc×nt Derived quantity Spatio-temporal variation in p in year t Ep(t) ns×nc Derived quantity Spatio-temporal variation in r in year t Er(t) ns×nc Derived quantity Spatial correlation in ɛps,c,t Rp ns×ns Derived quantity Spatial correlation in ɛrs,c,t Rr ns×ns Derived quantity Correlation among species in ɛps,c,t Vɛp nc×nc Derived quantity Correlation among species in ɛrs,c,t Vɛr nc×nc Derived quantity Correlation among species in δp(c,v) Vδp nc×nc Derived quantity Correlation among species in δr(c,v) Vδr nc×nc Derived quantity Index of abundance I(c,t) nc×nt Derived quantity Center of spatial distribution X(c,t) nc×nt Derived quantity Name . Symbol . Dimension . Type . Sample i — Index Taxon c — Index Location s — Index Year t — Index Taxon (species or species-group) c — Index Vessel v — Index Probability distribution for catches B — Random variable Catch data bi ni Data Area-swept for each sample wi ni Data Area associated with each location a(s) ns Data Statistic associated with each location x(s) ns Data Variance in positive catch rates σ2 — Fixed effect Intercept for p γpc,t nc×nt Fixed effect Intercept for r γrc,t nc×nt Fixed effect Decorrelation distance in ɛps,c,t κp — Fixed effect Decorrelation distance in ɛrs,c,t κr — Fixed effect Geometric anisotropy H 2×2 Fixed effect Factor approximation to Vɛp Lɛp nc×nf Fixed effect Factor approximation to Vɛr Lɛr nc×nf Fixed effect Factor approximation to Vδp Lδp nc×nf Fixed effect Factor approximation to Vδr Lδr nc×nf Fixed effect Spatio-temporal variation in p ɛps,c,t ns×nc×nt Random effect Spatio-temporal variation in r ɛrs,c,t ns×nc×nt Random effect Vessel effect for p δp(c,v) nc×nv Random effect Vessel effect for r δr(c,v) nc×nv Random effect Encounter probability p(s,c,t) ns×nc×nt Derived quantity Positive catch rates r(s,c,t) ns×nc×nt Derived quantity Local density d(s,c,t) ns×nc×nt Derived quantity Spatio-temporal variation in p in year t Ep(t) ns×nc Derived quantity Spatio-temporal variation in r in year t Er(t) ns×nc Derived quantity Spatial correlation in ɛps,c,t Rp ns×ns Derived quantity Spatial correlation in ɛrs,c,t Rr ns×ns Derived quantity Correlation among species in ɛps,c,t Vɛp nc×nc Derived quantity Correlation among species in ɛrs,c,t Vɛr nc×nc Derived quantity Correlation among species in δp(c,v) Vδp nc×nc Derived quantity Correlation among species in δr(c,v) Vδr nc×nc Derived quantity Index of abundance I(c,t) nc×nt Derived quantity Center of spatial distribution X(c,t) nc×nt Derived quantity Throughout, we use subscripts to indicate either properties of the ith sample (e.g. area-swept wi ⁠, location si ⁠), or parameters associated with encounter probability or positive catch rates (e.g. κp vs. κr ⁠). Indices or scalars have dimension of zero, indicated using a “—” symbol. Open in new tab Table 2 List of symbols representing indices, data, fixed effects, random effects, and derived quantities defined in the main text. Name . Symbol . Dimension . Type . Sample i — Index Taxon c — Index Location s — Index Year t — Index Taxon (species or species-group) c — Index Vessel v — Index Probability distribution for catches B — Random variable Catch data bi ni Data Area-swept for each sample wi ni Data Area associated with each location a(s) ns Data Statistic associated with each location x(s) ns Data Variance in positive catch rates σ2 — Fixed effect Intercept for p γpc,t nc×nt Fixed effect Intercept for r γrc,t nc×nt Fixed effect Decorrelation distance in ɛps,c,t κp — Fixed effect Decorrelation distance in ɛrs,c,t κr — Fixed effect Geometric anisotropy H 2×2 Fixed effect Factor approximation to Vɛp Lɛp nc×nf Fixed effect Factor approximation to Vɛr Lɛr nc×nf Fixed effect Factor approximation to Vδp Lδp nc×nf Fixed effect Factor approximation to Vδr Lδr nc×nf Fixed effect Spatio-temporal variation in p ɛps,c,t ns×nc×nt Random effect Spatio-temporal variation in r ɛrs,c,t ns×nc×nt Random effect Vessel effect for p δp(c,v) nc×nv Random effect Vessel effect for r δr(c,v) nc×nv Random effect Encounter probability p(s,c,t) ns×nc×nt Derived quantity Positive catch rates r(s,c,t) ns×nc×nt Derived quantity Local density d(s,c,t) ns×nc×nt Derived quantity Spatio-temporal variation in p in year t Ep(t) ns×nc Derived quantity Spatio-temporal variation in r in year t Er(t) ns×nc Derived quantity Spatial correlation in ɛps,c,t Rp ns×ns Derived quantity Spatial correlation in ɛrs,c,t Rr ns×ns Derived quantity Correlation among species in ɛps,c,t Vɛp nc×nc Derived quantity Correlation among species in ɛrs,c,t Vɛr nc×nc Derived quantity Correlation among species in δp(c,v) Vδp nc×nc Derived quantity Correlation among species in δr(c,v) Vδr nc×nc Derived quantity Index of abundance I(c,t) nc×nt Derived quantity Center of spatial distribution X(c,t) nc×nt Derived quantity Name . Symbol . Dimension . Type . Sample i — Index Taxon c — Index Location s — Index Year t — Index Taxon (species or species-group) c — Index Vessel v — Index Probability distribution for catches B — Random variable Catch data bi ni Data Area-swept for each sample wi ni Data Area associated with each location a(s) ns Data Statistic associated with each location x(s) ns Data Variance in positive catch rates σ2 — Fixed effect Intercept for p γpc,t nc×nt Fixed effect Intercept for r γrc,t nc×nt Fixed effect Decorrelation distance in ɛps,c,t κp — Fixed effect Decorrelation distance in ɛrs,c,t κr — Fixed effect Geometric anisotropy H 2×2 Fixed effect Factor approximation to Vɛp Lɛp nc×nf Fixed effect Factor approximation to Vɛr Lɛr nc×nf Fixed effect Factor approximation to Vδp Lδp nc×nf Fixed effect Factor approximation to Vδr Lδr nc×nf Fixed effect Spatio-temporal variation in p ɛps,c,t ns×nc×nt Random effect Spatio-temporal variation in r ɛrs,c,t ns×nc×nt Random effect Vessel effect for p δp(c,v) nc×nv Random effect Vessel effect for r δr(c,v) nc×nv Random effect Encounter probability p(s,c,t) ns×nc×nt Derived quantity Positive catch rates r(s,c,t) ns×nc×nt Derived quantity Local density d(s,c,t) ns×nc×nt Derived quantity Spatio-temporal variation in p in year t Ep(t) ns×nc Derived quantity Spatio-temporal variation in r in year t Er(t) ns×nc Derived quantity Spatial correlation in ɛps,c,t Rp ns×ns Derived quantity Spatial correlation in ɛrs,c,t Rr ns×ns Derived quantity Correlation among species in ɛps,c,t Vɛp nc×nc Derived quantity Correlation among species in ɛrs,c,t Vɛr nc×nc Derived quantity Correlation among species in δp(c,v) Vδp nc×nc Derived quantity Correlation among species in δr(c,v) Vδr nc×nc Derived quantity Index of abundance I(c,t) nc×nt Derived quantity Center of spatial distribution X(c,t) nc×nt Derived quantity Throughout, we use subscripts to indicate either properties of the ith sample (e.g. area-swept wi ⁠, location si ⁠), or parameters associated with encounter probability or positive catch rates (e.g. κp vs. κr ⁠). Indices or scalars have dimension of zero, indicated using a “—” symbol. Open in new tab Using this delta model, we separately develop a spatio-temporal model for encounter probabilities p and positive catch rates r ⁠. We approximate spatio-temporal variation in encounter probability psi,ci,ti using a logit-linked linear predictor : logit⁡[psi,ci,ti]=γpci,ti+ɛpsi,ci,ti+δp(ci,vi)(2) where γpci,ti is an intercept for encounter probability for each taxon c and time t ⁠, ɛpsi,ci,ti approximates spatio-temporal variation in encounter probability (in logit-space), and δp(ci,vi) is a “vessel effect” for the vessel vi conducting the ith sample when catching taxon ci ⁠. Vessel effects are included because the WCGBTS is obtained using 3–4 different vessels per year, and previous research indicates that vessels in each year have small, but important, variation in fishing behaviour and resulting catch rates (Helser et al., 2004; Thorson and Ward, 2014). Expected catch rates when a species is encountered rsi,ci,ti are similarly approximated using a log-linked linear predictor: log⁡[rsi,ci,ti]=γrci,ti+ɛrsi,ci,ti+δr(ci,vi)(3) where γrti,ci is an intercept, ɛrsi,ci,ti is the spatio-temporal variation, and δr(ci,vi) is a vessel effect for expected catch rates. ɛrsi,ci,ti ⁠, ɛpsi,ci,ti ⁠, δp(ci,vi) ⁠, and δr(ci,vi) approximate processes that affect density and catchability, respectively, but are not otherwise modelled explicitly (Thorson et al., 2016). Including these effects in the VAST model allows catch data bi for nearby locations to be correlated (via correlations in encounter probability p or positive catch rate r ⁠) and also improves density predictions at locations that otherwise have little data (Shelton et al., 2014). The VAST model involves specifying a probability distribution for spatio-temporal variation (⁠ ɛps,c,t and ɛr(s,c,t) ⁠) and vessel effects (⁠ δp(c,v) and δr(c,v) ⁠). For each modelled year, we therefore specify a three-dimensional Gaussian process for spatio-temporal variation: vec[Ep(t)]∼MVN(0,Rp⊗Vɛp)(4) where Ep(t) is a matrix composed of ɛps,c,t at every modelled location s and taxon c in a given year t ⁠, vec[Ept] is a vector composed of stacking every column of Ep(t) ⁠, Rp is the correlation matrix approximating similar encounter probability among locations, Vɛp is the covariance matrix approximating similar encounter probability among species, and ⊗ is the Kronecker product such that Rp⊗Vɛp is a covariance matrix between any taxon c and location s in year t [ Er(t) follows an identical distribution, but with Rr and Vɛr in place of Rp and Vɛp]. Spatial correlation Rp between location s and location s+h generally declines with an increased distance |h| between the two locations [sometimes termed Tobler's law of geography; Tobler, 1970]. We specify a Matérn function for this correlation, which includes a parameter κp governing the distance at which locations are essentially uncorrelated (increased κp leads to a decreased decorrelation distance) as well as a transformation matrix H representing geometric anisotropy (the tendency for correlations to decline faster in one direction than another): Rp(s,s+h)=12ν−1Γ(n)×(κp|hH|)n×Kν(κn|hH|)(5) where n is a smoothness parameter [fixed at 1.0; Simpson et al., 2012] and Kn is the Bessel function (⁠ Rr is defined identically but with κr in place of κp ⁠). Including geometric anisotropy is generally important for fishes along a narrow continental shelf like the US West Coast, where correlations decline faster moving onshore–offshore rather than moving alongshore (Thorson et al., 2015b). We do not know a priori which taxa are likely to be more or less correlated. Therefore, we model covariance Vɛp among species using a factor-analysis decomposition: Vɛp=LɛpLɛpT(6) where Lɛp is a nc by nf matrix defining the first nf columns of the Cholesky decomposition of covariance matrix Vɛp ⁠, LɛpT is the matrix-transpose of Lɛp ⁠, and Vɛr is defined identically but with Lɛr in place of Lɛp [see Thorson et al., 2016a,b and Warton et al., 2015 for details regarding this factor-analysis decomposition]. Similarly, we specify a factor-analysis decomposition for the covariance Vδp=LδpLδpT among vessel effects: δp(v)∼MVN(0,Vδp)(7) where δp(v) is the vector of vessel effects affecting encounter probability δpc,v for all taxa c and a given vessel v (and where Vδr is defined identically, but with Lδr in place of Lδp ⁠). This factor-analysis decomposition allows the analyst to select an appropriate number of factors nf for approximating spatio-temporal covariation or covariation among vessels, where 0<nf≤nc ⁠. Specifying a reduced number of factors (⁠ nf<nc ⁠) decreases the number of estimated parameters and, therefore, may result in smaller standard errors for other parameters or more precise predictions of local density (Thorson et al., 2015a). However, reducing the number of factors could also result in biased estimates of abundance trends (e.g. by shrinking dynamics for all species onto a small number of dimensions). We leave exploration of this bias-variance trade-off as a topic for future research and instead specify full rank for each covariance (⁠ nf=nc ⁠) to eliminate this source of potential bias. Parameters are estimated for the VAST model by maximizing the marginal likelihood of fixed effects given available data. We treat the intercept parameters for each species (⁠ γp(c,t) and γrc,t ⁠), the spatial scale of spatio-temporal variation (⁠ κλ and κr ⁠), the shape of geometric anisotropy (two parameters in H ⁠), the covariation among species (⁠ Lɛp and Lɛr ⁠), the covariation among vessels (⁠ Lδp and Lδr ⁠), and the magnitude of residual variation in positive catch rates for each species (⁠ σc2 ⁠) as fixed effects. We treat spatio-temporal variation [ ɛps,c,t and ɛrs,c,t] and catchability variation [ δp(c,v) and δr(c,v)]as random effects (Thorson and Minto, 2015). We define the joint likelihood as the product of the probability of random effects (given fixed effects) and the probability of the data (given random and fixed effects). We then calculate the marginal likelihood of fixed effects while integrating the joint likelihood with respect to random effects. We specifically use the Laplace approximation to approximate the multidimensional integral required to calculate the marginal likelihood (Skaug and Fournier, 2006). The Laplace approximation is implemented using Template Model Builder (Kristensen et al., 2016), and Template Model Builder also provides the gradient of the approximated marginal likelihood with respect to all fixed effects. We use a gradient-based nonlinear minimizer within the R statistical environment (R Core Team, 2015) to identify maximum-likelihood estimate (MLE) of fixed effects. To improve computational efficiency, we use Revolution Open R for low-level parallelization of matrix computations (http://www.revolutionanalytics.com/revolution-r-open), a stochastic partial differential equation (SPDE) approximation for all spatial processes (Lindgren et al., 2011), and the R-INLA software (Lindgren, 2012) to compute the triangulated mesh used in the SPDE approximation. We distribute code for applying the VAST model to other datasets as an R package on the author's website (www.github.com/james-thorson/VAST) and have confirmed that the VAST model provides identical parameter estimates to a previous spatio-temporal index standardization model (SpatialDeltaGLMM, Thorson et al., 2015b) when applied to data for a single species. However, the VAST model also incorporates most capabilities of spatial dynamic factor analysis for monitoring trends in community abundance or conducting species ordination (Thorson et al., 2016a,b). After parameters are estimated, we predict the value of random effects by identifying their values that maximize the joint likelihood, given the data and maximum likelihood estimates of fixed effects. We then use the predicted values for random effects to estimate total biomass I(c,t) for each taxon in each year [an “index of abundance”; Thorson et al., 2015b]: I(c,t)=∑s=1nsa(s)× logit ⁡−1[γp(c,t)+ɛp(s,c,t)]× exp ⁡[γr(c,t)+ɛr(s,c,t)](8) where a(s) is the area associated with location s [ I(c,t) does not include vessel effects δp or δr ⁠, because these are interpreted as representing variation in catchability]. We also estimate the centroid of the distribution X(c,t) for each species in each year [termed “center of gravity”; Thorson et al., 2016b]: X(c,t)=∑s=1nsx(s){a(s) logit ⁡−1[γp(c,t)+ɛp(s,c,t)]× exp ⁡[γr(c,t)+ɛr(s,c,t)]I(c,t)}(9) where x(s) can be any statistic used to summarize distribution. We are particularly interested in estimating shifts in fish distribution northward or southward along the US West Coast, so we define x(s) as the distance north of the equator (in kilometres) for location s ⁠. Standard errors for I(c,t) and X(c,t) are then calculated by Template Model Builder using a generalization of the delta method. We compare model performance when fitting all species simultaneously (the “multi-species analysis”) to a conventional “single-species analysis”, where each species is fitted individually using the VAST model. To compare performance between single- and multispecies models, we compute the Akaike information criterion (AIC) (Akaike, 1974). The AIC is a measure of model “parsimony” (see Figure 1.3 from Burnham and Anderson, 2002), which we use to identify the level of complexity that likely minimizes the combination of bias (from an overly simple model) and imprecision (from an overly complex model). We compute the “single-species” AIC as the sum of the AIC for the VAST model fitted to each individual species; this comparison is justified because the aggregate of single-species models is fitted to the same dataset as the multispecies model. We also conduct a tenfold cross-validation analysis to determine whether multi- or single-species analyses have greater predictive ability. To do so, we divide the data into 10 similarly sized partitions. For the first cross-validation, we estimate model parameters only using data in partitions 2–10, and calculate the probability of data in partition 1 using the predictive distribution, given estimated parameters. This process is repeated for all 10 partitions for the multispecies model. For the single-species model, we conduct this tenfold cross-validation for each species individually, and then sum the resulting log-predictive probabilities for each species. Results Inspection of density estimates for eight fishes and structure-forming invertebrates using the multispecies VAST model shows that species are unevenly distributed throughout the California Current (Figure 1, bottom row). By distribution, the fishes can be broadly classified into three groups: coastwide Sebastes spp. (splitnose, stripetail, and greenspotted), northern Sebastes spp. (POP, sharpchin, and darkblotched), and Sebastolobus spp. (longspine and shortspine thornyheads). The thornyheads are distinguished by having increased densities in the deepest waters farthest from the US coast. Structure-forming invertebrates are found at highest densities offshore near northern Oregon, close to the Oregon–California border and offshore from the south of Monterey Bay through the Southern California Bight. Figure 1 Open in new tabDownload slide Comparison of estimated density functions (averaged across all years for each taxon) for structure-forming invertebrates, six Sebastes, and two Sebastolobus species using the VAST multispecies model (bottom row) compared with a single-species estimate for each taxon (top row; see Table 1 for plotting code for taxa; inset colour bar shows average log-density in kg km-2 ⁠, and colours are defined identically for single- and multispecies models for each taxon). Figure 1 Open in new tabDownload slide Comparison of estimated density functions (averaged across all years for each taxon) for structure-forming invertebrates, six Sebastes, and two Sebastolobus species using the VAST multispecies model (bottom row) compared with a single-species estimate for each taxon (top row; see Table 1 for plotting code for taxa; inset colour bar shows average log-density in kg km-2 ⁠, and colours are defined identically for single- and multispecies models for each taxon). Our species classifications are supported by the estimated covariance matrices (Figure 2), where longspine and shortspine thornyheads have high pairwise correlations in both encounter probability and positive catch rates (0.5–0.7). Encounter probability of longspine is negatively associated with encounter probability of Sebastes spp., while encounter probability of shortspine has both positive and negative associations with different Sebastes spp. This difference between thornyhead species is also apparent in distribution maps (Figure 1), where longspine has the deepest distribution of any fish in our analysis, whereas shortspine occupies a more shoreward distribution that overlaps with the spatial distribution for several Sebastes (e.g. darkblotched, POP, and splitnose). The Sebastes spp. all generally have high correlations (0.5–0.9) with one another for encounter probability (⁠ Ep ⁠; Figure 2 top panel), but the northern vs. coastwide groups are strongly distinguished by correlations in positive catch rates (⁠ Er ⁠; Figure 2 bottom panel), where darkblotched, POP, and sharpchin have higher correlations with one another (0.9) than with splitnose, stripetail, and greenspotted (0.2–0.7). At this coastwide spatial scale, Sebastolobus spp. generally have increased encounter probability when structure-forming invertebrates are found, whereas coastwide and northern rockfish groups have somewhat decreased encounter probability in these cases. When fishes and SFI are encountered, however, an increased catch of SFI is associated with increased catch for all fishes except stripetail and sharpchin rockfish. Figure 2 Open in new tabDownload slide Analytic estimate of correlation among fishes and SFI using the VAST model (see Table 1 for taxa abbreviations) for encounter probabilities (⁠ Rεr ⁠, top panel) or positive catch rates (⁠ Rεp ⁠, bottom panel). Numbered columns correspond to the species groups indicated by the row labels, ordered from top to bottom. Figure 2 Open in new tabDownload slide Analytic estimate of correlation among fishes and SFI using the VAST model (see Table 1 for taxa abbreviations) for encounter probabilities (⁠ Rεr ⁠, top panel) or positive catch rates (⁠ Rεp ⁠, bottom panel). Numbered columns correspond to the species groups indicated by the row labels, ordered from top to bottom. We next compare estimates of biomass trends using multispecies and single-species estimates (Figure 3). Biomass trends are broadly similar between models, and particularly for SFI, which shows a trend of increased biomass since 2008. Biomass trend estimates are most different between multi- and single-species models for the group of northern Sebastes spp. (darkblotched, sharpchin, and POP; Figure 3 middle row). For example, the multispecies model estimates lower abundance for POP in 2008 than the single-species model. This lower estimate for POP in 2008 using the multispecies model reflects a similar decrease in abundance for darkblotched rockfish in 2008 using either model – the estimate for POP in this year for the multispecies model is “shrunk” towards the estimate for darkblotched rockfish. Figure 3 Open in new tabDownload slide Relative log-biomass I(c,t) in units of log-metric tons (Equation 8) for each species using single-species (grey) or multispecies (red) modelling (note different y-axis ranges for each species and see Table 1 for taxa abbreviations; top row: SFI and Sebastolobus group; middle row: northern Sebastes group; bottom row: coastwide Sebastes group). Figure 3 Open in new tabDownload slide Relative log-biomass I(c,t) in units of log-metric tons (Equation 8) for each species using single-species (grey) or multispecies (red) modelling (note different y-axis ranges for each species and see Table 1 for taxa abbreviations; top row: SFI and Sebastolobus group; middle row: northern Sebastes group; bottom row: coastwide Sebastes group). Estimates of variation and trends in center-of-gravity (COG) are also generally similar between multi- and single-species model outputs (Figure 4). The notable exceptions are again the northern Sebastes spp., specifically POP and sharpchin rockfish, which both have relatively few encounters relative to other species (ca. 500 each, see Table 1). For POP and sharpchin, the single-species estimates of COG are nearly 100 km farther south than COG estimates from the multispecies model (Figure 4 middle row). By sharing information about positive catch rates (as shown in the lower panel of Figure 2), the multispecies model estimates greater variation in density for these species between different locations off Oregon and Washington and, therefore, estimates a more northward distribution than the single-species model for POP and sharpchin (Figure 1 top row). We again interpret this as a consequence of statistical “shrinkage” for these species, where the multispecies model is sharing information among northern Sebastes spp. to infer density hotspots. Figure 4 Open in new tabDownload slide Northward center-of-gravity (COG), X(c,t) ⁠, in units of kilometres north of the equator (Equation 9) for each species using single-species (grey) or multispecies (red) modelling (see Table 1 for taxa abbreviations, and Figure 3 caption for details). Figure 4 Open in new tabDownload slide Northward center-of-gravity (COG), X(c,t) ⁠, in units of kilometres north of the equator (Equation 9) for each species using single-species (grey) or multispecies (red) modelling (see Table 1 for taxa abbreviations, and Figure 3 caption for details). Finally, a comparison of standard errors (Figure 5) shows that the multispecies model generally has slightly larger standard error for estimating log-biomass (median 0.01 increase relative to single-species model) and center-of-gravity (median 2.4 km increase). This increased standard error presumably occurs because the multispecies model estimates greater spatial variation in density (Figure 5). For POP, for example, the single-species model estimates little spatial pattern except an increase in density moving northward along the coast, while the multispecies model estimates density hotspots in the same mid-depth areas off the Washington coast as it estimates as good habitat for splitnose and darkblotched rockfishes (Figure 1). Figure 5 Open in new tabDownload slide Comparison of standard error estimates for log-index of abundance, ln⁡[Ic,t] in units of log-metric tonnes (Equation 8, top row) and centre-of-gravity (COG) X(c,t) in units of kilometres north of the equator (Equation 9, bottom row) from single-species and multispecies VAST models, where the scatterplots compare standard error estimates for each year t and taxon c (first column; x-axis shows single-species and y-axis shows multispecies standard errors) and the histograms show the difference between these standard error estimates (second column; where a positive value indicates that the multispecies model had a wider confidence interval than the single-species model for that taxon and year). The dotted line in each histogram indicates the median increase in log-standard error (top row) or standard error (bottom row) for the multispecies relative to the single-species model, and the number in the top-right corner indicates the median increase. Figure 5 Open in new tabDownload slide Comparison of standard error estimates for log-index of abundance, ln⁡[Ic,t] in units of log-metric tonnes (Equation 8, top row) and centre-of-gravity (COG) X(c,t) in units of kilometres north of the equator (Equation 9, bottom row) from single-species and multispecies VAST models, where the scatterplots compare standard error estimates for each year t and taxon c (first column; x-axis shows single-species and y-axis shows multispecies standard errors) and the histograms show the difference between these standard error estimates (second column; where a positive value indicates that the multispecies model had a wider confidence interval than the single-species model for that taxon and year). The dotted line in each histogram indicates the median increase in log-standard error (top row) or standard error (bottom row) for the multispecies relative to the single-species model, and the number in the top-right corner indicates the median increase. Despite estimating wider standard errors for abundance indices and distribution shifts, the multispecies model provides a more fit to available data. The multispecies model has an AIC score that is 5692.0 better than the combined AIC for single-species models, despite the multispecies model estimating an additional 112 parameters (409 fixed effects for the multispecies vs. 297 total among all single-species models). The improvement in fit for the multispecies model is also supported by the tenfold cross-validation analysis, where the multispecies model has a 4–5% greater predictive probability than when analysing each species individually (Table 3). This improvement in predictive score presumably arises because the multispecies model identifies fine-scale differences in species density for all taxa (e.g. comparing Fig. 1 top and bottom rows for splitnose), and these fine-scale density estimates are on average a useful prediction of variation in catch rates. Table 3 Predictive negative log-likelihood for left-out samples (where a low number indicates better fit for the multispecies model) from a tenfold cross-validation experiment comparing single-species models to a multispecies VAST model that was estimated for all species simultaneously, as well as the ratio of predictive probability for the multispecies model relative to the single-species model (a value >1.0 indicates better predictive performance for the multispecies model than the single-species model). Partition . Number of cross-validation samples . Predictive negative log-likelihood . Ratio of predictive probability . Single-species VAST model . Multispecies VAST model . 1 6889 8459.71 8078.84 1.057 2 6832 8560.26 8209.90 1.053 3 6835 8768.53 8377.45 1.059 4 6890 8203.66 7815.19 1.058 5 6799 8335.34 8009.96 1.049 6 6828 8548.48 8154.91 1.059 7 6800 8426.61 8133.02 1.044 8 6997 8688.89 8335.62 1.052 9 6743 8439.27 8077.01 1.055 10 6859 8594.30 8271.42 1.048 Partition . Number of cross-validation samples . Predictive negative log-likelihood . Ratio of predictive probability . Single-species VAST model . Multispecies VAST model . 1 6889 8459.71 8078.84 1.057 2 6832 8560.26 8209.90 1.053 3 6835 8768.53 8377.45 1.059 4 6890 8203.66 7815.19 1.058 5 6799 8335.34 8009.96 1.049 6 6828 8548.48 8154.91 1.059 7 6800 8426.61 8133.02 1.044 8 6997 8688.89 8335.62 1.052 9 6743 8439.27 8077.01 1.055 10 6859 8594.30 8271.42 1.048 Open in new tab Table 3 Predictive negative log-likelihood for left-out samples (where a low number indicates better fit for the multispecies model) from a tenfold cross-validation experiment comparing single-species models to a multispecies VAST model that was estimated for all species simultaneously, as well as the ratio of predictive probability for the multispecies model relative to the single-species model (a value >1.0 indicates better predictive performance for the multispecies model than the single-species model). Partition . Number of cross-validation samples . Predictive negative log-likelihood . Ratio of predictive probability . Single-species VAST model . Multispecies VAST model . 1 6889 8459.71 8078.84 1.057 2 6832 8560.26 8209.90 1.053 3 6835 8768.53 8377.45 1.059 4 6890 8203.66 7815.19 1.058 5 6799 8335.34 8009.96 1.049 6 6828 8548.48 8154.91 1.059 7 6800 8426.61 8133.02 1.044 8 6997 8688.89 8335.62 1.052 9 6743 8439.27 8077.01 1.055 10 6859 8594.30 8271.42 1.048 Partition . Number of cross-validation samples . Predictive negative log-likelihood . Ratio of predictive probability . Single-species VAST model . Multispecies VAST model . 1 6889 8459.71 8078.84 1.057 2 6832 8560.26 8209.90 1.053 3 6835 8768.53 8377.45 1.059 4 6890 8203.66 7815.19 1.058 5 6799 8335.34 8009.96 1.049 6 6828 8548.48 8154.91 1.059 7 6800 8426.61 8133.02 1.044 8 6997 8688.89 8335.62 1.052 9 6743 8439.27 8077.01 1.055 10 6859 8594.30 8271.42 1.048 Open in new tab Discussion We have used a JDSDM to illustrate strong associations (both positive and negative) between deep-water demersal fishes and structure-forming invertebrates at broad spatial scales along the US portion of the California Current. These associations vary substantially between two genera Sebastolobus (thornyheads) and Sebastes, where Sebastes can be further divided into northern and coastwide species. Previous work has shown phylogenetic signals in covariation among fishes (Thorson et al., 2015a, 2016a,b) or other species (Ovaskainen et al., 2010), but ours is the first study to (i) use a spatio-temporal statistical model to estimate covariance between fishes and structure-forming invertebrates, and (ii) decompose this covariation into components representing encounter probabilities vs. positive catch rates [i.e. using the delta models that are conventional in fisheries science; Maunder and Punt, 2004]. Although the JDSDM was more parsimonious and had better predictive performance than single-species models (as shown by AIC and cross-validation analysis), the multispecies analysis resulted in slightly wider confidence interval estimates than analysing data for each species individually. At a coastwide spatial scale, we estimate an increased encounter probability for Sebastolobus and a decreased encounter probability for Sebastes species where SFIs are present. In contrast, alternative visual sampling at fine spatial scales often shows a large increase in Sebastes density, given the presence of SFIs, and Sebastolobus densities are less often reported to be associated with biogenic habitat (Brodeur, 2001; Tissot et al., 2008; Yoklavich and O’Connell, 2008; du Preez and Tunnicliffe, 2011). Recent research suggests that correlations in distribution among species will often differ when looking at small and large scales (Ovaskainen et al., 2016), and this may explain why our results differ from those from fine-scale visual sampling. Alternatively, differences in results may arise because visual sampling often occurs in rocky habitats, whereas our analysis relies on bottom trawl data that are primarily available in soft-sediment habitats. We recommend future research combining data from small and coastwide scales (and both hard- and soft-bottom habitat) within a single spatio-temporal statistical model, where density-variation at fine scales could be obtained by either fishery-dependent catch-rate data or direct observations (Jagielo et al., 2003; Shelton et al., 2014; Rooper et al., 2016; Thorson et al., 2016). We also recommend future research to include habitat variables and associations within size-structured spatio-temporal models (e.g. Kristensen et al., 2014; Nielsen et al., 2014). These models could then estimate separate habitat associations for juvenile and adult fishes and be used to target spatial management towards the more vulnerable or sensitive life stage for protected species. Based on our results, we find that simultaneously modelling fishes and SFI yields more parsimonious predictions of density and also facilitates estimating variation in density at finer spatial scales than single-species models, even for species with few encounters (e.g. POP and sharpchin rockfishes). However, incorporating these associations when estimating trends in abundance or distribution does not shrink confidence intervals. For an ecologist conducting a stock assessment, incorporating multispecies data may complicate their description of estimated abundance indices, thereby decreasing stakeholder trust in the stock assessment process. Therefore, we imagine that our results will encourage many assessment scientists to continue using single-species models for estimating abundance indices. From a broader perspective, however, the increased parsimony and out-of-sample predictive ability of the multispecies model indicate that estimates of local density are generally improved by jointly modelling multiple species (including both fishes and biogenic habitat). Precise predictions of local density for rare species might be particularly useful for ecosystem modellers, who often initialize spatial ecosystem models using sparse sampling data for rare species or ecosystem components. These estimates of local density could also be used to prioritize areas for spatial management that have a high density of structure-forming invertebrates and fishes. Therefore, we suggest further research regarding the association of fished species and biogenic habitat, including the likely impact of spatial management on fishery productivity in the West Coast groundfish fishery. Acknowledgements We thank the NWFSC FRAM Fisheries Research Survey Team and the crew on the US West Coast Groundfish Bottom Trawl Survey for collecting the data, and we thank Michelle McClure, Trevor Branch, and Tim Essington for helpful comments and discussion. We also thank Mary Yoklavich, Joe Bizzarro, Chris Rooper, and two anonymous reviewers for comments on an earlier draft. Funding L.A.K.B. gratefully acknowledges funding from the Joint Institute for the Study of the Atmosphere and Ocean (JISAO) under NOAA Cooperative Agreement No. NA15OAR4320063, Contribution No. 2016-01-38. References Akaike H. 1974 . New look at statistical-model identification . IEEE Transactions on Automatic Control , AC19 : 716 – 723 . Google Scholar Crossref Search ADS WorldCat Bradburn M. J. , Keller A. A., Horness B. H. 2011 . The 2003 to 2008 US West Coast bottom trawl surveys of groundfish resources off Washington, Oregon, and California: Estimates of distribution, abundance, length, and age composition . NOAA Technical Memorandum, NMFS-NWFSC-114 . 323 pp. Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Brodeur R. D. 2001 . Habitat-specific distribution of Pacific ocean perch (Sebastes alutus) in Pribilof Canyon, Bering Sea . Continental Shelf Research , 21 : 207 – 224 . Google Scholar Crossref Search ADS WorldCat Burnham K. P. , Anderson D. 2002 . Model Selection and Multimodel Inference: A Practical Information-Theoretic Approach , 2nd edn. Springer , New York . 488 pp. Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Clark J. S. , Gelfand A. E., Woodall C. W., Zhu K. 2014 . More than the sum of the parts: forest climate response from joint species distribution models . Ecological Applications , 24 : pp. 990 – 999 . Google Scholar Crossref Search ADS PubMed WorldCat Francis R. I. C. C. 2011 . Data weighting in statistical fisheries stock assessment models . Canadian Journal of Fisheries and Aquatic Sciences , 68 : 1124 – 1138 . Google Scholar Crossref Search ADS WorldCat Helser T. E. , Punt A. E., Methot R. D. 2004 . A generalized linear mixed model analysis of a multi-vessel fishery resource survey . Fisheries Research , 70 : 251 – 264 . Google Scholar Crossref Search ADS WorldCat Hyde J. R. , Vetter R. D. 2007 . The origin, evolution, and diversification of rockfishes of the genus Sebastes (Cuvier) . Molecular Phylogenetics and Evolution , 44 : 790 – 811 . Google Scholar Crossref Search ADS PubMed WorldCat Ingram T. , Shurin J. B. 2009 . Trait-based assembly and phylogenetic structure in northeast Pacific rockfish assemblages . Ecology , 90 : 2444 – 2453 . Google Scholar Crossref Search ADS PubMed WorldCat Jagielo T. , Hoffmann A., Tagart J., Zimmermann M. 2003 . Demersal groundfish densities in trawlable and untrawlable habitats off Washington: Implications for the estimation of habitat bias in trawl surveys . Fishery Bulletin , 101 : 545 – 565 . Google Scholar OpenURL Placeholder Text WorldCat Kristensen K. , Nielsen A., Berg C. W., Skaug H., Bell B. M. 2016 . TMB: Automatic differentiation and Laplace approximation . Journal of Statistical Software , 70 : 1 – 21 . Google Scholar Crossref Search ADS WorldCat Kristensen K. , Thygesen U. H., Andersen K. H., Beyer J. E. 2014 . Estimating spatio-temporal dynamics of size-structured populations . Canadian Journal of Fisheries and Aquatic Sciences , 71 : 326 – 336 . Google Scholar Crossref Search ADS WorldCat Latimer A. M. , Banerjee S., Sang H. Jr, Mosher E. S., Silander J. A. Jr, 2009 . Hierarchical models facilitate spatial analysis of large data sets: a case study on invasive plant species in the northeastern United States . Ecology Letters , 12 : 144 – 154 . Google Scholar Crossref Search ADS PubMed WorldCat Lindgren F. 2012 . Continuous domain spatial models in R-INLA . ISBA Bulletin , 19 : 14 – 20 . Google Scholar OpenURL Placeholder Text WorldCat Lindgren F. , Rue H., Lindström J. 2011 . An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach . Journal of the Royal Statistical Society: Series B (Statistical Methodology) , 73 : 423 – 498 . Google Scholar Crossref Search ADS WorldCat Love M. S. , Yoklavich M. M., Thorsteinson L. K. 2002 . The Rockfishes of the Northeast Pacific . University of California Press , Berkeley . 416 pp. Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Mangel M. , Kindsvater H. K., Bonsall M. B. 2007 . Evolutionary analysis of life span, competition, and adaptive radiation, motivated by the Pacific rockfishes (Sebastes) . Evolution , 61 : 1208 – 1224 . Google Scholar Crossref Search ADS PubMed WorldCat Martin T. G. , Wintle B. A., Rhodes J. R., Kuhnert P. M., Field S. A., Low-Choy S. J., Tyre A. J. et al. 2005 . Zero tolerance ecology: Improving ecological inference by modelling the source of zero observations . Ecology Letters , 8 : 1235 – 1246 . Google Scholar Crossref Search ADS PubMed WorldCat Maunder M. N. , Punt A. E. 2004 . Standardizing catch and effort data: a review of recent approaches . Fisheries Research , 70 : 141 – 159 . Google Scholar Crossref Search ADS WorldCat Maunder M. N. , Punt A. E. 2013 . A review of integrated analysis in fisheries stock assessment . Fisheries Research , 142 : 61 – 74 . Google Scholar Crossref Search ADS WorldCat Methot R. D. 2009 . Stock assessment: Operational models in support of fisheries management. In The Future of Fisheries Science in North America , pp. 137 – 165 . Ed. by Beamish R. J., Rothschild B. J.. Springer , Netherlands, Dordrecht . Google Scholar Crossref Search ADS Google Preview WorldCat COPAC Nielsen J. R. , Kristensen K., Lewy P., Bastardie F. 2014 . A statistical model for estimation of fish density including correlation in size, space, time and between species from research survey data . PLoS One , 9 : e99151. Google Scholar Crossref Search ADS PubMed WorldCat Ovaskainen O. , Abrego N., Halme P., Dunson D. 2016 . Using latent variable models to identify large networks of species-to-species associations at different spatial scales . Methods in Ecology and Evolution , 7 : 549 – 555 . Google Scholar Crossref Search ADS WorldCat Ovaskainen O. , Hottola J., Siitonen J. 2010 . Modeling species co-occurrence by multivariate logistic regression generates new hypotheses on fungal interactions . Ecology , 91 : 2514 – 2521 . Google Scholar Crossref Search ADS PubMed WorldCat Perry A. L. , Low P. J., Ellis J. R., Reynolds J. D. 2005 . Climate change and distribution shifts in marine fishes . Science , 308 : 1912 – 1915 . Google Scholar Crossref Search ADS PubMed WorldCat Pinsky M. L. , Worm B., Fogarty M. J., Sarmiento J. L., Levin S. A. 2013 . Marine taxa track local climate velocities . Science , 341 : 1239 – 1242 . Google Scholar Crossref Search ADS PubMed WorldCat du Preez C. , Tunnicliffe V. 2011 . Shortspine thornyhead and rockfish (Scorpaenidae) distribution in response to substratum, biogenic structures and trawling . Marine Ecology Progress Series , 425 : 217 – 231 . Google Scholar Crossref Search ADS WorldCat R Core Team . 2015 . R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. Available from https://www.R-project.org/. Rooper C. N. , Sigler M. F., Goddard P., Malecha P., Towler R., Williams K., Wilborn R. et al. 2016 . Validation and improvement of species distribution models for structure-forming invertebrates in the eastern Bering Sea with an independent survey . Marine Ecology Progress Series , 551 : 117 – 130 . Google Scholar Crossref Search ADS WorldCat Shelton A. O. , Thorson J. T., Ward E. J., Feist B. E. 2014 . Spatial semiparametric models improve estimates of species abundance and distribution . Canadian Journal of Fisheries and Aquatic Sciences , 71 : 1655 – 1666 . Google Scholar Crossref Search ADS WorldCat Simpson D. , Lindgren F., Rue H. 2012 . In order to make spatial statistics computationally feasible, we need to forget about the covariance function . Environmetrics , 23 : 65 – 74 . Google Scholar Crossref Search ADS WorldCat Skaug H. , Fournier D. 2006 . Automatic approximation of the marginal likelihood in non-Gaussian hierarchical models . Computational Statistics & Data Analysis , 51 : 699 – 709 . Google Scholar Crossref Search ADS WorldCat Thorson J. T. , Fonner R., Haltuch M., Ono K., Winker H. 2016 . Accounting for spatiotemporal variation and fisher targeting when estimating abundance from multispecies fishery data . Canadian Journal of Fisheries and Aquatic Sciences , doi: 10.1139/cjfas-2015-0598. Google Scholar OpenURL Placeholder Text WorldCat Thorson J. T. , Ianelli J. N., Larsen E. A., Ries L., Scheuerell M. D., Szuwalski C., Zipkin E. F. 2016a . Joint dynamic species distribution models: a tool for community ordination and spatio-temporal monitoring . Global Ecology and Biogeography , 25 : 1144 – 1158 . Google Scholar Crossref Search ADS WorldCat Thorson J. T. , Minto C. 2015 . Mixed effects: a unifying framework for statistical modelling in fisheries biology . ICES Journal of Marine Science , 72 : 1245 – 1256 . Google Scholar Crossref Search ADS WorldCat Thorson J. T. , Pinsky M. L., Ward E. J. 2016b . Model-based inference for estimating shifts in species distribution, area occupied and centre of gravity . Methods in Ecology and Evolution , 7 : 990 – 1002 . Google Scholar Crossref Search ADS WorldCat Thorson J. T. , Scheuerell M. D., Shelton A. O., See K. E., Skaug H. J., Kristensen K. 2015a . Spatial factor analysis: a new tool for estimating joint species distributions and correlations in species range . Methods in Ecology and Evolution , 6 : 627 – 637 . Google Scholar Crossref Search ADS WorldCat Thorson J. T. , Shelton A. O., Ward E. J., Skaug H. J. 2015b . Geostatistical delta-generalized linear mixed models improve precision for estimated abundance indices for West Coast groundfishes . ICES Journal of Marine Science , 72 : 1297 – 1310 . Google Scholar Crossref Search ADS WorldCat Thorson J. T. , Stewart I., Punt A. 2011 . Accounting for fish shoals in single- and multi-species survey data using mixture distribution models . Canadian Journal of Fisheries and Aquatic Sciences , 68 : 1681 – 1693 . Google Scholar Crossref Search ADS WorldCat Thorson J. T. , Ward E. J. 2014 . Accounting for vessel effects when standardizing catch rates from cooperative surveys . Fisheries Research , 155 : 168 – 176 . Google Scholar Crossref Search ADS WorldCat Tissot B. N. , Wakefield W. W., Hixon M. A., Clemons J. E. 2008 . Twenty years of fish-habitat studies on Heceta Bank, Oregon. In Marine Habitat Mapping Technology for Alaska , pp. 203 – 217 . Ed. by Reynolds J. R., Greene H. G.. Alaska Sea Grant College Program, University of Alaska Fairbanks . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Tobler W. R. 1970 . A computer movie simulating urban growth in the Detroit region . Economic Geography , 46 : 234 – 240 . Google Scholar Crossref Search ADS WorldCat Walters C. 2003 . Folly and fantasy in the analysis of spatial catch rate data . Canadian Journal of Fisheries and Aquatic Sciences , 60 : 1433 – 1436 . Google Scholar Crossref Search ADS WorldCat Warton D. I. , Blanchet F. G., O'Hara R. B., Ovaskainen O., Taskinen S., Walker S. C., Hui F. K. 2015 . So many variables: Joint modeling in community ecology . Trends in Ecology & Evolution , 30 : 766 – 779 . Google Scholar Crossref Search ADS PubMed WorldCat Yoklavich M. M. , O’Connell V. 2008 . Twenty years of research on demersal communities using the Delta submersible in the Northeast Pacific. In Marine Habitat Mapping Technology for Alaska , pp. 143 – 155 . Ed. by Reynolds J. R., Greene H. G.. Alaska Sea Grant College Program, University of Alaska Fairbanks . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Zipkin E. F. , Andrew Royle J., Dawson D. K., Bates S. 2010 . Multi-species occurrence models to evaluate the effects of conservation and management actions . Biological Conservation , 143 : 479 – 484 . Google Scholar Crossref Search ADS WorldCat Published by International Council for the Exploration of the Sea 2017. This work is written by US Government employees and is in the public domain in the US. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model) Published by International Council for the Exploration of the Sea 2017. This work is written by US Government employees and is in the public domain in the US.
TI - Comparing estimates of abundance trends and distribution shifts using single- and multispecies models of fishes and biogenic habitat
JF - ICES Journal of Marine Science: Journal du Conseil
DO - 10.1093/icesjms/fsw193
DA - 2017-05-01
UR - https://www.deepdyve.com/lp/oxford-university-press/comparing-estimates-of-abundance-trends-and-distribution-shifts-using-mBufxXQsem
SP - 1311
EP - 1321
VL - 74
IS - 5
DP - DeepDyve
ER -