Ecological forecasting of tree growth: Regional fusion of tree-ring and forest inventory data to quantify drivers and characterize uncertainty.

Robust ecological forecasting of tree growth under future climate conditions is critical to anticipate future forest carbon storage and flux. Here, we apply three ingredients of ecological forecasting that are key to improving forecast skill: data fusion, confronting model predictions with new data, and partitioning forecast uncertainty. Specifically, we present the first fusion of tree-ring and forest inventory data within a Bayesian state-space model at a multi-site, regional scale, focusing on Pinus ponderosa var. brachyptera in the southwestern US. Leveraging the complementarity of these two data sources, we parsed the ecological complexity of tree growth into the effects of climate, tree size, stand density, site quality, and their interactions, and quantified uncertainties associated with these effects. New measurements of trees, an ongoing process in forest inventories, were used to confront forecasts of tree diameter with observations, and evaluate alternative tree growth models. We forecasted tree diameter and increment in response to an ensemble of climate change projections, and separated forecast uncertainty into four different causes: initial conditions, parameters, climate drivers, and process error. We found a strong negative effect of fall-spring maximum temperature, and a positive effect of water-year precipitation on tree growth. Furthermore, tree vulnerability to climate stress increases with greater competition, with tree size, and at poor sites. Under future climate scenarios, we forecast increment declines of 22%-117%, while the combined effect of climate and size-related trends results in a 56%-91% decline. Partitioning of forecast uncertainty showed that diameter forecast uncertainty is primarily caused by parameter and initial conditions uncertainty, but increment forecast uncertainty is mostly caused by process error and climate driver uncertainty. This fusion of tree-ring and forest inventory data lays the foundation for robust ecological forecasting of aboveground biomass and carbon accounting at tree, plot, and regional scales, including iterative improvement of model skill.

Trends in bird abundance differ among protected forests but not bird guilds.

Improved monitoring and associated inferential tools to efficiently identify declining bird populations, particularly of rare or sparsely distributed species, is key to informed conservation and management across large spatiotemporal regions. We assess abundance trends for 106 bird species in a network of eight forested national parks located within the northeast United States from 2006 to 2019 using a novel hierarchical model. We develop a multispecies, multiregion, removal-sampling model that shares information across species and parks to enable inference on rare species and sparsely sampled parks and to evaluate the effects of local forest structure. Trends in bird abundance over time varied widely across parks, but species showed similar trends within parks. Three parks (Acadia National Park and Marsh-Billings-Rockefeller and Morristown National Historical Parks [NHP]) decreased in bird abundance across all species, while three parks (Saratoga NHP and Roosevelt-Vanderbilt and Weir-Farm National Historic Sites) increased in abundance. Bird abundance peaked at medium levels of basal area and high levels of percent forest and forest regeneration, with percent forest having the largest effect. Variation in these effects across parks could be a result of differences in forest structural stage and diversity. By sharing information across both communities and parks, our novel hierarchical model enables uncertainty-quantified estimates of abundance across multiple geographical (i.e., network, park) and taxonomic (i.e., community, guild, species) levels over a large spatiotemporal region. We found large variation in abundance trends across parks but not across bird guilds, suggesting that local forest condition might have a broad and consistent effect on the entire bird community within a given park. Research should target the three parks with overall decreasing trends in bird abundance to further identify what specific factors are driving observed declines across the bird community. Understanding how bird communities respond to local forest structure and other stressors (e.g., pest outbreaks, climate change) is crucial for informed and lasting management.

Towards connecting biodiversity and geodiversity across scales with satellite remote sensing.

