Robust estimation of the pseudo cross‐variogram for cokriging soil properties

The pseudo cross‐variogram can be used for cokriging two or more soil properties when few or none of the sampling locations have values recorded for all of them. The usual estimator of the pseudo cross‐variogram is susceptible to the effects of extreme data (outliers). This will lead to overestimation of the error variance of predictions obtained by cokriging. A solution to this problem is to use robust estimators of the pseudo cross‐variogram, and three such estimators are proposed in this paper. The robust estimators were demonstrated on simulated data in the presence of different numbers of outlying data drawn from different contaminating distributions. The robust estimators were less sensitive to the outliers than the non‐robust one, but they had larger variances. Outliers tend to obscure the spatial structure of the cross‐correlation of the simulated variables as described by the non‐robust estimator. The several estimators of the pseudo cross‐variogram were applied to a multitemporal data set on soil water content. Since these were obtained non‐destructively, direct measurements of temporal change can be made. A prediction subset of the data was subsampled as if obtained by destructive analysis and the remainder used for validation. Estimators of the auto‐variogram and pseudo cross‐variogram were applied to the prediction data, then used to predict the change in water content at the validation sites by cokriging. The estimation variances of these predictions were best calculated with a robustly estimated model of coregionalization, although the validation set was too small to conclude that the non‐robust estimators were unsuitable in this instance.     

A comparison of some robust estimators of the variogram for use in soil survey

The standard estimator of the variogram is sensitive to outlying data, a few of which can cause overestimation of the variogram. This will result in incorrect variances when estimating the value of a soil property by kriging or when designing a sampling grid to map the property to a required precision. Several robust estimators of the variogram, based on location and scale estimation, have been proposed as improvements. They seem to be suitable for analysis of soil data in circumstances where the standard estimator is likely to be affected by outliers. Robust estimators are based on assumptions about the distribution of the data which will not always hold and which need not be made in kriging or in estimating the variogram by the standard estimator. The estimators are reviewed. Simulation studies show that the robust estimators vary in their susceptibility to moderate skew in the underlying distribution, but that the effects of outliers are generally greater. The estimators are applied to some soil data, and the resulting variograms used for ordinary kriging at sites in a separate validation data set. In most cases the variograms derived from the standard estimator gave kriging variances which appeared to overestimate the mean squared error of prediction (MSEP). Kriging with variograms based on robust estimators sometimes gave kriging variances which underestimated the MSEP or did not differ significantly from it. Estimates of kriging variance and the MSEP derived from the validation data were generally close to estimates from cross‐validation on the prediction set used to derive the variograms. This indicates that variogram models derived from different estimators could be compared by cross‐validation.     

Estimating variograms of soil properties by the method-of-moments and maximum likelihood

Variograms of soil properties are usually obtained by estimating the variogram for distinct lag classes by the method‐of‐moments and fitting an appropriate model to the estimates. An alternative is to fit a model by maximum likelihood to data on the assumption that they are a realization of a multivariate Gaussian process. This paper compares the two using both simulation and real data. The method‐of‐moments and maximum likelihood were used to estimate the variograms of data simulated from stationary Gaussian processes. In one example, where the simulated field was sampled at different intensities, maximum likelihood estimation was consistently more efficient than the method‐of‐moments, but this result was not general and the relative performance of the methods depends on the form of the variogram. Where the nugget variance was relatively small and the correlation range of the data was large the method‐of‐moments was at an advantage and likewise in the presence of data from a contaminating distribution. When fields were simulated with positive skew this affected the results of both the method‐of‐moments and maximum likelihood. The two methods were used to estimate variograms from actual metal concentrations in topsoil in the Swiss Jura, and the variograms were used for kriging. Both estimators were susceptible to sampling problems which resulted in over‐ or underestimation of the variance of three of the metals by kriging. For four other metals the results for kriging using the variogram obtained by maximum likelihood were consistently closer to the theoretical expectation than the results for kriging with the variogram obtained by the method‐of‐moments, although the differences between the results using the two approaches were not significantly different from each other or from expectation. Soil scientists should use both procedures in their analysis and compare the results.     

On spatial prediction of soil properties in the presence of a spatial trend: the empirical best linear unbiased predictor (E-BLUP) with REML

