Implementation Pipeline: A Spatial Data Example
Loading Required Libraries
The following libraries provide the computational stack required for spatial data processing, structural variogram modeling, BME interpolation, and geospatial visualization:
library(BMEmapping)
library(ggplot2)
library(sf)
library(gstat)
library(dplyr)
library(tidyr)
library(scales)
library(knitr)
library(gridExtra)
data("utsnowload")To demonstrate the practical application of the package, this
vignette utilizes the utsnowload dataset, which contains
detrended reliability-targeted design snow load (RTDSL) measurements
collected from 232 observation sites across Utah.
The dataset is partitioned to reflect typical heterogeneous spatial data collections:
- Observations 1 to 65: Configured as precise, exact measurements treated as hard data.
- Observations 66 to 232: Formatted as interval-valued bounds representing uncertain soft data.
- Spatial Coordinates: Defined explicitly by the
longitudeandlatitudefields. - Uncertainty Bounds: Quantified by the
lowerandupperinterval parameters defining the soft data limits.
The package example data is loaded and summarized below:
# Load the sample spatial dataset
data("utsnowload")
# Display baseline distributions and structural summaries
summary(utsnowload)
#> latitude longitude hard lower
#> Min. :37.02 Min. :-114.0 Min. :-1.408481 Min. :-3.3063
#> 1st Qu.:38.58 1st Qu.:-112.3 1st Qu.:-0.777562 1st Qu.:-2.4155
#> Median :40.00 Median :-111.7 Median :-0.426066 Median :-1.8106
#> Mean :39.72 Mean :-111.6 Mean :-0.381925 Mean :-1.7760
#> 3rd Qu.:40.77 3rd Qu.:-111.1 3rd Qu.: 0.005383 3rd Qu.:-1.1578
#> Max. :41.97 Max. :-109.1 Max. : 1.252338 Max. : 0.1011
#> NA's :165 NA's :67
#> upper
#> Min. :-0.6965
#> 1st Qu.: 0.0194
#> Median : 0.2961
#> Mean : 0.3423
#> 3rd Qu.: 0.6516
#> Max. : 1.4599
#> NA's :67Detailed documentation for the dataset, including variable descriptions and geographical metadata, can be accessed using the standard R help command:
Organizing Hard and Soft Data
A critical step in BME modeling involves separating the observations
into hard and soft components. The spatial coordinates associated with
hard and soft observations are extracted independently, together with
the corresponding observed values and interval bounds. For illustrative
purposes, a subset of the utsnowload dataset consisting 30
hard-data locations and 100 soft-data locations is selected. This
reduced dataset facilitates demonstration of the data organization and
visualization workflow while limiting computational demands and leaves
10 soft-data locations as prediction locations.
# hard data locations
ch <- utsnowload[1:30, c("longitude", "latitude")]
# hard data values
zh <- utsnowload[1:30, c("hard")]
# soft data locations
cs <- utsnowload[68:167, c("longitude", "latitude")]
# lower and upper bounds of soft data (intervals)
a <- utsnowload[68:167, c("lower")]
b <- utsnowload[68:167, c("upper")]The extracted information is subsequently organized into a BMEmapping object, which serves as the central data structure for posterior density estimation and BME prediction.
Visualizing Hard and Soft Data
Exploratory visualization provides important insight into the spatial coverage and uncertainty structure of the dataset prior to variogram estimation and spatial interpolation.
This visualization facilitates identification of spatial clustering, regions with sparse observations, and areas associated with substantial interval uncertainty. Such preliminary assessment is important because the quality and spatial distribution of both hard and soft observations strongly influence posterior density estimation and predictive performance.
BME prediction is demonstrated using a subset of the soft-data
locations from the utsnowload dataset as prediction
locations. The prediction locations, denoted by x_k, are
obtained by extracting the geographic coordinates of the selected target
locations. These coordinates are subsequently supplied to the BME
prediction functions for posterior density estimation and computation of
summary statistics such as the posterior mean and posterior mode. The
prediction locations used in this illustration are defined below:
BME Estimation Methods
Two BME prediction approaches are considered throughout this tutorial.
Conventional Bayesian Maximum Entropy (CBME)
The conventional Bayesian Maximum Entropy framework assumes that the soft observations follow a uniform distribution over the specified interval bounds. Under this assumption, analytical expressions for the interval mean and variance are incorporated into the covariance structure used during posterior estimation.
The CBME approach is computationally efficient and has been widely used in applications involving interval uncertainty. However, the uniform assumption may oversimplify uncertainty structures when the true distribution within the interval is asymmetric or nonuniform.
Quantile-Based Bayesian Maximum Entropy (QBME)
Spatial continuity is a foundational concept in geostatistics, capturing the degree to which values of a spatial variable are correlated across space to ensure accurate interpolation and simulation. Within the CBME framework, this spatial dependence is typically represented by a theoretical variogram model derived from precise physical measurements (hard data). However, in settings where hard data are sparse, practitioners frequently rely on a hybrid dataset that combines these precise measurements with the midpoints of interval-valued (soft data). While convenient, this approach introduces an important methodological gap: it implicitly assumes that a single midpoint adequately summarizes the underlying uncertainty of an entire interval. By ignoring the full span of the interval, this simplification can introduce systemic bias into the spatial covariance structure, ultimately degrading prediction precision and compromising the integrity of the localized uncertainty characterization.
To systematically address this limitation, we introduce the QBME framework, which models spatial continuity by explicitly accounting for the uncertainty inherent in soft-interval data. This approach discretizes each soft interval into a finite set of representative quantiles, which are individually paired with the available hard data to construct a sequence of distinct hybrid datasets. The framework fits a unique variogram model to each hybrid configuration to capture spatial dependencies across the full spectrum of plausible values. The resulting quantile-specific covariance matrices are then averaged into a single integrated covariance matrix that explicitly propagates the soft-data uncertainty throughout the entire domain. By shifting from a single midpoint approximation to an integrated ensemble of conditional variograms, the QBME implementation preserves spatial coherence in a principled, computationally tractable manner that remains strictly consistent with its information-theoretic foundations. A detailed description of the QBME approach and its applications is provided in https://doi.org/10.1016/j.spasta.2026.100974
CBME Implementation
Variogram Estimation
Specification of an appropriate variogram model is a fundamental component of CBME interpolation because the covariance structure determines the spatial dependence relationships used during posterior estimation. Because the number of hard observations is relatively limited, a combined dataset consisting of hard observations and representative values derived from the soft intervals is commonly used for variogram estimation. In this tutorial, the midpoints of the soft intervals are combined with the hard observations to provide a practical approximation of the underlying spatial variability.
The empirical semivariogram and its corresponding theoretical
covariance model are estimated using the gstat package.
This step converts the combined hard observations and soft-interval
midpoints into an explicit model of spatial dependence:
df <- data.frame(rbind(ch, cs), z = c(zh, (a + b) / 2))
sf_data <- sf::st_as_sf(df, coords = c("longitude", "latitude"))
vg <- variogram(z ~ 1, data = sf_data)
vg_model <- fit.variogram(vg, model = vgm(c("Exp", "Sph")))
vg_model
#> model psill range
#> 1 Nug 0.05502877 0.0000000
#> 2 Exp 0.33581773 0.6471111Posterior Density Estimation
One of the principal advantages of the BME framework is the ability to estimate complete posterior probability density functions rather than producing only single-point predictions. Posterior densities can be estimated at selected prediction locations using the BME prediction functions provided in the package.
Before performing full-scale BME prediction, it is often instructive
to examine the posterior distributions at a few selected locations.
Although this step is not strictly necessary and can be skipped, it
provides valuable insights into the shape, spread, and asymmetry of the
posterior distributions. Moreover, it allows users to verify that the
chosen domain for evaluation (zk_range) sufficiently
encompasses the possible values of the target variable and to make
adjustments if necessary.
The posterior density functions for selected prediction locations can be evaluated and visualized using the following code:
# Extract spatial parameters
model <- as.character(vg_model[2, 1])
nugget <- vg_model[1, 2]
sill <- vg_model[2, 2]
range <- vg_model[2, 3]
# default zk_range
p_1 <- prob_zk(xk[1,], data_object, model, nugget, sill, range)
p_2 <- prob_zk(xk[2,], data_object, model, nugget, sill, range)
p_3 <- prob_zk(xk[3,], data_object, model, nugget, sill, range)
p_4 <- prob_zk(xk[4,], data_object, model, nugget, sill, range)
p_df <- cbind.data.frame(p_1, p_2[, 2], p_3[, 2], p_4[, 2])
names(p_df) <- c("zk_i", "p1", "p2", "p3", "p4")
# Function to generate ggplot for a given column name
plot_prob_curve <- function(df, pi_col) {
ggplot(df, aes(x = zk_i, y = .data[[pi_col]])) +
geom_line(color = "darkblue", linewidth = 0.7) +
labs(x = "Spatial Field Value (z)", y = "Posterior Density f(z)") +
theme_minimal() +
theme(
panel.background = ggplot2::element_rect(fill = "white", color = "black")
)
}
# Generate individual plots
p1 <- plot_prob_curve(p_df, "p1")
p2 <- plot_prob_curve(p_df, "p2")
p3 <- plot_prob_curve(p_df, "p3")
p4 <- plot_prob_curve(p_df, "p4")
# Arrange in a 2x2 grid
grid.arrange(p1, p2, p3, p4, ncol = 2)
Inspection of the posterior density plots above reveals that the
estimated density functions reach zero at certain sub-regions within the
default zk_range, for example \(z
< -2\) and \(z > 1.5\).
This indicates that the true support of the posterior distribution is
likely much narrower than the initially chosen default range. To improve
the accuracy and computational efficiency of the BME predictions, it is
advisable to refine the posterior domain by excluding regions where
\(f(z) = 0\), focusing on the interval
where the posterior density is strictly positive. In this example, the
refined range can be set to [-2, 2], which better captures
the plausible values of the target variable. Using the updated range,
the posterior densities at the selected prediction locations are
recomputed as shown below:
# updated zk_range: [-2, 2]
q_1 <- prob_zk(xk[1,], data_object, model, nugget, sill, range, zk_range = c(-2, 2))
q_2 <- prob_zk(xk[2,], data_object, model, nugget, sill, range, zk_range = c(-2, 2))
q_3 <- prob_zk(xk[3,], data_object, model, nugget, sill, range, zk_range = c(-2, 2))
q_4 <- prob_zk(xk[4,], data_object, model, nugget, sill, range, zk_range = c(-2, 2))
q_df <- cbind.data.frame(q_1, q_2[, 2], q_3[, 2], q_4[, 2])
names(q_df) <- c("zk_i", "q1", "q2", "q3", "q4")
# Generate individual plots
q1 <- plot_prob_curve(q_df, "q1")
q2 <- plot_prob_curve(q_df, "q2")
q3 <- plot_prob_curve(q_df, "q3")
q4 <- plot_prob_curve(q_df, "q4")
grid.arrange(q1, q2, q3, q4, ncol = 2)
Spatial Prediction
Predictions within the BME framework can be computed using the posterior mean, posterior median, or posterior mode, depending on the desired interpretation and application:
- Posterior Mean: Represents the expected value of the posterior probability density function and provides an optimal point estimate that minimizes mean squared prediction error (MSPE).
- Posterior Median: Identifies the 50th percentile of the localized distribution, providing an alternative central tendency metric that is robust to extreme outliers and heavy tails.
- Posterior Mode: Corresponds to the value associated with the highest localized probability density, tracking the most likely spatial state when the conditional distribution is asymmetric or multimodal.
In addition to point predictions, the framework quantifies prediction reliability and spatial confidence through full uncertainty propagation:
- Posterior Variance: Quantifies the spread of the posterior density at each location, reflecting the proximity and information quality of neighboring hard and soft observations.
- Credible Intervals: Computed directly from the integrated posterior probability density function to define explicit probabilistic bounds (e.g., a 90% credible interval) that contain the true value with a specified probability, capturing localized non-Gaussian asymmetries that standard kriging variance intervals overlook.
Together, these point estimators, variance fields, and credible intervals offer a comprehensive probabilistic characterization of the spatial prediction process. The posterior point estimators across the unobserved target locations are executed below.
Posterior Mode Point Predictions
The posterior mode captures the peak of the localized probability density function:
CBME_mode <- bme_predict(xk, data_object, model, nugget, sill, range,
n = 100, zk_range = c(-2, 2), type = "mode")
head(CBME_mode)
#> longitude latitude mode
#> 201 -111.1472 40.9383 -0.1818
#> 202 -111.7836 40.6189 -0.2222
#> 203 -111.5289 40.4044 -0.2222
#> 204 -111.4336 40.9656 -0.0202
#> 205 -111.4164 39.7483 -0.1414Posterior Mean Point Predictions
The posterior mean minimizes the mean squared prediction error across the field:
CBME_mean <- bme_predict(xk, data_object, model, nugget, sill, range,
n = 100, zk_range = c(-2, 2), type = "mean")
head(CBME_mean)
#> longitude latitude mean variance
#> 201 -111.1472 40.9383 -0.1708 0.3590
#> 202 -111.7836 40.6189 -0.2031 0.3160
#> 203 -111.5289 40.4044 -0.2050 0.3499
#> 204 -111.4336 40.9656 -0.0161 0.3693
#> 205 -111.4164 39.7483 -0.1516 0.3737Posterior Median Point Predictions
The posterior median provides an alternative central tendency metric robust to distribution skewness:
CBME_median <- bme_predict(xk, data_object, model, nugget, sill, range,
n = 100, zk_range = c(-2, 2), type = "median")
head(CBME_median)
#> longitude latitude median
#> 201 -111.1472 40.9383 -0.1924
#> 202 -111.7836 40.6189 -0.2238
#> 203 -111.5289 40.4044 -0.2267
#> 204 -111.4336 40.9656 -0.0364
#> 205 -111.4164 39.7483 -0.1733Probabilistic Uncertainty Bounds (90% Credible Intervals)
To construct explicit probability bounds around our predictions, a 90% credible interval is integrated directly from the localized conditional probability density fields:
CBME_interval <- bme_predict_ci(xk, data_object, model, nugget, sill, range,
n = 100, zk_range = c(-2, 2), level = 0.90)
head(CBME_interval)
#> longitude latitude lower_90 upper_90
#> 201 -111.1472 40.9383 -1.1801 0.7996
#> 202 -111.7836 40.6189 -1.1505 0.7034
#> 203 -111.5289 40.4044 -1.2016 0.7522
#> 204 -111.4336 40.9656 -1.0396 0.9673
#> 205 -111.4164 39.7483 -1.1814 0.8393The QBME Implementation
Unlike the classical approach, the QBME implementation completely
bypasses manual variogram modeling. Spatial continuity is internally and
adaptively captured through an ensemble of quantile-specific variogram
estimations across different levels of the data distribution.
Consequently, the user-specified spatial parameters (model,
nugget, sill, and range) are
omitted entirely from the function calls. Instead, the user provides a
key configuration argument: nq, which defines the number of
discrete quantile levels used to partition and approximate the
soft-interval bounds.
Despite these distinct core structural changes in structural parameter estimation, the QBME framework produces the identical suite of localized posterior characterizations, point estimators, and uncertainty measures as the CBME workflow.
Posterior Distribution Estimation at Individual Locations
Just as with the CBME approach, it is highly recommended to inspect
localized posterior distribution curves at targeted test locations to
pick an optimal support window (zk_range). The conditional
density field can be extracted using the quantile-specific utility
function.
(Note: The following chunks illustrate the implementation syntax and are not evaluated here).
Spatial Prediction and Uncertainty Propagation
Point predictions and probabilistic interval bounds are generated
using the q_bme_predict() and
q_bme_predict_ci() engines. The implementation syntax
mirrors the classical execution, requiring only the removal of the
explicit variogram parameters and the inclusion of the nq
argument.
(Note: The code chunks below demonstrate the syntax configuration for the point estimators and credible intervals and are not evaluated here).
# Compute posterior mode point predictions
QBME_mode <- q_bme_predict(xk, data_object, n = 100, type = "mode", nq = 8)
# Compute posterior mean point predictions
QBME_mean <- q_bme_predict(xk, data_object, n = 100, type = "mean", nq = 8)
# Compute posterior median point predictions
QBME_median <- q_bme_predict(xk, data_object, n = 100, type = "median", nq = 8)
# Integrated 90% probabilistic credible interval bounds
QBME_int <- q_bme_predict_ci(xk, data_object, n = 100, level = 0.90, nq = 8)Cross-Validation Systems & Diagnostic Testing
To rigorously evaluate the predictive performance and structural
integrity of both interpolation frameworks, the package provides
dedicated cross-validation engines: bme_cv() for the
classical framework (CBME) and q_bme_cv() for the
quantile-based framework (QBME). These functions execute internal
validation routines exclusively at hard data locations, supporting both
flexible K-fold data partitioning and exact Leave-One-Out
Cross-Validation (LOOCV). By setting the fold parameter k
equal to the total number of hard observations, the engine automatically
executes a full LOOCV routine. During each validation step, the targeted
hard observation is systematically omitted from the conditioning
neighborhood, a localized estimation is generated using the remaining
hard and soft spatial datasets, and the result is verified against the
true, measured baseline value.
This resampling approach maximizes data utility, making it highly effective for small or heterogeneous spatial datasets where separating data into fixed training and validation subsets would severely degrade the sample size. Predictive performance across both validation routines is comprehensively quantified using standard geostatistical error metrics:
- Mean Error (ME): Evaluates systemic model bias; values near zero indicate a lack of global over- or under-estimation.
- Mean Absolute Error (MAE): Measures the average magnitude of absolute prediction residuals, providing a robust baseline of average error absolute limits.
- Root Mean Square Error (RMSE): Quantifies prediction precision, penalizing larger outliers more heavily to capture severe spatial variance deviations.
- Coefficient of Determination (\(R^2\)): Represents the proportion of total spatial variance successfully captured and accounted for by the BME model.
The resulting validation objects inherit customized S3 methods
(summary() and plot()) allowing users to
instantaneously compute standard geostatistical error metrics or render
diagnostic plots to evaluate model adequacy through visual assessment of
residual behavior. These plots facilitate identification of deviations
from assumptions including spatial independence, homoscedasticity, and
stationarity.
There is a small discrepancy in the original snippet: the text
mentions a 5-fold cross-validation, but the code argument is set to
k = 10 (10-fold). The version below fixes this parameter
discrepancy to ensure the prose matches the code execution exactly.
Executing a 5-fold cross-validation (\(k = 5\)) using the CBME framework is performed as follows:
# Execute 5-fold cross-validation across the hard data locations
CBME_cv <- bme_cv(data_object, model, nugget, sill, range, n = 100,
zk_range = c(-2, 2), type = "mean", k = 5)
CBME_cv
#> longitude latitude observed mean variance residual fold
#> 1 -112.24 40.44 0.09696012 -0.1989 0.3419 0.2959 5
#> 2 -112.41 39.94 0.12258678 -0.2824 0.3225 0.4050 4
#> 3 -113.40 37.51 -0.02302358 -0.0446 0.3435 0.0216 4
#> 4 -113.85 37.49 0.50354362 -0.1112 0.3571 0.6147 3
#> 5 -109.53 39.31 -0.68611327 -0.0224 0.3846 -0.6637 5
#> 6 -109.54 40.72 -0.53000397 -0.5104 0.2943 -0.0196 3
#> 7 -109.89 40.61 -0.71923519 -0.6573 0.3030 -0.0619 2
#> 8 -109.96 40.91 -1.31503404 -0.8351 0.2640 -0.4799 1
#> 9 -109.67 40.74 -0.94879597 -0.5384 0.2770 -0.4104 5
#> 10 -110.19 40.92 -1.39798035 -0.7088 0.3051 -0.6892 5
#> 11 -110.48 40.95 -1.21900906 -0.8769 0.1880 -0.3421 3
#> 12 -110.43 40.60 -1.24787225 -0.7845 0.2526 -0.4634 4
#> 13 -110.69 40.55 -0.55027484 -0.6149 0.2517 0.0646 4
#> 14 -110.50 40.91 -1.06708711 -1.0017 0.1752 -0.0654 2
#> 15 -110.47 40.72 -1.14044998 -0.9098 0.2328 -0.2306 3
#> 16 -110.59 40.58 -0.94551554 -0.7337 0.2316 -0.2118 1
#> 17 -110.80 40.86 -0.83840015 -0.5612 0.2612 -0.2772 2
#> 18 -110.01 40.77 -1.24671792 -0.7743 0.2663 -0.4724 5
#> 19 -110.88 40.80 -0.65036211 -0.5307 0.2355 -0.1197 5
#> 20 -110.95 40.68 -0.37127802 -0.3648 0.2984 -0.0065 2
#> 21 -110.75 39.89 -0.80367306 -0.2290 0.3538 -0.5747 4
#> 22 -110.99 39.96 -0.54230365 -0.2712 0.3529 -0.2711 1
#> 23 -111.94 41.38 0.94099563 0.1522 0.2600 0.7888 3
#> 24 -111.45 41.31 0.24796667 0.0401 0.3035 0.2079 1
#> 25 -111.83 41.41 0.47642403 0.6251 0.2475 -0.1487 4
#> 26 -111.92 41.38 1.25233814 0.1854 0.2473 1.0669 3
#> 27 -111.63 41.90 0.61655171 0.0162 0.3383 0.6004 1
#> 28 -111.42 41.68 0.18443361 -0.0237 0.3176 0.2081 2
#> 29 -111.54 41.41 0.11223798 0.1657 0.2170 -0.0535 2
#> 30 -111.50 41.47 0.10561343 0.0817 0.2297 0.0239 1Key summary performance metrics can be extracted using the S3 summary method:
# Extract comprehensive cross-validation error metrics
summary(CBME_cv)
#> Summary of BME Cross-Validation
#> for hard data locations
#> --------------------------------
#> Number of predictions: 30
#>
#> ME -0.0421
#> MAE 0.3287
#> RMSE 0.4218
#> Rsquared 0.6506Diagnostic plots are generated via the S3 plot method to visually inspect residual patterns:
Similarly, diagnostic cross-validation can be executed within the
QBME framework using q_bme_cv(). The code chunk below
demonstrates how to perform a full Leave-One-Out Cross-Validation
(LOOCV) by setting the fold parameter k equal to the total
number of hard observations (nrow(ch)).
(Note: The following chunks illustrate the implementation syntax and are not evaluated here).
# Execute Leave-One-Out Cross-Validation (LOOCV) under the QBME framework
QBME_cv <- q_bme_cv(data_object, n = 100, nq = 8, type = "mean", k = nrow(ch))
QBME_cvJust as with the classical framework, the summary performance statistics and diagnostic residual plots are easily extracted using standard S3 methods:
Spatial Visualization of Prediction Fields
To complement the numerical predictions, the package provides
comprehensive spatial visualization tools built directly into its
architecture. The standard S3 plot() method has been
extended to accept prediction collection objects from both the CBME and
QBME workflows (e.g., outputs from bme_predict() or
q_bme_predict()). Leveraging this function allows users to
instantly transform tabular spatial predictions into continuous
geographic maps. This visual diagnostic tool is crucial for analyzing
macro-scale geographic prediction trends, identifying spatial anomalies,
and mapping localized hotspots across your targeted study domain.
Furthermore, because the BME framework tracks full uncertainty propagation, the plotting engine can map both the estimated point prediction surfaces (mean, median, or mode) and their corresponding spatial uncertainty metrics, such as the posterior variance. This dual visualization provides a complete spatial perspective on both the predicted state of the field and the localized reliability of those predictions.
The continuous spatial prediction surface can be rendered using the following implementation syntax:
