GeoCV: Repeated holdout kriging cross-validation

In GeoModels: Procedures for Gaussian and Non Gaussian Geostatistical (Large) Data Analysis

Repeated holdout kriging cross-validation

Description

The procedure uses the GeoKrig or GeoKrigloc function to compute repeated holdout kriging cross-validation using information from a GeoFit object. The function returns prediction scores.

Usage

GeoCV(
  fit, K = 100, estimation = TRUE,
  optimizer = NULL, lower = NULL, upper = NULL,
  n.fold = 0.05, local = FALSE,
  neighb = NULL, maxdist = NULL, maxtime = NULL,
  sparse = FALSE, type_krig = "Simple", which = 1,
  parallel = FALSE, ncores = 6, progress = TRUE,
  seed = NULL, predictor = "linear", conditional_nrep = 1000L,
  conditional_n_iter = 25L, conditional_mcmc_thin = 1L
)

Arguments

fit	An object of class GeoFit.
K	Positive integer greater than or equal to 2 giving the number of cross-validation iterations. Non-integer values are rejected rather than rounded.
estimation	Logical; if TRUE, the model is re-estimated on the training observations at each iteration and the fold-specific estimates are used for prediction. If FALSE, the parameter estimates in fit are reused; this is therefore a conditional predictive assessment given parameters estimated from the complete dataset, not a fully refitted out-of-sample cross-validation.
optimizer	The type of optimization algorithm if estimation is TRUE. See GeoFit for details. If NULL, then the optimization algorithm stored in fit is used.
lower	An optional named list giving the values for the lower bounds of the parameters when bounded optimization is used and estimation is TRUE.
upper	An optional named list giving the values for the upper bounds of the parameters when bounded optimization is used and estimation is TRUE.
n.fold	Numeric; the fraction of observations randomly deleted and predicted in each cross-validation iteration. At least one training and one prediction observation are retained. In the space-time case, sampling also retains at least one training observation at every observed time.
local	Logical; if TRUE, local kriging is performed. The default is FALSE.
neighb	Numeric; an optional positive integer indicating the order of neighborhood if local kriging is performed.
maxdist	Numeric; an optional positive value indicating the spatial neighborhood distance if local kriging is performed.
maxtime	Numeric; an optional non-negative temporal-distance threshold, expressed in the same units as the fitted coordt, when local kriging is performed.
sparse	Logical; if TRUE, kriging and simulation are computed with sparse matrix algorithms using the spam package. The default is FALSE. It should be used with compactly supported covariance models.
type_krig	String; the type of kriging. If "Simple", the default, the fitted mean coefficients are treated as plug-in values. If "Universal", the prediction MSE is additionally corrected using the covariance matrix of the estimated mean coefficients. With estimation = TRUE, Universal cross-validation is supported for full-likelihood refits, for which varest = TRUE is requested automatically in each fold. It is deliberately rejected for composite-likelihood refits because each fold would require its own Godambe covariance matrix. With estimation = FALSE, a stored fit$varcov may be used; for composite likelihood it can be obtained with GeoVarest.
which	Numeric; in the case of bivariate cokriging, it indicates which variable to predict. It can be 1 or 2.
predictor	Character vector selecting the prediction method. The default "linear" preserves the historical behaviour based on GeoKrig or GeoKrigloc. For direct spatial Binomial and BinomialNeg random fields without a copula, "conditional" uses the conditional mean estimated by GeoSimcond. Supplying c("linear", "conditional") evaluates both predictors on exactly the same held-out observations and, when estimation=TRUE, from the same fold-specific refit. Conditional prediction is currently restricted to purely spatial, univariate, non-copula direct count models with local=FALSE and sparse=FALSE.
conditional_nrep	Positive integer; number of retained conditional simulations used by GeoSimcond to estimate the conditional mean within each cross-validation iteration. Used only when predictor includes "conditional". Default is 1000.
conditional_n_iter	Positive integer; burn-in length, in full Gibbs sweeps, for the direct count conditional sampler used within each cross-validation iteration. Used only when predictor includes "conditional". Default is 25.
conditional_mcmc_thin	Positive integer; thinning interval for retained direct-count conditional simulations within each cross-validation iteration. Used only when predictor includes "conditional". Default is 1.
parallel	Logical; default FALSE. Set TRUE to evaluate the cross-validation iterations in parallel.
ncores	Positive integer or NULL; default 6. With parallel=TRUE, an explicit integer requests that many workers, capped by detected cores and the number of cross-validation jobs. Set ncores=NULL for automatic selection, capped at six workers and normally leaving one detected core free. Non-integer values are rejected.
progress	Logical; if TRUE, a progress bar is shown.
seed	Optional finite integer seed used to make the random selection of folds and any stochastic refitting step reproducible. Non-integer values are rejected. If NULL, the current random number generator state is left unchanged. When supplied, the previous .Random.seed is restored on exit.