Geostatistical estimates of a soil property by kriging are equivalent to the best linear unbiased predictions (BLUPs). Universal kriging is BLUP with a fixed‐effect model that is some linear function of spatial co‐ordinates, or more generally a linear function of some other secondary predictor variable when it is called kriging with external drift. A problem in universal kriging is to find a spatial variance model for the random variation, since empirical variograms estimated from the data by method‐of‐moments will be affected by both the random variation and that variation represented by the fixed effects. The geostatistical model of spatial variation is a special case of the linear mixed model where our data are modelled as the additive combination of fixed effects (e.g. the unknown mean, coefficients of a trend model), random effects (the spatially dependent random variation in the geostatistical context) and independent random error (nugget variation in geostatistics). Statisticians use residual maximum likelihood (REML) to estimate variance parameters, i.e. to obtain the variogram in a geostatistical context. REML estimates are consistent (they converge in probability to the parameters that are estimated) with less bias than both maximum likelihood estimates and method‐of‐moment estimates obtained from residuals of a fitted trend. If the estimate of the random effects variance model is inserted into the BLUP we have the empirical BLUP or E‐BLUP. Despite representing the state of the art for prediction from a linear mixed model in statistics, the REML–E‐BLUP has not been widely used in soil science, and in most studies reported in the soils literature the variogram is estimated with methods that are seriously biased if the fixed‐effect structure is more complex than just an unknown constant mean (ordinary kriging). In this paper we describe the REML–E‐BLUP and illustrate the method with some data on soil water content that exhibit a pronounced spatial trend.     

The assessment of point and diffuse metal pollution of soils from an urban geochemical survey of Sheffield, England

Abstract. A model of soil variability as a continuous background process with superimposed point contamination was applied to 569 measurements of metal concentrations (Cr, Ni and Pb) in the topsoils of Sheffield, England. Robust estimators of the variogram were shown to be required to describe spatial variation of the metal concentrations at most sampled locations. This is diagnostic of the presence of a contaminant process. Values of the standardized kriging error from the cross‐validation of each datum were used to identify spatial outliers for each metal. The ordinary kriged estimates of Cr, Ni and Pb were mapped after removing the outliers to estimate the background variation. Each of the 35 spatial outliers that occured in gardens have concentrations exceeding their Soil Guideline Value for residential land use with plant uptake, highlighting a potentially significant exposure pathway. The frequent observation of coal and furnace waste at these sites suggests that their dispersal, following domestic use and industrial processes, respectively, represents a significant point contaminant process. There was no evidence for spatial clustering of the point process. However, the spatial outliers of Cr and Ni showed a significant association with disturbed sites identified from historical land use maps, in part due to their prevalence in areas of historical steel manufacture. The magnitude of diffuse pollution for each metal in the urban soil was estimated by removing the spatial outliers and comparing robust measures of location with those from a survey of soils developed over the same parent materials in adjacent rural and peri‐urban environments. The Winsorized mean Pb concentrations in urban topsoil (203 mg kg^?1) were twice the value in the rural environment (101 mg kg^?1), highlighting a very substantial diffuse Pb load to urban soils. The equivalent estimated diffuse components in urban soils for Cr and Ni were, respectively, 25% and 14% higher than the rural soils.     

Estimating temporal change in soil monitoring: I. Statistical theory

Detecting small temporal change of spatially varying soil properties demands precise estimation. Design– and model–based methods are compared for estimating temporal change of soil properties over finite areas. Analytical expressions for the estimators and their variances arc derived for the two approaches, and formulae for the expectations of the variances under the random–process model are developed. Among the randomized designs simple, stratified, and systematic random sampling using the arithmetic mean as estimator have been studied. Pairing the sampling positions on the different occasions increases the precision of design–based estimation if the observations are positively cross–correlated. The relative precisions of the means of stratified and systematic samples depends on the spatial correlation. Neither is more precise than the other in all circumstances. The stratified design provides an unbiased estimator for the sampling error, which is not available from systematic samples. Theoretically, the geostatistical global estimator is more precise than the estimates derived from any of the classical designs when many realizations arc repeatedly sampled at random. In practice, with only a single realization of the process, this is no longer relevant. Moreover, errors in estimating the variograms add to the total error of the method. It seems that only by sampling from large auto–correlated random fields can the precisions of the methods be compared in practice.     

