Mapping depth to Pleistocene sand with Bayesian generalized linear geostatistical models

Spatial soil applications frequently involve binomial variables. If relevant environmental covariates are available, using a Bayesian generalized linear model (BGLM) might be a solution for mapping such discrete soil properties. The geostatistical extension, a Bayesian generalized linear geostatistical model (BGLGM), adds spatial dependence and is thus potentially better equipped. The objective of this work was to evaluate whether it pays off to extend from a BGLM to a BGLGM for mapping binary soil properties, evaluated in terms of prediction accuracy and modelling complexity. As motivating example, we mapped the presence/absence of the Pleistocene sand layer within 120 cm from the land surface in the Dutch province of Flevoland, using the BGLGM implementation in the R-package geoRglm. We found that BGLGM yields considerably better statistical validation metrics compared to a BGLM, especially with – as in our case – a large (n = 1,000) observation sample and few relevant covariates available. Also, the inferred posterior BGLGM parameters enable the quantification of spatial relationships. However, calibrating and applying a BGLGM is quite demanding with respect to the minimal required sample size, tuning the algorithm, and computational costs. We replaced manual tuning by an automated tuning algorithm (which eases implementing applications) and found a sample composition that delivers meaningful results within 50 h calculation time. With the gained insights and shared scripts spatial soil practitioners and researchers can – for their specific cases – evaluate if using BGLGM is feasible and if the extra gain is worth the extra effort.