Compositional Bayesian indicator estimation

Indicator kriging is widely used for mapping spatial binary variables and for estimating the global and local spatial distributions of variables in geosciences. For continuous random variables, indicator kriging gives an estimate of the cumulative distribution function, for a given threshold, which is then the estimate of a probability. Like any other kriging procedure, indicator kriging provides an estimation variance that, although not often used in applications, should be taken into account as it assesses the uncertainty of the estimate. An alternative approach to indicator estimation is proposed in this paper. In this alternative approach the complete probability density function of the indicator estimate is evaluated. The procedure is described in a Bayesian framework, using a multivariate Gaussian likelihood and an a priori distribution which are both combined according to Bayes theorem in order to obtain a posterior distribution for the indicator estimate. From this posterior distribution, point estimates, interval estimates and uncertainty measures can be obtained. Among the point estimates, the median of the posterior distribution is the maximum entropy estimate because there is a fifty-fifty chance of the unknown value of the estimate being larger or smaller than the median; that is, there is maximum uncertainty in the choice between two alternatives. Thus in some sense, the latter is an indicator estimator, alternative to the kriging estimator, that includes its own uncertainty. On the other hand, the mode of the posterior distribution estimator, assuming a uniform prior, is coincidental with the simple kriging estimator. Additionally, because the indicator estimate can be considered as a two-part composition which domain of definition is the simplex, the method is extended to compositional Bayesian indicator estimation. Bayesian indicator estimation and compositional Bayesian indicator estimation are illustrated with an environmental case study in which the probability of the content of a geochemical element in soil being over a particular threshold is of interest. The computer codes and its user guides are public domain and freely available.     

Estimation and prediction of the spatial occurrence of fish species using Bayesian latent Gaussian models

A methodological approach for modelling the occurrence patterns of species for the purpose of fisheries management is proposed here. The presence/absence of the species is modelled with a hierarchical Bayesian spatial model using the geographical and environmental characteristics of each fishing location. Maps of predicted probabilities of presence are generated using Bayesian kriging. Bayesian inference on the parameters and prediction of presence/absence in new locations (Bayesian kriging) are made by considering the model as a latent Gaussian model, which allows the use of the integrated nested Laplace approximation ( INLA ) software (which has been seen to be quite a bit faster than the well-known MCMC methods). In particular, the spatial effect has been implemented with the stochastic partial differential equation (SPDE) approach. The methodology is evaluated on Mediterranean horse mackerel (Trachurus mediterraneus) in the Western Mediterranean. The analysis shows that environmental and geographical factors can play an important role in directing local distribution and variability in the occurrence of species. Although this approach is used to recognize the habitat of mackerel, it could also be for other different species and life stages in order to improve knowledge of fish populations and communities.     

Kriging with external drift in a Birnbaum–Saunders geostatistical model

Spatial models to describe dependent georeferenced data are applied in different fields and, particularly, are used to analyze earth and environmental data. Most of these applications are carried out under Gaussian spatial models. The Birnbaum–Saunders distribution is a unimodal and positively skewed model which has received considerable attention in several areas, including earth and environmental sciences. In addition, theoretical arguments have been provided to justify its usage in the data modeling from these sciences, at least in the same settings where the lognormal distribution can be employed. This paper presents kriging with external drift based on a Birnbaum–Saunders spatial model. The maximum likelihood method is considered to estimate its parameters. The results obtained in the paper are illustrated by an experimental data set related to agricultural management. Specifically, in this illustration, the spatial variability of magnesium content in the soil as a function of calcium content is analyzed.     

Copula-based geostatistical modeling of continuous and discrete data including covariates

It is common in geostatistics to use the variogram to describe the spatial dependence structure and to use kriging as the spatial prediction methodology. Both methods are sensitive to outlying observations and are strongly influenced by the marginal distribution of the underlying random field. Hence, they lead to unreliable results when applied to extreme value or multimodal data. As an alternative to traditional spatial modeling and interpolation we consider the use of copula functions. This paper extends existing copula-based geostatistical models. We show how location dependent covariates e.g. a spatial trend can be accounted for in spatial copula models. Furthermore, we introduce geostatistical copula-based models that are able to deal with random fields having discrete marginal distributions. We propose three different copula-based spatial interpolation methods. By exploiting the relationship between bivariate copulas and indicator covariances, we present indicator kriging and disjunctive kriging. As a second method we present simple kriging of the rank-transformed data. The third method is a plug-in prediction and generalizes the frequently applied trans-Gaussian kriging. Finally, we report on the results obtained for the so-called Helicopter data set which contains extreme radioactivity measurements.     

