Slow climate velocities of mountain streams portend their role as refugia for cold-water biodiversity.

The imminent demise of montane species is a recurrent theme in the climate change literature, particularly for aquatic species that are constrained to networks and elevational rather than latitudinal retreat as temperatures increase. Predictions of widespread species losses, however, have yet to be fulfilled despite decades of climate change, suggesting that trends are much weaker than anticipated and may be too subtle for detection given the widespread use of sparse water temperature datasets or imprecise surrogates like elevation and air temperature. Through application of large water-temperature databases evaluated for sensitivity to historical air-temperature variability and computationally interpolated to provide high-resolution thermal habitat information for a 222,000-km network, we estimate a less dire thermal plight for cold-water species within mountains of the northwestern United States. Stream warming rates and climate velocities were both relatively low for 1968-2011 (average warming rate = 0.101 °C/decade; median velocity = 1.07 km/decade) when air temperatures warmed at 0.21 °C/decade. Many cold-water vertebrate species occurred in a subset of the network characterized by low climate velocities, and three native species of conservation concern occurred in extremely cold, slow velocity environments (0.33-0.48 km/decade). Examination of aggressive warming scenarios indicated that although network climate velocities could increase, they remain low in headwaters because of strong local temperature gradients associated with topographic controls. Better information about changing hydrology and disturbance regimes is needed to complement these results, but rather than being climatic cul-de-sacs, many mountain streams appear poised to be redoubts for cold-water biodiversity this century.

Last call: Passive acoustic monitoring shows continued rapid decline of critically endangered vaquita.

The vaquita is a critically endangered species of porpoise. It produces echolocation clicks, making it a good candidate for passive acoustic monitoring. A systematic grid of sensors has been deployed for 3 months annually since 2011; results from 2016 are reported here. Statistical models (to compensate for non-uniform data loss) show an overall decline in the acoustic detection rate between 2015 and 2016 of 49% (95% credible interval 82% decline to 8% increase), and total decline between 2011 and 2016 of over 90%. Assuming the acoustic detection rate is proportional to population size, approximately 30 vaquita (95% credible interval 8-96) remained in November 2016.

Modelling dendritic ecological networks in space: an integrated network perspective.

Dendritic ecological networks (DENs) are a unique form of ecological networks that exhibit a dendritic network topology (e.g. stream and cave networks or plant architecture). DENs have a dual spatial representation; as points within the network and as points in geographical space. Consequently, some analytical methods used to quantify relationships in other types of ecological networks, or in 2-D space, may be inadequate for studying the influence of structure and connectivity on ecological processes within DENs. We propose a conceptual taxonomy of network analysis methods that account for DEN characteristics to varying degrees and provide a synthesis of the different approaches within the context of stream ecology. Within this context, we summarise the key innovations of a new family of spatial statistical models that describe spatial relationships in DENs. Finally, we discuss how different network analyses may be combined to address more complex and novel research questions. While our main focus is streams, the taxonomy of network analyses is also relevant anywhere spatial patterns in both network and 2-D space can be used to explore the influence of multi-scale processes on biota and their habitat (e.g. plant morphology and pest infestation, or preferential migration along stream or road corridors).

spmodel: Spatial statistical modeling and prediction in [Formula: see text].

spmodel is an [Formula: see text] package used to fit, summarize, and predict for a variety spatial statistical models applied to point-referenced or areal (lattice) data. Parameters are estimated using various methods, including likelihood-based optimization and weighted least squares based on variograms. Additional modeling features include anisotropy, non-spatial random effects, partition factors, big data approaches, and more. Model-fit statistics are used to summarize, visualize, and compare models. Predictions at unobserved locations are readily obtainable.

Indexing and partitioning the spatial linear model for large data sets.

