Prediction of spatial soil property information from ancillary sensor data using ordinary linear regression: Model derivations, residual assumptions and model validation tests

Geospatial measurements of ancillary sensor data, such as bulk soil electrical conductivity or remotely sensed imagery data, are commonly used to characterize spatial variation in soil or crop properties. Geostatistical techniques like kriging with external drift or regression kriging are often used to calibrate geospatial sensor data to specific soil or crop properties. More traditional statistical methods such as ordinary linear regression models are also commonly used. Unfortunately, some soil scientists see these as competing and unrelated modeling approaches and are unaware of their relationship. In this article we review the connection between the ordinary linear regression model and the more comprehensive geostatistical mixed linear model and describe when and under what conditions ordinary linear regression models represent valid spatial prediction models. The formulas for the ordinary linear regression model parameter estimates and best linear unbiased predictions are derived from the geostatistical mixed linear model under two different residual error assumptions; i.e., strictly uncorrelated (SU) residuals and effectively uncorrelated (EU) residuals. The theoretically optimal (best linear unbiased) and computable (linear unbiased) predictions and variance estimates derived under the EU error assumption are examined in detail. Statistical tests for detecting spatial correlation in LR model residuals are also reviewed, in addition to three LR model validation tests derived from classical linear modeling theory. Two case studies are presented that highlight and demonstrate the various parameter estimation, response variable prediction and model validation techniques discussed in this article.     

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 accuracy between ordinary kriging and regression kriging was compared based on the combined consideration of sample size, spatial structure, and auxiliary variables (terrain indices and electromagnetic induction surveys) for a variety of soil properties in two contrasting landscapes (agricultural vs. forested). When spatial structure could not be well captured by point-based observations (e. g., when the ratio of sample spacing over correlation range was > 0.5), or when a strong relationship existed between target soil properties and auxiliary variables (e. g., their R² was > 0.6), regression kriging (RK) was more accurate for interpolating soil properties in both landscapes studied. Otherwise, ordinary kriging (OK) was better. Soil depth and wetness condition did not appear to affect the selection of kriging for soil moisture interpolation, because they did not significantly change the ratio of sample spacing over correlation range and the relationship with the auxiliary variables. Because of a smaller ratio of elevation change over total study area (E/A = 1.2) and multiple parent materials in the agricultural land, OK was generally more accurate in that landscape. In contrast, a larger E/A ratio of 6.8 and a single parent material led to RK being preferable in the steep-sloped forested catchment. The results from this study can be useful for selecting kriging for various soil properties and landscapes.     

Kriging a soil variable with a simple nonstationary variance model

Kriging is a standard tool in the environmental sciences for spatial prediction from limited sample data, subject to the assumption of intrinsic stationarity, made about the underlying spatially correlated random function. It is generally well understood how the assumption of stationarity in the mean can be relaxed within the linear mixed model framework, using residual maximum likelihood to estimate variance parameters for the random effects. The Best Linear Unbiased Predictor (BLUP) is equivalent to the kriging predictor in these circumstances. However, nonstationarity in the variance is a harder problem to solve. Stationarity assumptions are necessary if the spatial covariance of a random process is to be estimated from the single realization which nature provides. However, they are not always plausible for variables arising from processes in complex landscapes across contrasting topography and geology.     

Spatial designs and properties of spatial correlation: Effects on covariance estimation

In a spatial regression context, scientists are often interested in a physical interpretation of components of the parametric covariance function. For example, spatial covariance parameter estimates in ecological settings have been interpreted to describe spatial heterogeneity or “patchiness” in a landscape that cannot be explained by measured covariates. In this article, we investigate the influence of the strength of spatial dependence on maximum likelihood (ML) and restricted maximum likelihood (REML) estimates of covariance parameters in an exponential-with-nugget model, and we also examine these influences under different sampling designs—specifically, lattice designs and more realistic random and cluster designs—at differing intensities of sampling (n=144 and 361). We find that neither ML nor REML estimates perform well when the range parameter and/or the nugget-to-sill ratio is large—ML tends to underestimate the autocorrelation function and REML produces highly variable estimates of the autocorrelation function. The best estimates of both the covariance parameters and the autocorrelation function come under the cluster sampling design and large sample sizes. As a motivating example, we consider a spatial model for stream sulfate concentration.     