Application of the BME approach to soil texture mapping

In order to derive accurate space/time maps of soil properties, soil scientists need tools that combine the usually scarce hard data sets with the more easily accessible soft data sets. In the field of modern geostatistics, the Bayesian maximum entropy (BME) approach provides new and powerful means for incorporating various forms of physical knowledge (including hard and soft data, soil classification charts, land cover data from satellite pictures, and digital elevation models) into the space/time mapping process. BME produces the complete probability distribution at each estimation point, thus allowing the calculation of elaborate statistics (even when the distribution is not Gaussian). It also offers a more rigorous and systematic method than kriging for integrating uncertain information into space/time mapping. In this work, BME is used to estimate the three textural fractions involved in a texture map. The first case study focuses on the estimation of the clay fraction, whereas the second one considers the three textural fractions (sand, silt and clay) simultaneously. The BME maps obtained are informative (important soil characteristics are identified, natural variations are well reproduced, etc.). Furthermore, in both case studies, the estimates obtained by BME were more accurate than the simple kriging (SK) estimates, thus offering a better picture of soil reality. In the multivariate case, classification error rate analysis in terms of BME performs considerably better than in terms of kriging. Analysis in terms of BME can offer valuable information to be used in sampling design, in optimizing the hard to soft data ratio, etc.     

Kriging with external drift for functional data for air quality monitoring

Functional data featured by a spatial dependence structure occur in many environmental sciences when curves are observed, for example, along time or along depth. Recently, some methods allowing for the prediction of a curve at an unmonitored site have been developed. However, the existing methods do not allow to include in a model exogenous variables that, for example, bring meteorology information in modeling air pollutant concentrations. In order to introduce exogenous variables, potentially observed as curves as well, we propose to extend the so-called kriging with external drift—or regression kriging—to the case of functional data by means of a three-step procedure involving functional modeling for the trend and spatial interpolation of functional residuals. A cross-validation analysis allows to choose smoothing parameters and a preferable kriging predictor for the functional residuals. Our case study considers daily PM₁₀ concentrations measured from October 2005 to March 2006 by the monitoring network of Piemonte region (Italy), with the trend defined by meteorological time-varying covariates and orographical constant-in-time variables. The performance of the proposed methodology is evaluated by predicting PM₁₀ concentration curves on 10 validation sites, even with simulated realistic datasets on a larger number of spatial sites. In this application the proposed methodology represents an alternative to spatio-temporal modeling but it can be applied more generally to spatially dependent functional data whose domain is not a time interval.     

Regularized kriging as a generalization of simple, universal, and bayesian kriging

In this article the properties of regularized kriging (RK) are studied. RK is obtained as a result of relaxing the universal kriging (UK) non-bias condition by using the support vectors methodology. More specifically, we demonstrate how RK is a continuum of solutions in function of the regularizing parameter, which includes as particular and extreme cases, simple kriging (SK) and UK, and as an intermediate case, Bayesian kriging (BK). Likewise, expressions are obtained for the mean, variance and mean squared error (MSE), as also the expression for the corresponding estimator of the coefficients of the mean. Finally, we investigate the relationship between RK and the support vector machines. By means of simulations we compare the MSE for RK with those for BK and UK, for different association models, for different levels of noise, and for differently sized mean coefficients. The RK results prove to be an improvement on the UK and BK results, and, moreover, these improvements are proportionally greater for greater levels of noise.     

Spatial sampling design and covariance-robust minimax prediction based on convex design ideas

