Improving estimates of species distribution change by incorporating local trends

flounder Atheresthes stomias as an example, the observed southward shift over time in the center of gravity is not reflective of a uniform shift in densities but local trends of decreasing density in the northern region and rapidly increasing density at the southern edge of the species’ range. Thus, estimating local trends with spatiotemporal models improves interpretation of species distribution change.


Introduction
In the fields of natural resource conservation, management and global change biology, demand for and implementation of tools for assessing species distribution shifts has grown dramatically in recent decades (Elith and Leathwick 2009).These approaches are widely applicable, from studies of plants (Lenoir et al. 2008), terrestrial vertebrates (Hitch and Leberg 2007) and marine fishes (Pinsky et al. 2013).However, the way distribution shifts are quantified has changed relatively little (Elith et al. 2010).At the simplest level, researchers often use existing tools to estimate occurrence probability, present maps of how the extent and distribution of suitable habitat is expected to change, and sometimes present descriptive statistics on the mean change throughout a region (Yackulic et al. 2013).Yet, when reliable abundance data are available, distribution shifts are more robustly quantified by spatial predictions of population size because these are a richer form of data that are less sensitive to detection issues and anomalous observations of single individuals (Tingley and Beissinger 2009), while being more likely to detect persistent distribution shifts caused by more nuanced factors than absolute physiological limits.For example, while much species distribution modeling focuses on how drivers such as climate change may predict change in species' range limits, the core of a species' distribution may shift due to the influence of multiple drivers on the geography of abundance via movement, dispersal and heterogeneity in demographic rates (e.g.age-or size-specific fecundity, somatic growth and mortality; Sagarin et al. 2006).Shifting distributions of abundance or population density may also be qualitatively conveyed through maps, but quantitative spatial indicators can also be provided, such as the mean location weighted by population density (also termed the 'center of gravity', COG; Thorson et al. 2015).
Spatial distributions of population density are often complex and heterogeneous (Sagarin andGaines 2002, Sagarin et al. 2006).Heterogeneity may be present in the distribution of a species throughout its range, but the change in a species' population density over time may also have a spatially varying component.Mechanisms responsible for spatial variability in change might be biological (e.g.variation in birth and death rates) or forced by spatially structured pressures (such as environmental or habitat change, and anthropogenic disturbances, Barnett et al. 2019 and references therein).Consequently, attempting to describe a uniform shift in distribution across a broad geographic range can be misleading (Sagarin et al. 2006), particularly when different regions exhibit contrasting trends.For example, if densities increase at opposing range boundaries at an equivalent rate, there may be no trend in the range-wide COG, masking finer-scale shifts.Thus, when using spatial indicators to describe species distribution shifts, the spatial scale of aggregation can affect inference (Connor et al. 2019), as in the classic problem of pattern and scale in ecology (Levin 1992).Therefore, there is a general need to develop objective methods for defining appropriate scales to evaluate changes in species distributions.
Techniques for estimating how populations vary over space and time evolved rapidly with increases in computational power and the development of novel methods and applications of tools such as hierarchical statistical models.Some of the largest methodological changes have been advances in spatiotemporal analyses that model space continuously and explicitly account for spatial autocorrelation between spatially referenced observations that are proximate in both space and time (Banerjee et al. 2008, Finley et al. 2009, Latimer et al. 2009, Cressie and Wikle 2011, Shelton et al. 2014, Thorson et al. 2015).There are a number of advantages of estimating a species' density in a framework that accounts for spatial or spatiotemporal variation.First, explicitly accounting for spatial variation in density has been shown to increase precision of estimated temporal trends (Thorson et al. 2015).Second, modeling spatial or spatiotemporal variation in population density can be performed within flexible and established frameworks such as mixedeffect models where the spatial or spatiotemporal components are estimated as random effects (Latimer et al. 2009, Cressie and Wikle 2011, Shelton et al. 2014, Finley et al. 2015).Similar to their non-spatial predecessors, these spatiotemporal modeling approaches often treat time as a discrete factor to allow for unbiased estimates of temporal trends.With such approaches, the spatial distribution of density can either be constant (modeled as a single spatial field) or timevarying (with variability modeled either as independent over time, or as an autoregressive process).
These, and similar spatial model-based estimators, have in many applications replaced conventional design-based estimators of population density (e.g.area-weighted mean of observations), which assume that density is homogenous within sampling strata (Chen et al. 2004).In addition to being used for estimating population density or spatial distributions, output from these modeling approaches has been used to generate model-based summaries to track change in species distributions, including the COG or area occupied, with more robust estimation than that provided by design-based estimates (Thorson et al. 2016a).As tools to implement these methods have become accessible in open source software, such as INLA (Rue et al. 2009), VAST (Thorson 2019b) or sdm-TMB (Anderson et al. 2019(Anderson et al. , 2020)), these approaches have seen broad application to populations in diverse ecosystems around the world, including terrestrial plants (Banerjee et al. 2008, Finley et al. 2009, Latimer et al. 2009) and animals (Thorson et al. 2016b), freshwater (Hocking et al. 2018) and marine communities (Shelton et al. 2014, Thorson et al. 2015, 2016a, Thorson and Barnett 2017, Anderson and Ward 2019).
Spatiotemporal modeling of population density has particularly flourished in the field of marine fisheries, where practical constraints limit the use of spatial capture-recapture, distance sampling and other survey methods typically applied in terrestrial and freshwater systems that have led to parallel development of models for similar purposes (Efford 2004, 2011, Royle and Young 2008, Royle et al. 2013).The most reliable estimates of marine fish densities at broad scales are generally derived from fishery-independent surveys where observations of population density are taken to be proportional to the catch-per-unit-effort of a fishing gear, often implemented using some form of stratified random sampling design.In addition to providing the data on population size and structure that is needed for managing individual fish populations, fishery-independent survey data are used for purposes such as informing ecosystem-based management (Link et al. 2002, Nicholson and Jennings 2004, Harvey et al. 2018), evaluating the impacts of harvesting non-target species (Stock et al. 2019), and quantifying shifts in species distributions (Pinsky et al. 2013, Thorson et al. 2016a).Changes in the spatial distribution of marine fishes have significant implications for community structure and national food security (Rice and Garcia 2011).Thus, robust predictions of marine fish distribution shifts are needed, yet difficult to obtain given sparse and often uneven sampling effort, short monitoring time series, limited ability to repeatsample individuals, and the complexities of surveying open populations with rich spatial structure.Spatial heterogeneity is particularly strong in these marine ecosystems, where complex coastline, bathymetric topography and geology interact with physical oceanographic drivers (Levin et al. 2010).Therefore, it is critical to determine to what degree such heterogeneity must be accounted for to adequately characterize marine fish distribution shifts.
Here, we introduce a new approach to address how estimates of change in species distributions are dependent on the spatial scale of quantitative indicators of species distribution.We describe the development of a modeling technique that explicitly accounts for spatial variability in how population densities change through time to estimate finer-scale indicators of species distribution shifts (local trends).While widely applicable to a wide range of data types (presence-absence, discrete counts or continuous measures such as biomass density), we focus on an application to changes in the distribution of commercially fished marine species.These represent 19 species from a 15-yr publicly available trawl survey dataset.We illustrate how our new approach may be used to infer changes over time, and also how output from this modeling approach may be useful in identifying spatial regions where change is greater or lesser than average.Specifically, we compare interpretations of species distribution shifts along a spectrum of indicators from coarse-scales (global COG trends calculated over an entire survey domain), to moderate-(regional COG trends) and fine-scales (local trend).