Fluctuations in method‐of‐moments variograms caused by clustered sampling and their elimination by declustering and residual maximum likelihood estimation

Soil data accumulated in national and regional archives derive from many sources and tend to be concentrated in zones of particular interest. Experimental variograms computed from such data by the usual method of moments can appear highly erratic, and therefore models fitted to them are likely to be unreliable. We have explored two methods of avoiding the effects, one by computing declustering weights and incorporating them into the method of moments, the other using residual maximum likelihood. The methods are illustrated with data on bulk density, exchangeable magnesium, cation exchange capacity and organic carbon of 4182 samples of soil from numerous soil surveys in the whole of Australia and stored in the CSIRO's national archive. The experimental variograms of all four variables are erratic. Cell declustering produced much smoother sequences of estimates to which plausible models could be fitted with confidence. The residual maximum likelihood models closely matched those models over several hundred km. Finally values were simulated at the same sampling points from the residual maximum likelihood models, reproducing ‘spiky’ experimental variograms such as those computed from the data. The simulation showed that clustered design of sampling causes spiky artefacts. We conclude that where data are clustered experimental variograms should be computed with declustered weighting or variogram models be fitted by residual maximum likelihood.     

Optimal interpolation and isarithmic mapping of soil properties: I The semi‐variogram and punctual kriging

Kriging is a means of spatial prediction that can be used for soil properties. It is a form of weighted local averaging. It is optimal in the sense that it provides estimates of values at unrecorded places without bias and with minimum and known variance. Isarithmic maps made by kriging are alternatives to conventional soil maps where properties can be measured at close spacings. Kriging depends on first computing an accurate semi‐variogram, which measures the nature of spatial dependence for the property. Estimates of semi‐variance are then used to determine the weights applied to the data when computing the averages, and are presented in the kriging equations. The method is applied to three sets of data from detailed soil surveys in Central Wales and Norfolk. Sodium content at Plas Gogerddan was shown to vary isotropically with a linear semi‐variogram. Ordinary punctual kriging produced a map with intricate isarithms and fairly large estimation variance, attributed to a large nugget effect. Stoniness on the same land varied anisotropically with a linear semi‐variogram, and again the estimation error of punctual kriging was fairly large. At Hole Farm, Norfolk, the thickness of cover loam varied isotropically, but with a spherical semi‐variogram. Its parameters were estimated and used to krige point values and produce a map showing substantial short‐range variation.     

Bayesian Calibration of Blue Crab (<Emphasis Type="Italic">Callinectes sapidus</Emphasis>) Abundance Indices Based on Probability Surveys

Abundance and standard error estimates in surveys of fishery resources typically employ classical design-based approaches, ignoring the influences of non-design factors such as varying catchability. We developed a Bayesian approach for estimating abundance and associated errors in a fishery survey by incorporating sampling and non-sampling variabilities. First, a zero-inflated spatial model was used to quantify variance components due to non-sampling factors; second, the model was used to calibrate the estimated abundance index and its variance using pseudo empirical likelihood. The approach was applied to a winter dredge survey conducted to estimate the abundance of blue crabs (Callinectes sapidus) in the Chesapeake Bay. We explored the properties of the calibration estimators through a limited simulation study. The variance estimator calibrated on posterior sample performed well, and the mean estimator had comparable performance to design-based approach with slightly higher bias and lower (about 15% reduction) mean squared error. The results suggest that application of this approach can improve estimation of abundance indices using data from design-based fishery surveys.     

Micro-Spatial Variability of Soil Nitrate Following Nitrogen Fertilization and Drip Irrigation