This paper presents new ideas on sampling design and minimax prediction in a geostatistical model setting. Both presented methodologies are based on regression design ideas. For this reason the appendix of this paper gives an introduction to optimum Bayesian experimental design theory for linear regression models with uncorrelated errors. The presented methodologies and algorithms are then applied to the spatial setting of correlated random fields. To be specific, in Sect. 1 we will approximate an isotropic random field by means of a regression model with a large number of regression functions with random amplitudes, similarly to Fedorov and Flanagan (J Combat Inf Syst Sci: 23, 1997). These authors make use of the Karhunen Loeve approximation of the isotropic random field. We use the so-called polar spectral approximation instead; i.e. we approximate the isotropic random field by means of a regression model with sine-cosine-Bessel surface harmonics with random amplitudes and then, in accordance with Fedorov and Flanagan (J Combat Inf Syst Sci: 23, 1997), apply standard Bayesian experimental design algorithms to the resulting Bayesian regression model. Section 2 deals with minimax prediction when the covariance function is known to vary in some set of a priori plausible covariance functions. Using a minimax theorem due to Sion (Pac J Math 8:171–176, 1958) we are able to formulate the minimax problem as being equivalent to an optimum experimental design problem, too. This makes the whole experimental design apparatus available for finding minimax kriging predictors. Furthermore some hints are given, how the approach to spatial sampling design with one a priori fixed covariance function may be extended by means of minimax kriging to a whole set of a priori plausible covariance functions such that the resulting designs are robust. The theoretical developments are illustrated with two examples taken from radiological monitoring and soil science.     

BME prediction of continuous geographical properties using auxiliary variables

Using auxiliary information to improve the prediction accuracy of soil properties in a physically meaningful and technically efficient manner has been widely recognized in pedometrics. In this paper, we explored a novel technique to effectively integrate sampling data and auxiliary environmental information, including continuous and categorical variables, within the framework of the Bayesian maximum entropy (BME) theory. Soil samples and observed auxiliary variables were combined to generate probability distributions of the predicted soil variable at unsampled points. These probability distributions served as soft data of the BME theory at the unsampled locations, and, together with the hard data (sample points) were used in spatial BME prediction. To gain practical insight, the proposed approach was implemented in a real-world case study involving a dataset of soil total nitrogen (TN) contents in the Shayang County of the Hubei Province (China). Five terrain indices, soil types, and soil texture were used as auxiliary variables to generate soft data. Spatial distribution of soil total nitrogen was predicted by BME, regression kriging (RK) with auxiliary variables, and ordinary kriging (OK). The results of the prediction techniques were compared in terms of the Pearson correlation coefficient (r), mean error (ME), and root mean squared error (RMSE). These results showed that the BME predictions were less biased and more accurate than those of the kriging techniques. In sum, the present work extended the BME approach to implement certain kinds of auxiliary information in a rigorous and efficient manner. Our findings showed that the BME prediction technique involving the transformation of variables into soft data can improve prediction accuracy considerably, compared to other techniques currently in use, like RK and OK.     

Why do we need and how should we implement Bayesian kriging methods

The spatial prediction methodology that has become known under the heading of kriging is largely based on the assumptions that the underlying random field is Gaussian and the covariance function is exactly known. In practical applications, however, these assumptions will not hold. Beyond Gaussianity of the random field, lognormal kriging, disjunctive kriging, (generalized linear) model-based kriging and trans-Gaussian kriging have been proposed in the literature. The latter approach makes use of the Box–Cox-transform of the data. Still, all the alternatives mentioned do not take into account the uncertainty with respect to the distribution (or transformation) and the estimated covariance function of the data. The Bayesian trans-Gaussian kriging methodology proposed in the present paper is in the spirit of the “Bayesian bootstrap” idea advocated by Rubin (Ann Stat 9:130–134, 1981) and avoids the unusual specification of noninformative priors often made in the literature and is entirely based on the sample distribution of the estimators of the covariance function and of the Box–Cox parameter. After some notes on Bayesian spatial prediction, noninformative priors and developing our new methodology finally we will present an example illustrating our pragmatic approach to Bayesian prediction by means of a simulated data set.     

Selecting best mapping strategies for storm runoff modeling in a mountainous semi‐arid area