We consider four main goals when fitting spatial linear models: 1) estimating covariance parameters, 2) estimating fixed effects, 3) kriging (making point predictions), and 4) block-kriging (predicting the average value over a region). Each of these goals can present different challenges when analyzing large spatial data sets. Current research uses a variety of methods, including spatial basis functions (reduced rank), covariance tapering, etc, to achieve these goals. However, spatial indexing, which is very similar to composite likelihood, offers some advantages. We develop a simple framework for all four goals listed above by using indexing to create a block covariance structure and nearest-neighbor predictions while maintaining a coherent linear model. We show exact inference for fixed effects under this block covariance construction. Spatial indexing is very fast, and simulations are used to validate methods and compare to another popular method. We study various sample designs for indexing and our simulations showed that indexing leading to spatially compact partitions are best over a range of sample sizes, autocorrelation values, and generating processes. Partitions can be kept small, on the order of 50 samples per partition. We use nearest-neighbors for kriging and block kriging, finding that 50 nearest-neighbors is sufficient. In all cases, confidence intervals for fixed effects, and prediction intervals for (block) kriging, have appropriate coverage. Some advantages of spatial indexing are that it is available for any valid covariance matrix, can take advantage of parallel computing, and easily extends to non-Euclidean topologies, such as stream networks. We use stream networks to show how spatial indexing can achieve all four goals, listed above, for very large data sets, in a matter of minutes, rather than days, for an example data set.

Discretized and aggregated: modeling dive depth of harbor seals from ordered categorical data with temporal autocorrelation.

Ordered categorical data are pervasive in environmental and ecological data, and often arise from constraints that require discretizing a continuous variable into ordered categories. A great deal of data have been collected toward the study of marine mammal dive behavior using satellite depth recorders (SDRs), which often discretize a continuous variable such as depth. Additionally, data storage or transmission constraints may also necessitate the aggregation of data over time intervals of a specified length. The categorization and aggregation create a time series of ordered multicategory counts for each animal, which present challenges in terms of statistical modeling and practical interpretation. We describe an intuitive strategy for modeling such aggregated, ordered categorical data allowing for inference regarding the category probabilities and a measure of central tendency on the original scale of the data (e.g., meters), along with incorporation of temporal correlation and overdispersion. The strategy extends covariate-specific cutpoint models for ordinal data. We demonstrate the method in an analysis of SDR dive-depth data collected on harbor seals in Alaska. The primary goal of the analysis is to assess the relationship of covariates, such as time of day, with number of dives and maximum depth of dives. We also predict missing values and introduce novel graphical summaries of the data and results.

A Bayesian hierarchical model of Antarctic fur seal foraging and pup growth related to sea ice and prey abundance.

We created a Bayesian hierarchical model (BHM) to investigate ecosystem relationships between the physical ecosystem (sea ice extent), a prey measure (krill density), predator behaviors (diving and foraging effort of female Antarctic fur seals, Arctocephalus gazella, with pups) and predator characteristics (mass of maternal fur seals and pups). We collected data on Antarctic fur seals from 1987/1988 to 1994/1995 at Seal Island, Antarctica. The BHM allowed us to link together predators and prey into a model that uses all the data efficiently and accounts for major sources of uncertainty. Based on the literature, we made hypotheses about the relationships in the model, which we compared with the model outcome after fitting the BHM. For each BHM parameter, we calculated the mean of the posterior density and the 95% credible interval. Our model confirmed others' findings that increased sea ice was related to increased krill density. Higher krill density led to reduced dive intensity of maternal fur seals, as measured by dive depth and duration, and to less time spent foraging by maternal fur seals. Heavier maternal fur seals and lower maternal foraging effort resulted in heavier pups at 22 d. No relationship was found between krill density and maternal mass, or between maternal mass and foraging effort on pup growth rates between 22 and 85 days of age. Maternal mass may have reflected environmental conditions prior to the pup provisioning season, rather than summer prey densities. Maternal mass and foraging effort were not related to pup growth rates between 22 and 85 d, possibly indicating that food was not limiting, food sources other than krill were being used, or differences occurred before pups reached age 22 d.

A mixed-model moving-average approach to geostatistical modeling in stream networks.

