| Title: | Singularity Regression Kriging for Spatial Prediction |
|---|---|
| Description: | Implements the Singularity Regression Kriging ('SRK') model for spatial prediction by integrating covariate singularity feature construction, nonlinear trend estimation via random forest, and geostatistical interpolation of residuals using ordinary kriging. Singularity-based anomaly indices are computed from environmental covariates at multiple spatial scales to capture local multiscale heterogeneity and augment the random forest feature set for trend estimation. The resulting residuals are interpolated using ordinary kriging to generate final spatial predictions with uncertainty quantification. Tools for spatial block cross-validation, parameter sensitivity analysis, and diagnostic visualization are also provided. Methods are based on Ren, Song, Chen, and Yu (2026) <doi:10.1080/15481603.2026.2690341>, with singularity theory from Cheng (2012) <doi:10.1016/j.gexplo.2012.07.007> and Cheng (2017) <doi:10.1016/j.gr.2017.07.011>, random forest methodology from Breiman (2001) <doi:10.1023/A:1010933404324>, and regression kriging framework from Hengl, Heuvelink, and Rossiter (2007) <doi:10.1016/j.cageo.2007.05.001>. |
| Authors: | Shikhar Tyagi [aut, cre] (ORCID: <https://orcid.org/0000-0003-1606-0844>), Arvind Pandey [aut], Bhupendra Singh [aut], Vrijesh Tripathi [aut] |
| Maintainer: | Shikhar Tyagi <shikhar1093tyagi@gmail.com> |
| License: | GPL (>= 3) |
| Version: | 0.1.0 |
| Built: | 2026-10-08 10:14:29 UTC |
| Source: | https://github.com/cran/SingRegKrig |
Computes singularity-based anomaly indices from environmental covariates at multiple spatial scales. For each location, the singularity index characterises how the local covariate intensity changes with spatial scale, following the singularity theory of Cheng (2012, 2017).
compute_singularity( coords, values, scales = NULL, min_neighbours = 3L, min_scales = 2L )compute_singularity( coords, values, scales = NULL, min_neighbours = 3L, min_scales = 2L )
coords |
A numeric matrix or data frame with two columns representing
spatial coordinates (e.g., easting and northing). Coordinate units must
match the units used in |
values |
A numeric vector of covariate values at the locations
specified by |
scales |
A numeric vector of spatial scales at which to compute the
singularity index. Scales should be positive and in the same units as
|
min_neighbours |
Integer. Minimum number of valid neighbouring samples required within a scale window for that scale to be used. Default is 3. |
min_scales |
Integer. Minimum number of valid scales required to estimate the singularity index at a location. Locations with fewer valid scales are assigned the neutral value of 2. Default is 2. |
The singularity index is computed following the method described in
Ren et al. (2026) and based on the singularity theory of
Cheng (2012, 2017, 2018). For each location and spatial
scale , the local covariate intensity is
estimated as the mean of the absolute covariate values within a square
window of half-width centred at . The singularity index
is then estimated as the slope of the ordinary
least-squares regression of on , plus 2:
where is the OLS slope from regressing
on across the selected scales.
A value of corresponds to a spatially uniform field;
indicates local enrichment; indicates
local depletion.
A numeric vector of singularity indices with the same length as the
number of rows in coords. A value of 2 indicates a spatially
uniform field (neutral reference); values less than 2 indicate local
enrichment (positive anomaly); values greater than 2 indicate local
depletion.
Cheng, Q. (2012). Singularity theory and methods for mapping geochemical anomalies caused by buried sources and for predicting undiscovered mineral deposits in covered areas. Journal of Geochemical Exploration, 122, 55–70. doi:10.1016/j.gexplo.2012.07.007
Cheng, Q. (2017). Singularity analysis of global zircon U-Pb age series and implication of continental crust evolution. Gondwana Research, 51, 51–63. doi:10.1016/j.gr.2017.07.011
Ren, K., Song, Y., Chen, M., and Yu, Q. (2026). A singularity regression kriging for spatial prediction. GIScience & Remote Sensing, 63(1), 2690341. doi:10.1080/15481603.2026.2690341
set.seed(42) coords <- as.matrix(expand.grid(x = 1:10, y = 1:10)) values <- rnorm(100) alpha <- compute_singularity(coords, values, scales = c(1, 2, 3, 4, 5)) summary(alpha) # Uniform field returns singularity = 2 alpha_uniform <- compute_singularity(coords, rep(1, 100), scales = c(1, 2, 3, 4, 5)) all(abs(alpha_uniform - 2) < 1e-10)set.seed(42) coords <- as.matrix(expand.grid(x = 1:10, y = 1:10)) values <- rnorm(100) alpha <- compute_singularity(coords, values, scales = c(1, 2, 3, 4, 5)) summary(alpha) # Uniform field returns singularity = 2 alpha_uniform <- compute_singularity(coords, rep(1, 100), scales = c(1, 2, 3, 4, 5)) all(abs(alpha_uniform - 2) < 1e-10)
Produces diagnostic plots for a fitted SRK model, including residual histogram, variable importance, observed vs fitted, and variogram plots.
## S3 method for class 'srk' plot(x, which = c(1L, 2L, 3L, 4L), ...)## S3 method for class 'srk' plot(x, which = c(1L, 2L, 3L, 4L), ...)
x |
An object of class |
which |
Integer vector specifying which plots to produce.
|
... |
Additional arguments passed to plotting functions. |
Invisibly returns NULL.
sim_data <- sim_srk_data(n_side = 10, scenario = "normal", seed = 42) model <- srk(z ~ cov1, data = sim_data, coords = ~x + y, scales = c(1, 2, 3, 4, 5), ntree = 100) plot(model, which = 3)sim_data <- sim_srk_data(n_side = 10, scenario = "normal", seed = 42) model <- srk(z ~ cov1, data = sim_data, coords = ~x + y, scales = c(1, 2, 3, 4, 5), ntree = 100) plot(model, which = 3)
Produces an observed vs predicted scatter plot from spatial block cross-validation.
## S3 method for class 'srk_cv' plot(x, ...)## S3 method for class 'srk_cv' plot(x, ...)
x |
An object of class |
... |
Additional arguments passed to |
Invisibly returns NULL.
Produces a heatmap of performance metrics across parameter combinations.
## S3 method for class 'srk_sensitivity' plot(x, metric = "R2", ...)## S3 method for class 'srk_sensitivity' plot(x, metric = "R2", ...)
x |
An object of class |
metric |
Character. Metric to plot. One of |
... |
Additional arguments passed to |
Invisibly returns NULL.
Generates spatial predictions at new locations using a fitted SRK model. Computes singularity features at prediction locations, evaluates the random forest trend, and performs ordinary kriging of residuals.
## S3 method for class 'srk' predict(object, newdata, ...)## S3 method for class 'srk' predict(object, newdata, ...)
object |
An object of class |
newdata |
A data frame containing covariate and coordinate columns at
prediction locations. Column names must match those used in the original
|
... |
Additional arguments (currently unused). |
The prediction process follows three steps:
Singularity features are computed at prediction locations using combined training and prediction covariate data for neighbourhood completeness.
The random forest trend model is evaluated at prediction locations using the augmented feature set.
Residuals from training locations are interpolated to prediction locations via ordinary kriging, providing both predictions and uncertainty estimates.
A data frame with columns:
First coordinate of prediction location.
Second coordinate of prediction location.
Final SRK prediction (trend plus kriged residual).
Random forest trend estimate.
Ordinary kriging interpolation of residuals.
Kriging standard error (prediction uncertainty).
Ren, K., Song, Y., Chen, M., and Yu, Q. (2026). A singularity regression kriging for spatial prediction. GIScience & Remote Sensing, 63(1), 2690341. doi:10.1080/15481603.2026.2690341
set.seed(42) train <- sim_srk_data(n_side = 10, scenario = "normal", seed = 42) model <- srk(z ~ cov1, data = train, coords = ~x + y, scales = c(1, 2, 3, 4, 5), ntree = 100) # Create prediction locations pred_locs <- expand.grid(x = seq(2, 9, by = 1), y = seq(2, 9, by = 1)) pred_locs$cov1 <- rnorm(nrow(pred_locs), mean = mean(train$cov1)) preds <- predict(model, newdata = pred_locs) head(preds)set.seed(42) train <- sim_srk_data(n_side = 10, scenario = "normal", seed = 42) model <- srk(z ~ cov1, data = train, coords = ~x + y, scales = c(1, 2, 3, 4, 5), ntree = 100) # Create prediction locations pred_locs <- expand.grid(x = seq(2, 9, by = 1), y = seq(2, 9, by = 1)) pred_locs$cov1 <- rnorm(nrow(pred_locs), mean = mean(train$cov1)) preds <- predict(model, newdata = pred_locs) head(preds)
Print an SRK Model Object
## S3 method for class 'srk' print(x, ...)## S3 method for class 'srk' print(x, ...)
x |
An object of class |
... |
Additional arguments (currently unused). |
Invisibly returns x.
sim_data <- sim_srk_data(n_side = 10, scenario = "normal", seed = 42) model <- srk(z ~ cov1, data = sim_data, coords = ~x + y, scales = c(1, 2, 3, 4, 5), ntree = 100) print(model)sim_data <- sim_srk_data(n_side = 10, scenario = "normal", seed = 42) model <- srk(z ~ cov1, data = sim_data, coords = ~x + y, scales = c(1, 2, 3, 4, 5), ntree = 100) print(model)
Print Cross-Validation Results
## S3 method for class 'srk_cv' print(x, ...)## S3 method for class 'srk_cv' print(x, ...)
x |
An object of class |
... |
Additional arguments (currently unused). |
Invisibly returns x.
Print Sensitivity Analysis Results
## S3 method for class 'srk_sensitivity' print(x, ...)## S3 method for class 'srk_sensitivity' print(x, ...)
x |
An object of class |
... |
Additional arguments (currently unused). |
Invisibly returns x.
Print Summary of an SRK Model
## S3 method for class 'summary.srk' print(x, ...)## S3 method for class 'summary.srk' print(x, ...)
x |
An object of class |
... |
Additional arguments (currently unused). |
Invisibly returns x.
Generates simulated spatial datasets on a regular grid with spatially autocorrelated fields and an associated covariate, following the simulation design described in Ren et al. (2026, Section 2.3).
sim_srk_data( n_side = 20L, scenario = c("normal", "skewed", "longtail"), seed = 42L )sim_srk_data( n_side = 20L, scenario = c("normal", "skewed", "longtail"), seed = 42L )
n_side |
Integer. Number of grid points per side.
Total sample size is |
scenario |
Character. Distribution scenario for the response variable.
One of |
seed |
Integer. Random seed for reproducibility. Default is 42. |
The data generation process follows the simulation design in Section 2.3.1 of Ren et al. (2026):
A spatially autocorrelated base field is generated using a spherical covariance model on a regular grid.
A covariate field is constructed to maintain a strong but imperfect spatial association with the base field.
The response variable is derived from the base field under
one of three distributional scenarios:
"normal" (linear shift),
"skewed" (exponential transformation),
or "longtail" (stronger exponential transformation).
A data frame with columns:
Horizontal grid coordinate.
Vertical grid coordinate.
Response variable.
Spatially autocorrelated covariate.
Ren, K., Song, Y., Chen, M., and Yu, Q. (2026). A singularity regression kriging for spatial prediction. GIScience & Remote Sensing, 63(1), 2690341. doi:10.1080/15481603.2026.2690341
# Normal distribution scenario sim_normal <- sim_srk_data(n_side = 10, scenario = "normal", seed = 42) head(sim_normal) # Skewed distribution scenario sim_skew <- sim_srk_data(n_side = 10, scenario = "skewed", seed = 42) summary(sim_skew$z)# Normal distribution scenario sim_normal <- sim_srk_data(n_side = 10, scenario = "normal", seed = 42) head(sim_normal) # Skewed distribution scenario sim_skew <- sim_srk_data(n_side = 10, scenario = "skewed", seed = 42) summary(sim_skew$z)
Fits the Singularity Regression Kriging (SRK) model for spatial prediction. SRK integrates covariate singularity feature construction, nonlinear trend estimation via random forest, and geostatistical interpolation of residuals using ordinary kriging.
srk( formula, data, coords, pred_data = NULL, scales = NULL, ntree = 500L, sd_threshold = 0.5, min_neighbours = 3L, min_scales = 2L, variogram_model = "auto" )srk( formula, data, coords, pred_data = NULL, scales = NULL, ntree = 500L, sd_threshold = 0.5, min_neighbours = 3L, min_scales = 2L, variogram_model = "auto" )
formula |
A formula specifying the response variable and covariates,
e.g., |
data |
A data frame containing the response, covariates, and coordinate
columns referenced in |
coords |
A one-sided formula specifying the coordinate columns in
|
pred_data |
An optional data frame containing covariate and coordinate
columns at prediction locations. Column names must match those in
|
scales |
A numeric vector of spatial scales for singularity computation.
If |
ntree |
Integer. Number of trees in the random forest. Default is 500. |
sd_threshold |
Numeric. Minimum standard deviation for a singularity feature to be retained. Features with lower standard deviation are excluded. Default is 0.5. |
min_neighbours |
Integer. Minimum neighbours for singularity
computation. Default is 3. See |
min_scales |
Integer. Minimum scales for singularity computation.
Default is 2. See |
variogram_model |
Character. Variogram model type for residual kriging.
Options are |
The SRK model follows three main steps as described in Ren et al. (2026):
Singularity feature construction: For each covariate,
singularity indices are computed at multiple spatial scales using
compute_singularity. Features with standard deviation
below sd_threshold are discarded.
Trend estimation: A random forest (Breiman, 2001) is trained on the augmented feature set (original covariates plus retained singularity features) to estimate the spatial trend.
Residual kriging: The residuals from the random forest trend are interpolated using ordinary kriging with an automatically fitted variogram model to generate final predictions with uncertainty quantification.
An object of class "srk" containing:
The matched function call.
The model formula.
Names of coordinate columns.
Names of covariate columns.
Matrix of training coordinates.
Data frame of training covariates.
Numeric vector of observed response values.
Numeric vector of fitted values at training locations.
Numeric vector of random forest trend estimates at training locations.
Numeric vector of residuals (observed minus trend).
The fitted randomForest object.
The empirical variogram of residuals.
The fitted variogram model.
A data frame of singularity indices at training locations.
Character vector of retained singularity feature names.
Named numeric vector of variable importance scores from the random forest.
The spatial scales used.
The standard deviation threshold used.
Number of training observations.
A named list with training metrics: R2,
RMSE, and MAE.
If pred_data was provided, a data frame
with columns x, y, prediction, trend,
kriged_residual, and kriging_se.
Ren, K., Song, Y., Chen, M., and Yu, Q. (2026). A singularity regression kriging for spatial prediction. GIScience & Remote Sensing, 63(1), 2690341. doi:10.1080/15481603.2026.2690341
Breiman, L. (2001). Random forests. Machine Learning, 45(1), 5–32. doi:10.1023/A:1010933404324
Hengl, T., Heuvelink, G. B., and Rossiter, D. G. (2007). About regression-kriging: from equations to case studies. Computers & Geosciences, 33(10), 1301–1315. doi:10.1016/j.cageo.2007.05.001
# Generate simulated data sim_data <- sim_srk_data(n_side = 10, scenario = "normal", seed = 42) # Fit SRK model model <- srk(z ~ cov1, data = sim_data, coords = ~x + y, scales = c(1, 2, 3, 4, 5), ntree = 100) print(model) summary(model)# Generate simulated data sim_data <- sim_srk_data(n_side = 10, scenario = "normal", seed = 42) # Fit SRK model model <- srk(z ~ cov1, data = sim_data, coords = ~x + y, scales = c(1, 2, 3, 4, 5), ntree = 100) print(model) summary(model)
Performs spatial block cross-validation to evaluate the predictive accuracy of the SRK model. The study area is divided into spatial blocks, which are then assigned to folds for cross-validation, following the evaluation framework described in Ren et al. (2026, Section 3.2.3).
srk_block_cv( formula, data, coords, nfold = 5L, block_size, scales = NULL, ntree = 500L, sd_threshold = 0.5, min_neighbours = 3L, min_scales = 2L, variogram_model = "auto" )srk_block_cv( formula, data, coords, nfold = 5L, block_size, scales = NULL, ntree = 500L, sd_threshold = 0.5, min_neighbours = 3L, min_scales = 2L, variogram_model = "auto" )
formula |
A formula specifying the response and covariates. |
data |
A data frame containing the variables. |
coords |
A one-sided formula specifying coordinate columns. |
nfold |
Integer. Number of cross-validation folds. Default is 5. |
block_size |
Numeric. Size of spatial blocks (in coordinate units). |
scales |
Numeric vector of spatial scales for singularity computation.
See |
ntree |
Integer. Number of trees in the random forest. Default is 500. |
sd_threshold |
Numeric. Singularity feature standard deviation threshold. Default is 0.5. |
min_neighbours |
Integer. Minimum neighbours for singularity computation. Default is 3. |
min_scales |
Integer. Minimum scales for singularity computation. Default is 2. |
variogram_model |
Character. Variogram model type. Default
|
The study area is divided into non-overlapping rectangular spatial blocks
of the specified size. Each block is assigned to one of nfold folds.
For each fold, the SRK model is fitted on all observations outside the
fold, and predictions are generated for observations within the fold.
Performance metrics (R-squared, RMSE, MAE) are computed per fold and
overall.
An object of class "srk_cv" containing:
Numeric vector of observed values.
Numeric vector of cross-validated predictions.
Integer vector of fold assignments.
A data frame with per-fold R2, RMSE, and MAE.
A named list with overall R2, RMSE, and MAE.
Number of folds used.
Block size used.
Ren, K., Song, Y., Chen, M., and Yu, Q. (2026). A singularity regression kriging for spatial prediction. GIScience & Remote Sensing, 63(1), 2690341. doi:10.1080/15481603.2026.2690341
set.seed(42) sim_data <- sim_srk_data(n_side = 15, scenario = "normal", seed = 42) cv_result <- srk_block_cv(z ~ cov1, data = sim_data, coords = ~x + y, nfold = 3, block_size = 5, scales = c(1, 2, 3, 4, 5), ntree = 100) print(cv_result)set.seed(42) sim_data <- sim_srk_data(n_side = 15, scenario = "normal", seed = 42) cv_result <- srk_block_cv(z ~ cov1, data = sim_data, coords = ~x + y, nfold = 3, block_size = 5, scales = c(1, 2, 3, 4, 5), ntree = 100) print(cv_result)
Evaluates the robustness of the SRK model under different combinations of maximum singularity scale and singularity feature selection threshold, following the sensitivity analysis framework described in Ren et al. (2026, Section 4.4).
srk_sensitivity( formula, data, coords, scale_range, threshold_range, nfold = 5L, block_size, ntree = 500L, min_neighbours = 3L, min_scales = 2L, variogram_model = "auto", scale_step = 1 )srk_sensitivity( formula, data, coords, scale_range, threshold_range, nfold = 5L, block_size, ntree = 500L, min_neighbours = 3L, min_scales = 2L, variogram_model = "auto", scale_step = 1 )
formula |
A formula specifying the response and covariates. |
data |
A data frame containing the variables. |
coords |
A one-sided formula specifying coordinate columns. |
scale_range |
Numeric vector of maximum singularity scales to test.
For each value, a scale sequence from |
threshold_range |
Numeric vector of standard deviation thresholds to test for singularity feature selection. |
nfold |
Integer. Number of cross-validation folds. Default is 5. |
block_size |
Numeric. Size of spatial blocks for cross-validation. |
ntree |
Integer. Number of trees in the random forest. Default is 500. |
min_neighbours |
Integer. Minimum neighbours for singularity computation. Default is 3. |
min_scales |
Integer. Minimum scales for singularity computation. Default is 2. |
variogram_model |
Character. Variogram model type. Default
|
scale_step |
Numeric. Step size for generating scale sequences. Default is 1. |
For each combination of maximum singularity scale and feature selection
threshold, the function performs spatial block cross-validation using
srk_block_cv and records the resulting performance metrics.
This allows assessment of model robustness to parameter choices as
described in the paper.
An object of class "srk_sensitivity" containing:
A data frame with columns max_scale,
threshold, R2, RMSE, and MAE for each
parameter combination.
The scale range tested.
The threshold range tested.
Ren, K., Song, Y., Chen, M., and Yu, Q. (2026). A singularity regression kriging for spatial prediction. GIScience & Remote Sensing, 63(1), 2690341. doi:10.1080/15481603.2026.2690341
set.seed(42) sim_data <- sim_srk_data(n_side = 10, scenario = "normal", seed = 42) sens <- srk_sensitivity(z ~ cov1, data = sim_data, coords = ~x + y, scale_range = c(3, 5), threshold_range = c(0.3, 0.5), nfold = 3, block_size = 3, ntree = 50) print(sens)set.seed(42) sim_data <- sim_srk_data(n_side = 10, scenario = "normal", seed = 42) sens <- srk_sensitivity(z ~ cov1, data = sim_data, coords = ~x + y, scale_range = c(3, 5), threshold_range = c(0.3, 0.5), nfold = 3, block_size = 3, ntree = 50) print(sens)
A simulated spatial dataset on a 20 by 20 regular grid with a spatially autocorrelated response variable and covariate, generated following the simulation design described in Ren et al. (2026, Section 2.3.1).
srk_simsrk_sim
A data frame with 400 rows and 4 variables:
Horizontal grid coordinate (integer, 1 to 20).
Vertical grid coordinate (integer, 1 to 20).
Response variable (spatially autocorrelated, normal scenario).
Spatially autocorrelated covariate.
Generated using sim_srk_data with n_side = 20,
scenario = "normal", and seed = 42.
Ren, K., Song, Y., Chen, M., and Yu, Q. (2026). A singularity regression kriging for spatial prediction. GIScience & Remote Sensing, 63(1), 2690341. doi:10.1080/15481603.2026.2690341
data(srk_sim) head(srk_sim) summary(srk_sim)data(srk_sim) head(srk_sim) summary(srk_sim)
Provides a detailed summary of a fitted SRK model including variable importance, variogram parameters, and accuracy metrics.
## S3 method for class 'srk' summary(object, ...)## S3 method for class 'srk' summary(object, ...)
object |
An object of class |
... |
Additional arguments (currently unused). |
An object of class "summary.srk".
sim_data <- sim_srk_data(n_side = 10, scenario = "normal", seed = 42) model <- srk(z ~ cov1, data = sim_data, coords = ~x + y, scales = c(1, 2, 3, 4, 5), ntree = 100) summary(model)sim_data <- sim_srk_data(n_side = 10, scenario = "normal", seed = 42) model <- srk(z ~ cov1, data = sim_data, coords = ~x + y, scales = c(1, 2, 3, 4, 5), ntree = 100) summary(model)