Accurate runoff and soil erosion modeling is constrained by data availability, particularly for physically based models such as OpenLISEM that are data demanding, as the processes are calculated on a cell‐by‐cell basis. The first decision when using such models is to select mapping units that best reflect the spatial variability of the soil and hydraulic properties in the catchment. In environments with limited data, available maps are usually generic, with large units that may lump together the values of the soil properties, affecting the spatial patterns of the predictions and output values in the outlet. Conversely, the output results may be equally acceptable, following the principle of equifinality. To studyhow the mapping method selected affects the model outputs, four types of input maps with different degrees of complexity were created: average values allocated to general soil map units (ASG1), average values allocated to detailed map units (ASG2), values interpolated by ordinary kriging (OK) and interpolated by kriging with external drift (KED). The study area was Ribeira Seca, a 90 km² catchment located in Santiago Island, Cape Verde (West Africa), a semi‐arid country subject to scarce but extreme rainfall during the short tropical summer monsoon. To evaluate the influence of rainfall on runoff and erosion, two storm events with different intensity and duration were considered. OK and KED inputs produced similar results, with the latter being closer to the observed hydrographs. The highest soil losses were obtained with KED (43 ton ha^{? 1} for the strongest event). To improve the results of soil loss predictions, higher accurate spatial information on the processes is needed; however, spatial information of input soil properties alone is not enough in complex landscapes. The results demonstrate the importance of selecting the appropriate mapping strategy to obtain reliable runoff and erosion estimates. Copyright © 2013 John Wiley & Sons, Ltd.     

Spatial prediction of categorical variables: the Bayesian maximum entropy approach

 Being a non-linear method based on a rigorous formalism and an efficient processing of various information sources, the Bayesian maximum entropy (BME) approach has proven to be a very powerful method in the context of continuous spatial random fields, providing much more satisfactory estimates than those obtained from traditional linear geostatistics (i.e., the various kriging techniques). This paper aims at presenting an extension of the BME formalism in the context of categorical spatial random fields. In the first part of the paper, the indicator kriging and cokriging methods are briefly presented and discussed. A special emphasis is put on their inherent limitations, both from the theoretical and practical point of view. The second part aims at presenting the theoretical developments of the BME approach for the case of categorical variables. The three-stage procedure is explained and the formulations for obtaining prior joint distributions and computing posterior conditional distributions are given for various typical cases. The last part of the paper consists in a simulation study for assessing the performance of BME over the traditional indicator (co)kriging techniques. The results of these simulations highlight the theoretical limitations of the indicator approach (negative probability estimates, probability distributions that do not sum up to one, etc.) as well as the much better performance of the BME approach. Estimates are very close to the theoretical conditional probabilities, that can be computed according to the stated simulation hypotheses.     

Variability in surface moisture content along a hillslope transect: Rattlesnake Hill, Texas

Surface soil moisture content exhibits a high degree of spatial and temporal variability. The purpose of this study was (a) to characterize variations in moisture content in the 0–5 cm surface soil layer along a hillslope transect by means of intensive sampling in both space and time; and (b) to make inferences regarding the environmental factors that influence this variability. Over a period of seven months, soil moisture content was measured (gravimetric method) on a near-daily basis at 10 m intervals along a 200 m downslope transect at the Rattlesnake Hill field site in Austin, Texas. Results indicate that significant variability in soil moisture content exists along the length of the transect; that variability decreases with decreasing transect-mean moisture content as the hillslope dries down following rain events; and that the dominant influences on moisture content variability are dependent upon the moisture conditions on the hillslope. While topographic and soil attributes operate jointly to redistribute soil water following storm events, under wet conditions, variability in surface moisture content is most strongly influenced by porosity and hydraulic conductivity, and under dry conditions, correlations are strongest to relative elevation, aspect and clay content. Consequently, the dominant influence on soil moisture variability gradually changes from soil heterogeneity to joint control by topographic and soil properties as the transect dries following significant rain events.     

Use of variograms for field magnetometry analysis in Upper Silesia Industrial Region

Measurements of topsoil magnetic susceptibility are often used for quick assessment of soil contamination of anthropogenic origin, with heavy metals or other pollutants. However, because of complicated correlations between low-field magnetic susceptibility (shortened to magnetic susceptibility) of topsoil and soil pollution, the outcome of a field magnetometry survey can not be related directly to soil pollution. For each case study, the results should be interpreted on their own taking into account not only the type of pollution but also pedogenic, biogenic and environmental factors. In practice, it is very difficult to measure and consider all these factors. Here we illustrate the merit of geostatistical methods, which are focused on the spatial variability of a phenomenon, in the interpretation of soil magnetometry results. This article presents the analysis of spatial variability of top soil layers magnetic susceptibility-within the Upper Silesia Industrial Region (USIR)-using semivariance analysis. It also explains how to adjust the sampling density of field magnetometry measurements to spatial variability of the soil pollution as well as to the spatial scale of the investigated area. For this purpose, the values of magnetic susceptibility have been measured by using various sampling densities at areas of different size located within USIR. This enabled to determine the main scales of magnetic susceptibility spatial variability of soils within USIR using semivariance. A few distinct scales of variability were found from the site scale to a more regional scale. Variability ranges of 30 km, 12 km, and 5 km refer to the large regional scale, whereas smaller ranges of few hundreds down to a few tens of meters, can be attributed to the local (site) scale. In addition, the precision of the measuring campaigns, performed within USIR with different sampling densities, was compared through the analysis of the spatial variability of the soil magnetic susceptibility signal by using ordinary kriging. jarek97@yahoo.com, piotr.fabijanczyk@is.pw.edu.pl     

