GeoCV: 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
)

Value

Returns a list containing the following components:

rmse: The vector of root mean squared errors.
mae: The vector of mean absolute errors.
mad: The vector of median absolute errors.
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.
seed: The seed used for reproducibility, or NULL.

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.
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.

Author

Moreno Bevilacqua, moreno.bevilacqua89@gmail.com, https://sites.google.com/view/moreno-bevilacqua/home, Víctor Morales Oñate, victor.morales@uv.cl, https://sites.google.com/site/moralesonatevictor/, Christian Caamaño-Carrillo, chcaaman@ubiobio.cl, https://www.researchgate.net/profile/Christian-Caamano

Details

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 using kriging, and computes a set of predictive scores.

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 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.

Examples

Run this code

# \donttest{
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)
# }

Run the code above in your browser using DataLab