Spatial GLMM overview
The majority of recent applications of species distribution models (SDMs) to marine fish survey data have been implemented in a GLMM (generalized linear mixed-effects model) framework, where random effects are used to describe spatial or spatiotemporal components.Spatial components are differentiated from spatiotemporal components in that the former are constant, whereas the latter vary through time.Examples include applications to Gaussian predictive process models (Shelton et al. 2014, Anderson andWard 2019), and predictive modeling using integrated nested Laplace approximations (INLA; Rue et al. 2009, Ruiz-Cárdenas et al. 2012, Thorson et al. 2015).The latter approach has been particularly useful for large datasets, where substantial gains in computational efficiency are accomplished by taking advantage of sparse matrix approximations to the variance-covariance matrix (Thorson and Barnett 2017).Regardless of the estimation approach used, the general formulation of these models uses a link function g(•) to relate the observed response to covariates and a latent spatial process.For example, where u s,t is the expectation at location s and time t, X s,t are covariates, b represents a vector of estimated coefficients, w s is the mean spatial component at location s (constant through time) and ε s,t is the spatiotemporal process at location s and time t.The spatiotemporal process describing ε is flexible in that it can be removed from the model (leaving a model with a spatial but no spatiotemporal component), may be independent for each time slice, or modeled with an autoregressive process (allowing hotspots to persist through time; Ward et al. 2015, Thorson et al. 2015, Anderson and Ward 2019).Previous applications to marine fishes have either used a delta-GLMM framework to model presence-absence and population density (e.g.biomass per unit area) when present separately (Thorson et al. 2015) or a Tweedie distribution to model total variation in density (Anderson et al. 2019).Within this GLMM framework, non-stationary changes in the spatial predictions through time can only be modeled with inclusion of dynamic covariates, or by modeling spatiotemporal variability as an autoregressive spatial process through time.While inclusion of covariates can improve predictive performance in some cases (Shelton et al. 2014, Johnson et al. 2019), this requires additional data and can introduce new challenges associated with finding the most appropriate form of the covariate effect, thus for generality and simplicity we focus here primarily on a latent variable approach for describing patterns in spatially explicit temporal trends (hereafter local trends) rather than directly inferring their drivers.Estimates of local trends may be derived from spatial and spatiotemporal fields post-hoc; however, such post-hoc estimation results in biases (Supporting information), specifically a low bias caused by partial pooling, which effectively pulls the intercept deviations toward the mean.To explicitly account for non-stationary trends in densities, we extend the above framework to include a trend parameter as an additional spatial random field for the slopes over time (in the simplest case, each value in the field represents the spatially-explicit linear trend of the response over the modeled time period).Extending the model above, this becomes where z s represents the spatially varying temporal trend, or local trend.This local trend field can be thought of as the spatial variability in how a species' density changes through time, which differentiates such trends from time-independent spatiotemporal random fields (Fig. 1).Note that this framework could also be extended to model systems in which most spatially explicit responses are non-linear by either modifying the model structure (e.g.…+ z s t 2 or a new and truly nonlinear term in the strict sense) or by fitting separate models to each stanza during which a linear trend is suspected.