Defining probability of saturation with indicator kriging on hard and soft data

In humid, well-vegetated areas, such as in the northeastern US, runoff is most commonly generated from relatively small portions of the landscape becoming completely saturated, however, little is known about the spatial and temporal behavior of these saturated regions. Indicator kriging provides a way to use traditional water table data to quantify probability of saturation to evaluate predicted spatial distributions of runoff generation risk, especially for the new generation of water quality models incorporating saturation excess runoff theory. When spatial measurements of a variable are transformed to binary indicators (i.e., 1 if above a given threshold value and 0 if below) and the resulting indicator semivariogram is modeled, indicator kriging produces the probability of the measured variable to exceed the threshold value. Indicator kriging gives quantified probability of saturation or, consistent with saturation excess runoff theory, runoff generation risk with depth to water table as the variable and the threshold set near the soil surface. The probability of saturation for a 120 m × 180 m hillslope based upon 43 measurements of depth to water table is investigated with indicator semivariograms for six storm events. The indicator semivariograms show high spatial structure in saturated regions with large antecedent rainfall conditions. The temporal structure of the data is used to generate interpolated (soft) data to supplement measured (hard) data. This improved the spatial structure of the indicator semivariograms for lower antecedent rainfall conditions. Probability of saturation was evaluated through indicator kriging incorporating soft data showing, based on this preliminary study, highly connected regions of saturation as expected for the wet season (April through May) in the Catskill Mountain region of New York State. Supplementation of hard data with soft data incorporates physical hydrology of the hillslope to capture significant patterns not available when using hard data alone for indicator kriging. With the need for water quality models incorporating appropriate runoff generation risk estimates on the rise, this manner of data will lay the groundwork for future model evaluation and development.     

Test of statistical means for the extrapolation of soil depth point information using overlays of spatial environmental data and bootstrapping techniques

Hydrological modelling depends highly on the accuracy and uncertainty of model input parameters such as soil properties. Since most of these data are field surveyed, geostatistical techniques such as kriging, classification and regression trees or more sophisticated soil‐landscape models need to be applied to interpolate point information to the area. Most of the existing interpolation techniques require a random or regular distribution of points within the study area but are not adequate to satisfactorily interpolate soil catena or transect data. The soil landscape model presented in this study is predicting soil information from transect or catena point data using a statistical mean (arithmetic, geometric and harmonic mean) to calculate the soil information based on class means of merged spatial explanatory variables. A data set of 226 soil depth measurements covering a range of 0–6·5 m was used to test the model. The point data were sampled along four transects in the Stubbetorp catchment, SE‐Sweden. We overlaid a geomorphology map (8 classes) with digital elevation model‐derived topographic index maps (2–9 classes) to estimate the range of error the model produces with changing sample size and input maps. The accuracy of the soil depth predictions was estimated with the root mean square error (RMSE) based on a testing and training data set. RMSE ranged generally between 0·73 and 0·83 m ± 0·013 m depending on the amount of classes the merged layers had, but were smallest for a map combination with a low number of classes predicted with the harmonic mean (RMSE = 0·46 m). The results show that the prediction accuracy of this method depends on the number of point values in the sample, the value range of the measured attribute and the initial correlations between point values and explanatory variables, but suggests that the model approach is in general scale invariant. Copyright © 2009 John Wiley & Sons, Ltd.     

A non-stationary geostatistical approach to multigaussian kriging for local reserve estimation