应用土壤质地预测干旱区葡萄园土壤饱和导水率空间分布

田间表层土壤饱和导水率的空间变异性是影响灌溉水分入渗和土壤水分再分布的主要因素之一,研究土壤饱和导水率的空间变化规律,有助于定量估计土壤水分的空间分布和设计农田的精准灌溉管理制度。为了探究应用其他土壤性质如质地、容重、有机质预测土壤饱和导水率空间分布的可行性,试验在7.6 hm2的葡萄园内,采用均匀网格25 m×25 m与随机取样相结合的方式,测定了表层(0～10 cm)土壤饱和导水率、粘粒、粉粒、砂粒、容重和有机质含量,借助经典统计学和地统计学,分析了表层土壤饱和导水率的空间分布规律、与土壤属性的空间相关性,并对普通克里格法、回归法和回归克里格法预测土壤饱和导水率空间分布的结果进行了对比。结果表明:1)土壤饱和导水率具有较强的变异性,平均值为1.64 cm/d,变异系数为1.17;2)表层土壤饱和导水率60%的空间变化是由随机性或小于取样尺度的空间变异造成;3)土壤饱和导水率与粘粒、粉粒、砂粒和有机质含量具有一定空间相关性,而与土壤容重几乎没有空间相关性;4)在中值区以土壤属性辅助的回归克里格法对土壤饱和导水率的预测精度较好,在低值和高值区其与普通克里格法表现类似。研究结果将为更好地描述土壤饱和导水率空间变异结构及更准确地预测其空间分布提供参考。     

Spatial prediction of soil properties using EBLUP with the Matérn covariance function

Spatial prediction with the presence of spatially dense ancillary variables has attracted research in pedometrics. While soil survey and analysis of soil properties are still expensive and time consuming, the secondary data can be made available on a dense grid for the whole area of interest. The main aim of using the ancillary data is to enhance prediction of soil properties by making use of the ancillary variables as covariates. Methods that can be used for this purpose are kriging with external drift, cokriging, regression kriging, and REML-EBLUP (Residual Maximum Likelihood-Empirical Best Linear Unbiased Predictor). Regression kriging is a sub-optimal method that has been utilised extensively because it is easy to use and has been shown empirically to perform as well as other methods. A statically sound method is REML-EBLUP. This paper examines the use of REML-EBLUP in combination with the Matérn covariance function for spatial prediction of soil properties. Methods for estimating parameters of the Matérn variogram using REML, and prediction with EBLUP are described. The prediction capability of REML-EBLUP, regression kriging, and ordinary kriging is compared for four datasets. Results show that although REML-EBLUP generally improves the prediction, the improvement is small compared with regression kriging. Thus, for practical applications regression kriging appears to be a robust method. REML-EBLUP is useful when the trend is strong, and the number of observations is small (< 200). We concluded that improvement in the prediction of soil properties does not rely on more sophisticated statistical methods, but rather on gathering more useful and higher quality data.     

Geostatistical modelling and mapping of nematode-based soil ecological quality indices in a polluted nature reserve

Nematodes are indicators of soil quality and soil health. Knowledge of the relationships between nematode-based soil quality indices and environmental properties is beneficial for assessing environmental threats on soil biota. This study evaluated the spatial distribution of nematode-based soil quality indices in a 23-ha heavy metal-polluted nature reserve using geostatistical methods. We expected that a selection of abiotic soil properties (pH and moisture, clay, organic matter, cadmium (Cd), and zinc (Zn) contents) could explain a significant portion of the spatial variation of the indices and that regression kriging could more accurately model their spatial distribution than ordinary kriging. A stratified simple random sampling scheme was used to select 80 locations where soil samples were taken to extract nematodes and derive the indices. The area had a distinct gradient in soil properties with Cd and Zn content ranging from 0.07 to 68.9 and 5.3 to 1329 mg kg^-1, respectively. Linear regression models were fitted to describe the relationships between the indices and soil properties. By also modelling the spatial correlation structure of regression residuals using spherical semivariograms, regression kriging was used to produce maps of the indices. The regression models explained between 21% and 44% of the total original variance in the indices. Soil pH was a significant explanatory variable in almost all cases, while heavy metal conent had a remarkably low effect. In some cases, the regression residuals had spatial structure. Independent validation indicated that in all cases, regression kriging performed slightly better because of having lower values of the root mean square prediction error and a mean prediction error closer to zero than ordinary kriging. This study showed the importance of soil properties in explaining the spatial distribution of biological soil quality indices in ecological risk assessment.     