Spatial autocorrelation is an intrinsic characteristic in freshwater stream environments where nested watersheds and flow connectivity may produce patterns that are not captured by Euclidean distance. Yet, many common autocovariance functions used in geostatistical models are statistically invalid when Euclidean distance is replaced with hydrologic distance. We use simple worked examples to illustrate a recently developed moving-average approach used to construct two types of valid autocovariance models that are based on hydrologic distances. These models were designed to represent the spatial configuration, longitudinal connectivity, discharge, and flow direction in a stream network. They also exhibit a different covariance structure than Euclidean models and represent a true difference in the way that spatial relationships are represented. Nevertheless, the multi-scale complexities of stream environments may not be fully captured using a model based on one covariance structure. We advocate using a variance component approach, which allows a mixture of autocovariance models (Euclidean and stream models) to be incorporated into a single geostatistical model. As an example, we fit and compare "mixed models," based on multiple covariance structures, for a biological indicator. The mixed model proves to be a flexible approach because many sources of information can be incorporated into a single model.

A model-based approach for making ecological inference from distance sampling data.

We consider a fully model-based approach for the analysis of distance sampling data. Distance sampling has been widely used to estimate abundance (or density) of animals or plants in a spatially explicit study area. There is, however, no readily available method of making statistical inference on the relationships between abundance and environmental covariates. Spatial Poisson process likelihoods can be used to simultaneously estimate detection and intensity parameters by modeling distance sampling data as a thinned spatial point process. A model-based spatial approach to distance sampling data has three main benefits: it allows complex and opportunistic transect designs to be employed, it allows estimation of abundance in small subregions, and it provides a framework to assess the effects of habitat or experimental manipulation on density. We demonstrate the model-based methodology with a small simulation study and analysis of the Dubbo weed data set. In addition, a simple ad hoc method for handling overdispersion is also proposed. The simulation study showed that the model-based approach compared favorably to conventional distance sampling methods for abundance estimation. In addition, the overdispersion correction performed adequately when the number of transects was high. Analysis of the Dubbo data set indicated a transect effect on abundance via Akaike's information criterion model selection. Further goodness-of-fit analysis, however, indicated some potential confounding of intensity with the detection function.

Long-term recovery patterns of arctic tundra after winter seismic exploration.

In response to the increasing global demand for energy, oil exploration and development are expanding into frontier areas of the Arctic, where slow-growing tundra vegetation and the underlying permafrost soils are very sensitive to disturbance. The creation of vehicle trails on the tundra from seismic exploration for oil has accelerated in the past decade, and the cumulative impact represents a geographic footprint that covers a greater extent of Alaska's North Slope tundra than all other direct human impacts combined. Seismic exploration for oil and gas was conducted on the coastal plain of the Arctic National Wildlife Refuge, Alaska, USA, in the winters of 1984 and 1985. This study documents recovery of vegetation and permafrost soils over a two-decade period after vehicle traffic on snow-covered tundra. Paired permanent vegetation plots (disturbed vs. reference) were monitored six times from 1984 to 2002. Data were collected on percent vegetative cover by plant species and on soil and ground ice characteristics. We developed Bayesian hierarchical models, with temporally and spatially autocorrelated errors, to analyze the effects of vegetation type and initial disturbance levels on recovery patterns of the different plant growth forms as well as soil thaw depth. Plant community composition was altered on the trails by species-specific responses to initial disturbance and subsequent changes in substrate. Long-term changes included increased cover of graminoids and decreased cover of evergreen shrubs and mosses. Trails with low levels of initial disturbance usually improved well over time, whereas those with medium to high levels of initial disturbance recovered slowly. Trails on ice-poor, gravel substrates of riparian areas recovered better than those on ice-rich loamy soils of the uplands, even after severe initial damage. Recovery to pre-disturbance communities was not possible where trail subsidence occurred due to thawing of ground ice. Previous studies of disturbance from winter seismic vehicles in the Arctic predicted short-term and mostly aesthetic impacts, but we found that severe impacts to tundra vegetation persisted for two decades after disturbance under some conditions. We recommend management approaches that should be used to prevent persistent tundra damage.

Comparing spatial regression to random forests for large environmental data sets.

Environmental data may be "large" due to number of records, number of covariates, or both. Random forests has a reputation for good predictive performance when using many covariates with nonlinear relationships, whereas spatial regression, when using reduced rank methods, has a reputation for good predictive performance when using many records that are spatially autocorrelated. In this study, we compare these two techniques using a data set containing the macroinvertebrate multimetric index (MMI) at 1859 stream sites with over 200 landscape covariates. A primary application is mapping MMI predictions and prediction errors at 1.1 million perennial stream reaches across the conterminous United States. For the spatial regression model, we develop a novel transformation procedure that estimates Box-Cox transformations to linearize covariate relationships and handles possibly zero-inflated covariates. We find that the spatial regression model with transformations, and a subsequent selection of significant covariates, has cross-validation performance comparable to random forests. We also find that prediction interval coverage is close to nominal for each method, but that spatial regression prediction intervals tend to be narrower and have less variability than quantile regression forest prediction intervals. A simulation study is used to generalize results and clarify advantages of each modeling approach.

Robustness of close-kin mark-recapture estimators to dispersal limitation and spatially varying sampling probabilities.

Close-kin mark-recapture (CKMR) is a method for estimating abundance and vital rates from kinship relationships observed in genetic samples. CKMR inference only requires animals to be sampled once (e.g., lethally), potentially widening the scope of population-level inference relative to traditional monitoring programs.One assumption of CKMR is that, conditional on individual covariates like age, all animals have an equal probability of being sampled. However, if genetic data are collected opportunistically (e.g., via hunters or fishers), there is potential for spatial variation in sampling probability that can bias CKMR estimators, particularly when genetically related individuals stay in close proximity.We used individual-based simulation to investigate consequences of dispersal limitation and spatially biased sampling on performance of naive (nonspatial) CKMR estimators of abundance, fecundity, and adult survival. Population dynamics approximated that of a long-lived mammal species subject to lethal sampling.Naive CKMR abundance estimators were relatively unbiased when dispersal was unconstrained (i.e., complete mixing) or when sampling was random or subject to moderate levels of spatial variation. When dispersal was limited, extreme variation in spatial sampling probabilities negatively biased abundance estimates. Reproductive schedules and survival were well estimated, except for survival when adults could emigrate out of the sampled area. Incomplete mixing was readily detected using Kolmogorov-Smirnov tests.Although CKMR appears promising for estimating abundance and vital rates with opportunistically collected genetic data, care is needed when dispersal limitation is coupled with spatially biased sampling. Fortunately, incomplete mixing is easily detected with adequate sample sizes. In principle, it is possible to devise and fit spatially explicit CKMR models to avoid bias under dispersal limitation, but development of such models necessitates additional complexity (and possibly additional data). We suggest using simulation studies to examine potential bias and precision of proposed modeling approaches prior to implementing a CKMR program.

SSNdesign-An R package for pseudo-Bayesian optimal and adaptive sampling designs on stream networks.

Streams and rivers are biodiverse and provide valuable ecosystem services. Maintaining these ecosystems is an important task, so organisations often monitor the status and trends in stream condition and biodiversity using field sampling and, more recently, autonomous in-situ sensors. However, data collection is often costly, so effective and efficient survey designs are crucial to maximise information while minimising costs. Geostatistics and optimal and adaptive design theory can be used to optimise the placement of sampling sites in freshwater studies and aquatic monitoring programs. Geostatistical modelling and experimental design on stream networks pose statistical challenges due to the branching structure of the network, flow connectivity and directionality, and differences in flow volume. Geostatistical models for stream network data and their unique features already exist. Some basic theory for experimental design in stream environments has also previously been described. However, open source software that makes these design methods available for aquatic scientists does not yet exist. To address this need, we present SSNdesign, an R package for solving optimal and adaptive design problems on stream networks that integrates with existing open-source software. We demonstrate the mathematical foundations of our approach, and illustrate the functionality of SSNdesign using two case studies involving real data from Queensland, Australia. In both case studies we demonstrate that the optimal or adaptive designs outperform random and spatially balanced survey designs implemented in existing open-source software packages. The SSNdesign package has the potential to boost the efficiency of freshwater monitoring efforts and provide much-needed information for freshwater conservation and management.

Accounting for uncertainty in ecological analysis: the strengths and limitations of hierarchical statistical modeling.