Testing the ability to recover local trends
We conducted a simulation analysis to evaluate our ability to recover an added spatial field representing the true local trend.Given results from previous work on similar classes of models (Finley et al. 2015, Auger-Méthé et al. 2016, Finley and Banerjee 2020), we focused our simulations on understanding how the magnitude of spatiotemporal variation and observation error variation affect our ability to recover the local trend (details in Supporting information, or in the publicly available code used to generate these analyses).We also performed similar sensitivity analyses to verify that the magnitude of spatial variance and local trend would affect our ability to recover the local trend in predictable ways.All simulations were conducted following these five steps: 1) for each evaluated (time-invariant) value of spatiotemporal variation and observation error, we simulated a random spatial field.2) We then simulated a latent spatiotemporal process over 10 time steps, using spatial and spatiotemporal components (modeled as independent from year to year) along with the local trend field.3) To include measurement or observation error, we simulated normally distributed observations from this spatiotemporal process.4) Fit a spatial GLMM to the simulated data and assumed the model structure to be known.5) Compared estimated values of the local trend at the locations of the data with known values to generate statistical summaries (bias [expectation of difference], variance [sample variance of difference] and Pearson correlations between predicted and observed values).For each combination of parameter values, we simulated 100 random datasets.

West coast groundfish application
As This survey represents an ideal case study because it has been used extensively in testing new index standardization methods for stock assessments (Thorson et al. 2015), is publicly available (<www.nwfsc.noaa.gov/data/map>),and has been used to develop coast-wide indicators, including shifts in center of gravity (Thorson et al. 2016a).We selected 19 groundfish species to model in this analysis based on a combination of high commercial landings, market value, conservation concern and prevalence in the survey data (Supporting information).It is important to note that the distributions of many of these species extend farther to the north and south.Therefore, conclusions from these analyses only describe the dynamics of their density distribution within the survey area, and not their entire range.While the method does allow for making predictions outside the spatial domain of the data, we chose not extrapolate in this case because doing so would reduce the precision of our estimates.We fit spatial GLMMs with and without a local trend ( z s ) to each species to evaluate whether the local trend may be appropriate for modeling how these 19 species change over time.We allowed both models to include spatial and spatiotemporal components (independent by year, because preliminary testing indicated that including temporal structure was not typically supported, as the 95% confidence interval around the estimate of the first-order autoregressive correlation parameter included 0), depth modeled as a quadratic effect (Thorson et al. 2015), and year as a factor.Below we describe in detail the full model including the local trend.
Because of the positive continuous nature of the recorded fish densities combined with many zeros, we modeled the response y s,t (catch per unit effort [CPUE] at point in space s and time t) with a Tweedie distribution because it performs well with modeling such data (Tweedie 1984, Dunn and Smyth 2005, Anderson et al. 2019), and a log link: where µ represents the mean, p represents the power parameter and φ represents the dispersion parameter.Thorson et al. 2015).
We approximated the continuous random fields using a triangulated mesh (Lindgren et al. 2011) with vertices at 350 'knots' (a representative set of locations sensu Latimer et al. 2009) as calculated with the INLA R package (Rue et al. 2009).We found the minimum log likelihood using the R nlminb optimization routine with Template Model Builder (TMB; Kristensen et al. 2016) implementing the Laplace approximation to the marginal likelihood.TMB uses the generalized delta-method to calculate standard errors.Specifically, we fit all models in R ver.3.5.3(R Core Team) using the package sdmTMB (Anderson et al. 2019(Anderson et al. , 2020) ) which interfaces automatic differentiation in Template Model Builder (Kristensen et al. 2016) with INLA (Rue et al. 2009).
To compare models with different random effect structures (with and without the local trend field), we used restricted maximum likelihood (REML, Zuur et al. 2009) to generate Akaike's information criterion values for each model (AIC, Akaike 1973).AIC is a relative measure of goodnessof-fit that is penalized by the number of model parameters.
Using AIC as a model screening tool, we found broad support for the inclusion of the local trend for these 19 species, with the trend model generating lower AIC values in 17 of the 19 cases, and AIC scores differing by less than two in the remaining two cases (Supporting information).To verify that AIC was effective at selecting the model most consistent with the data-generating process, we performed parallel contrasts (between models with and without the local trend) using simulated data.
Given the evidence supporting the local trend model as the most parsimonious model, we used this model structure to evaluate changes in species' density distributions over time.
To obtain a smooth surface of predicted density across the footprint of the survey area (Fig. 2), we predicted density at each location on a grid, where covariate values were defined using a composite of depth layers defined by NOAA bathymetry data (<www.ngdc.noaa.gov/mgg/coastal/crm.html>).These data were spatially aggregated using bilinear interpolation to match the resolution of the survey sampling grid (~ 2.8 × 3.7 km), which is the spatial resolution we used for all analyses.We implemented a number of diagnostics using spot checks on these predictions and model fits to further analyze whether a local trend was appropriate (e.g.examining spatial patterns in residuals and the estimated spatiotemporal component).