Do model choice and sample ratios separately or simultaneously influence soil organic matter prediction?

This study was performed to examine the separate and simultaneous influence of predictive models’ choice alongside sample ratios selection in soil organic matter (SOM). The research was carried out in northern Morocco, characterized by relatively cold weather and diverse geological conditions. The dataset herein used accounted for 1591 soil samples, which were randomly split into the following ratios: 10% (～150 sample ratio), 20% (～250 sample ratio), 35% (～450 sample ratio), 50% (～600 sample ratio) and 95% (～1200 sample ratio). Models herein involved were ordinary kriging (OK), regression kriging (RK), multiple linear regression (MLR), random forest (RF), quantile regression forest (QRF), Gaussian process regression (GPR) and an ensemble model. The findings in the study showed that the accuracy of SOM prediction is sensitive to both predictive models and sample ratios. OK combined with 95% sample ratio performed equally to RF in conjunction with all the sample ratios, as the latter did not show much sensitivity to sample ratios. ANOVA results revealed that RF with a ～10% sample ratio could also be optimum for predicting SOM in the study area. In conclusion, the findings herein reported could be instrumental for producing cost-effective detailed and accurate spatial estimation of SOM in other sites. Furthermore, they could serve as a baseline study for future research in the region or elsewhere. Therefore, we recommend conducting series of simulation of all possible combinations between various predictive models and sample ratios as a preliminary step in soil organic matter prediction.     

The Bayesian Group Lasso for Confounded Spatial Data

Generalized linear mixed models for spatial processes are widely used in applied statistics. In many applications of the spatial generalized linear mixed model (SGLMM), the goal is to obtain inference about regression coefficients while achieving optimal predictive ability. When implementing the SGLMM, multicollinearity among covariates and the spatial random effects can make computation challenging and influence inference. We present a Bayesian group lasso prior with a single tuning parameter that can be chosen to optimize predictive ability of the SGLMM and jointly regularize the regression coefficients and spatial random effect. We implement the group lasso SGLMM using efficient Markov chain Monte Carlo (MCMC) algorithms and demonstrate how multicollinearity among covariates and the spatial random effect can be monitored as a derived quantity. To test our method, we compared several parameterizations of the SGLMM using simulated data and two examples from plant ecology and disease ecology. In all examples, problematic levels multicollinearity occurred and influenced sampling efficiency and inference. We found that the group lasso prior resulted in roughly twice the effective sample size for MCMC samples of regression coefficients and can have higher and less variable predictive accuracy based on out-of-sample data when compared to the standard SGLMM.Supplementary materials accompanying this paper appear online.     

A practical two‐step approach for mixed model‐based kriging,with an application to the prediction of soil organic carbon concentration

Soil scientists often use prediction models to obtain values at unsampled locations. The spatial variation in the soil is best captured by using the empirical best linear unbiased predictor (EBLUP) based on a restricted maximum likelihood (REML) approach that efficiently exploits available data on both mean trends and correlation structures. We proposed a practical two‐step implementation of the REML approach for model‐based kriging, exemplified by predicting soil organic carbon (SOC) concentrations in mineral soils in Estonia from the large‐scale digital soil map information and a previously established prediction model. The prediction model was a linear mixed model with soil type, physical clay content (particle size < 0.01 mm) and A‐horizon thickness as fixed effects and site, transect, plot, year, year‐transect random intercepts and site‐specific random slopes for clay content. We used only the site‐specific intercept EBLUPs for estimating spatial correlation parameters as they described most of the variation in the random effects (86.8%). Fitting an exponential correlation model to these EBLUPs resulted in an estimated range of 10.5 km and the estimated proportion of the variance from the nugget effect was 0.23. The results of a simulation study showed a downwards bias that decreased with sample size. The results were validated through an external dataset, resulting in root mean square errors (RMSE) of 1.06 and 1.07% for the two‐step approach for kriging and the model with only fixed effects (no kriging), respectively. These results indicate that using the two‐step approach for kriging may improve prediction.     

Selection of a Geostatistical Method to Interpolate Soil Properties of the State Crop Testing Fields using Attributes of a Digital Terrain Model