Analyses of ecological data should account for the uncertainty in the process(es) that generated the data. However, accounting for these uncertainties is a difficult task, since ecology is known for its complexity. Measurement and/or process errors are often the only sources of uncertainty modeled when addressing complex ecological problems, yet analyses should also account for uncertainty in sampling design, in model specification, in parameters governing the specified model, and in initial and boundary conditions. Only then can we be confident in the scientific inferences and forecasts made from an analysis. Probability and statistics provide a framework that accounts for multiple sources of uncertainty. Given the complexities of ecological studies, the hierarchical statistical model is an invaluable tool. This approach is not new in ecology, and there are many examples (both Bayesian and non-Bayesian) in the literature illustrating the benefits of this approach. In this article, we provide a baseline for concepts, notation, and methods, from which discussion on hierarchical statistical modeling in ecology can proceed. We have also planted some seeds for discussion and tried to show where the practical difficulties lie. Our thesis is that hierarchical statistical modeling is a powerful way of approaching ecological analysis in the presence of inevitable but quantifiable uncertainties, even if practical issues sometimes require pragmatic compromises.

Decline towards extinction of Mexico's vaquita porpoise (Phocoena sinus).

The vaquita (Phocoena sinus) is a small porpoise endemic to Mexico. It is listed by IUCN as Critically Endangered because of unsustainable levels of bycatch in gillnets. The population has been monitored with passive acoustic detectors every summer from 2011 to 2018; here we report results for 2017 and 2018. We combine the acoustic trends with an independent estimate of population size from 2015, and visual observations of at least seven animals in 2017 and six in 2018. Despite adoption of an emergency gillnet ban in May 2015, the estimated rate of decline remains extremely high: 48% decline in 2017 (95% Bayesian credible interval (CRI) 78% decline to 9% increase) and 47% in 2018 (95% CRI 80% decline to 13% increase). Estimated total population decline since 2011 is 98.6%, with greater than 99% probability the decline is greater than 33% yr-1. We estimate fewer than 19 vaquitas remained as of summer 2018 (posterior mean 9, median 8, 95% CRI 6-19). From March 2016 to March 2019, 10 dead vaquitas killed in gillnets were found. The ongoing presence of illegal gillnets despite the emergency ban continues to drive the vaquita towards extinction. Immediate management action is required if the species is to be saved.

A general framework for the analysis of animal resource selection from telemetry data.

We propose a general framework for the analysis of animal telemetry data through the use of weighted distributions. It is shown that several interpretations of resource selection functions arise when constructed from the ratio of a use and availability distribution. Through the proposed general framework, several popular resource selection models are shown to be special cases of the general model by making assumptions about animal movement and behavior. The weighted distribution framework is shown to be easily extended to readily account for telemetry data that are highly autocorrelated; as is typical with use of new technology such as global positioning systems animal relocations. An analysis of simulated data using several models constructed within the proposed framework is also presented to illustrate the possible gains from the flexible modeling framework. The proposed model is applied to a brown bear data set from southeast Alaska.

Habitat selection and seasonal movements of young bearded seals (Erignathus barbatus) in the Bering Sea.

The first year of life is typically the most critical to a pinniped's survival, especially for Arctic phocids which are weaned at only a few weeks of age and left to locate and capture prey on their own. Their seasonal movements and habitat selection are therefore important factors in their survival. During a cooperative effort between scientists and subsistence hunters in October 2004, 2005, and 2006, 13 female and 13 male young (i.e., age <2) bearded seals (Erignathus barbatus) were tagged with satellite-linked dive recorders (SDRs) in Kotzebue Sound, Alaska. Shortly after being released, most seals moved south with the advancing sea-ice through the Bering Strait and into the Bering Sea where they spent the winter and early spring. The SDRs of 17 (8 female and 9 male) seals provided frequent high-quality positions in the Bering Sea; their data were used in our analysis. To investigate habitat selection, we simulated 20 tracks per seal by randomly selecting from the pooled distributions of the absolute bearings and swim speeds of the tagged seals. For each point in the observed and simulated tracks, we obtained the depth, sea-ice concentration, and the distances to sea-ice, open water, the shelf break and coastline. Using logistic regression with a stepwise model selection procedure, we compared the simulated tracks to those of the tagged seals and obtained a model for describing habitat selection. The regression coefficients indicated that the bearded seals in our study selected locations near the ice edge. In contrast, aerial surveys of the bearded seal population, predominantly composed of adults, indicated higher abundances in areas farther north and in heavier pack ice. We hypothesize that this discrepancy is the result of behavioral differences related to age. Ice concentration was also shown to be a statistically significant variable in our model. All else being equal, areas of higher ice concentration are selected for up to about 80%. The effects of sex and bathymetry were not statistically significant. The close association of young bearded seals to the ice edge in the Bering Sea is important given the likely effects of climate warming on the extent of sea-ice and subsequent changes in ice edge habitat.