Details

The option GeoModels.cv_strict = TRUE rejects any excluded observation, non-finite prediction, or invalid Gaussian predictive MSE. The default is FALSE: exclusions are counted explicitly, and a fold with no usable predictions returns NA scores rather than an unreported change of denominator. Warnings raised by iterations are collected and summarized in the calling R session; options(warn = 2) continues to treat warnings as errors.

For a spatio-temporal GeoFit object, the stored temporal coordinates are reused without imposing equal spacing. Local temporal neighborhoods interpret maxtime as a distance threshold in the same units as those coordinates.

The function randomly removes a fraction n.fold of the observations at each iteration, predicts the removed observations, and computes predictive scores. With the default predictor="linear", prediction uses kriging as in previous versions. For supported direct binomial and negative-binomial random fields, predictor="conditional" instead estimates the conditional mean from repeated calls to the count-constrained conditional simulator. When both predictors are requested, the fold is constructed and the model is refitted only once, so the comparison is paired within each cross-validation iteration.

If estimation = TRUE, the model is re-estimated at each cross-validation iteration before prediction. Refits preserve the fitted likelihood settings, including pair weighting, pair thinning, distance-memory handling and anisotropy. For misspecified estimators, the response model stored in fit$model remains distinct from the working likelihood stored in fit$estimation_model; the latter is used for each refit. Estimated anisotropy parameters are passed once through anisopars, avoiding duplicate angle/ratio entries in the starting or fixed parameter lists.

For conditional prediction of direct Binomial and BinomialNeg random fields, each fold-specific fitted object is passed to GeoSimcond. The conditional predictor is the empirical mean of conditional_nrep retained conditional simulations after conditional_n_iter burn-in sweeps, with thinning conditional_mcmc_thin. Parallelisation, when requested, is performed over cross-validation iterations; the nested GeoSimcond call is deliberately run serially within each worker to avoid nested parallel execution.

For Universal kriging with fold-specific re-estimation, full-likelihood fits request varest = TRUE within each fold. Composite-likelihood refits are not assigned the Godambe covariance matrix from the full dataset: GeoCV instead rejects estimation = TRUE, type_krig = "Universal" for composite likelihood, because a statistically coherent analysis would require a fold-specific GeoVarest calculation.

If estimation = FALSE, the parameter estimates from the original fit are reused after removing the validation observations from the conditioning data. Consequently, the reported errors assess prediction conditional on parameters that were estimated using the complete dataset and may be more optimistic than a fully refitted out-of-sample cross-validation.

For Poisson, Binomial, and BinomialNeg with copula="SkewGaussian", global cross-validation (local=FALSE) supports both fold refitting and fixed-parameter prediction. Discrete Clayton-like copula models still require estimation=FALSE until their pairwise likelihood is implemented.

Regular-grid fits are converted to explicit point coordinates after observations are removed, because each training sample is no longer a complete Cartesian grid. Space-time covariates may be stored either as one matrix in observation order or as a list containing one matrix per time.

For Gaussian fitted models all documented scores are computed from the Gaussian predictive mean and MSE. For non-Gaussian fitted models, only RMSE, MAE and MAD are returned; brie, crps, lscore, pit, intscore and coverage are set to NA, because an exact non-Gaussian predictive distribution is not available from only a mean and MSE.