Multigaussian kriging technique has many applications in mining, soil science, environmental science and other fields. Particularly, in the local reserve estimation of a mineral deposit, multigaussian kriging is employed to derive panel-wise tonnages by predicting conditional probability of block grades. Additionally, integration of a suitable change of support model is also required to estimate the functions of the variables with larger support than that of the samples. However, under the assumption of strict stationarity, the grade distributions and important recovery functions are estimated by multigaussian kriging using samples within a supposedly spatial homogeneous domain. Conventionally, the underlying random function model is required to be stationary in order to carry out the inference on ore grade distribution and relevant statistics. In reality, conventional stationary model often fails to represent complicated geological structure. Traditionally, the simple stationary model neither considers the obvious changes in local means and variances, nor is it able to replicate spatial continuity of the deposit and hence produces unreliable outcomes. This study deals with the theoretical design of a non-stationary multigaussian kriging model allowing change of support and its application in the mineral reserve estimation scenario. Local multivariate distributions are assumed here to be strictly stationary in the neighborhood of the panels. The local cumulative distribution function and related statistics with respect to the panels are estimated using a distance kernel approach. A rigorous investigation through simulation experiments is performed to analyze the relevance of the developed model followed by a case study on a copper deposit.     

Relative model score: a scoring rule for evaluating ensemble simulations with application to microbial soil respiration modeling

This paper defines a new scoring rule, namely relative model score (RMS), for evaluating ensemble simulations of environmental models. RMS implicitly incorporates the measures of ensemble mean accuracy, prediction interval precision, and prediction interval reliability for evaluating the overall model predictive performance. RMS is numerically evaluated from the probability density functions of ensemble simulations given by individual models or several models via model averaging. We demonstrate the advantages of using RMS through an example of soil respiration modeling. The example considers two alternative models with different fidelity, and for each model Bayesian inverse modeling is conducted using two different likelihood functions. This gives four single-model ensembles of model simulations. For each likelihood function, Bayesian model averaging is applied to the ensemble simulations of the two models, resulting in two multi-model prediction ensembles. Predictive performance for these ensembles is evaluated using various scoring rules. Results show that RMS outperforms the commonly used scoring rules of log-score, pseudo Bayes factor based on Bayesian model evidence (BME), and continuous ranked probability score (CRPS). RMS avoids the problem of rounding error specific to log-score. Being applicable to any likelihood functions, RMS has broader applicability than BME that is only applicable to the same likelihood function of multiple models. By directly considering the relative score of candidate models at each cross-validation datum, RMS results in more plausible model ranking than CRPS. Therefore, RMS is considered as a robust scoring rule for evaluating predictive performance of single-model and multi-model prediction ensembles.     

Model-based geostatistical interpolation of the annual number of ozone exceedance days in the Netherlands

This paper discusses two model-based geostatistical methods for spatial interpolation of the number of days that ground level ozone exceeds a threshold level. The first method assumes counts to approximately follow a Poisson distribution, while the second method assumes a log-Normal distribution. First, these methods were compared using an extensive data set covering the Netherlands, Belgium and Germany. Second, the focus was placed on only the Netherlands, where only a small data set was used. Bayesian techniques were used for parameter estimation and interpolation. Parameter estimates are comparable due to the log-link in both models. Incorporating data from adjacent countries improves parameter estimation. The Poisson model predicts more accurately (maximum kriging standard deviation of 2.16 compared to 2.69) but shows smoother surfaces than the log-Normal model. The log-Normal approach ensures a better representation of the observations and gives more realistic patterns (an RMSE of 2.26 compared to 2.44). Model-based geostatistical procedures are useful to interpolate limited data sets of counts of ozone exceedance days. Spatial risk estimates using existing prior information can be made relating health effects to environmental thresholds.     

Modeling locally varying anisotropy of CO2 emissions in the United States

In a spatial property modeling context, the variables of interest to be modeled often display complex nonlinear features. Techniques to incorporate these nonlinear features, such as multiple point statistics or cummulants, are often complex with input parameters that are difficult to infer. The methodology proposed in this paper uses a classical vector-based definition of locally varying anisotropy to characterize nonlinear features and incorporate locally varying anisotropy into numerical property models. The required input is an exhaustive field of anisotropy orientation and magnitude. The methodology consists of (1) using the shortest path distance between locations to define the covariance between points in space (2) multidimensional scaling of the domain to ensure positive definite kriging equations and (3) estimation or simulation with kriging or sequential Gaussian simulation. The only additional parameter required when kriging or simulating with locally varying anisotropy is the number of dimensions to retain in multidimensional scaling. The methodology is demonstrated on a CO₂ emissions data set for the United States in 2002 and shows an improvement in cross validation results as well as a visual reproduction of nonlinear features.