The three most common techniques to interpolate soil properties at a field scale—ordinary kriging (OK), regression kriging with multiple linear regression drift model (RK + MLR), and regression kriging with principal component regression drift model (RK + PCR)—were examined. The results of the performed study were compiled into an algorithm of choosing the most appropriate soil mapping technique. Relief attributes were used as the auxiliary variables. When spatial dependence of a target variable was strong, the OK method showed more accurate interpolation results, and the inclusion of the auxiliary data resulted in an insignificant improvement in prediction accuracy. According to the algorithm, the RK + PCR method effectively eliminates multicollinearity of explanatory variables. However, if the number of predictors is less than ten, the probability of multicollinearity is reduced, and application of the PCR becomes irrational. In that case, the multiple linear regression should be used instead.     

Prediction of soil salinity in the Yellow River Delta using geographically weighted regression

It is essential to determine the content and spatial distribution of soil salinity in a timely manner because soil salinization can cause land degradation on a regional scale. Geographically weighted regression (GWR) is a local regression method that can achieve the spatial extension of dependent variables based on the relationships between the dependent variables and environment variables and the spatial distances between the sample points and predicted locations. This study aimed to explore the feasibility of GWR in predicting soil salinity because the existing interpolation methods for soil salinity in the Yellow River Delta are still of low precision. Additionally, multiple linear regressions, cokriging and regression kriging were added to compare the accuracy of GWRs. The results showed that GWR predicted soil salinity with high accuracy. Furthermore, the accuracy was improved when compared to other methods. The root mean square error, correlation coefficient, regression coefficient and adjustment coefficients between the observed values and predicted values of the validation points were 0.31, 0.65, 0.57 and 0.42, respectively, which were better than that of other methods, indicating that GWR is an optimal method.     

基于磁感式土壤表观电导率空间变异性的插值方法比较

针对目前黄河三角洲地区存在的土壤盐渍障碍问题,以该地区典型地块为研究对象,运用磁感大地电导率仪(EM38),结合GIS与地统计学研究了不同样点密度下土壤表观电导率的空间变异特征,确定了最佳的空间插值模式,并采用偏差指数对各分布图层的空间相似性进行了评价。结果表明：不同样点密度下的土壤表观电导率均呈中等变异强度,并服从对数正态分布;各样点密度下的土壤表观电导率均表现为强空间相关性,其空间变异主要表现在小于10 m的田间尺度上,且对于预测精度,泛克立格>普通克立格>反距离权重>局部多项式。偏差指数法分析表明,各插值方法分布图的空间相似性均随样点密度的降低而下降;对于相同的样点密度减小比例下的空间信息保留度,泛克立格>普通克立格>反距离权重>局部多项式,采用泛克立格法可以在保证预测精度的基础上合理降低样点密度。该研究为黄河三角洲地区盐渍土地磁感式田间数据采集的合理密度确定以及优化空间插值模式选取提供了理论依据与技术参考。     

An Efficient Approach to Spatiotemporal Analysis and Modeling of Air Pollution Data

A statistically efficient approach is adopted for modeling spatial time-series of large data sets. The estimation of the main diagnostic tool such as the likelihood function in Gaussian spatiotemporal models is a cumbersome task when using extended spatial time-series such as air pollution. Here, using the Innovation Algorithm, we manage to compute it for many spatiotemporal specifications. These specifications refer to the spatial periodic-trend, the spatial autoregressive moving average, the spatial autoregressive integrated and fractionally integrated moving average Gaussian models. Our method is applied to daily pollutants over a large metropolitan area like Athens. In the applied part of our paper, we first diagnose temporal and spatial structures of data using non-likelihood based criteria, such as the empirical autocorrelation and covariance functions. Second, we use likelihood and non-likelihood based criteria to select a spatiotemporal model among various specifications. Finally, using kriging we regionalize the resulting parameter estimates of the best-fitted model in space at any unmonitored location in the Athens region. The results show that a specific autoregressive integrated moving average spatiotemporal model can optimally perform in within and out of spatial sample estimation. Supplemental materials for this article are available from the journal website.     

结合高光谱信息的土壤有机碳密度地统计模型