When seed is supplied, the cross-validation samples are reproducible. The function preserves and restores the user's random number generator state.

Value

Returns a list containing predictive score vectors. With the default predictor="linear", the historical top-level structure is unchanged. When conditional prediction is requested, components linear and/or conditional contain the corresponding score lists. For the conditional predictor, RMSE, MAE and MAD are populated and the Gaussian-only probability scores remain NA. For backward compatibility, the top-level score aliases refer to the linear predictor when it is requested, and otherwise to the conditional predictor.

linear	When requested, a list of score vectors for the linear predictor.
conditional	When requested, a list of score vectors for the conditional-mean predictor.
rmse	The vector of root mean squared errors for the primary predictor.
mae	The vector of mean absolute errors for the primary predictor.
mad	The vector of median absolute errors for the primary predictor.
brie	The vector of Brier scores, or NA for non-Gaussian fits.
crps	The vector of continuous ranked probability scores, or NA for non-Gaussian fits.
lscore	The vector of log-scores, or NA for non-Gaussian fits.
pit_mean	The vector of mean probability integral transform values within each holdout fold, or NA for non-Gaussian fits. This is a descriptive PIT summary, not a calibration score: a mean near 0.5 does not by itself imply a uniform PIT distribution.
pit	Backward-compatible alias of pit_mean.
intscore	The vector of interval scores, or NA for non-Gaussian fits.
coverage	The vector of empirical coverage values, or NA for non-Gaussian fits.
n_requested	Number of requested predictions in each iteration.
n_valid	Number of finite observation/prediction pairs used for error scores.
n_failed	Number of excluded pairs; n_requested - n_valid.
n_invalid_mse	Among valid pairs, the number with non-finite or nonpositive predictive MSE. NA when Gaussian probability scores are not applicable.
n_valid_prob	Number of pairs used for Gaussian probability scores; NA when not applicable.
warnings	A frequency table of warning messages from all iterations, or NULL when none were raised.
seed	The seed used for reproducibility, or NULL.
predictor	The prediction method or methods requested.
conditional_nrep	When conditional prediction is requested, the number of retained conditional simulations per fold.
conditional_n_iter	When conditional prediction is requested, the Gibbs burn-in length per fold.
conditional_mcmc_thin	When conditional prediction is requested, the Gibbs thinning interval.

Author(s)

Moreno Bevilacqua, [e-mail], https://sites.google.com/view/moreno-bevilacqua/home, Víctor Morales Oñate, [e-mail], https://sites.google.com/site/moralesonatevictor/, Christian Caamaño-Carrillo, [e-mail], https://www.researchgate.net/profile/Christian-Caamano

See Also

GeoKrig, GeoKrigloc, GeoSimcond, GeoFit.

Examples

library(GeoModels)

###############################################################
### Example of spatial kriging cross-validation
###############################################################

model <- "Gaussian"
set.seed(79)
x <- runif(400, 0, 1)
y <- runif(400, 0, 1)
coords <- cbind(x, y)

corrmodel <- "GenWend"
mean <- 0
sill <- 5
nugget <- 0
scale <- 0.2
smooth <- 0
power2 <- 4

param <- list(
  mean = mean, sill = sill, nugget = nugget,
  scale = scale, smooth = smooth, power2 = power2
)

data <- GeoSim(coordx = coords, corrmodel = corrmodel,
               param = param)$data

fixed <- list(nugget = nugget, smooth = 0, power2 = power2)
start <- list(mean = 0, scale = scale, sill = 1)
I <- Inf
lower <- list(mean = -I, scale = 0, sill = 0)
upper <- list(mean = I, scale = I, sill = I)

fit <- GeoFit(
  data, coordx = coords, corrmodel = corrmodel,
  model = model, likelihood = "Marginal", type = "Pairwise",
  neighb = 3, optimizer = "nlminb", lower = lower,
  upper = upper, start = start, fixed = fixed
)

#a <- GeoCV(fit, K = 100, estimation = TRUE,
#           parallel = TRUE, seed = 123)
#mean(a$rmse)

GeoModels documentation built on Sept. 23, 2026, 5:07 p.m.