Micro-spatial analysis of nitrate (NO₃), an environmental contaminant partially attributed to nitrogen fertilization, can be useful for estimating its distribution in soils. A study was conducted to determine the micro-spatial distribution of soil NO₃ using kriging and cokriging in a drip-irrigated and nitrogen-fertilized field. One hundred soil samples were collected in a regular grid pattern from a 10 m × 20 m plot, and analyzed for soil NO₃ and pH. The effect of reduced sample size on NO₃ estimation was also evaluated. The pH data indicated the soils were slightly acidic to neutral with log[NO₃] values ranging from 1.66 to 2.95. These parameters were inversely related; which was probably an attribute of soil nitrification process. Sample variograms and cross-variograms suggested that the spatial distribution of pH and log[NO₃] could be described by linear models in the area studied, as indicated by small MSE (mean sum error), and RKV (reduced kriging variance) values close to 1. Contour maps based on kriging and cokriging estimates indicated greater homogeneity of the variables in the south-north direction than the east-west, except for zones of high NO₃ and low pH in the north-central edge and north-east corner of the grid area. Cokriging of log[NO₃] estimation, using pH data, improved MSE, MSSE (mean sum square error), MKV (mean kriging variance), RKV, CEE (correlation between estimated data and error), CEM (correlation between estimated and measured data) by 46, 31, 30, 22, 96, and 98%, respectively, as compared to kriging. Lower cokriging variance for any estimated log[NO₃] value, as compared to the kriging analysis, indicated that cokriging provided more accurate estimates. With reduced sample observations (n) for NO₃ similar conclusions were obtained; and the estimation accuracy was maintained up to n >70. Cokriging analysis with reduced n also curtailed the analytical cost, and facilitated NO₃ estimation by means of pH, which was measured at a cheaper cost.     

Incorporating parameter uncertainty into prediction intervals for spatial data modeled via a parametric variogram

In spatial predictions, researchers usually treat the estimated theoretical variogram parameters as known without error and ignore the variability of the parameter estimators. Although the prediction is still unbiased, the prediction error is usually underestimated. Therefore, the coverage probability of the prediction interval usually is lower than the nominal probability. A simulation study is performed to show how the coverage probability for prediction relates to the true range and sill of an exponential variogram. This article proposes two parametric bootstrap methods to incorporate the variability of the corresponding parameter estimators. A simulation study is performed to evaluate the coverage probability of these proposed methods. Finally, we apply the parametric bootstrap methods to a real dataset and compare the results with those from naive (i.e., treating estimated parameters as known) and Bayesian methods.     

Comparing sampling needs for variograms of soil properties computed by the method of moments and residual maximum likelihood

It has been generally accepted that the method of moments (MoM) variogram, which has been widely applied in soil science, requires about 100 sites at an appropriate interval apart to describe the variation adequately. This sample size is often larger than can be afforded for soil surveys of agricultural fields or contaminated sites. Furthermore, it might be a much larger sample size than is needed where the scale of variation is large. A possible alternative in such situations is the residual maximum likelihood (REML) variogram because fewer data appear to be required. The REML method is parametric and is considered reliable where there is trend in the data because it is based on generalized increments that filter trend out and only the covariance parameters are estimated. Previous research has suggested that fewer data are needed to compute a reliable variogram using a maximum likelihood approach such as REML, however, the results can vary according to the nature of the spatial variation. There remain issues to examine: how many fewer data can be used, how should the sampling sites be distributed over the site of interest, and how do different degrees of spatial variation affect the data requirements? The soil of four field sites of different size, physiography, parent material and soil type was sampled intensively, and MoM and REML variograms were calculated for clay content. The data were then sub-sampled to give different sample sizes and distributions of sites and the variograms were computed again. The model parameters for the sets of variograms for each site were used for cross-validation. Predictions based on REML variograms were generally more accurate than those from MoM variograms with fewer than 100 sampling sites. A sample size of around 50 sites at an appropriate distance apart, possibly determined from variograms of ancillary data, appears adequate to compute REML variograms for kriging soil properties for precision agriculture and contaminated sites.     

Estimating Distance Sampling Detection Functions When Distances Are Measured With Errors

Distance sampling methods assume that distances are known but in practice there are often errors in measuring them. These can have substantial impact on the bias and precision of distance sampling estimators. In this paper we develop methods that accommodate both systematic and stochastic measurement errors. We use the methods to estimate detection probability in two surveys with substantial measurement error. The first is a shipboard line transect survey in the North Sea in which information on measurement error comes from photographically measured distances to a subset of detections. The second is an aerial cue-counting survey off Iceland in which information on measurement error comes from pairs of independently estimated distances to a subset of detections. Different methods are required for measurement error estimation in the two cases. We investigate by simulation the properties of the new estimators and compare them to conventional estimators. They are found to perform better than conventional estimators in the presence of measurement error, more so in the case of cue-counting and point transect estimators than line transect estimators. An appendix on the asymptotic distributions of conditional and full likelihood estimators is available online.     