传统线性回归模型在借助光谱信息进行土壤属性预测时,通常忽略了土壤自身所具有的空间异质性和依赖性,并且未考虑模型残差的空间结构。针对以上不足,该文以江汉平原232个土壤样本为研究对象,以土壤反射光谱为辅助变量,采用偏最小二乘回归、普通克里格、协同克里格以及回归克里格分别构建土壤有机碳密度预测模型,选取决定系数(R~2)、均方根误差、标准差与预测均方根误差比(ratio of performance to deviation,RPD)对模型预测精度进行对比评价。结果显示,结合高光谱信息,且同时考虑残差空间结构的回归克里格模型表现优于其他模型,预测决定系数R~2为0.617,RPD为1.614。鉴于土壤光谱信息同时还具有测定简单、省时、无损等优点,因此土壤光谱是土壤有机碳密度空间插值的理想辅助因子。     

县域土壤质量数字制图方法比较

土壤质量研究几乎涵盖土壤研究的所有领域,土壤质量制图理论与方法是土壤质量研究的一项重要研究内容。该研究以北京市密云县为研究区,基于土壤质量评价最小数据集和指数和法计算的土壤质量指数,探究了在地学模型支持下区域土壤质量数字制图方法。研究设计了5种区域土壤质量数字制图方法,并比较了不同方法的空间数字制图精度。结果显示,目前广泛使用的基于参评指标空间插值结果的土壤质量数字制图方法精度最低、工序较繁琐,且无法反映研究区景观高度异质的特点;而基于计算后的土壤质量指数（soil quality index,SQI）,借助于地统计学方法的土壤质量数字制图方法相对比较科学合理,其中又以基于计算后的SQI和回归克里格法预测效果最好,均方根误差最小,仅为0.01897,相对于基于参评指标空间插值结果的土壤质量数字制图方法,精度相对提高率最大,达到50%以上。综合考虑空间制图精度、工序的繁简程度,在该研究设计的5种方法中基于计算的SQI和回归克里格法最佳,该法避免了地统计插值在景观高度异质区的应用局限性,预测结果与实际最为相符。     

Analysis of Genebank Evaluation Data by using Geostatistical Methods

Genebanks often characterize accessions based on evaluation trials. This paper evaluates geostatistical methods as a tool to increase the utility of evaluation data. These methods were selected to overcome limitations resulting from a relative lack of replication and the scarcity of standards or check varieties. The data employed in the present study comprise nine characteristics of spring and winter barley, evaluated mostly as ratings. Ratings with quasi-metric scales were transformed by using the folded exponential transformation. To estimate the genetic component of the total effect, we compared two methods: Method 1 whereby a variogram is fitted by non-linear regression, and subsequently the implied spatial correlation is embedded into a mixed model analysis, which estimates the genetic effect by Best Linear Unbiased Prediction (BLUP); and Method 2 where each data value is re-estimated by kriging to correct for spatial effects and then the corrected data are submitted to a mixed model analysis. For practical application we propose Method 1 (though occasionally we met convergence problems): Fit the short range of the empirical variogram, visually choose the suitable covariance model. Use this and the initial values from non-linear regression fit with the mixed model, fixing the spatial parts at their starting values from non-linear regression, and estimate genetic effects by BLUP by using the fitted mixed model. To improve performance, we recommend that more standard or check varieties be used and, wherever possible, replace rating scales with metric scales or free-percentage scales (without categories).     

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.     

The use of region partitioning to improve the representation of geo statistically mapped soil attributes

An attempt to improve the representation of a geo statistically mapped soil attribute, clay content of the surface soil, through partitioning of the study area into two new regions was made. A topographic boundary divided the study area into hill and plain regions. Possible global non-stationarity or non-stationarity within the two newly defined regions was dealt with through the use of intrinsic random functions (IRF) of order k. Cross-validation of generalized covariance functions suggested that ordinary kriging might also have been appropriate. Exponential variogram models were subsequently fitted to the experimental variograms for each region. IRF-k block kriging and ordinary block kriging were then used as the primary methods of estimation. Both IRF-k and ordinary kriging performed badly in the vicinity of the topographic boundary when global models were used. This discontinuity was removed, at the expense of the introduction of some additional edge effects, when the hill and plain regions were kriged using models appropriate to each region. Independent zero-order generalized covariance functions with nugget and linear terms and exponential variogram models produced similar representations of clay content within each region, when used with their respective estimators. Splitting the region resulted in a 6% reduction in mean absolute deviation and a 14% reduction in mean squared deviation of predicted clay contents compared with a global model.