ISSUE: Geodiversity (i.e., the variation in Earth's abiotic processes and features) has strong effects on biodiversity patterns. However, major gaps remain in our understanding of how relationships between biodiversity and geodiversity vary over space and time. Biodiversity data are globally sparse and concentrated in particular regions. In contrast, many forms of geodiversity can be measured continuously across the globe with satellite remote sensing. Satellite remote sensing directly measures environmental variables with grain sizes as small as tens of metres and can therefore elucidate biodiversity-geodiversity relationships across scales. EVIDENCE: We show how one important geodiversity variable, elevation, relates to alpha, beta and gamma taxonomic diversity of trees across spatial scales. We use elevation from NASA's Shuttle Radar Topography Mission (SRTM) and c. 16,000 Forest Inventory and Analysis plots to quantify spatial scaling relationships between biodiversity and geodiversity with generalized linear models (for alpha and gamma diversity) and beta regression (for beta diversity) across five spatial grains ranging from 5 to 100 km. We illustrate different relationships depending on the form of diversity; beta and gamma diversity show the strongest relationship with variation in elevation. CONCLUSION: With the onset of climate change, it is more important than ever to examine geodiversity for its potential to foster biodiversity. Widely available satellite remotely sensed geodiversity data offer an important and expanding suite of measurements for understanding and predicting changes in different forms of biodiversity across scales. Interdisciplinary research teams spanning biodiversity, geoscience and remote sensing are well poised to advance understanding of biodiversity-geodiversity relationships across scales and guide the conservation of nature.

Spatial Factor Models for High-Dimensional and Large Spatial Data: An Application in Forest Variable Mapping.

Gathering information about forest variables is an expensive and arduous activity. As such, directly collecting the data required to produce high-resolution maps over large spatial domains is infeasible. Next generation collection initiatives of remotely sensed Light Detection and Ranging (LiDAR) data are specifically aimed at producing complete-coverage maps over large spatial domains. Given that LiDAR data and forest characteristics are often strongly correlated, it is possible to make use of the former to model, predict, and map forest variables over regions of interest. This entails dealing with the high-dimensional (~102) spatially dependent LiDAR outcomes over a large number of locations (~105-106). With this in mind, we develop the Spatial Factor Nearest Neighbor Gaussian Process (SF-NNGP) model, and embed it in a two-stage approach that connects the spatial structure found in LiDAR signals with forest variables. We provide a simulation experiment that demonstrates inferential and predictive performance of the SF-NNGP, and use the two-stage modeling strategy to generate complete-coverage maps of forest variables with associated uncertainty over a large region of boreal forests in interior Alaska.

Variable effects of climate on forest growth in relation to climate extremes, disturbance, and forest dynamics.

Changes in the frequency, duration, and severity of climate extremes are forecast to occur under global climate change. The impacts of climate extremes on forest productivity and health remain difficult to predict due to potential interactions with disturbance events and forest dynamics-changes in forest stand composition, density, size and age structure over time. Such interactions may lead to non-linear forest growth responses to climate involving thresholds and lag effects. Understanding how forest dynamics influence growth responses to climate is particularly important given stand structure and composition can be modified through management to increase forest resistance and resilience to climate change. To inform such adaptive management, we develop a hierarchical Bayesian state space model in which climate effects on tree growth are allowed to vary over time and in relation to past climate extremes, disturbance events, and forest dynamics. The model is an important step toward integrating disturbance and forest dynamics into predictions of forest growth responses to climate extremes. We apply the model to a dendrochronology data set from forest stands of varying composition, structure, and development stage in northeastern Minnesota that have experienced extreme climate years and forest tent caterpillar defoliation events. Mean forest growth was most sensitive to water balance variables representing climatic water deficit. Forest growth responses to water deficit were partitioned into responses driven by climatic threshold exceedances and interactions with insect defoliation. Forest growth was both resistant and resilient to climate extremes with the majority of forest growth responses occurring after multiple climatic threshold exceedances across seasons and years. Interactions between climate and disturbance were observed in a subset of years with insect defoliation increasing forest growth sensitivity to water availability. Forest growth was particularly sensitive to climate extremes during periods of high stem density following major regeneration events when average inter-tree competition was high. Results suggest the resistance and resilience of forest growth to climate extremes can be increased through management steps such as thinning to reduce competition during early stages of stand development and small-group selection harvests to maintain forest structures characteristic of older, mature stands.

Seedling survival responses to conspecific density, soil nutrients, and irradiance vary with age in a tropical forest.

Understanding processes that promote species coexistence is integral to diversity maintenance. In hyperdiverse tropical forests, local conspecific density (LCD) and light are influential to woody seedling recruitment and soil nutrients are often limiting, yet the simultaneous effects of these factors on seedling survival across time remain unknown. We fit species- and age-specific models to census and resource data of seedlings of 68 woody species from a Costa Rican wet tropical forest. In decreasing order of prevalence, seedling survivorship was related to LCD, soil base cations, irradiance, nitrogen, and phosphorus. Species-specific responses to factors did not covary, providing evidence that species life history strategies have not converged to one continuum of high-surviving stress tolerant to low-surviving stress intolerant species. Survival responses to all factors varied over the average seedling's lifetime, indicating seedling requirements change with age and conclusions drawn about processes important to species coexistence depend on temporal resolution.

Predicting tree biomass growth in the temperate-boreal ecotone: Is tree size, age, competition, or climate response most important?

As global temperatures rise, variation in annual climate is also changing, with unknown consequences for forest biomes. Growing forests have the ability to capture atmospheric CO2 and thereby slow rising CO2 concentrations. Forests' ongoing ability to sequester C depends on how tree communities respond to changes in climate variation. Much of what we know about tree and forest response to climate variation comes from tree-ring records. Yet typical tree-ring datasets and models do not capture the diversity of climate responses that exist within and among trees and species. We address this issue using a model that estimates individual tree response to climate variables while accounting for variation in individuals' size, age, competitive status, and spatially structured latent covariates. Our model allows for inference about variance within and among species. We quantify how variables influence aboveground biomass growth of individual trees from a representative sample of 15 northern or southern tree species growing in a transition zone between boreal and temperate biomes. Individual trees varied in their growth response to fluctuating mean annual temperature and summer moisture stress. The variation among individuals within a species was wider than mean differences among species. The effects of mean temperature and summer moisture stress interacted, such that warm years produced positive responses to summer moisture availability and cool years produced negative responses. As climate models project significant increases in annual temperatures, growth of species like Acer saccharum, Quercus rubra, and Picea glauca will vary more in response to summer moisture stress than in the past. The magnitude of biomass growth variation in response to annual climate was 92-95% smaller than responses to tree size and age. This means that measuring or predicting the physical structure of current and future forests could tell us more about future C dynamics than growth responses related to climate change alone.

Landscape fragmentation affects responses of avian communities to climate change.

Forecasting the consequences of climate change is contingent upon our understanding of the relationship between biodiversity patterns and climatic variability. While the impacts of climate change on individual species have been well-documented, there is a paucity of studies on climate-mediated changes in community dynamics. Our objectives were to investigate the relationship between temporal turnover in avian biodiversity and changes in climatic conditions and to assess the role of landscape fragmentation in affecting this relationship. We hypothesized that community turnover would be highest in regions experiencing the most pronounced changes in climate and that these patterns would be reduced in human-dominated landscapes. To test this hypothesis, we quantified temporal turnover in avian communities over a 20-year period using data from the New York State Breeding Atlases collected during 1980-1985 and 2000-2005. We applied Bayesian spatially varying intercept models to evaluate the relationship between temporal turnover and temporal trends in climatic conditions and landscape fragmentation. We found that models including interaction terms between climate change and landscape fragmentation were superior to models without the interaction terms, suggesting that the relationship between avian community turnover and changes in climatic conditions was affected by the level of landscape fragmentation. Specifically, we found weaker associations between temporal turnover and climatic change in regions with prevalent habitat fragmentation. We suggest that avian communities in fragmented landscapes are more robust to climate change than communities found in contiguous habitats because they are comprised of species with wider thermal niches and thus are less susceptible to shifts in climatic variability. We conclude that highly fragmented regions are likely to undergo less pronounced changes in composition and structure of faunal communities as a result of climate change, whereas those changes are likely to be greater in contiguous and unfragmented habitats.

Joint species distribution models with imperfect detection for high-dimensional spatial data.

Determining the spatial distributions of species and communities is a key task in ecology and conservation efforts. Joint species distribution models are a fundamental tool in community ecology that use multi-species detection-nondetection data to estimate species distributions and biodiversity metrics. The analysis of such data is complicated by residual correlations between species, imperfect detection, and spatial autocorrelation. While many methods exist to accommodate each of these complexities, there are few examples in the literature that address and explore all three complexities simultaneously. Here we developed a spatial factor multi-species occupancy model to explicitly account for species correlations, imperfect detection, and spatial autocorrelation. The proposed model uses a spatial factor dimension reduction approach and Nearest Neighbor Gaussian Processes to ensure computational efficiency for data sets with both a large number of species (e.g., >100) and spatial locations (e.g., 100,000). We compared the proposed model performance to five alternative models, each addressing a subset of the three complexities. We implemented the proposed and alternative models in the spOccupancy software, designed to facilitate application via an accessible, well documented, and open-source R package. Using simulations, we found that ignoring the three complexities when present leads to inferior model predictive performance, and the impacts of failing to account for one or more complexities will depend on the objectives of a given study. Using a case study on 98 bird species across the continental US, the spatial factor multi-species occupancy model had the highest predictive performance among the alternative models. Our proposed framework, together with its implementation in spOccupancy, serves as a user-friendly tool to understand spatial variation in species distributions and biodiversity while addressing common complexities in multi-species detection-nondetection data.

Highly Scalable Bayesian Geostatistical Modeling via Meshed Gaussian Processes on Partitioned Domains.

We introduce a class of scalable Bayesian hierarchical models for the analysis of massive geostatistical datasets. The underlying idea combines ideas on high-dimensional geostatistics by partitioning the spatial domain and modeling the regions in the partition using a sparsity-inducing directed acyclic graph (DAG). We extend the model over the DAG to a well-defined spatial process, which we call the Meshed Gaussian Process (MGP). A major contribution is the development of a MGPs on tessellated domains, accompanied by a Gibbs sampler for the efficient recovery of spatial random effects. In particular, the cubic MGP (Q-MGP) can harness high-performance computing resources by executing all large-scale operations in parallel within the Gibbs sampler, improving mixing and computing time compared to sequential updating schemes. Unlike some existing models for large spatial data, a Q-MGP facilitates massive caching of expensive matrix operations, making it particularly apt in dealing with spatiotemporal remote-sensing data. We compare Q-MGPs with large synthetic and real world data against state-of-the-art methods. We also illustrate using Normalized Difference Vegetation Index (NDVI) data from the Serengeti park region to recover latent multivariate spatiotemporal random effects at millions of locations. The source code is available at github.com/mkln/meshgp.

Adaptive Gaussian Predictive Process Models for Large Spatial Datasets.

Large point referenced datasets occur frequently in the environmental and natural sciences. Use of Bayesian hierarchical spatial models for analyzing these datasets is undermined by onerous computational burdens associated with parameter estimation. Low-rank spatial process models attempt to resolve this problem by projecting spatial effects to a lower-dimensional subspace. This subspace is determined by a judicious choice of "knots" or locations that are fixed a priori. One such representation yields a class of predictive process models (e.g., Banerjee et al., 2008) for spatial and spatial-temporal data. Our contribution here expands upon predictive process models with fixed knots to models that accommodate stochastic modeling of the knots. We view the knots as emerging from a point pattern and investigate how such adaptive specifications can yield more flexible hierarchical frameworks that lead to automated knot selection and substantial computational benefits.

Over half of western United States' most abundant tree species in decline.

Changing forest disturbance regimes and climate are driving accelerated tree mortality across temperate forests. However, it remains unknown if elevated mortality has induced decline of tree populations and the ecological, economic, and social benefits they provide. Here, we develop a standardized forest demographic index and use it to quantify trends in tree population dynamics over the last two decades in the western United States. The rate and pattern of change we observe across species and tree size-distributions is alarming and often undesirable. We observe significant population decline in a majority of species examined, show decline was particularly severe, albeit size-dependent, among subalpine tree species, and provide evidence of widespread shifts in the size-structure of montane forests. Our findings offer a stark warning of changing forest composition and structure across the western US, and suggest that sustained anthropogenic and natural stress will likely result in broad-scale transformation of temperate forests globally.

Hierarchical spatial modeling of additive and dominance genetic variance for large spatial trial datasets.

SUMMARY: This article expands upon recent interest in Bayesian hierarchical models in quantitative genetics by developing spatial process models for inference on additive and dominance genetic variance within the context of large spatially referenced trial datasets. Direct application of such models to large spatial datasets are, however, computationally infeasible because of cubic-order matrix algorithms involved in estimation. The situation is even worse in Markov chain Monte Carlo (MCMC) contexts where such computations are performed for several iterations. Here, we discuss approaches that help obviate these hurdles without sacrificing the richness in modeling. For genetic effects, we demonstrate how an initial spectral decomposition of the relationship matrices negate the expensive matrix inversions required in previously proposed MCMC methods. For spatial effects, we outline two approaches for circumventing the prohibitively expensive matrix decompositions: the first leverages analytical results from Ornstein-Uhlenbeck processes that yield computationally efficient tridiagonal structures, whereas the second derives a modified predictive process model from the original model by projecting its realizations to a lower-dimensional subspace, thereby reducing the computational burden. We illustrate the proposed methods using a synthetic dataset with additive, dominance, genetic effects and anisotropic spatial residuals, and a large dataset from a Scots pine (Pinus sylvestris L.) progeny study conducted in northern Sweden. Our approaches enable us to provide a comprehensive analysis of this large trial, which amply demonstrates that, in addition to violating basic assumptions of the linear model, ignoring spatial effects can result in downwardly biased measures of heritability.

Improving the performance of predictive process modeling for large datasets.

Advances in Geographical Information Systems (GIS) and Global Positioning Systems (GPS) enable accurate geocoding of locations where scientific data are collected. This has encouraged collection of large spatial datasets in many fields and has generated considerable interest in statistical modeling for location-referenced spatial data. The setting where the number of locations yielding observations is too large to fit the desired hierarchical spatial random effects models using Markov chain Monte Carlo methods is considered. This problem is exacerbated in spatial-temporal and multivariate settings where many observations occur at each location. The recently proposed predictive process, motivated by kriging ideas, aims to maintain the richness of desired hierarchical spatial modeling specifications in the presence of large datasets. A shortcoming of the original formulation of the predictive process is that it induces a positive bias in the non-spatial error term of the models. A modified predictive process is proposed to address this problem. The predictive process approach is knot-based leading to questions regarding knot design. An algorithm is designed to achieve approximately optimal spatial placement of knots. Detailed illustrations of the modified predictive process using multivariate spatial regression with both a simulated and a real dataset are offered.

Efficient algorithms for Bayesian Nearest Neighbor Gaussian Processes.

We consider alternate formulations of recently proposed hierarchical Nearest Neighbor Gaussian Process (NNGP) models (Datta et al., 2016a) for improved convergence, faster computing time, and more robust and reproducible Bayesian inference. Algorithms are defined that improve CPU memory management and exploit existing high-performance numerical linear algebra libraries. Computational and inferential benefits are assessed for alternate NNGP specifications using simulated datasets and remotely sensed light detection and ranging (LiDAR) data collected over the US Forest Service Tanana Inventory Unit (TIU) in a remote portion of Interior Alaska. The resulting data product is the first statistically robust map of forest canopy for the TIU.

A Case Study Competition Among Methods for Analyzing Large Spatial Data.

The Gaussian process is an indispensable tool for spatial data analysts. The onset of the "big data" era, however, has lead to the traditional Gaussian process being computationally infeasible for modern spatial data. As such, various alternatives to the full Gaussian process that are more amenable to handling big spatial data have been proposed. These modern methods often exploit low-rank structures and/or multi-core and multi-threaded computing environments to facilitate computation. This study provides, first, an introductory overview of several methods for analyzing large spatial data. Second, this study describes the results of a predictive competition among the described methods as implemented by different groups with strong expertise in the methodology. Specifically, each research group was provided with two training datasets (one simulated and one observed) along with a set of prediction locations. Each group then wrote their own implementation of their method to produce predictions at the given location and each was subsequently run on a common computing environment. The methods were then compared in terms of various predictive diagnostics. Supplementary materials regarding implementation details of the methods and code are available for this article online. ELECTRONIC SUPPLEMENTARY MATERIAL: Supplementary materials for this article are available at 10.1007/s13253-018-00348-w.

Assessing impact of exogenous features on biotic phenomena in the presence of strong spatial dependence: A lake sturgeon case study in natural stream settings.

Modeling spatially explicit data provides a powerful approach to identify the effects of exogenous features associated with biological processes, including recruitment of stream fishes. However, the complex spatial and temporal dynamics of the stream and the species' reproductive and early life stage behaviors present challenges to drawing valid inference using traditional regression models. In these settings it is often difficult to ensure the spatial independence among model residuals-a key assumption that must be met to ensure valid inference. We present statistical models capable of capturing complex residual anisotropic patterns through the addition of spatial random effects within an inferential framework that acknowledges uncertainty in the data and parameters. Proposed models are used to explore the impact of environmental variables on Lake sturgeon (Acipenser fulvescens) reproduction, particularly questions about patterns in egg deposition. Our results demonstrate the need to apply valid statistical methods to identify relationships between response variables, e.g., egg counts, across locations, and environmental covariates in the presence of strong and anisotropic autocorrelation in stream systems. The models may be applied to other settings where gamete distribution or, more generally, other biotic phenomena may be associated with spatially dynamic and anisotropic processes.

spBayes: An R Package for Univariate and Multivariate Hierarchical Point-referenced Spatial Models.

Scientists and investigators in such diverse fields as geological and environmental sciences, ecology, forestry, disease mapping, and economics often encounter spatially referenced data collected over a fixed set of locations with coordinates (latitude-longitude, Easting-Northing etc.) in a region of study. Such point-referenced or geostatistical data are often best analyzed with Bayesian hierarchical models. Unfortunately, fitting such models involves computationally intensive Markov chain Monte Carlo (MCMC) methods whose efficiency depends upon the specific problem at hand. This requires extensive coding on the part of the user and the situation is not helped by the lack of available software for such algorithms. Here, we introduce a statistical software package, spBayes, built upon the R statistical computing platform that implements a generalized template encompassing a wide variety of Gaussian spatial process models for univariate as well as multivariate point-referenced data. We discuss the algorithms behind our package and illustrate its use with a synthetic and real data example.

Spatial Variation in Nutrient and Water Color Effects on Lake Chlorophyll at Macroscales.

The nutrient-water color paradigm is a framework to characterize lake trophic status by relating lake primary productivity to both nutrients and water color, the colored component of dissolved organic carbon. Total phosphorus (TP), a limiting nutrient, and water color, a strong light attenuator, influence lake chlorophyll a concentrations (CHL). But, these relationships have been shown in previous studies to be highly variable, which may be related to differences in lake and catchment geomorphology, the forms of nutrients and carbon entering the system, and lake community composition. Because many of these factors vary across space it is likely that lake nutrient and water color relationships with CHL exhibit spatial autocorrelation, such that lakes near one another have similar relationships compared to lakes further away. Including this spatial dependency in models may improve CHL predictions and clarify how well the nutrient-water color paradigm applies to lakes distributed across diverse landscape settings. However, few studies have explicitly examined spatial heterogeneity in the effects of TP and water color together on lake CHL. In this study, we examined spatial variation in TP and water color relationships with CHL in over 800 north temperate lakes using spatially-varying coefficient models (SVC), a robust statistical method that applies a Bayesian framework to explore space-varying and scale-dependent relationships. We found that TP and water color relationships were spatially autocorrelated and that allowing for these relationships to vary by individual lakes over space improved the model fit and predictive performance as compared to models that did not vary over space. The magnitudes of TP effects on CHL differed across lakes such that a 1 µg/L increase in TP resulted in increased CHL ranging from 2-24 µg/L across lake locations. Water color was not related to CHL for the majority of lakes, but there were some locations where water color had a positive effect such that a unit increase in water color resulted in a 2 µg/L increase in CHL and other locations where it had a negative effect such that a unit increase in water color resulted in a 2 µg/L decrease in CHL. In addition, the spatial scales that captured variation in TP and water color effects were different for our study lakes. Variation in TP-CHL relationships was observed at intermediate distances (~20 km) compared to variation in water color-CHL relationships that was observed at regional distances (~200 km). These results demonstrate that there are lake-to-lake differences in the effects of TP and water color on lake CHL and that this variation is spatially structured. Quantifying spatial structure in these relationships furthers our understanding of the variability in these relationships at macroscales and would improve model prediction of chlorophyll a to better meet lake management goals.

On nearest-neighbor Gaussian process models for massive spatial data.

Gaussian Process (GP) models provide a very flexible nonparametric approach to modeling location-and-time indexed datasets. However, the storage and computational requirements for GP models are infeasible for large spatial datasets. Nearest Neighbor Gaussian Processes (Datta A, Banerjee S, Finley AO, Gelfand AE. Hierarchical nearest-neighbor gaussian process models for large geostatistical datasets. J Am Stat Assoc 2016., JASA) provide a scalable alternative by using local information from few nearest neighbors. Scalability is achieved by using the neighbor sets in a conditional specification of the model. We show how this is equivalent to sparse modeling of Cholesky factors of large covariance matrices. We also discuss a general approach to construct scalable Gaussian Processes using sparse local kriging. We present a multivariate data analysis which demonstrates how the nearest neighbor approach yields inference indistinguishable from the full rank GP despite being several times faster. Finally, we also propose a variant of the NNGP model for automating the selection of the neighbor set size.