Wavelet analysis of soil variation at nanometre‐ to micrometre‐scales: an example of organic carbon content in a micro‐aggregate

This paper demonstrates the potential of wavelet analysis to investigate fine‐scale spatial variation in soil without statistical assumptions that are generally implausible. We analysed the optical densities of different forms of carbon which were measured at intervals of 50 nm along a 16‐µm transect on a soil micro‐aggregate using near‐edge X‐ray fine‐structure spectroscopy (NEXAFS). We found different patterns of scale‐dependent variation between the carbon forms, which could be represented by pair‐wise wavelet correlations at the different scales, and by principal components analysis of all the correlations at each scale. These results represent only one small soil micro‐aggregate and are not presented as general findings about soil carbon, but they do indicate that fine‐scale variation of soil carbon can be complex in ways that the wavelet analysis can accommodate but alternative spatial statistics such as variograms cannot. Among the patterns of variation that the analysis could identify were scale‐dependent correlations of the different forms of carbon. In some cases, positive correlations were found at coarser scales and negative at the finest scales, suggesting a multi‐scale pattern in which contrasting forms of carbon are deposited in common clumps but at finer scales either one or the other form dominates. Aromatic and carboxylic carbon varied jointly in this way. Other forms of carbon, such as carboxylic and aliphatic carbon, were strongly correlated at the finest scales but not the coarser scales. We found evidence for changes in the variance and correlation of forms of carbon along the transect, indicating that the spatial distribution of carbon at these fine scales may be very complex in ways that are inconsistent with the assumptions of geostatistics. This quantitative analysis of the spatial patterns of different soil components at micro‐scales offers a basis for formulating and testing specific hypotheses on replicated samples.     

Adaptive sampling and reconnaissance surveys for geostatistical mapping of the soil

The effort required to survey a soil variable depends upon the acceptable uncertainty of estimates and the variogram of the variable. The variogram is unknown prior to sampling, so it must be inferred from a reconnaissance survey before an efficient survey can be designed. The results of reconnaissance surveys are subject to uncertainty, which depends upon the variogram and the number and location of observations. Here, we develop an adaptive approach for optimizing reconnaissance surveys. The observations within these reconnaissance surveys are collected in distinct phases. After each phase, a probability density function of the required sampling density of the main survey is calculated within a Bayesian framework. The number and location of observations within further phases are selected to reduce efficiently the uncertainty of the estimate of the required sampling density. In simulation studies, the number and location of observations in Bayesian adaptive reconnaissance surveys vary according to the variogram of the property of interest. For variograms with a short range, the reconnaissance surveys are intensive with a large proportion of clustered locations. Fewer, more evenly spread locations are required for variables with a long range. Bayesian adaptive reconnaissance surveys lead to more efficient surveys than conventional approaches because the reconnaissance survey is specifically designed for the variable of interest. A hand‐held field system is implemented and tested in a survey of soil moisture content over a field.     

基于最大似然法与矩法的黄土高原小流域土壤碳氮空间变异分析

以黄土高原寺底沟小流域为研究对象,根据不同土地利用方式采集46个样点的土壤样品,通过地统计方法对土壤有机碳和全氮的空间变异特征进行了分析。采用受限最大似然法(REML)和矩法(MOM)两种方法分别对变异函数进行了估计,通过交叉检验选择克里金预测效果较好的变异函数进行地统计插值。(1)与矩法(MOM)相比,在多数情况下受限最大似然法(REML)估计的变异函数进行克里金插值更加准确。(2)土层深度对土壤全氮空间变异影响较小,对土壤有机碳影响较大,表层土壤有机碳含量及变异程度明显高于下层土壤。(3)土地利用方式对土壤有机碳和全氮的空间分布有重要影响,灌木林和天然草地土壤有机碳和全氮水平最高,弃耕地其次,梯田、果园、人工草地最低,表明退耕还林对提高土壤碳氮水平有重要贡献。     

A new variance estimator for parameters of semiparametric generalized additive models