Quasi-Poisson vs. negative binomial regression: how should we model overdispersed count data?

Quasi-Poisson and negative binomial regression models have equal numbers of parameters, and either could be used for overdispersed count data. While they often give similar results, there can be striking differences in estimating the effects of covariates. We explain when and why such differences occur. The variance of a quasi-Poisson model is a linear function of the mean while the variance of a negative binomial model is a quadratic function of the mean. These variance relationships affect the weights in the iteratively weighted least-squares algorithm of fitting models to data. Because the variance is a function of the mean, large and small counts get weighted differently in quasi-Poisson and negative binomial regression. We provide an example using harbor seal counts from aerial surveys. These counts are affected by date, time of day, and time relative to low tide. We present results on a data set that showed a dramatic difference on estimating abundance of harbor seals when using quasi-Poisson vs. negative binomial regression. This difference is described and explained in light of the different weighting used in each regression method. A general understanding of weighting can help ecologists choose between these two methods.

Trends in spatial patterns of heavy metal deposition on national park service lands along the Red Dog Mine haul road, Alaska, 2001-2006.

Spatial patterns of Zn, Pb and Cd deposition in Cape Krusenstern National Monument (CAKR), Alaska, adjacent to the Red Dog Mine haul road, were characterized in 2001 and 2006 using Hylocomium moss tissue as a biomonitor. Elevated concentrations of Cd, Pb, and Zn in moss tissue decreased logarithmically away from the haul road and the marine port. The metals concentrations in the two years were compared using Bayesian posterior predictions on a new sampling grid to which both data sets were fit. Posterior predictions were simulated 200 times both on a coarse grid of 2,357 points and by distance-based strata including subsets of these points. Compared to 2001, Zn and Pb concentrations in 2006 were 31 to 54% lower in the 3 sampling strata closest to the haul road (0-100, 100-2000 and 2000-4000 m). Pb decreased by 40% in the stratum 4,000-5,000 m from the haul road. Cd decreased significantly by 38% immediately adjacent to the road (0-100m), had an 89% probability of a small decrease 100-2000 m from the road, and showed moderate probabilities (56-71%) for increase at greater distances. There was no significant change over time (with probabilities all ≤ 85%) for any of the 3 elements in more distant reference areas (40-60 km). As in 2001, elemental concentrations in 2006 were higher on the north side of the road. Reductions in deposition have followed a large investment in infrastructure to control fugitive dust escapement at the mine and port sites, operational controls, and road dust mitigation. Fugitive dust escapement, while much reduced, is still resulting in elevated concentrations of Zn, Pb and Cd out to 5,000 m from the haul road. Zn and Pb levels were slightly above arctic baseline values in southern CAKR reference areas.

Geographic pattern of serum antibody prevalence for Brucella spp. in caribou, grizzly bears, and wolves from Alaska, 1975-1998.

Blood samples were collected from 2,635 caribou (Rangifer tarandus), 1,238 grizzly bears (Ursus arctos), and 930 wolves (Canis lupus) from throughout mainland Alaska during 1975-98. Sera were tested for evidence of exposure to Brucella spp. Serum antibody prevalences were highest in the northwestern region of the state. In any specific area, antibody prevalences for caribou and wolves were of a similar magnitude, whereas antibody prevalence for bears in these same areas were two to three times higher.