Using local trends as indicators of change
We compared inferences of changes in species' density distributions obtained from metrics calculated on a spectrum of spatial resolution to demonstrate the utility of understanding fine-scale temporal trends.Quantifying change at multiple spatial scales has implications for the management of marine fishes and has utility as a spatial indicator within the California Current ecosystem.Specifically, we compared the fine-scale interpretation of the local trend and mean predicted density over all years to coarse-scale interpretations of: 1) the local trend, 2) regional COGs and 3) coastwide COG calculated from predicted densities y for each location , where L s is the y coordinate of location s).We use the mean predicted density over time as a benchmark for describing how species' distributions change because this can be interpreted as a 'weight' on the local trend, where the product of the two defines the absolute magnitude of the change in density over time.We evaluated whether local trend estimates from our model can be used to identify discrete areas of change that may reflect stock structure.One approach among many possible options for doing this is to apply post-hoc cluster analyses to model outputs or covariates; for our groundfish application, we used the partitioning around medoids (PAM) algorithm with estimation of the number of clusters to demonstrate a possible framework for boundary detection to be used with other system information to define appropriate spatial scales for summarizing monitoring data (implemented with R packages 'fpc' and 'cluster', Hennig 2019, Maechler et al. 2019).PAM is a robust clustering algorithm that minimizes the sum of Euclidean dissimilarities (root of sum-of-squares of differences) between observations and cluster values (Reynolds et al. 2006, Kaufman andRousseeuw 2009).We used latitude and the predicted local trends as clustering variables given that the majority of the contrast in dynamics along the US west coast is in the latitudinal direction.For other applications, additional metrics could also be included in clustering including longitude, habitat features, environmental covariates or human impacts such as fisheries removals.We chose the number of clusters (constrained between 1 and 10) that maximized the average silhouette width across all predictions for a given species (Kaufman and Rousseeuw 2009).Code and data necessary to replicate all above analyses are included in the repository for this project (<https://github.com/fate-spatialindicators/spatial-trend>).

Simulation testing
Results from our simulation indicated that, as expected, both observation error variation and spatiotemporal variation degraded our ability to estimate the true local trend (Fig. 3).When both variance parameters were small, estimates were precisely estimated and unbiased (i.e. the distribution of differences between estimated and true values was centered at approximately zero, as was its greatest mass; Fig. 3a, d); however, large values of either (observation error σ or spatiotemporal σ ≥ 0.5) limited the ability to recover the trend by increasing the standard deviation (Fig. 3b, e) and decreasing correlation between true and estimated values of ζ s (Fig. 3c, f ).Results of further sensitivity analysis were also as expected (Supporting information), with spatial variation having no effect on local trend estimates, while estimates of the local trend were only poor when the variation of the local trend field was extremely low (i.e. the signal was barely present and obscured by variation in other components, causing low correlation between estimated and true local trends; Supporting information).Furthermore, we found that the correct model (the model including the local trend) was easily distinguished by model selection using AIC except when observation error or spatiotemporal variation was extremely high, or when the local trend variance was extremely low (Supporting information).
Figure 3. Simulation testing the effects of observation error and spatiotemporal variation on the ability to recover the local trend.The symbols θ and q refer to the local trend random effect values at each location ζ s and their estimate, respectively.Accurately recovering the local trend would be indicated by a bias of approximately zero, low standard deviation (SD) and high correlation between true and estimated values of ζ s .Each violin represents the distribution of location by location comparisons from 100 simulations and the dots represent the median value.In all cases, the standard deviation σ of the non-varying parameter is held at 0.01, while σ of the parameter subject to the sensitivity test varies along {0.01, 0.25, 0.5, 0.75}.Note that these results were also computed for σ = 1 (Supporting information), yet are omitted here as they were very similar to results from σ = 0.75.

West coast groundfish application
Predictions of the spatially explicit temporal trend from the local trend model revealed intricate fine-scale spatial structure in rates of change of species in the west coast groundfish community.Our cluster analysis of the estimated local trend and latitude helped to delineate areas with the greatest relative rate of change in density over time.Over 89% of species had non-trivial spatial heterogeneity in local trends (Fig. 4), where Pacific halibut and spotted ratfish were the only species that had similar trends across space.For the majority (58%) of species, at least one of the breaks between local trend clusters (unique columns of points in Fig. 4) occurred at a latitude within 100 km of the two predominant biogeographic breaks in this ecosystem: Point Conception in southern California and Cape Mendocino in northern California (Fig. 4).Furthermore, all species had at least one cluster break falling at latitudes between these two biogeographic boundaries and 79% of species had at least half of their cluster breaks in this area.However, there was variability among species in the precise location of the boundaries of the local trend clusters.Results were similar, yet clusters were often less spatially contiguous (e.g.compare bocaccio between Fig. 4 and Supporting information for an extreme example), when clustering on only the local trend without latitude (Supporting information).Given the general proximity between trend cluster breaks and the established biogeographic boundaries, we chose to evaluate the latitudinal center of gravity (COG) within each biogeographic region (rather than within each species-specific local trend cluster) to compare with the coastwide COG.Early exploration indicated that results were qualitatively similar between these two approaches for defining the unit on which to compute COGs; however, we chose not to compute COGs by local trend cluster over concerns that this could minimize and obscure distribution shifts (e.g. if density is changing uniformly over time within a cluster, one would expect the COG of that cluster to remain relatively constant).
We highlight results for six groundfish species with unique distributional responses (Fig. 5; see Supporting information for results from additional species and for predicted density distributions for all 19 species).Within each of the six species, there was support for 2-3 trends (Fig. 5; second column).Comparison of the local trend predictions and clusters (Fig. 5; first two columns) and the mean density from the full model (Fig. 5; third column) revealed how several unique patterns of regional relationships can contribute to nuanced and difficult to detect broad-scale distributional changes Figure 4. Strip plot showing each unique cluster of latitude and local trend (slope over time) by species.Within a species, each set of points associated with a given cluster are represented by a different column and colored by their deviation from the mean coast-wide trend for a given species.Grey points in the vertical plane represent clusters from which the local trend did not differ from (was within 0.01 of ) the mean coast-wide trend.Horizontal black lines represent approximate positions of known biogeographic breaks: Cape Mendocino, California, in the north; Point Conception, California, in the south (see map in Fig. 2).Horizontal gray shading represents a buffer of 100 km around each biogeographic boundary, which provide a benchmark for the proximity statistics described in the main text.California (south); and between these two landmarks (central).The third column represents the mean density over all years.The fourth column shows the time series of the center of gravity (COG, or latitude weighted by density) with 95% confidence intervals along the same y-axis scale as the maps in adjacent columns.The black line with grey interval represents the COG calculated from predicted densities coastwide, whereas the colored lines represent the COGs for each unique biogeographic region (corresponding to the range of northings indicated by the map panels).Line color represents the proportion of a species' relative biomass in a given region.
including northward, southward and bi-directional (convergent or divergent) density shifts (examples below), in addition to localized offshore shifts (e.g.sablefish off Oregon and northern California -red near shore and blue further offshore).Furthermore, the interpretation of the distributional change often varied between spatial scales of metrics.Typically, inference differed the most between the fine-scale map-based interpretation of the local trend and the coastwide COG.The map of estimated mean density allows one to visually weight the local trend map to better understand where absolute changes in density are greatest.
Examining the predictions of the local trend and density indicated that arrowtooth flounder Atheresthes stomias had a southward density shift and shortspine thornyhead Sebastolobus alascanus had a northward shift, yet the COG inferences differed to some degree between species.The predicted mean density (column 3) indicated that the majority of arrowtooth flounder (Fig. 5, first row) was in the northern region, yet the local trend pattern indicates that their density is increasing at the highest rate in the central region (column 1).Combined, these regional results suggest a southward shift driven by increases at the southern range edge, similar to the traveling wave pattern demonstrated by many species invasions.The time series of the coast-wide COG (black line in column 4, Fig. 5) is in agreement of a southward shift, yet the interpretation is not as clear because the coast-wide pattern is heavily weighted by the high densities in the far northern portion of the study area.When the COG from each biogeographic region is calculated (colored lines in column 4, Fig. 5), we can see that coast-wide COG has been driven further south in the latter half of the time series by decreases in the COG in the central region while the northern COG had almost no trend, providing additional support for the possibility that the change is due to increased density or southward shifts in the central region.
For other species in our analysis, even the region-specific COG does not accurately capture the nuanced spatial changes described by the local trend field.For example, shortspine thornyhead is distributed coast-wide, yet their density is increasing fastest in the north-central area and decreasing in the south and within some isolated patches in the far northern end of the region (Fig. 5, last row, column 1).In this case, the coast-wide COG indicates a northward distribution shift, yet the region-specific COG indicates converging trends, perhaps indicative of contraction of the core range: slightly southward shifting in the northern region and slight northward shift in the central region.Thus, the interpretation from the COGs at both scales are relatively consistent with the fine-scale interpretation of the local trend (Fig. 5 column 1), yet these coarse-scale metrics still mask underlying patterns, in this case the decreased density in the southern region.
Other species demonstrated additional patterns of changes in spatial distribution of density and contrasting inference among metrics, including bocaccio rockfish Sebastes paucispinis, English sole Parophrys vetulus, petrale sole Eopsetta jordani and sablefish Anoplopoma fimbria.Bocaccio were typically more abundant in the southern and central areas (Fig. 5 column 3) yet were experiencing the fastest increases in density in the north (Fig. 5 column 1), indicating a northward density shift.These observations contrast with those from the COG, where the coast-wide COG for bocaccio was highly variable with either no trend or a very slight southward trend, and the COG of the northern region indicates a southward shift in some years (Fig. 5 column 4).Divergent density shifts were observed for English sole and, to some extent, petrale sole.English sole were typically present in relatively similar densities coast-wide, yet the local trend indicated that densities were increasing fastest at the northern and southern ends of the region.However, the coast-wide COG reveals only a slight southward shift, while the region-specific COGs show only a slight northward shift in the northern region.Petrale sole had a complex local trend field, increasing fastest in the north with the exception of isolated declining patches on the inshore side.These changes are somewhat consistent with the coast-wide COG indicating a slight northward trend amidst moderate interannual variability.However, COGs of the northern and central regions -where petrale sole are typically most prevalent -indicated a divergent pattern, in which densities were shifting northward in the northern region and slightly southward in the central region.Finally, no obvious directional shift in density was apparent in the local trend for sablefish, yet the coast-wide COG time series indicated that a northward shift had occurred in the most recent 5-6 yr.The region-specific COG indicates that this was driven by density changes in the northern and to some extent central regions.Thus, depending on the evidence used, one could either conclude that there was a recent northward density shift, or simply an increase in productivity or aggregation near the core of the range within the north-central area.

Discussion
The complex spatial distribution of biotic and abiotic drivers of population productivity and habitat suitability in ecosystems suggests that fine-scale descriptors may provide a more accurate representation of changes in species distributions than global indicators calculated across an entire region.Here, we introduced a new approach to modeling and summarizing spatially referenced time series data on species population densities to calculate area-specific trends in population size.Our approach was able to recover local trends in simulated data and reveal nuanced local trends in the dynamics of 19 marine fishes off the west coast of the USA that often differed from conventional descriptors of larger scale distribution shifts (Woillez et al. 2009, Pinsky et al. 2013, Thorson et al. 2016a).Furthermore, the ability of our models to detect geographic boundaries between regions with different trends was supported as these boundaries were often congruent with known biogeographic breaks (yet we acknowledge that this may be influenced by assumptions affecting the optimization of the number of clusters and other factors that require deeper investigation to strengthen this initial finding).Therefore, boundary detection techniques applied to a local trend field may be valuable for helping to define appropriate spatial scales for summarizing monitoring products (such as abundance time series) or regulating the spatial allocation of human impacts (e.g.allowable take of animals or plants), especially in cases where little other information on spatial population and community structure is available.Furthermore, we note that the local trend model described here has potential for broad applications to a variety of other data types beyond population density (e.g.spatial patterns of temperature variability) and other systems (as an extension of existing applications, Banerjee et al. 2008, Finley et al. 2009, Latimer et al. 2009, Cressie and Wikle 2011, Thorson et al. 2016b, Hocking et al. 2018) Our simulations and application of the local trend model indicate that our proposed approach can improve estimation and communication of spatially varying temporal trends in population density.In particular, our application to marine fish survey data indicated that models including a local trend field were more parsimonious than those without a local trend.This result is consistent with a recent study incorporating a spatially varying influence of an oceanographic index on groundfish distributions in the eastern Bering Sea (Thorson 2019a).Furthermore, according to our simulations the estimated local trends were less biased than those estimated post-hoc from predictions of a model without the local trend field.However, the local trend model is somewhat sensitive to observation error and spatiotemporal variation.Such sources of variation can obscure the local trend, yet this is to be expected in the same way that any trend is less detectable given noisier data (Weatherhead et al. 1998).Therefore, our method is likely most skillful at detecting spatial structure in population or community dynamics from observations with precise measurement within systems with low temporal variation in spatial structure (e.g.those consisting of species with higher longevity, generation time and site fidelity, and lower rates of movement and variation in dispersal paths).We expect that the predictions in our example application in this study are robust to the sensitivity of the method to spatiotemporal variation because the estimated spatiotemporal variance is much lower than the spatial variance for groundfish species in this system.Observation error in trawl surveys can include a wide range of values as a result of variance in sampling efficiency (Kotwicki and Ono 2019), but relating such values to the observation error scale parameter evaluated in our simulations requires additional research.Additional ways to constrain the variance parameters, such as developing informative Bayesian priors from similar surveys might extend the detectability of local trend structure over the models used here.
We showed how inference about shifts in species' population density depend on the spatial scale at which they are summarized.When we applied the local trend model to marine fishes, the resulting maps of the local trend and density from the model revealed nuanced patterns of heterogeneity and directional change in groundfish density.Taking the predicted density to represent the underlying spatial heterogeneity, the local trend random field conveyed finescale information about potential range dynamics that were masked when evaluating coast-wide COG time series.The disparity of inference was greatest in cases where density was increasing fastest at opposing ends of a range, density was spatially diverging, or where density among patches were converging toward the center of the distribution.Furthermore, when examining only the coast-wide COG, one is unable to differentiate between shifts due to an increase in density in one area or a decrease in density in another area.
For complex ecosystems such as the US west coast, and other coastal upwelling systems where physical variables like temperature do not follow a simple monotonic gradient over broad geographic scales, it may be too simplistic to expect clear coast-wide trends in COG across multiple species as a result of climate change.These coast-wide patterns are observed in systems with broader continental shelves such as the northeast US (Pinsky et al. 2013, Kleisner et al. 2016) where the major boundary currents are far from the coast.However, along coastlines with narrower continental margins, such as the US west coast, fish may be able to find equivalent temperatures by moving much shorter distances perpendicular to the shelf break (Li et al. 2019).Furthermore, population and community density distributions are inherently patchy, particularly for species associated with patchy reef habitats, meaning that detecting a redistribution over time may require careful examination of the microstructure of density distribution rather than a region-wide shift in mean density distribution.We encourage future research on species distribution shifts that begins with more specific and nuanced hypotheses regarding the expected response at shorter and perhaps longer time scales than those explored here, as spatially explicit trends are likely to differ between intra-annual, inter-annual and inter-decadal time scales.For example, event-scale analyses of the local trend could help test how different species population density distributions respond as a result of movement or demography to disease outbreaks, intensive harvesting or extreme climate events such as marine heat waves.If climate change causes a global intensification of upwelling over longer time scales as some researchers predict (Bakun et al. 2010), one could hypothesize that density distributions will become patchier over time in response to increasing contrast in local physical conditions, or that distributions will shift deeper as larvae are transported further offshore before settling.
The future of environmental conservation and natural resource management relies on greater incorporation of spatial information into models that inform such policies and into the decision-making process itself (Berger et al. 2017, Lowerre-Barbieri et al. 2019).By defining the geography of population trends and the breaks between clusters of locations with similar trends, our modeling framework provides a data-driven method to objectively define the spatial scale and boundaries for summarizing monitoring data and structuring these inputs to resource management models.This is an important advancement over non-spatial resource assessments or the reliance on the use of jurisdictional boundaries to structure resource assessments.Our vision is that these and subsequent methods for boundary detection will aid the development of spatial resource assessment models and stimulate further applications of such approaches to more disparate management solutions such as invasive species management.Furthermore, extensions of the methods presented here may lead to the creation and improvement of spatial indicators for monitoring factors affecting emergent ecological properties (Barnett et al. 2019 and references therein).Novel indicators of ecological stability could arise from metrics of the spatial structure of temporal trends or oscillations in population density (Kéfi et al. 2014, Walter et al. 2017), by drawing on the evidence that spatial heterogeneity can increase population and community stability by disrupting synchrony across space or among species (Huffaker 1958, Tilman and Kareiva 1997, Hassell 2000).
an example of how the local trend model can be applied to improve the interpretation of changes in spatial distribution, we fit the local trend model to groundfish data collected from a fishery-independent survey along the US west coast: the NOAA Fisheries, Northwest Fisheries Science Center, US West Coast Groundfish Bottom Trawl Survey (Keller et al. 2017) from 2003 to 2018.The annual survey uses a stratified random sampling design, with strata defined by depth and latitude, to estimate population density (in terms of catch per area swept by the net) along the continental shelf and upper slope (from 55 to 1280 m depth) of California, Oregon and Washington state.Roughly 650 tows (the unit of observation) are performed during two passes from north to south, typically occurring between late May and the end of October.

Figure 1 .
Figure 1.Visualization of how the spatial distribution of population density changes over time when the temporal response differs among locations.Predictions are shown from the spatial and temporal random effects of a GLMM with (top row) and without (bottom row) a spatially varying temporal trend (i.e.local trend).Each panel shows a field representing the spatial variation in population density at a given time, where time progresses among panels from left to right.When a spatially varying temporal trend is present, some regions develop systematically higher (red) or lower (blue) density over time.In contrast, when a spatially varying temporal trend is absent, spatial deviations from year-to-year are independent.For this example, we have omitted all other sources of variability and error for simplicity.

s
and time t ( COG t

Figure 2 .
Figure 2. Map of the bathymetry within the US West Coast Groundfish Bottom Trawl Survey footprint.Cape Mendocino and Point Conception are labeled to represent the latitudinal boundaries between known predominant biogeographic regions (N.B. there is also more limited evidence for additional or alternative biogeographic breaks, e.g.Cape Blanco, the westernmost point in Oregon).

Figure 5 .
Figure 5. Spatial and temporal patterns of predicted density for selected species along the US west coast.Panels depict model outputs within the study area (Fig. 2), where the thin black polygon in maps is the land boundary.The first column shows maps of the predicted local trend (slope of log density across years).The second shows how each spatial location groups with a unique cluster of latitude and local trend, where each color represents a cluster.The first panel also shows how the latitudinal boundaries of established biogeographic breaks (thick black horizontal lines) separate the study area into three regions: north of Cape Mendocino, California (north); south of Point Conception,California (south); and between these two landmarks (central).The third column represents the mean density over all years.The fourth column shows the time series of the center of gravity (COG, or latitude weighted by density) with 95% confidence intervals along the same y-axis scale as the maps in adjacent columns.The black line with grey interval represents the COG calculated from predicted densities coastwide, whereas the colored lines represent the COGs for each unique biogeographic region (corresponding to the range of northings indicated by the map panels).Line color represents the proportion of a species' relative biomass in a given region.
(Cressie and Wikle 2011,nts the spatially varying coefficients that represent local trends through time, also drawn from Gaussian Markov random fields.Time, t, is entered into the model for multiplication with z s after centering it by its mean value.All three random fields have covariance matrices constrained by anisotropic Matérn covariance functions with independent scales but shared κ parameters controlling the rate of decay of spatial correlation with distance(Cressie and Wikle 2011, (Cressie and Wikle 2011)sent independent means estimated for each year, and β 1 and β 2 represent coefficients for log depth (D) and log depth squared (D 2 ).The symbols w s and ε s,t represent spatial and spatiotemporal random effects (respectively) drawn from Gaussian Markov random fields(Cressie and Wikle 2011)with covariance matrices S S ε and S S ω