Generalized additive models (GAMs) have become popular in the air pollution epidemiology literature. Two problems, recently surfaced, concern implementation of these semiparametric models. The first problem, easily corrected, was laxity of the default convergence criteria. The other, noted independently by Klein, Flanders, and Tolbert, and Ramsay, Burnett, and Krewski concerned variance estimates produced by commercially available software. In simulations, they were as much as 50% too small. We derive an expression for a variance estimator for the parametric component of generalized additive models that can include up to three smoothing splines, and show how the standard error (SE) estimated by this method differs from the corresponding SE estimated with error in a study of air pollution and emergency room admissions for cardiorespiratory disease. The derivation is based on asymptotic linearity. Using Monte Carlo experiments, we evaluated performance of the estimator in finite samples. The estimator performed well in Monte Carlo experiments, in the situations considered. However, more work is needed to address performance in additional situations. Using data from our study of air pollution and cardiovascular disease, the standard error estimated using the new method was about 10% to 20% larger than the biased, commercially available standard error estimate.     

Model‐based analysis using REML for inference from systematically sampled data on soil

The general linear model encompasses statistical methods such as regression and analysis of variance (anova ) which are commonly used by soil scientists. The standard ordinary least squares (OLS) method for estimating the parameters of the general linear model is a design‐based method that requires that the data have been collected according to an appropriate randomized sample design. Soil data are often obtained by systematic sampling on transects or grids, so OLS methods are not appropriate. Parameters of the general linear model can be estimated from systematically sampled data by model‐based methods. Parameters of a model of the covariance structure of the error are estimated, then used to estimate the remaining parameters of the model with known variance. Residual maximum likelihood (REML) is the best way to estimate the variance parameters since it is unbiased. We present the REML solution to this problem. We then demonstrate how REML can be used to estimate parameters for regression and anova ‐type models using data from two systematic surveys of soil. We compare an efficient, gradient‐based implementation of REML (ASReml) with an implementation that uses simulated annealing. In general the results were very similar; where they differed the error covariance model had a spherical variogram function which can have local optima in its likelihood function. The simulated annealing results were better than the gradient method in this case because simulated annealing is good at escaping local optima.     

On the Akaike Information Criterion for choosing models for variograms of soil properties

A problem in the application of geostatistics to soil is to find satisfactory models for variograms of soil properties. It is usually solved by fitting plausible models to the sample variogram by weighted least squares approximation. The residual sum of squares can always be diminished, and the fit improved in that sense, by adding parameters to the model. A satisfactory compromise between goodness of fit and parsimony can be achieved by applying the Akaike Information Criterion (AIC). For a given set of data the variable part of the AIC is estimated by where n is the number of experimental points on the variogram, R is the residual sum of squares and p is the number of parameters in the model. The model to choose is the one for which  is least.
The AIC is closely related to Akaike's earlier final prediction error and the Schwarz criterion. It is also equivalent to an F test when adding parameters in nested models.     

Cost-effective abundance estimation of rare animals: Testing performance of small-boat surveys for killer whales in British Columbia

Top predators are often rare, subject to anthropogenic mortality, and possess life-history traits that make them inherently vulnerable to extinction. IUCN criteria recognise populations as Critically Endangered when abundance is <250 mature individuals, but estimating abundance of rare species can be more challenging than for common ones. Cost-effective methods are needed to provide robust abundance estimates. In marine environments, small boats are more widely accessible than large ships for researchers conducting sightings surveys with limited funds, but studies are needed into efficacy of small-boat surveys. This study compares line transect and mark-recapture estimates from small-boat surveys in summer 2004 and 2005 for ‘northern resident’ killer whales in British Columbia to true population size, known from censuses conducted by Fisheries and Oceans Canada. The line transect estimate of 195 animals (95% CI 27-559) used model averaging to incorporate uncertainty in the detection function, while the mark-recapture estimate of 239 animals (CI 154-370) used a simple two-sample Chapman estimator. Both methods produced estimates close to the true population size, which numbered 219 animals in 2004 and 235 in 2006, but both suffered from the small sample sizes and violations of some model assumptions that will vex most pilot studies of rare species. Initial abundance estimates from relatively low-cost surveys can be thought of as hypotheses to be tested as new data are collected. For species of conservation concern, any cost-effective attempt to estimate absolute abundance will assist status assessments, as long as estimates are presented with appropriate caveats.