Package 'spatialkit'

Description

prep_model_data() and assign_features_to_polygons() return a layer with an attribute recording what happened to its rows ("dropped" and "ties" respectively). Those records describe the rows the layer was built with, and [ on an sf object copies attributes through unchanged, which would leave a subset reporting its parent's numbers for a different set of rows. Subsetting therefore returns a plain layer with the record removed, and so do the dplyr verbs that select or reorder rows (filter(), slice(), arrange(), distinct()); read the record from the layer the function returned, before subsetting it. Binding such layers (rbind(), dplyr::bind_rows()) likewise returns a plain sf layer.

Usage

## S3 method for class 'spatialkit_rows'
x[...]

Arguments

Details

Each record carries n_rows, the number of rows it was made for. sf::st_drop_geometry() keeps the rows, and with them the record: it returns a data frame of class c("spatialkit_rows", "data.frame"). Binding such data frames with rbind() or dplyr::bind_rows() keeps the first one's record and class, so the record then describes only the first input's rows: its n_rows no longer equals nrow() of the result. The package's own readers ignore a record whose n_rows does not match; when reading attr(x, "dropped") or attr(x, "ties") yourself from a layer that has been through such steps, check it the same way.

Value

The subset, without the row records and without this class. A subset that is not a data frame (a single column taken with drop = TRUE) is returned unchanged.

See Also

prep_model_data() and assign_features_to_polygons(), the two functions that attach the records this method removes.

Examples

library(sf)
dat <- st_as_sf(
  data.frame(x = 1:5, y = 5:1,
             resp = c(1, 2, NA, 4, 5), pred = c(1, 2, 3, 4, Inf)),
  coords = c("x", "y"), crs = 32632
)
clean <- prep_model_data(dat, "resp", "pred")
attr(clean, "dropped")$n          # 2
attr(clean[1:2, ], "dropped")     # NULL -- the record does not follow

Description

Computes the dissimilarity index (DI) of Meyer & Pebesma (2021) for a set of prediction locations and flags those that fall inside the model's area of applicability (AOA), the region of predictor space where the model's cross-validated performance estimate can be expected to hold.

Usage

area_of_applicability(
  newdata,
  model = NULL,
  train_sf = NULL,
  predictor_vars = NULL,
  weights = NULL,
  folds = NULL,
  threshold = NULL,
  normalizer_max_n = 5000L,
  seed = 123L,
  chunk_size = NULL,
  use_fnn = requireNamespace("FNN", quietly = TRUE)
)

Arguments

Value

An object of class aoa: a list with

• aoa: newdata with a numeric DI column and a logical AOA column added. This is the object the computation ran on, which for a coordinate-using model is newdata after pointizing, CRS reconciliation and the addition of the "..x" and "..y" columns. A row whose predictors are not all finite gets NA in both columns, and a row outside the training range of a predictor in dropped_vars gets DI = Inf and AOA = FALSE (see Limitations).

• threshold: the DI cut-off used.

• train_DI: the training points' own DI values.

• normalizer: the mean pairwise training distance.

• weights: the weight vector actually applied, named by predictor_vars.

• predictor_vars: the predictors used, including "..x"/"..y" when the model uses coordinates and excluding dropped_vars.

• dropped_vars: predictors dropped for negligible variance.

• scaling: a list with center and scale, each named by predictor_vars: the training means and standard deviations the index is computed in, so a location's DI can be traced to the predictor that put it outside.

• n_outliers: the number of training DI values above the Q3 + 1.5 * IQR fence, which the default threshold rule sets aside (the "outlier-removed" in its name); computed whether or not threshold was supplied.

• n_train, n_new, n_inside, n_outside, n_na: row counts. n_train counts the training rows that survived the finite-value filter; n_new is every row of newdata, so n_new = n_inside + n_outside + n_na.

• params records the call: folds_supplied, n_folds, folds_method, threshold_supplied, normalizer_max_n, normalizer_n_used, normalizer_subsampled, weights_supplied and seed. folds_method is the method of a make_folds() result ("block_kfold", "random_kfold", ...), "labels" for a vector of fold labels, "splits" for a bare list of train/test splits, and NA when no folds were supplied. It is printed with the object, because the threshold's meaning depends on it: the cross-validated DI that sets it comes from the same kind of hold-out as the CV error it should be quoted beside, so an AOA built on random_kfold folds is logged as a caution and does not belong next to a blocked cv_*() result.

Why a map alone is not enough

A fitted model will return a number for any location you hand it, including locations whose predictor values look nothing like anything it was trained on. Those predictions are extrapolations dressed as interpolations, and a cross-validation score says nothing about them, because the held-out folds were drawn from the same predictor distribution as the training data. The AOA marks where the score applies.

How it is computed

Predictors are centred and scaled using the training data's own means and standard deviations, then optionally weighted by variable importance. For a prediction point $p$, the DI is the distance to its nearest training point in that space, divided by the mean pairwise distance among training points. The same quantity is computed for the training data itself, using each point's nearest neighbour among the training rows of the fold that holds it out. That means everything outside its own fold for random and block folds, and the smaller training set that buffered and NNDM folds actually leave (see the next section). The threshold is then the largest training DI that is not an upper outlier, i.e. not above the fence Q3 + 1.5 * IQR of the training DI, with the quartiles of stats::quantile()'s default type 7. Prediction points at or below that threshold are inside the AOA.

That is the paper's "outlier-removed maximum". CAST, the reference implementation, computes the same fence with the same quartiles but uses the fence itself as the threshold, capped at the largest training DI. The two agree whenever no training DI lies above the fence (n_outliers is 0); otherwise CAST's threshold is the larger, and so is its AOA. (Earlier CAST releases used grDevices::boxplot.stats(), which gives the rule used here but with Tukey's hinges as the quartiles, so they can also differ when the number of training points is even.) To apply the current CAST rule to the same training DI, pass threshold = min(quantile(res$train_DI, 0.75) + 1.5 * IQR(res$train_DI), max(res$train_DI)) for an earlier result res.

The DI is invariant to the overall scale of weights: the numerator and the normaliser carry the same factor. Importance values can be passed as-is.

Each training point's reference is its nearest other training row, so an exact duplicate in predictor space (repeat visits to a site with static covariates, or covariates read off a raster coarser than the sampling) has a training DI of 0. Once about three quarters of the rows have a twin among their reference rows the threshold is 0, and only exact copies of a training row count as inside. That is logged as a caution; folds that keep the duplicates together (make_folds(method = "leave_location_out", group_var = ...)), or removing them, give the threshold its meaning back.

The fold scheme changes the answer, and should

With folds = NULL the training reference is each point's nearest neighbour anywhere in the training data, which for clustered data is very close, giving a small threshold and a conservative AOA. Passing the folds you actually validated with makes the reference distances larger and the AOA correspondingly wider. That is not a loophole. The AOA is defined relative to a performance estimate, and a spatially blocked estimate is a claim about predicting further away. Pass the same make_folds() result you passed to cv_spatial. Buffered and NNDM folds use the training set they actually left available, not merely "everything outside the fold".

Limitations

Predictors must be numeric; categorical variables are refused and never silently dummy-coded. Predictors whose variance is negligible relative to their own magnitude (the test is sd < sqrt(.Machine$double.eps) * max(abs(x)), so the same variable in metres and in gigametres is treated identically) are dropped from the distance. A prediction point taking a value there outside the training range (with the same relative tolerance) is extrapolation along a direction the training data never varied in: its scaled distance along it is infinite, so it gets DI = Inf, is outside the AOA, and a warning gives the count. A point missing that value is judged on the other predictors. Without weights every predictor counts equally, which overstates dissimilarity along directions the model barely uses.

Models fitted with the coordinates as predictors

When model was fitted with include_coords = TRUE the model splits on location, so the dissimilarity index has to measure location too: the coordinates are added to both sides as the predictors "..x" and "..y" and are then centred, scaled and weighted like any other column. Without this a prediction point far outside the training extent but with ordinary covariate values reads as inside the area of applicability. That is exactly the extrapolation this index exists to catch.

This path needs geometry on both sides, so train_sf and newdata must both be sf objects; a data.frame is refused rather than quietly measured without location. Non-POINT newdata (grid polygons, say) is reduced to representative points first, as coerce_to_points would. If exactly one side carries a CRS the other is brought into it: reprojected when its coordinates look like longitude/latitude, stamped otherwise, with a warning either way. This is done because degrees fed into a metre-space index silently understate the distances.

References

Meyer, H. and Pebesma, E. (2021). Predicting into unknown space? Estimating the area of applicability of spatial prediction models. Methods in Ecology and Evolution 12(9), 1620–1633. doi:10.1111/2041-210X.13650

See Also

predict_surface to build the grid, make_folds for the fold scheme.

Other cross-validation: cv_bayes(), cv_block_size_sweep(), cv_gwr(), cv_rf(), cv_spatial(), estimate_sac_range(), fold_separation(), gwr_model_selection(), make_folds(), sac_nugget(), select_features_forward()

Other prediction: predict_surface()

Examples

library(sf)
set.seed(1)
n <- 120
train <- st_as_sf(
  data.frame(x = 5e5 + runif(n, 0, 1000), y = 5e6 + runif(n, 0, 1000),
             a = rnorm(n), b = rnorm(n)),
  coords = c("x", "y"), crs = 32632
)
train$z <- 2 * train$a - train$b + rnorm(n, 0, 0.3)

# Prediction points, some of them well outside the training predictor range
newpts <- st_as_sf(
  data.frame(x = 5e5 + runif(50, 0, 1000), y = 5e6 + runif(50, 0, 1000),
             a = c(rnorm(40), rnorm(10, 8)), b = rnorm(50)),
  coords = c("x", "y"), crs = 32632
)

res <- area_of_applicability(newpts, train_sf = train,
                             predictor_vars = c("a", "b"))
res
table(res$aoa$AOA)

Description

Joins an sf layer of input features to a polygon layer via spatial join.

Usage

assign_features_to_polygons(
  features_sf,
  polygons_sf,
  polygon_id_col = "poly_id",
  keep_unassigned = FALSE,
  predicate = sf::st_intersects,
  largest = TRUE,
  tie_break = c("smallest_area", "first")
)

Arguments

Details

This is the second step of the package's pipeline: it labels every observation with the cell it falls in, which is what summarize_by_cell() then aggregates over. Prefer it to sf::st_join() when the join has to be unambiguous. It resolves features matching several polygons by an explicit tie_break rule instead of silently duplicating rows, so the assigned layer keeps one row per input feature and cell-level counts mean what they say.

The join runs in the CRS of polygons_sf whenever that CRS is projected, so cell edges are the straight lines the cells were drawn with and overlap areas are planar. A copy of features_sf is transformed for it, and the features come back with the coordinates they arrived with. Otherwise (the polygons are in lon/lat, or carry no CRS) the join runs in the CRS of features_sf.

Value

An sf object with polygon_id_col attached, one row per input feature (fewer if keep_unassigned = FALSE dropped unmatched ones), in the CRS features_sf arrived in. A column of features_sf already called polygon_id_col is dropped before the spatial join (with a warning), so re-assigning an already-assigned layer replaces the old IDs and does not fail. Every other column is kept, including one named like the polygons' own ID column when that is read from a fallback such as "id" (a site id joined to cells keyed by id). If no feature falls inside any polygon the result is empty (or all-NA with keep_unassigned = TRUE) and a warning is raised, since the usual cause is two layers in different places (a CRS that could only be stamped, not reprojected). The attribute "ties" records how many features matched more than one polygon (with largest, overlapped several by exactly the same area) and had the tie_break rule decide for them: a list with n, which (their row positions in features_sf), rule and n_rows (the number of rows returned, which the record was made for). A large n means the polygon layer overlaps, and per-cell counts built from the result depend on the rule. The record describes the rows this call returned and does not survive subsetting: joined[i, ], like dplyr::filter(), slice() or arrange() of it, is a plain layer with no "ties" attribute, so nothing reports the parent's count for a different set of rows. sf::st_drop_geometry() keeps the record, since the rows are the same; see [.spatialkit_rows for what binding such data frames does.

See Also

build_tessellation() to build the polygon layer; summarize_by_cell() for the aggregation step that consumes the result.

Other aggregation: determine_optimal_levels(), kriging_adequacy(), resolution_profile(), select_resolution(), summarize_by_cell(), summary.resolution_profile()

Examples

library(sf)
set.seed(1)
pts <- st_as_sf(
  data.frame(x = 5e5 + runif(20, 0, 100), y = 5e6 + runif(20, 0, 100),
             val = rnorm(20)),
  coords = c("x", "y"), crs = 32632
)
bnd <- st_sf(geometry = st_sfc(st_polygon(list(rbind(
  c(5e5, 5e6), c(5e5 + 100, 5e6), c(5e5 + 100, 5e6 + 100),
  c(5e5, 5e6 + 100), c(5e5, 5e6)
))), crs = 32632))
grid <- create_grid_polygons(bnd, target_cells = 9, type = "square")
assigned <- assign_features_to_polygons(pts, grid)
table(assigned$poly_id)

Description

The single entry point for turning a point pattern into analysis regions, and the first step of the package's pipeline. It wraps the four tessellation methods behind one interface that handles CRS projection, clipping and stable cell identifiers consistently, and returns the cell layer together with the point-to-cell index that assign_features_to_polygons() and summarize_by_cell() consume. Prefer it to the individual constructors whenever you might want to compare methods: the return shape does not change with method, so swapping "voronoi" for "hex" costs one argument.

Usage

build_tessellation(
  points_sf,
  boundary = NULL,
  method = c("voronoi", "triangles", "hex", "square"),
  approx_n_cells = NULL,
  cellsize = NULL,
  expand = 0,
  clip = TRUE,
  keep_duplicates = FALSE,
  crs = NULL,
  quiet = FALSE
)

Arguments

Details

Which method to reach for. "voronoi" gives one cell per point, so resolution follows sampling density. That is the choice when the observations themselves define the regions. "hex" and "square" give equal-area cells on a fixed grid, so cell size is a decision you make, and the data does not make it for you; hexagons avoid the axis-aligned artefacts of squares and have uniform neighbour distances. "triangles" returns the Delaunay triangulation, useful for interpolation and adjacency work; it is not meant as an aggregation unit. determine_optimal_levels() will suggest a cell count from the spatial structure of the data.

"voronoi" and "triangles" are built on the points' vertices: a MULTIPOINT feature with several vertices gets one cell (or triangle corner) per vertex, and its index entry is the smallest cell_id among the cells it touches. See create_voronoi_polygons().

Value

A list with components:

cells

An sf polygon layer, one row per cell. It always carries a cell_id column; the "hex" and "square" methods additionally carry poly_id, which holds the same values.

index

Integer vector of cell_id values, one per row of points_sf, and NA for a point that falls inside no cell, that is, one outside the study area. Only a point within a thousandth of the median cell width of a cell is snapped to it. That covers points sitting exactly on a shared edge, and leaves points outside the study area as NA. A summary built from index therefore counts only the points the tessellation actually covers.

boundary

The boundary used (possibly derived and/or reprojected, and for "voronoi" grown by expand).

method

The method actually used.

params

The parameters the tessellation was built with, plus snapped, the record of that nearest-cell repair: a list with n, which (row positions in points_sf) and distance (how far outside every cell each sat, in CRS units). For "hex" and "square" also cells_occupied and cells_empty, the number of cells that hold at least one point and that hold none.

See Also

Other tessellation: create_grid_polygons(), create_grid_polygons_cached(), create_voronoi_polygons(), ensure_stable_poly_id(), get_voronoi_seeds(), plot_tessellation_map(), voronoi_seeds_kmeans(), voronoi_seeds_random()

Examples

library(sf)
set.seed(1)
pts <- st_as_sf(
  data.frame(x = 5e5 + runif(20, 0, 100), y = 5e6 + runif(20, 0, 100)),
  coords = c("x", "y"), crs = 32632
)
tess <- build_tessellation(pts, method = "voronoi", quiet = TRUE)
tess$cells

Description

Removes the lazily-cached fitted() result so that the next call recomputes from the posterior. This is only necessary if the underlying brmsfit engine has been manually mutated after fitting. A change to data_sf invalidates the entry on its own, because the cached value carries a digest of the data it was computed from (see fitted.bayesian_fit). Normal usage never requires it.

Usage

clear_fitted_cache(object)

Arguments

Details

The cache environment is shared by every copy of a fit, so clearing it through one copy clears it for all of them. That is harmless: the others recompute.

Value

object, invisibly (called for side effect).

See Also

Other package options and caches: clear_grid_cache(), create_grid_polygons_cached(), spatialkit_quiet()

Examples

# Only a bayesian_fit carries the cache; on any other fit this is a no-op.
if (requireNamespace("ranger", quietly = TRUE)) {
  library(sf)
  set.seed(1)
  pts <- st_as_sf(
    data.frame(x = 5e5 + runif(60, 0, 1000), y = 5e6 + runif(60, 0, 1000),
               a = rnorm(60)),
    coords = c("x", "y"), crs = 32632
  )
  pts$z <- 2 * pts$a + rnorm(60, 0, 0.3)
  fit <- fit_rf_model(pts, "z", "a", num_trees = 50, seed = 1)
  clear_fitted_cache(fit)
}

Description

Removes all memoized grid results from the internal cache environment.

Usage

clear_grid_cache(cache_env = .gmt_cache)

Arguments

Value

Invisibly, the number of entries removed.

See Also

Other package options and caches: clear_fitted_cache(), create_grid_polygons_cached(), spatialkit_quiet()

Examples

library(sf)
bnd <- st_sf(geometry = st_sfc(st_polygon(list(rbind(
  c(0, 0), c(100, 0), c(100, 100), c(0, 100), c(0, 0)
))), crs = 32632))
g1 <- create_grid_polygons_cached(bnd, target_cells = 9)
g2 <- create_grid_polygons_cached(bnd, target_cells = 9)   # cache hit
clear_grid_cache()   # returns the number of entries removed, invisibly
print(clear_grid_cache())                                   # 0: already empty

Description

Resolves a single polygon to tessellate within. With a boundary it is that boundary (optionally buffered by expand); without one it is the axis-aligned bounding box of points_sf (the rectangle in the working CRS), again optionally buffered, or a small buffer around the points when they all share one x or one y, or nearly so (the short side of their bounding box below a millionth of the long side). Reach for it to build the boundary that method = "hex" and "square" require, to check that a study-area polygon actually contains the observations before tessellating, or to pass the same envelope to create_voronoi_polygons() and create_grid_polygons() so that two tessellations of one dataset cover identical ground.

Usage

clip_target_for(points_sf, boundary = NULL, expand = 0, quiet = FALSE)

Arguments

Details

It is not the target build_tessellation() derives on its own: method = "voronoi" without a boundary clips to the convex hull of the points buffered by 2 percent of its diagonal, and there expand is always a distance. A bounding box over a non-rectangular point cloud includes corners with no data, so a hex or square grid laid over it has cells that hold no points; pass the study-area polygon when there is one.

Value

An sf polygon layer representing the clip target. For lon/lat input the layer is returned in the automatically selected local projected CRS, not the input CRS; a message reports this unless quiet = TRUE.

See Also

Other spatial data preparation: coerce_to_points(), ensure_projected(), harmonize_crs(), prep_model_data()

Examples

library(sf)
set.seed(1)
pts <- st_as_sf(
  data.frame(x = 5e5 + runif(30, 0, 100), y = 5e6 + runif(30, 0, 100)),
  coords = c("x", "y"), crs = 32632
)
# No boundary: the bounding box, expanded by 10% of the extent
box <- clip_target_for(pts, expand = 0.1, quiet = TRUE)
st_area(box)

Description

Returns the posterior summary of the global (non-spatial) regression terms: estimate, error and credible interval per predictor, as brms::fixef() reports them. Reach for it to read the average effect of a predictor with its uncertainty attached (the Bayesian counterpart to a coefficient table), remembering that the Gaussian-process term has already absorbed the spatially structured part of the signal, so these are effects net of location.

Usage

## S3 method for class 'bayesian_fit'
coef(object, ...)

Arguments

Value

A matrix of fixed-effect posterior summaries, as returned by brms::fixef(), on the fitted scale (per standard deviation of each predictor under standardize_predictors = TRUE; see above). Never NULL: a missing 'brms' or a failing fixef() call errors, following the coef() contract described in new_spatial_fit.

Standardised predictors

The summaries are on the scale the model was fitted on. A fit made with standardize_predictors = TRUE was fitted on centred and scaled numeric predictors, so each slope is the change in the linear predictor per standard deviation of its predictor and the intercept is its value at the predictor means, not the raw-unit numbers stats::lm() reports on the same formula. Nothing on the returned matrix says so; print() on the fit does, and the centre and scale of each predictor are in object$info$predictor_scaling. To put a slope back in raw units divide its Estimate, Est.Error and interval bounds by that predictor's scale. The intercept's Estimate follows by linearity (subtract each raw-unit slope times its predictor's center), but its Est.Error and interval depend on the posterior covariance of the coefficients: transform the draws from brms::as_draws_df(object$engine) for those, or refit without standardising.

See Also

Other methods on a fitted model: coef.gwr_fit(), coef.rf_fit(), fitted.bayesian_fit(), fitted.gwr_fit(), fitted.rf_fit(), model_metrics(), predict.bayesian_fit(), predict.gwr_fit(), predict.rf_fit(), print.rf_fit(), print.spatial_fit(), residuals.bayesian_fit(), residuals.gwr_fit(), residuals.rf_fit(), summary.spatial_fit()

Description

Returns the whole surface of coefficients, one row per observation and one column per term, against the single global vector coef() returns for an lm. That table is the point of fitting a GWR at all: inspect the spread of a predictor's column to see where, and by how much, its relationship with the response changes across the study area, and join it back to object$data_sf to map it. Use plot.spatial_fit() for a quick look at that map.

Usage

## S3 method for class 'gwr_fit'
coef(object, ...)

Arguments

Value

A data.frame of local coefficient estimates: one row per observation, one column per model term. Never NULL: when the engine carries no SDF component this errors, following the coef() contract described in new_spatial_fit.

What is and is not returned

Only the model terms: the intercept and one column per predictor. GWmodel's SDF data slot carries a good deal more alongside them (standard errors, t-values, the observed response, the fitted values, the residuals, Local_R2): 15 columns for a two-predictor fit, of which 3 are coefficients. Returning the whole slot would have made coef(fit)$Local_R2 and coef(fit)$a_SE read like coefficients and ncol(coef(fit)) a meaningless number. Reach for object$engine$SDF when you want the rest; it is the unmodified GWmodel object.

If the model terms cannot be located in the SDF (a GWmodel that names its coefficient columns differently), the whole slot is returned with a warning saying so. The call neither errors nor returns a silently short table.

See Also

Other methods on a fitted model: coef.bayesian_fit(), coef.rf_fit(), fitted.bayesian_fit(), fitted.gwr_fit(), fitted.rf_fit(), model_metrics(), predict.bayesian_fit(), predict.gwr_fit(), predict.rf_fit(), print.rf_fit(), print.spatial_fit(), residuals.bayesian_fit(), residuals.gwr_fit(), residuals.rf_fit(), summary.spatial_fit()

Description

Consistent with coef.gwr_fit() and coef.bayesian_fit(), which also error, instead of returning NULL, when they cannot supply coefficients. See the coef() contract in new_spatial_fit.

Usage

## S3 method for class 'rf_fit'
coef(object, ...)

Arguments

Value

Never returns; always signals an error.

See Also

Other methods on a fitted model: coef.bayesian_fit(), coef.gwr_fit(), fitted.bayesian_fit(), fitted.gwr_fit(), fitted.rf_fit(), model_metrics(), predict.bayesian_fit(), predict.gwr_fit(), predict.rf_fit(), print.rf_fit(), print.spatial_fit(), residuals.bayesian_fit(), residuals.gwr_fit(), residuals.rf_fit(), summary.spatial_fit()

Description

Converts the geometry column of an sf object to POINTs using one of several strategies.

Usage

coerce_to_points(
  x,
  mode = c("auto", "centroid", "point_on_surface", "surface", "line_midpoint",
    "bbox_center"),
  tmp_project = TRUE
)

Arguments

Details

The result has one row per row of x, in the same order. An EMPTY geometry of any type, lines included, becomes an EMPTY POINT in its own row; prep_model_data() and make_folds() then drop such rows, as they drop any other empty geometry. Empty lines are never handed to sf::st_line_sample(): it yields no midpoint for them, which would misalign the result, and with sf 1.0.x an empty MULTILINESTRING (or an empty part of one) crashed the R session. An empty part inside a non-empty feature is ignored, so the feature gets the point its other parts give; GEOS's interior point, used by "point_on_surface" and by the temporary projection's choice of CRS, segfaulted on an empty line part too.

Value

An sf object with geometry coerced to POINTs, row for row with x; an empty input geometry gives an empty POINT.

See Also

Other spatial data preparation: clip_target_for(), ensure_projected(), harmonize_crs(), prep_model_data()

Examples

library(sf)
poly <- st_sf(
  id = 1,
  geometry = st_sfc(st_polygon(list(rbind(
    c(0, 0), c(2, 0), c(2, 2), c(0, 2), c(0, 0)
  ))), crs = 32632)
)
coerce_to_points(poly, "auto")  # interior representative point

Description

Takes a named list of already-fit spatial_fit objects and produces a tidy comparison table including in-sample (for a forest, out-of-bag) metrics and model-specific information criteria (AICc, LOOIC).

Usage

compare_models(fits, newdata = NULL, ...)

Arguments

Value

A data.frame comparing all models. Its metric_basis column says what each row's metrics were computed on (see evaluate_insample); a table that mixes "out-of-bag" and "in-sample" rows does not rank the models, and says so in the log. AICc (GWR) and LOOIC (Bayesian) are sums over the rows a model was fitted to, so each column is set to NA, with a warning, when the fits carrying it were fitted to different rows (a predictor with missing values drops rows, for example) or to different responses (a transformed response on the same rows), and the warning says which. convergence_ok is TRUE or FALSE for a Bayesian fit whose convergence was checked (see fit_bayesian_spatial_model) and NA otherwise; a fit that did not converge is ranked like the others, so it also raises a warning. Alongside the metrics it carries resid_morans_I, resid_morans_z, resid_morans_p and resid_morans_null, the last of which names the null residual_morans_i scored each model against, since that choice is per-fit and governs how much the p-value is worth. A significant p-value is noted in the log (not raised as an R warning): positive autocorrelation as structure the model may have missed, negative (resid_morans_z < 0) as the alternating residuals of a model that tracks its data closely. Read the caveats in ?residual_morans_i before treating silence as evidence of no residual structure.

What the metrics are computed on

With newdata = NULL the metrics come from fitted(object). That is in-sample for a gwr_fit or a bayesian_fit, but out-of-bag for an rf_fit, whose fitted() method returns out-of-bag predictions (see fit_rf_model). The data.frame model_metrics() returns carries no label distinguishing the two, so check object$info$fitted_are_oob before comparing numbers across backends; evaluate_insample() and compare_models() record it per model in a metric_basis column. compare_models_cv scores every backend the same way.

$R^2$ is $1 - RSS/TSS$ with the total sum of squares taken about the mean of the response the model was fitted to. In sample that is the ordinary $R^2$. With newdata it is out-of-sample $R^2$, the convention every cv_*() function uses: the model is measured against the prediction it had to beat, the training mean, not against the new rows' own mean, which it could not have known. It is below 0 when the model predicts the new rows worse than the training mean does, and it is NA when the response does not vary about that baseline by more than rounding error (100 machine epsilons of its magnitude, whatever its units).

Percentage errors on responses with zeros

MAPE divides by the observed value and SMAPE by $|y| + |\hat{y}|$, so neither is defined where its denominator is zero. Neither returns Inf or NaN. Both are averaged over the rows whose denominator is non-zero, and are NA when no row qualifies. Non-zero is judged at the scale of the data: a denominator no larger than 100 machine epsilons times the largest one counts as zero, so the rule does not depend on the units of the response. The n_MAPE and n_SMAPE columns record how many rows that was; the n column counts finite observation/prediction pairs. Read a percentage error next to its count: when n_MAPE < n, MAPE is an average over a subset of the data, whatever its value.

This bites on any response taking exact zeros: counts, rainfall, abundance, claim amounts. On a zero-inflated response with 62 zeros out of 120, MAPE is an average over the 58 non-zero rows, which n_MAPE = 58 now says. SMAPE fails differently and more subtly: it drops the rows where observation and prediction are both near zero (which on a well-fitted zero-inflated model are the rows it got right), so it averages the harder rows only and reads worse than the fit deserves; n_SMAPE shows how many rows it kept, and the count is only a label, not a repair.

RMSE, MAE and $R^2$ use every finite row and are unaffected; prefer them whenever the response can be zero. For a Bayesian fit, cv_bayes() additionally reports CRPS and interval coverage, which are proper scoring rules and have no such failure mode.

See Also

Other model evaluation: compare_models_cv(), evaluate_insample(), model_metrics(), residual_morans_i()

Examples

if (requireNamespace("ranger", quietly = TRUE)) {
  library(sf)
  set.seed(1)
  pts <- st_as_sf(
    data.frame(x = 5e5 + runif(60, 0, 1000), y = 5e6 + runif(60, 0, 1000),
               a = rnorm(60), b = rnorm(60)),
    coords = c("x", "y"), crs = 32632
  )
  # An offset keeps the response away from zero, where MAPE is undefined.
  pts$z <- 10 + 2 * pts$a - pts$b + rnorm(60, 0, 0.3)
  # Two predictor sets for the same learner: the model that leaves `b` out
  # against the one that has it.  These are in-sample (out-of-bag) numbers;
  # compare_models_cv() scores the same question on spatial folds.
  fits <- list(a_only = fit_rf_model(pts, "z", "a", num_trees = 100, seed = 1),
               a_and_b = fit_rf_model(pts, "z", c("a", "b"), num_trees = 100, seed = 1))
  compare_models(fits)
}

Description

Fits and cross-validates one or more model types, returning a unified comparison table. Unlike compare_models(), this function does perform fitting (inside CV folds), because CV inherently requires repeated fitting.

Usage

compare_models_cv(
  data_sf,
  response_var,
  predictor_vars,
  models = c("GWR", "Bayesian"),
  k = 5,
  seed = 123,
  folds = NULL,
  boundary = NULL,
  pointize = "auto",
  gwr_args = list(),
  bayes_args = list(),
  rf_args = list(),
  summary = c("mean", "median"),
  quiet = FALSE,
  block_size = NULL,
  auto_range = FALSE,
  metrics = NULL
)

Arguments

Value

A list with overall, by_fold, and per-model cv_results (gwr_cv, bayes_cv, rf_cv for the models that ran). overall has one row per model with the pooled metrics, n_pred (the rows they are computed on), the coverage and CRPS columns described above when a Bayesian model ran, and model as its last column. Shared folds do not guarantee shared rows: a model that fails on a fold (GWR with a fixed bandwidth across a gap in the data, say) or predicts NA for some rows pools fewer rows, usually without the hardest ones. When the models predicted different rows, the function warns and recomputes every model's pooled metrics, your own metrics included, on the rows all of them predicted, so n_pred is the same on every row that has predictions; a model that predicted nothing stays an NA row. Each model's metrics over all the rows it predicted stay in its *_cv element and in attr(overall, "all_rows"), a table of the same shape. by_fold and the Bayesian coverage and CRPS columns are per fold and are not recomputed. Only the models that actually ran appear, so check which names are present; there is not always one entry per requested model, because a backend whose package is missing is dropped with a message. When no requested backend is available there is nothing to return and the function errors with "no viable models." instead of returning an empty comparison.

Coverage and CRPS in the overall table

A model can predict well on average and still be badly calibrated, so a comparison read from RMSE alone can prefer the model whose uncertainty is wrong. Heaton et al. (2019) found good point prediction routinely alongside poor interval coverage. When "Bayesian" is among the models that ran, overall therefore also carries the columns of cv_bayes()$predictive_coverage: coverage_50, coverage_80, coverage_95 (the share of held-out rows inside the posterior predictive interval at each level, averaged across folds weighted by each fold's n_pred) and mean_CRPS (the continuous ranked probability score, lower is better). They are NA on the GWR and RF rows, which produce point predictions and no draws, and absent when no Bayesian model ran. Read coverage against its nominal level: 0.95 at coverage_95 is calibrated, well below it is overconfident, well above it is wider than it needs to be.

Percentage errors on responses with zeros

MAPE divides by the observed value and SMAPE by $|y| + |\hat{y}|$, so neither is defined where its denominator is zero. Neither returns Inf or NaN. Both are averaged over the rows whose denominator is non-zero, and are NA when no row qualifies. Non-zero is judged at the scale of the data: a denominator no larger than 100 machine epsilons times the largest one counts as zero, so the rule does not depend on the units of the response. The n_MAPE and n_SMAPE columns record how many rows that was; the n column counts finite observation/prediction pairs. Read a percentage error next to its count: when n_MAPE < n, MAPE is an average over a subset of the data, whatever its value.

This bites on any response taking exact zeros: counts, rainfall, abundance, claim amounts. On a zero-inflated response with 62 zeros out of 120, MAPE is an average over the 58 non-zero rows, which n_MAPE = 58 now says. SMAPE fails differently and more subtly: it drops the rows where observation and prediction are both near zero (which on a well-fitted zero-inflated model are the rows it got right), so it averages the harder rows only and reads worse than the fit deserves; n_SMAPE shows how many rows it kept, and the count is only a label, not a repair.

RMSE, MAE and $R^2$ use every finite row and are unaffected; prefer them whenever the response can be zero. For a Bayesian fit, cv_bayes() additionally reports CRPS and interval coverage, which are proper scoring rules and have no such failure mode.

Which metrics survive a non-Gaussian response

RMSE and MAE are defined for any numeric response and are what to read for a count, a rate or a bounded outcome. MAPE and SMAPE assume a response that is rarely zero (see the previous section), and R-squared and adjusted R-squared compare residual variance to total variance, which is the right comparison for a Gaussian response and a loose one for anything whose variance tracks its mean. None of the four is wrong to compute; each is Gaussian-shaped thinking, and on a Poisson or zero-inflated response should be read as a rough summary rather than a score.

For the Bayesian backend, cv_bayes() additionally reports CRPS and interval coverage at 50, 80 and 95 percent. Both are proper scoring rules computed from posterior draws, so they are meaningful for any family that predicts one number per row (a count, a rate, a binary or bounded outcome), and they are the numbers to compare when the response is not Gaussian. A categorical or ordinal family predicts a probability per response category instead, so cv_bayes() refuses one before fitting anything. When every fold fails, the fold_metrics frame cv_bayes() returns carries the CRPS column but not the coverage_* columns, so code that reads those columns must tolerate their absence.

References

Heaton, M. J., Datta, A., Finley, A. O., Furrer, R., Guinness, J., Guhaniyogi, R., Gerber, F., Gramacy, R. B., Hammerling, D., Katzfuss, M., Lindgren, F., Nychka, D. W., Sun, F. and Zammit-Mangion, A. (2019). A case study competition among methods for analyzing large spatial data. Journal of Agricultural, Biological and Environmental Statistics, 24(3), 398–425. doi:10.1007/s13253-018-00348-w

See Also

Other model evaluation: compare_models(), evaluate_insample(), model_metrics(), residual_morans_i()

Examples

if (requireNamespace("ranger", quietly = TRUE)) {
  library(sf)
  set.seed(1)
  n <- 120
  dat <- st_as_sf(
    data.frame(x = 5e5 + runif(n, 0, 1000), y = 5e6 + runif(n, 0, 1000),
               elev = rnorm(n)),
    coords = c("x", "y"), crs = 32632
  )
  dat$price <- 10 + 0.01 * (st_coordinates(dat)[, 1] - 5e5) +
    2 * dat$elev + rnorm(n)
  # The random forest, and GWR beside it when GWmodel is installed, on the
  # same folds.  "Bayesian" is left out here because a Stan fit per fold
  # takes minutes; add it to `models` for a run to report.
  models <- c("RF", if (requireNamespace("GWmodel", quietly = TRUE) &&
                       requireNamespace("sp", quietly = TRUE)) "GWR")
  cmp <- compare_models_cv(dat, "price", "elev", models = models, k = 3,
                           rf_args = list(num_trees = 100),
                           gwr_args = list(bandwidth = 30))
  cmp$overall
}

Description

Lays a regular grid of equal-area cells over boundary and clips it to that boundary. Reach for this rather than create_voronoi_polygons() when cell size should be a decision you make, instead of one dictated by where the observations happen to be. That is the case when you need per-cell rates comparable across the map, or a resolution that stays fixed as the sample grows. Hexagons (type = "hex") avoid the axis-aligned artefacts of squares and give every cell the same distance to all six neighbours, which matters for anything that reads neighbourhoods.

Usage

create_grid_polygons(
  boundary,
  target_cells = NULL,
  type = c("square", "hex"),
  cellsize = NULL,
  n = NULL,
  clip = TRUE,
  crs = NULL,
  quiet = FALSE,
  max_cells = 1e+06
)

Arguments

Details

Size the grid with exactly one of target_cells (roughly how many cells you want, the package derives the rest), cellsize (a fixed edge length in CRS units) or n (a fixed number of columns and rows). See ⁠@param cellsize⁠ for what happens when more than one is given.

Value

An sf polygon layer with poly_id column.

See Also

Other tessellation: build_tessellation(), create_grid_polygons_cached(), create_voronoi_polygons(), ensure_stable_poly_id(), get_voronoi_seeds(), plot_tessellation_map(), voronoi_seeds_kmeans(), voronoi_seeds_random()

Examples

library(sf)
bnd <- st_sf(geometry = st_sfc(st_polygon(list(rbind(
  c(0, 0), c(100, 0), c(100, 100), c(0, 100), c(0, 0)
))), crs = 32632))
grid_sq  <- create_grid_polygons(bnd, target_cells = 100, type = "square")
grid_hex <- create_grid_polygons(bnd, target_cells = 100, type = "hex")
nrow(grid_sq)
# Hex counts run above target because clipping keeps every hexagon that
# merely overhangs the boundary; the inflation is proportionally larger
# at small target_cells.
nrow(grid_hex)

Description

Builds a grid via create_grid_polygons() and memoizes the result so repeated calls with the same inputs return instantly.

Usage

create_grid_polygons_cached(
  boundary,
  target_cells = NULL,
  type = c("square", "hex"),
  ...,
  cache_env = .gmt_cache,
  max_entries = 50L
)

Arguments

Value

An sf data frame with a stable poly_id column. The rows are re-ordered and re-numbered by ensure_stable_poly_id, which create_grid_polygons does not do: the same cell therefore carries a different poly_id depending on which of the two builders produced it. Use one builder throughout an analysis; joining a summary keyed on IDs from one onto geometries from the other draws the values on the wrong polygons.

See Also

Other tessellation: build_tessellation(), create_grid_polygons(), create_voronoi_polygons(), ensure_stable_poly_id(), get_voronoi_seeds(), plot_tessellation_map(), voronoi_seeds_kmeans(), voronoi_seeds_random()

Other package options and caches: clear_fitted_cache(), clear_grid_cache(), spatialkit_quiet()

Examples

library(sf)
bnd <- st_sf(geometry = st_sfc(st_polygon(list(rbind(
  c(0, 0), c(100, 0), c(100, 100), c(0, 100), c(0, 0)
))), crs = 32632))
g <- create_grid_polygons_cached(bnd, target_cells = 16, type = "hex")
nrow(g)
# The IDs come from ensure_stable_poly_id(), so they follow the geometry:
# the same request with the boundary's vertices in another order gives the
# same ID to the same cell.
bnd2 <- st_sf(geometry = st_sfc(st_polygon(list(rbind(
  c(100, 100), c(0, 100), c(0, 0), c(100, 0), c(100, 100)
))), crs = 32632))
g2 <- create_grid_polygons_cached(bnd2, target_cells = 16, type = "hex")
same_cell <- match(st_as_text(st_geometry(g2)), st_as_text(st_geometry(g)))
all(g2$poly_id == g$poly_id[same_cell])

Description

Assigns every location in the study area to its nearest input point, giving one cell per point. This is the tessellation to reach for when the observations themselves define the regions of interest (sampling sites, monitoring stations, service points), because cell size then adapts to sampling density instead of being imposed by a fixed grid: dense areas get small cells and sparse areas large ones. Prefer create_grid_polygons() instead when you need equal-area cells or a resolution independent of where the data happen to be.

Usage

create_voronoi_polygons(
  points_sf,
  boundary = NULL,
  expand = 0,
  clip = TRUE,
  keep_duplicates = FALSE,
  crs = NULL,
  quiet = FALSE
)

Arguments

Details

The heavy lifting is sf::st_voronoi(). What this adds is the surrounding bookkeeping: projecting lon/lat input, building and buffering an envelope so edge cells are bounded, clipping to boundary, restoring the point-to-cell correspondence that st_voronoi() scrambles, and stamping stable cell_id values.

The generators are the points' vertices, not the features. A MULTIPOINT feature with several vertices therefore gets one cell per vertex, and its index entry is the smallest cell_id among the cells it touches; the others are referenced by no feature. Other functions in the package (prep_model_data(), make_folds()) reduce such a feature to its centroid instead, so cast to POINT, or take centroids, first if one cell per feature is what you want.

Value

A list with cells, index, boundary, method and params. index holds one cell_id per row of points_sf, and NA for a point that falls outside every cell, which means outside the study area (grown by expand when it is positive), so a summary built from it counts only the points the tessellation actually covers. boundary is the boundary the cells were built in: the one supplied or derived, grown by expand.

See Also

Other tessellation: build_tessellation(), create_grid_polygons(), create_grid_polygons_cached(), ensure_stable_poly_id(), get_voronoi_seeds(), plot_tessellation_map(), voronoi_seeds_kmeans(), voronoi_seeds_random()

Examples

library(sf)
set.seed(1)
pts <- st_as_sf(
  data.frame(x = 5e5 + runif(15, 0, 100), y = 5e6 + runif(15, 0, 100)),
  coords = c("x", "y"), crs = 32632
)
res <- create_voronoi_polygons(pts, quiet = TRUE)
res$cells   # one polygon per unique point, with stable cell_id
res$index   # cell_id assignment for each input point

Description

Refits the Gaussian-process model of fit_bayesian_spatial_model() on each training fold and scores it on the held-out fold. Beyond the point-prediction metrics the other CV wrappers report, this one scores the whole predictive distribution: predictive_coverage says what fraction of held-out observations fell inside the 50/80/95% intervals, and mean_CRPS rates sharpness and calibration together. That is the reason to reach for it. A Bayesian model is usually chosen for its uncertainty, and only held-out coverage shows whether those intervals are honest at locations the model has not seen.

Usage

cv_bayes(
  data_sf,
  response_var,
  predictor_vars,
  folds = NULL,
  k = 5,
  seed = 123,
  boundary = NULL,
  pointize = "auto",
  fit_args = list(),
  summary = c("mean", "median"),
  compute_pred_intervals = TRUE,
  coverage_levels = c(0.5, 0.8, 0.95),
  block_size = NULL,
  auto_range = FALSE,
  parallel = FALSE,
  metrics = NULL
)

Arguments

Details

It is the most expensive wrapper in the package by a wide margin: every fold is a full MCMC run. Use few folds, and parallel = TRUE if you have the cores. For a cheap first pass on the same question, cross-validate a forest with cv_rf() and come back here once the predictor set has settled.

Value

A list with overall, fold_metrics, predictions, folds, n_folds_attempted, n_folds_succeeded, fold_status, orphan_rows, n_unknown_ids, n_dropped, formula, predictive_coverage and coverage_levels (the nominal levels, named by their coverage_* column). The two fold counts make a run where every fold failed visible in the return value itself, beyond the warning, and fold_status (one row per fold: fold, status, message) keeps the reason each missing fold is missing, including the error text of a fold whose sampler failed, where a long run's console output would not. See cv_spatial for the five statuses and for orphan_rows. predictions carries, beyond the columns its siblings share, yhat_sd: the posterior predictive standard deviation of each held-out row, from the same draws that give the coverage below (NA when compute_pred_intervals = FALSE or the draws failed for that fold). overall$Adj_R2 is always NA, as for every cv_*(): see cv_spatial. fold_metrics carries, beyond the columns its siblings share, gp_k, gp_n_basis, n_draws, CRPS, the coverage_* columns and convergence_ok: TRUE or FALSE as fit_bayesian_spatial_model() judged that fold's sampler (R-hat, effective sample size, divergences), NA when fit_args sets check_convergence = FALSE. A fold that did not converge is scored like the others, so a run with any FALSE raises one warning naming those folds. The predictive_coverage entries (one per coverage_levels value, plus mean_CRPS) are averages across folds weighted by each fold's n_pred, because the per-fold values in fold_metrics are themselves means over that fold's test rows; an unweighted average would not be the pooled quantity when fold sizes differ, which for spatially blocked folds they routinely do.

Percentage errors on responses with zeros

MAPE divides by the observed value and SMAPE by $|y| + |\hat{y}|$, so neither is defined where its denominator is zero. Neither returns Inf or NaN. Both are averaged over the rows whose denominator is non-zero, and are NA when no row qualifies. Non-zero is judged at the scale of the data: a denominator no larger than 100 machine epsilons times the largest one counts as zero, so the rule does not depend on the units of the response. The n_MAPE and n_SMAPE columns record how many rows that was; the n column counts finite observation/prediction pairs. Read a percentage error next to its count: when n_MAPE < n, MAPE is an average over a subset of the data, whatever its value.

This bites on any response taking exact zeros: counts, rainfall, abundance, claim amounts. On a zero-inflated response with 62 zeros out of 120, MAPE is an average over the 58 non-zero rows, which n_MAPE = 58 now says. SMAPE fails differently and more subtly: it drops the rows where observation and prediction are both near zero (which on a well-fitted zero-inflated model are the rows it got right), so it averages the harder rows only and reads worse than the fit deserves; n_SMAPE shows how many rows it kept, and the count is only a label, not a repair.

RMSE, MAE and $R^2$ use every finite row and are unaffected; prefer them whenever the response can be zero. For a Bayesian fit, cv_bayes() additionally reports CRPS and interval coverage, which are proper scoring rules and have no such failure mode.

Which metrics survive a non-Gaussian response

RMSE and MAE are defined for any numeric response and are what to read for a count, a rate or a bounded outcome. MAPE and SMAPE assume a response that is rarely zero (see the previous section), and R-squared and adjusted R-squared compare residual variance to total variance, which is the right comparison for a Gaussian response and a loose one for anything whose variance tracks its mean. None of the four is wrong to compute; each is Gaussian-shaped thinking, and on a Poisson or zero-inflated response should be read as a rough summary rather than a score.

For the Bayesian backend, cv_bayes() additionally reports CRPS and interval coverage at 50, 80 and 95 percent. Both are proper scoring rules computed from posterior draws, so they are meaningful for any family that predicts one number per row (a count, a rate, a binary or bounded outcome), and they are the numbers to compare when the response is not Gaussian. A categorical or ordinal family predicts a probability per response category instead, so cv_bayes() refuses one before fitting anything. When every fold fails, the fold_metrics frame cv_bayes() returns carries the CRPS column but not the coverage_* columns, so code that reads those columns must tolerate their absence.

See Also

Other cross-validation: area_of_applicability(), cv_block_size_sweep(), cv_gwr(), cv_rf(), cv_spatial(), estimate_sac_range(), fold_separation(), gwr_model_selection(), make_folds(), sac_nugget(), select_features_forward()

Examples

## Not run: 
# Not run: fits with Stan, which needs a working C++ toolchain and takes
# minutes of MCMC -- both outside what an example may assume.
if (requireNamespace("brms", quietly = TRUE)) {
  library(sf)
  set.seed(1)
  n <- 60
  dat <- st_as_sf(
    data.frame(x = 5e5 + runif(n, 0, 1000), y = 5e6 + runif(n, 0, 1000),
               elev = rnorm(n)),
    coords = c("x", "y"), crs = 32632
  )
  dat$price <- 10 + 0.01 * (st_coordinates(dat)[, 1] - 5e5) +
    2 * dat$elev + rnorm(n)
  # Two short chains per fold keep this to a few minutes and leave the
  # posterior rough, so the intervals below will run narrow; a run to
  # report uses brms's defaults (chains = 4, iter = 2000).
  cv <- cv_bayes(dat, "price", "elev", k = 2,
                 fit_args = list(chains = 2, iter = 1000))
  print(cv$overall)
  cv$predictive_coverage  # coverage at 50/80/95% plus mean CRPS
}

## End(Not run)

Description

A blocked cross-validation with blocks smaller than the autocorrelation range leaks: every held-out point has a near-identical neighbour in the training set, and the score is optimistic in proportion. The single number a cv_*() call returns cannot show this. This runs the same cross-validation at a ladder of block sizes, and by default once with random folds as the fully leaky reference. It returns the metric at each, with the estimated autocorrelation range alongside so that the curve can be read against it: it rises as the blocks pass the range and then plateaus, and the height of the rise is how much the random-fold number overstated the model.

Usage

cv_block_size_sweep(
  data_sf,
  response_var,
  predictor_vars,
  fit_fn,
  block_sizes = NULL,
  n_sizes = 6L,
  k = 5L,
  metric = "RMSE",
  include_random = TRUE,
  max_fits = 60L,
  sac = NULL,
  seed = 123L,
  quiet = FALSE,
  ...
)

Arguments

Value

A data.frame of class "block_size_sweep" with one row per cross-validation: block_size (NA for the random reference), method, blocks_used, k (the folds actually built), n_folds_succeeded, value (the pooled metric), fold_min, fold_max and fold_sd (its spread across folds). Attributes: metric, sac_range (the effective range, or NA), crs, k (the folds requested, which print() and plot() report), n_fits, response_var, and results, the full cv_spatial() result at every size. plot() draws it.

The fit budget

Each block size is a full cross-validation, so the cost is length(block_sizes) * k fits, plus k for the random reference. max_fits caps that (default 60: room for up to eleven sizes at k = 5 plus the reference; the default six-size ladder at k = 5 needs at most 35). A sweep that would run past the cap refuses to start, naming the number of fits it would have needed. Raise max_fits deliberately; the RF example below takes seconds, a Bayesian fit_fn takes minutes per fit.

The ladder

When block_sizes is NULL, n_sizes values are log-spaced from a twenty-fifth of the shorter side of the data's extent (of the longer side when the points lie on one line parallel to an axis) to the largest size, at most half that side, at which the grid still holds k blocks. Half the side gives a grid two blocks across, enough for k up to 4 and, for larger k, on an extent long enough in the other direction; on a roughly square extent at the default k = 5 the top is about a third of the side. Any size at which the grid would hold more than the 1,000,000 make_folds() will build is dropped. The count is of grid cells: on clustered data a grid of k or more cells can have fewer than k that hold points, and make_folds() then lowers k at that size, which the k column shows. On an extent much longer than it is wide every rung can fall below the autocorrelation range while longer blocks would still fit k times along the longer side; the sweep warns when that happens, and block_sizes is then the way to reach past the range. Sizes are in the units of the CRS the folds are built in (make_folds()'s params$crs, metres for geographic input), and the returned table records that CRS. Each is the minimum block edge handed to make_folds(), which fits a whole number of cells across the extent, so the cells are somewhat longer than the size on the axis (up to twice as long).

See Also

plot.block_size_sweep().

Other cross-validation: area_of_applicability(), cv_bayes(), cv_gwr(), cv_rf(), cv_spatial(), estimate_sac_range(), fold_separation(), gwr_model_selection(), make_folds(), sac_nugget(), select_features_forward()

Examples

if (requireNamespace("ranger", quietly = TRUE) &&
    requireNamespace("gstat", quietly = TRUE)) {
  library(sf)
  set.seed(1)
  n <- 200
  x <- 5e5 + runif(n, 0, 1000); y <- 5e6 + runif(n, 0, 1000)
  d <- as.matrix(dist(cbind(x, y)))
  field <- as.numeric(t(chol(exp(-d / 100) + diag(1e-6, n))) %*% rnorm(n))
  dat <- st_as_sf(data.frame(x = x, y = y, a = rnorm(n)), coords = c("x", "y"),
                  crs = 32632)
  dat$z <- field + 0.5 * dat$a + rnorm(n, 0, 0.2)
  rf_fn <- function(train_sf)
    fit_rf_model(train_sf, "z", "a", include_coords = TRUE, num_trees = 100, seed = 1)
  sw <- cv_block_size_sweep(dat, "z", "a", fit_fn = rf_fn, k = 4, n_sizes = 4,
                            quiet = TRUE)
  print(sw)               # the curve as a table, with the random-fold reference
  if (requireNamespace("ggplot2", quietly = TRUE)) plot(sw)
}

Description

Refits a geographically weighted regression from scratch on each training fold and scores it on the held-out fold, so the reported error is what the model achieves at locations it did not see. Reach for it whenever you need a defensible accuracy figure for a GWR: the in-sample $R^2$ that model_metrics() reports on a gwr_fit is close to meaningless, because a local regression with a small bandwidth can track the training points almost exactly. Bandwidth is re-selected per fold unless you fix it with bandwidth, which keeps the selection itself inside the cross-validation, with no tuning on the full data first.

Usage

cv_gwr(
  data_sf,
  response_var,
  predictor_vars,
  folds = NULL,
  k = 5,
  seed = 123,
  adaptive = TRUE,
  bandwidth = NULL,
  kernel = c("bisquare", "gaussian", "tricube", "boxcar", "exponential"),
  boundary = NULL,
  pointize = "auto",
  block_size = NULL,
  auto_range = FALSE,
  parallel = FALSE,
  metrics = NULL
)

Arguments

Details

Folds default to spatial blocks (make_folds(method = "block_kfold")), not random ones. With autocorrelated data a random split leaves a held-out point's neighbours in the training set and the score comes back flattering. Use cv_bayes() for the same treatment of a Bayesian GP model, cv_rf() for a forest, and compare_models_cv() to score several backends on one set of folds.

Value

A list with overall, fold_metrics, predictions, folds, n_folds_attempted, n_folds_succeeded, fold_status, orphan_rows, n_unknown_ids, n_dropped, formula and adaptive. The two fold counts make a run where every fold failed visible in the return value itself, beyond the warning, since overall is a well-formed all-NA row either way, and fold_status (one row per fold: fold, status, message) says why each missing fold is missing. See cv_spatial for the five statuses and for orphan_rows. Adj_R2 is NA in both overall and fold_metrics: the pooled predictions have no single parameter count, and a GWR's effective parameter count is not its predictor count (see cv_spatial).

Percentage errors on responses with zeros

MAPE divides by the observed value and SMAPE by $|y| + |\hat{y}|$, so neither is defined where its denominator is zero. Neither returns Inf or NaN. Both are averaged over the rows whose denominator is non-zero, and are NA when no row qualifies. Non-zero is judged at the scale of the data: a denominator no larger than 100 machine epsilons times the largest one counts as zero, so the rule does not depend on the units of the response. The n_MAPE and n_SMAPE columns record how many rows that was; the n column counts finite observation/prediction pairs. Read a percentage error next to its count: when n_MAPE < n, MAPE is an average over a subset of the data, whatever its value.

This bites on any response taking exact zeros: counts, rainfall, abundance, claim amounts. On a zero-inflated response with 62 zeros out of 120, MAPE is an average over the 58 non-zero rows, which n_MAPE = 58 now says. SMAPE fails differently and more subtly: it drops the rows where observation and prediction are both near zero (which on a well-fitted zero-inflated model are the rows it got right), so it averages the harder rows only and reads worse than the fit deserves; n_SMAPE shows how many rows it kept, and the count is only a label, not a repair.

RMSE, MAE and $R^2$ use every finite row and are unaffected; prefer them whenever the response can be zero. For a Bayesian fit, cv_bayes() additionally reports CRPS and interval coverage, which are proper scoring rules and have no such failure mode.

See Also

Other cross-validation: area_of_applicability(), cv_bayes(), cv_block_size_sweep(), cv_rf(), cv_spatial(), estimate_sac_range(), fold_separation(), gwr_model_selection(), make_folds(), sac_nugget(), select_features_forward()

Examples

if (requireNamespace("GWmodel", quietly = TRUE) &&
    requireNamespace("sp", quietly = TRUE)) {
  library(sf)
  set.seed(1)
  n <- 60
  dat <- st_as_sf(
    data.frame(x = 5e5 + runif(n, 0, 1000), y = 5e6 + runif(n, 0, 1000),
               elev = rnorm(n)),
    coords = c("x", "y"), crs = 32632
  )
  dat$price <- 10 + 0.01 * (st_coordinates(dat)[, 1] - 5e5) +
    2 * dat$elev + rnorm(n)
  cv <- cv_gwr(dat, "price", "elev", k = 3, bandwidth = 30)
  print(cv$overall)       # pooled over the held-out rows of every fold
  cv$fold_metrics         # and fold by fold
}

Description

Thin wrapper over cv_spatial that refits a ranger forest on each training fold. Unlike the out-of-bag error, this holds out whole spatial blocks, so neighbours of a held-out point are not sitting in the training set.

Usage

cv_rf(
  data_sf,
  response_var,
  predictor_vars,
  folds = NULL,
  k = 5,
  seed = 123,
  parallel = FALSE,
  block_size = NULL,
  auto_range = FALSE,
  boundary = NULL,
  pointize = "auto",
  metrics = NULL,
  ...
)

Arguments

Value

The cv_spatial result.

Percentage errors on responses with zeros

MAPE divides by the observed value and SMAPE by $|y| + |\hat{y}|$, so neither is defined where its denominator is zero. Neither returns Inf or NaN. Both are averaged over the rows whose denominator is non-zero, and are NA when no row qualifies. Non-zero is judged at the scale of the data: a denominator no larger than 100 machine epsilons times the largest one counts as zero, so the rule does not depend on the units of the response. The n_MAPE and n_SMAPE columns record how many rows that was; the n column counts finite observation/prediction pairs. Read a percentage error next to its count: when n_MAPE < n, MAPE is an average over a subset of the data, whatever its value.

This bites on any response taking exact zeros: counts, rainfall, abundance, claim amounts. On a zero-inflated response with 62 zeros out of 120, MAPE is an average over the 58 non-zero rows, which n_MAPE = 58 now says. SMAPE fails differently and more subtly: it drops the rows where observation and prediction are both near zero (which on a well-fitted zero-inflated model are the rows it got right), so it averages the harder rows only and reads worse than the fit deserves; n_SMAPE shows how many rows it kept, and the count is only a label, not a repair.

RMSE, MAE and $R^2$ use every finite row and are unaffected; prefer them whenever the response can be zero. For a Bayesian fit, cv_bayes() additionally reports CRPS and interval coverage, which are proper scoring rules and have no such failure mode.

See Also

Other cross-validation: area_of_applicability(), cv_bayes(), cv_block_size_sweep(), cv_gwr(), cv_spatial(), estimate_sac_range(), fold_separation(), gwr_model_selection(), make_folds(), sac_nugget(), select_features_forward()

Examples

if (requireNamespace("ranger", quietly = TRUE)) {
  library(sf)
  set.seed(1)
  n <- 200
  dat <- st_as_sf(
    data.frame(x = 5e5 + runif(n, 0, 1000), y = 5e6 + runif(n, 0, 1000),
               a = rnorm(n)),
    coords = c("x", "y"), crs = 32632
  )
  dat$z <- 2 * dat$a + rnorm(n, 0, 0.3)
  cv_rf(dat, "z", "a", k = 4)$overall
}

Description

Run K-fold CV for any model that returns a spatial_fit object. This is the extensibility point: to plug in a new model type, supply a fit_fn(train_sf) that returns a spatial_fit.

Usage

cv_spatial(
  data_sf,
  response_var,
  predictor_vars,
  fit_fn,
  folds = NULL,
  k = 5,
  seed = 123,
  boundary = NULL,
  pointize = "auto",
  predict_args = list(),
  fold_info_fn = NULL,
  p = NULL,
  block_size = NULL,
  auto_range = FALSE,
  parallel = FALSE,
  metrics = NULL,
  .caller = "cv_spatial"
)

Arguments

Value

A list with overall, fold_metrics, predictions, folds, the two fold counts n_folds_attempted and n_folds_succeeded, and four elements that say what happened to the difference between them and to the rows: fold_status, orphan_rows, n_unknown_ids and n_dropped. The counts are reported deliberately: a fit_fn that fails on every fold otherwise looks like a successful run that happened to score NA, so compare them before trusting overall. fold_status is a data.frame with one row per fold supplied (fold, status and message), where status is "ok" (message is empty unless fold_info_fn failed on the fold); "error" (the fit or its predict() threw; message is the error text); "skipped" (nothing scorable: too few matched rows, a prediction of the wrong length, or no finite observed/predicted pair); "dropped" (an empty test set, or fewer than two training rows, once incomplete rows were removed, so the fold never reached the fitter); or "worker_error" (a parallel worker died, for example killed for lack of memory). Each fold runs in its own worker, so a failure costs that fold only, and an error that stops a sequential run (a metrics or fold_info_fn return value of the wrong shape) stops a parallel one too, naming the fold. Every fold missing from fold_metrics has its reason there, which matters most when the console output of a long run is gone. When some folds, but not all, end as "error", "skipped" or "worker_error", the function warns, naming them and how many rows overall covers: it is pooled over the folds that produced predictions, and the fold that fails is often the hardest to predict, so it may flatter the model. (A "dropped" fold has its own warning.) orphan_rows holds the ..row_ids of rows in the data that no fold names (they enter no training set and are never scored; non-empty only when the folds were built on a different or subsetted layer), and n_unknown_ids counts the distinct row IDs the folds name that the data does not have. Each such row is named by every fold, once as a test row and once in each other fold's training set, and this counts the row, not the mentions (expected when rows were removed for missing values). n_dropped is how many rows prep_model_data() removed for missing or non-finite values or a bad geometry before any fold was fitted. The fold column of fold_metrics, predictions and fold_status carries the fold's index in the folds object that was supplied, so it lines up with make_folds()$assignment$fold even when some folds were unusable and dropped. Splits that already carry a fold_id, as the folds of a cv_*() result do, keep it: handing one run's folds to another labels each fold as the first run and fold_separation() do, a dropped fold's gap included. R2 is out-of-sample $R^2$: the total sum of squares is taken about the mean of the training rows (in overall, each held-out row about its own fold's training mean), the null prediction available when the fold is predicted, and not about the held-out rows' own mean, a null model that would know the test data. On spatial blocks the two can differ widely; R2 is below 0 when the model predicts worse than the training mean. model_metrics(newdata = ) uses the same baseline. overall$Adj_R2 is always NA: the pooled out-of-sample predictions come from k separately fitted models and have no single parameter count to adjust for. The per-fold fold_metrics$Adj_R2 carries the adjusted value when p is supplied, and is NA otherwise.

Name collision with blockCV

The blockCV package exports a function of the same name that does the opposite job: blockCV::cv_spatial() builds spatial folds, where this function runs a cross-validation over folds it is given. With both packages attached, whichever was attached last masks the other; spatialkit::cv_spatial() always resolves to this one. The two cooperate: blockCV::cv_spatial() returns its fold assignment as $folds_ids, a vector of fold labels, and that vector is accepted directly as the folds argument here and in every other cv_*() function. Fold construction is blockCV's home ground and this package does not try to match its breadth there; what this package adds is the path from irregular points through data-drawn regions to a cross-validated, compared model.

Your own metrics

The built-in columns (RMSE, MAE, MAPE, SMAPE, R2, Adj_R2) are the Gaussian regression set, and the section below says which of them survive a count or a bounded response. metrics is the way to score what they cannot: a function(y, yhat) that returns a named numeric vector (a named list of scalars, or a one-row data frame, also serve), for example a Poisson deviance, a log score on a probability, or a loss with your own weights. It is applied twice, in the same way the built-in metrics are: once per fold, to that fold's held-out rows, so each name becomes a column of fold_metrics; and once to the pooled out-of-sample predictions of every fold, so each name becomes a column of overall. Only the pairs the built-in metrics use reach the function (both y and yhat finite), so its columns describe the same rows as RMSE, and n_pred counts them.

The contract: every element named, names unique and not one of the built-in column names, one number per name. The built-in names include the per-fold extras of the backend or of your fold_info_fn (bandwidth for cv_gwr(); CRPS, coverage_*, gp_k, gp_n_basis and n_draws for cv_bayes()), and mean_CRPS, which compare_models_cv() writes. Anything else is an error, because a scoring function that returns the wrong shape is a mistake to surface instead of a fold to skip. A function that throws on a fold is logged and its columns are NA for that fold (and for overall, if it throws on the pooled predictions); a fold is never dropped for it. When no fold produced a prediction the empty fold_metrics frame still carries the function's columns, typed, provided the function can be called on zero-length input.

fold_info_fn is the per-fold half of the same mechanism, with access to the fitted object and the held-out layer; metrics sees only the two vectors but is also pooled. compare_models_cv() hands one metrics to every backend, so the columns are comparable across the rows of its overall.

Percentage errors on responses with zeros

MAPE divides by the observed value and SMAPE by $|y| + |\hat{y}|$, so neither is defined where its denominator is zero. Neither returns Inf or NaN. Both are averaged over the rows whose denominator is non-zero, and are NA when no row qualifies. Non-zero is judged at the scale of the data: a denominator no larger than 100 machine epsilons times the largest one counts as zero, so the rule does not depend on the units of the response. The n_MAPE and n_SMAPE columns record how many rows that was; the n column counts finite observation/prediction pairs. Read a percentage error next to its count: when n_MAPE < n, MAPE is an average over a subset of the data, whatever its value.

This bites on any response taking exact zeros: counts, rainfall, abundance, claim amounts. On a zero-inflated response with 62 zeros out of 120, MAPE is an average over the 58 non-zero rows, which n_MAPE = 58 now says. SMAPE fails differently and more subtly: it drops the rows where observation and prediction are both near zero (which on a well-fitted zero-inflated model are the rows it got right), so it averages the harder rows only and reads worse than the fit deserves; n_SMAPE shows how many rows it kept, and the count is only a label, not a repair.

RMSE, MAE and $R^2$ use every finite row and are unaffected; prefer them whenever the response can be zero. For a Bayesian fit, cv_bayes() additionally reports CRPS and interval coverage, which are proper scoring rules and have no such failure mode.

See Also

new_spatial_fit() for the constructor a fit_fn must use; cv_gwr(), cv_bayes() and cv_rf() for the built-in backends, which are thin wrappers over this function.

Other cross-validation: area_of_applicability(), cv_bayes(), cv_block_size_sweep(), cv_gwr(), cv_rf(), estimate_sac_range(), fold_separation(), gwr_model_selection(), make_folds(), sac_nugget(), select_features_forward()

Examples

library(sf)
set.seed(1)
n <- 80
site <- st_as_sf(
  data.frame(x = 5e5 + runif(n, 0, 1000), y = 5e6 + runif(n, 0, 1000),
             elev = rnorm(n)),
  coords = c("x", "y"), crs = 32632
)
site$price <- 10 + 0.01 * (st_coordinates(site)[, 1] - 5e5) +
  2 * site$elev + rnorm(n)

# 1. A fit_fn returning a spatial_fit of your own subclass.
lm_fit <- function(train_sf) {
  new_spatial_fit(
    subclass       = "lm_fit",
    engine         = lm(price ~ elev, st_drop_geometry(train_sf)),
    formula        = price ~ elev,
    response_var   = "price",
    predictor_vars = "elev",
    data_sf        = train_sf
  )
}

# 2. The predict() method cv_spatial() scores each fold with. Without it
#    every fold fails and `overall` comes back all-NA.
predict.lm_fit <- function(object, newdata = NULL, ...) {
  if (is.null(newdata)) newdata <- object$data_sf
  as.numeric(stats::predict(object$engine, st_drop_geometry(newdata)))
}
registerS3method("predict", "lm_fit", predict.lm_fit)

cv <- cv_spatial(site, "price", "elev", fit_fn = lm_fit, k = 3, seed = 1)
cv$overall
# Compare these before trusting the metrics above; fold_status says why
# any fold is missing from fold_metrics.
c(attempted = cv$n_folds_attempted, succeeded = cv$n_folds_succeeded)
cv$fold_status

# 3. A metric of your own beside the built-in ones: per fold and pooled.
med_ae <- function(y, yhat) c(MedAE = stats::median(abs(y - yhat)))
cv2 <- cv_spatial(site, "price", "elev", fit_fn = lm_fit, k = 3, seed = 1,
                  metrics = med_ae)
cv2$overall$MedAE
cv2$fold_metrics[, c("fold", "RMSE", "MedAE")]

Description

Computes a WSS curve over k=1..K_max using k-means on projected feature coordinates and selects candidate k values around the elbow.

Usage

determine_optimal_levels(
  data_sf,
  max_levels = 12L,
  top_n = 3L,
  sample_n = 1500L,
  set_seed = 123L,
  response_var = NULL,
  predictor_vars = NULL,
  criterion = c("geometric", "morans_i", "combined"),
  select_on = c("all", "split")
)

Arguments

Details

The elbow is read on log-log axes, and there may be none. Points with no cluster structure have a WSS curve close to $c/k$, and the classical rule, the point furthest below the chord from the first to the last k, still finds a "knee" on it on linear axes, at about $\sqrt{K_{max}}$: 4 at the default max_levels = 12 and 13 at 160, whatever the data. On $\log k$ against $\log$ WSS that curve is a straight line, while separated clusters fall faster than it until there is one cell per cluster and like it after, a bend at the cluster count. The elbow is therefore the k whose $\log$ WSS sags furthest below the straight line joining k = 1 and $K_{max}$ on those axes, and it counts as one only when the sag is at least $\log 1.25$ (the WSS a fifth below the power law through the ends). Measured on 60 to 1500 points with ladders to 3–40, uniform layouts over squares, discs, triangles, an L-shape, density gradients and jittered lattices sagged at most 0.16, and two to ten separated clusters at least 0.6 once max_levels passed the cluster count. With no elbow the function warns and returns the linear-axis answer, which the ladder chose, not the data. An elongated extent also bends, at about its aspect ratio, because the first cuts go across its long axis (0.2–0.45 for a 4:1 rectangle); past the threshold that bend is reported as an elbow, and it describes the extent's shape rather than clusters in it. When locations repeat (stations visited many times), the sweep can reach one cell per distinct location, where the WSS is zero (to within $10^{-12}$ of the total). That k has no place on log axes and is left out of the line, but the fall to zero is the sharpest bend there is: when the rest of the curve has no elbow, the number of distinct locations is the elbow.

When response_var and predictor_vars are provided, the geometric WSS elbow is supplemented with Moran's I computed on OLS residuals at each candidate k. The Moran's I profile measures how much spatial autocorrelation in the response remains unexplained at a given tessellation resolution. It is a direct reflection of the spatial process being modeled instead of the mere geometric compactness of coordinates. The combined criterion selects the k that best balances geometric parsimony and residual spatial independence.

To keep memory use and runtime bounded for large max_levels, the initial k-means sweep records only within-cluster sum-of-squares (WSS) without retaining cluster assignments. Moran's I is then evaluated lazily: k-means is re-run only for a focused neighbourhood around the elbow (±4 by default, or ±top_n if larger), so that only the most promising candidate k values incur the cost of the full Moran's I computation.

The WSS curve is read for its shape, so it is fitted to be smooth. Every point on it is a k-means local optimum, and with a few random restarts the curve mixes level effects with optimisation noise: measured on clustered layouts, a sweep over k = 1..30 at stats::kmeans(nstart = 5) rose at one or two steps in three of five draws, and an earlier form of the elbow rule once selected such a bump. Each k is therefore fitted as the best of 25 restarts seeded by k-means++ (Arthur and Vassilvitskii 2007), the budget at which the gain from further restarts saturates (Fränti and Sieranoja 2019; Steinley 2003 on why the usual handful is not enough). The same sweeps then had no increase at all. A curve that still rises somewhere is reported but not refused: a bumpy curve is uncertain, not unidentified. A logged warning names the number of rising steps, and the model-aware paths return it as wss_bumps in the "diagnostics" attribute beside wss_spread, the relative spread of WSS across the restarts at each k. Because the optimiser changed, a selection made by an earlier version on a curve that had such a bump can differ from the one made now; where the earlier curve was clean, the selection is the same.

The model-aware criteria rank on the standardised deviate, not on |Moran's I|. Both $E[I]$ and $Var[I]$ depend on the number of cells, so $|I|$ falls as k grows whether or not the finer tessellation is capturing anything. Measured over 300 replicates of a response with no spatial structure, mean $|I|$ fell monotonically from 0.114 at k = 10 to 0.050 at k = 60 ($-56\%$), which made an $|I|$ ranking prefer the largest candidate for arithmetic reasons alone. Candidates are therefore ordered by $|z| = |I - E[I]| / \mathrm{sd}(I)$ using the Cliff & Ord regression residual moments. The cell-level residuals are OLS residuals by construction, which is the case those moments are derived for, but the derivation also assumes errors of equal variance, and a cell mean over $n_j$ points has a variance proportional to $1/n_j$. Over the same runs $z$ had mean $\approx 0$, $\mathrm{sd} \approx 1$ and a two-sided 5% rejection rate of 0.040–0.057 at every k, and it stayed calibrated on gradient and moderately clustered layouts; with single-point cells next to cells of 70 or more points its mean rose to 0.2–0.34 and its rejection rate to 7–8% at 20 cells. Where structure remains, $|z|$ mixes its size with the number of cells it is measured on, since $\mathrm{sd}(I)$ shrinks as cells are added. Both quantities are reported in the "diagnostics" attribute, as moran_i and moran_z.

Resolution floor on the model-aware criteria. Moran's I is computed on cell-level residuals with an 8-nearest-neighbour weight matrix, so it only carries information once there are more than nine cells. At nine or fewer, every cell is a neighbour of every other, the row-standardised weight matrix is complete, and Moran's I collapses to exactly $-1/(k - 1)$ for any residual vector (a function of k alone). The criterion ranks on $|z|$, not on $|I|$, and at the floor the residual moments give $E[I] = I$ and $\mathrm{Var}[I] = 0$ identically (the algebra holds to $10^{-16}$), so the standardised deviate is $0/0$: it carries no information about the tessellation, and whichever way rounding noise resolves it those candidates would rank first or last on nothing. They therefore return NA and are excluded from the model-aware ranking. When no candidate in the elbow neighbourhood clears the floor, the whole call falls back to the geometric ranking and logs a warning that says so. That is the usual outcome well past max_levels = 10: the neighbourhood is the elbow plus or minus max(4, top_n), and on points with no cluster structure the elbow sits near $\sqrt{K_{max}}$, so the neighbourhood reaches ten cells only from about max_levels = 40. Measured on 1000 uniform points, max_levels of 12, 20 and 30 all fell back, and 40 scored k = 10 and 11 alone. On such a layer resolution_profile(), which scores Moran's z at every level of its ladder, is the model-aware view. Under criterion = "combined", an elbow below ten cells is a count Moran's I cannot score, so it cannot weigh the response there: the geometric ranking is returned, with a logged warning, and the diagnostics record it (see Value). (Ranking the window anyway put the smallest count Moran's I scores, ten, first whatever the response did.) Otherwise the candidates below the floor that sit alongside candidates above it all take the last place on the Moran's I axis, after every candidate it scored, while still competing on the geometric axis. That axis is the elbow's own log-log sag at each candidate; when the WSS curve has no elbow it is flat, every candidate tied, and Moran's I alone orders them.

"combined" is not an estimate of the number of clusters. Where the elbow is at ten cells or more, Moran's I can move the pick away from it when the response is still spatially structured at the elbow's resolution. Use "geometric" when the cell count should follow the clustering of the points, and resolution_profile() when the response should drive it.

Value

An integer vector of candidate level counts, best first: under the geometric criterion the elbow, then its lower and upper neighbours (with a warning when the WSS curve has no elbow and the first is the ladder's choice; see Details); under the model-aware criteria the candidates in rank order. k[1] is therefore the top-ranked count on every path, and top_n = 1 returns it alone. When criterion != "geometric", an attribute "diagnostics" is attached with per-k Moran's I values (moran_i) and their standardised deviates (moran_z), the WSS curve (wss) with the relative between-restart spread at each k (wss_spread), the number of rising steps on it (wss_bumps), the restart budget (nstart), the geometric elbow the evaluated neighbourhood was drawn around (knee_k), the k at which k-means failed (failed_k; their wss entries are interpolated from the neighbours, not measured) and the k the model-aware pass actually scored (eval_ks, the elbow's neighbourhood). Under "combined" it also carries the WSS of the re-run clustering at those k (wss_eval), the rank average that ordered them (combined_rank, named by k) and criterion = "combined"; when "combined" returned the geometric ranking because the elbow is below ten cells, it carries criterion = "geometric" and fallback (the reason) in place of combined_rank. When the model-aware path itself falls back to the geometric result (no viable k in the elbow neighbourhood, or Moran's I could not be computed for any candidate), no diagnostics are available and the attribute is absent. Both fallbacks are logged as warnings. The geometric path returns a plain integer vector; a rising WSS curve is still logged there. With select_on = "split" every path adds a "split" attribute: a list with selection and estimation (integer row positions in data_sf), method and seed. For a full per-level table of criteria, cell support, restart spread and the flat region, see resolution_profile().

Post-selection inference

When the selection reads the response (here, whenever both response_var and predictor_vars are supplied), everything estimated afterwards on the chosen cells is estimated on data that already influenced the choice, and its standard errors are post-selection ones: descriptive, not at nominal coverage (Gao, Bien and Witten 2022; Chen and Witten 2023 give the exact selective test for two k-means clusters, the first link of this chain). The exposure is narrower than "the partition was chosen on the response": every partition here is k-means on the coordinates alone, and the response only decides which count is ranked first. It is not zero, because the count determines every cell the downstream standard errors are computed over.

select_on = "split" is sample splitting: the layer is cut into two spatially blocked halves (make_folds(k = 2, method = "block_kfold")), the criteria that read the response (Moran's I on the cell means) read the first half only, and the row positions of both halves come back in the "split" attribute (selection and estimation). The WSS curve and the k-means cells still use every point: they read coordinates alone, and the count is for a tessellation of every point, so it is chosen on that layer's extent and clusters rather than on half of them. Build the tessellation on every point (cells are geometry), but aggregate and fit on data_sf[attr(x, "split")$estimation, ], whose response the selection never saw. That keeps the selection's use of the response out of the estimates, but only as far as the two halves are independent: the split has no buffer, so points near the border between them are correlated with the selection half over the autocorrelation range, and coverage is nominal only when that range is short against the blocks. The price is precision: half the points estimate, and García Rasines and Young (2023) show a contiguous spatial half is less efficient than the exchangeable split the i.i.d. theory assumes, because the two halves are not interchangeable. Two alternatives keep the whole sample: data thinning for count responses (Neufeld et al. 2024) and data fission for Gaussian-like ones (Leiner et al. 2023). Neither is implemented here; the split needs no distributional assumption, which is why it comes first. Selection on coordinates alone ("geometric" with no response) is not exposed in this way, and "split" then changes nothing but the attribute.

The estimation rows cover only the estimation half's blocks, while the cells cover the whole layer, so the estimation rows do not fill every cell. A cell inside the selection half gets none and comes back NA; a cell across the border between the halves is estimated from its estimation-half points alone, and its standard error describes that part, not the cell (on 400 simulated fields a nominal 95% interval covered 0.97 in cells inside the estimation half and 0.88 in cells across the border). Count, per cell, how many of its points are estimation rows, and read inferential results only off cells whose points all are.

References

Arthur, D. and Vassilvitskii, S. (2007). k-means++: the advantages of careful seeding. Proceedings of the 18th Annual ACM-SIAM Symposium on Discrete Algorithms, 1027–1035.

Fränti, P. and Sieranoja, S. (2019). How much can k-means be improved by using better initialization and repeats? Pattern Recognition, 93, 95–112. doi:10.1016/j.patcog.2019.04.014

Steinley, D. (2003). Local optima in K-means clustering: what you don't know may hurt you. Psychological Methods, 8(3), 294–304. doi:10.1037/1082-989X.8.3.294

Gao, L. L., Bien, J. and Witten, D. (2022). Selective inference for hierarchical clustering. Journal of the American Statistical Association, 119, 332–342. doi:10.1080/01621459.2022.2116331

Chen, Y. T. and Witten, D. M. (2023). Selective inference for k-means clustering. Journal of Machine Learning Research, 24(152), 1–41. https://jmlr.org/papers/v24/22-0371.html

García Rasines, D. and Young, G. A. (2023). Splitting strategies for post-selection inference. Biometrika, 110(3), 597–614. doi:10.1093/biomet/asac070

Leiner, J., Duan, B., Wasserman, L. and Ramdas, A. (2023). Data fission: splitting a single data point. Journal of the American Statistical Association, 120(549), 135–146. doi:10.1080/01621459.2023.2270748

Neufeld, A., Dharamshi, A., Gao, L. L. and Witten, D. (2024). Data thinning for convolution-closed distributions. Journal of Machine Learning Research, 25(57), 1–35. https://jmlr.org/papers/v25/23-0446.html

See Also

build_tessellation() and get_voronoi_seeds(), which accept this function's result directly as approx_n_cells and n (the first candidate is used, and the output records that it came from here); assign_features_to_polygons() and summarize_by_cell() for the steps that follow.

Other aggregation: assign_features_to_polygons(), kriging_adequacy(), resolution_profile(), select_resolution(), summarize_by_cell(), summary.resolution_profile()

Examples

library(sf)
set.seed(1)
# Two clearly separated clusters: the elbow should sit near k = 2
pts <- st_as_sf(
  data.frame(x = 5e5 + c(runif(25, 0, 10), runif(25, 90, 100)),
             y = 5e6 + c(runif(25, 0, 10), runif(25, 90, 100))),
  coords = c("x", "y"), crs = 32632
)
determine_optimal_levels(pts, max_levels = 6)   # 2 1 3: the elbow first

Description

Coerces spatial objects to a projected coordinate reference system suitable for distance/area calculations.

Usage

ensure_projected(x, target_crs = NULL, purpose = c("distance", "area"))

Arguments

Details

An object that already has a projected CRS is returned untouched. Only geographic (lon/lat) input is transformed, and the CRS chosen depends on the extent of the data. It is not always UTM:

Local extents

The UTM zone containing the data's centre (EPSG:326xx north of the equator, EPSG:327xx south). Distances and areas are close to true over a few degrees of longitude, which is the case this package is usually in. The centre is the centroid on the sphere, computed with s2 whatever sf::sf_use_s2() is set to, so data near a zone edge get the same zone in every session.

Wide extents

Once the data reach well beyond the roughly 3 degrees a UTM zone is designed for, a single zone can distort distances by several percent, and that error propagates straight into variogram ranges, block sizes, GWR bandwidths and GP length-scales. Which projection is actually best is then measured, not assumed: the zone, a Lambert azimuthal equal-area centred on the data and (where its standard parallels do not degenerate) an Albers conic are each scored by projecting representative points of the data (a non-POINT layer is reduced to one point per feature, plus the vertices of its outline when it has fewer than 40 features, so a single study-area polygon is scored too) and comparing planar with geodesic pairwise distances, and the one that distorts least is used. The choice, both error figures and this argument are logged (see the logging note under spatialkit_quiet()); they are not R warnings, so tryCatch(warning = ) does not see them.

Antimeridian

Data straddling ±180° have a bounding box wider than a hemisphere. The wrap is detected from the coordinates (one very large gap in the sorted longitudes) and an equal-area projection centred on the true extent is used.

Around a pole

Data spanning more than 180 degrees of longitude with no such gap surround a pole. When every point also lies on one side of the equator (Antarctic stations, a pan-Arctic network), a Lambert azimuthal equal-area centred on that pole is used, provided it measures a smaller distance error than the global fallback. Web Mercator splits such a layer at +/-180 degrees and stretches it towards the pole: a ring of Antarctic stations measured a worst-case distance error near 20,000 percent in it, against about 2 percent in the polar projection. Only the coverage left over, spanning both hemispheres or a low-latitude belt the polar projection fits worse, falls back to EPSG:3857 (Equal Earth for purpose = "area").

Missing CRS

With no target_crs, a bounding box that looks like lon/lat means EPSG:4326 is assumed (a real warning) and the rules above then apply; coordinates the heuristic declines are left exactly as they are. With target_crs supplied there is no source CRS to reproject from, so the same heuristic decides between two outcomes: lon/lat-looking coordinates are read as EPSG:4326 and reprojected to the target (a real warning), and anything else has the target stamped on without reprojection. That is a relabel, logged only, so verify the coordinates really are in that CRS. Set the CRS explicitly to suppress either.

target_crs overrides all of this. Pass it whenever you need a specific, reproducible projection: comparing runs, matching an existing layer, or fixing the units that make_folds()'s block_size will be interpreted in.

Value

x, potentially with a new projected CRS. When a projection was chosen here (lon/lat input, no target_crs) the result carries attr(x, "crs_choice"): a data.frame with one row per projection considered, holding name, crs (its definition), the measured worst-case distance_error (relative, over sampled pairs; NA where it could not be measured) and chosen. The figure the log line quotes for the winner is therefore recoverable for every candidate, and is measured for the single candidate on the paths where no comparison runs (a UTM zone on a local extent, the equal-area projection chosen for a layer straddling the antimeridian). It is NULL exactly when no local projection was chosen here: input that already carried a projected CRS, a target_crs you supplied, or an extent no local projection fits, which falls back to Web Mercator or Equal Earth. CRS-less input additionally carries attr(x, "crs_assumed"): "EPSG:4326" when the lon/lat heuristic fired, "none" when it declined. That attribute is also read on the way IN. An object already carrying "none" is returned untouched, with the heuristic skipped, which is how a predict() method replays a fit's negative decision so that a subset of the training rows is not judged differently from the whole.

See Also

Other spatial data preparation: clip_target_for(), coerce_to_points(), harmonize_crs(), prep_model_data()

Examples

library(sf)
pts_ll <- st_as_sf(
  data.frame(lon = c(9.1, 9.2), lat = c(48.7, 48.8)),
  coords = c("lon", "lat"), crs = 4326
)
# A local extent gets the containing UTM zone; the zone's measured
# distance error over the extent travels with the result.
st_crs(ensure_projected(pts_ll))$epsg  # 32632
attr(ensure_projected(pts_ll), "crs_choice")

# A continental extent is scored against the zone and may get an equal-area
# projection instead; the choice and both error figures are LOGGED, not
# warned -- see Details and ?spatialkit_quiet.
wide <- st_as_sf(
  data.frame(lon = c(-120, -70), lat = c(30, 48)),
  coords = c("lon", "lat"), crs = 4326
)
st_crs(ensure_projected(wide))$proj4string

# target_crs overrides the choice entirely.
st_crs(ensure_projected(pts_ll, target_crs = 3035))$epsg  # 3035

# For densities per cell the CRS has to be equal-area.  Two candidates are
# scored, an Albers conic and a Lambert azimuthal, both centred on the data,
# and whichever distorts distance less over the extent is used: here Albers.
st_crs(ensure_projected(pts_ll, purpose = "area"))$proj4string

Description

Ensures that a polygon layer has a reproducible, deterministic identifier column by sorting features using representative point coordinates (and secondary tie-breakers) and then assigning sequential IDs.

Usage

ensure_stable_poly_id(
  polygons_sf,
  id_col = "poly_id",
  method = c("centroid", "surface_point", "bbox_center"),
  make_valid = TRUE,
  transform_for_sort = 4326
)

Arguments

Details

A feature that is valid in its own CRS is not always one s2 can measure. Only the vertices are transformed to the sort CRS, and on the sphere they are joined by great-circle arcs, so a long edge that is straight in the layer's CRS and passes about a metre from another vertex of the same ring can end up on the other side of that vertex. The ring then crosses itself, and s2 refuses it (1 of 291 Voronoi cells of Texas clipped to the state outline in EPSG:5070; "Loop 1 is not valid: Edge 36 crosses edge 52"). The repair of the sort copy does not split crossing edges, so a feature s2 still refuses is handled on its own. Its sort copy is repaired a second time with the crossing edges split. If s2 refuses that as well, or with make_valid = FALSE, which asks for no repair, the feature's centroid and area are measured as plane geometry in the layer's own CRS and the centroid is transformed to the sort CRS. With the methods "surface_point" and "bbox_center" only the area needs this, because those points are not taken with s2.

The function warns once, saying how many features took each route. Their keys are close to the spherical ones and not equal to them. The second repair moved the Texas cell's area by about 40 square metres in 2,949 square kilometres. Its centroid taken in the plane lay 20 m from the one the second repair gave on the sphere, and the two kinds of centroid lay 5 m apart at the median and 61 m at most for the other 290 cells. A feature measured in the plane can therefore sort on the other side of a neighbour whose centre is that close to its own in longitude. Every other feature's key is unaffected, and the geometry returned is never the sort copy.

Value

An sf polygon layer re-ordered with sequential IDs in id_col. Non-polygonal rows are dropped (with a warning), so the result can have fewer rows than the input; if no polygonal rows remain, an error is raised.

See Also

Other tessellation: build_tessellation(), create_grid_polygons(), create_grid_polygons_cached(), create_voronoi_polygons(), get_voronoi_seeds(), plot_tessellation_map(), voronoi_seeds_kmeans(), voronoi_seeds_random()

Examples

library(sf)
bnd <- st_sf(geometry = st_sfc(st_polygon(list(rbind(
  c(0, 0), c(100, 0), c(100, 100), c(0, 100), c(0, 0)
))), crs = 32632))
g <- create_grid_polygons(bnd, target_cells = 9)
# Reverse the rows and re-derive the IDs: the SAME cell gets the same ID,
# which is the property a row-position ID does not have.  The check joins
# the two layers on geometry, because comparing sorted ID vectors would
# pass for any two permutations of 1:9.
fwd <- ensure_stable_poly_id(g)
rev <- ensure_stable_poly_id(g[nrow(g):1, ])
same_cell <- match(st_as_text(st_geometry(rev)), st_as_text(st_geometry(fwd)))
all(rev$poly_id == fwd$poly_id[same_cell])

Description

Fits exponential (or spherical) variogram models and returns the effective range: for the exponential model, three times the fitted range parameter, which is where the semivariance reaches ~95 % of the sill; for the spherical model (fitted only when the exponential fit is singular or does not converge) the fitted range itself, which is where the spherical semivariance reaches its sill exactly. Both are the distance beyond which two observations are (near) uncorrelated, which is what a block or a buffer has to exceed.

Usage

estimate_sac_range(
  points_sf,
  response_var,
  predictor_vars = NULL,
  n_max = 5000L,
  cutoff = 0.5,
  range_frac = 1,
  seed = 123L,
  detrend = c("ols", "reml"),
  reml_max_n = 400L,
  keep_directional_fits = FALSE
)

Arguments

Details

The exponential model is kept whenever it converges, without comparing it with the spherical fit, and on fields smoother than exponential that makes the range long. Measured on simulated fields (n = 300 on a 1000 m square, 30 draws each): about 1.8–2.1 times the practical range of a Gaussian covariance, and 1.3–1.4 times the range of a spherical one, while an exponential field came back at 0.97 of its effective range. The error is on the safe side (blocks too large, cross-validation pessimistic), and it is kept on purpose: choosing the family by the smaller weighted sum of squares corrects the spherical case but sends exponential fields low, to about 0.82 of the truth, which is the direction that leaks.

The estimate is the omnidirectional (all-pairs) fit. Directional variograms are fitted as well, at 0° (N–S), 45°, 90° (E–W) and 135° azimuths with a ±22.5° tolerance. Those four windows tile all 180 distinct azimuths exactly once, and their ranges are returned in the directional attribute, with their largest-over-smallest ratio in anisotropy. They are a diagnostic, not the answer, for two reasons. Each direction sees about a quarter of the point pairs, and the maximum of four quarter-sample fits is biased upward: on simulated isotropic fields it came in about 40% above the truth, and no hurdle placed in front of it (all four directions fitted, ratio above 1.5, maximum above 1.5× the all-pairs fit) kept it out. One isotropic field rotated in 10° steps "established" anisotropy in 14 of 18 orientations. And the windows are fixed to the coordinate axes, so any answer built from them changes when the layer is rotated, which a property of the field must not do. The all-pairs fit is the best-powered estimate available and is invariant to rotation.

Where a field is known to be anisotropic, blocks must be at least as large as the longest autocorrelation range to avoid leakage, and the conservative choice is to size them from the longest directional range explicitly. Read it from directional_fitted, not directional: on a strongly anisotropic field the major axis is the direction most likely to run past the fitted lags, which leaves it NA in directional, so max() of that is NA, or with na.rm = TRUE the second-longest range. Check directional_status first: a major axis marked "over_cutoff" has no identified range at all, and a longer cutoff or make_folds(method = "nndm") is the way on. A ratio above 1.5 is written to the package log at INFO level with that advice, which reaches the session log file but not the console (the ratio passes 1.5 on most isotropic fields too); print() shows the directional ranges and the ratio, and attr(range, "anisotropy") holds it. Only when the omnidirectional fit is singular or did not converge is the directional maximum returned in its place, and anisotropy_used is TRUE in that case alone. An omnidirectional fit that converged to a range past the fitted lags is refused (see the Value section) whatever the directions found: the directions that reached a sill are the shorter ones, so their maximum is a lower bound, not an estimate.

A direction whose fit fails, does not converge, or reports a range beyond the longest fitted lag is excluded and recorded as NA in the directional attribute.

Every variogram model is fitted with a nugget. A nugget-free model forces the curve through the origin, and on any real measurement (which has one) gstat's default N/h² weights buy that constraint by collapsing the range: with a 50% nugget the fitted range came back at about 0.45 of the truth, so make_folds(auto_range = TRUE) built blocks less than half the correlation length it reported.

The lags are binned the way gstat bins them by default, 15 bins out to the cutoff, each cutoff * max_dist / 15 wide (about 47 m on a 1000 m square at the defaults). A range spanning only one or two bins is resolved coarsely and comes out long: exponential fields with an effective range of 60 m (n = 300 on a 1000 m square, 30 draws) returned a median of 89–102 m, where the same fields binned over a 200 m cutoff gave 65–68, and at a range of 300 m there was no bias. A range shorter than the first bin cannot be resolved at all and can come back several times too long: an effective range of 24 m (n = 1500 on a 1000 m square, 8 draws) returned 93–479 m, five of them as the directional maximum, against 19–32 m at cutoff = 0.1; the first bin's semivariance was 92–99 percent of the fitted sill in all eight. So whenever the empirical variogram is already at its sill in the first one or two bins (plot() the result), run it again with a smaller cutoff, whatever range was fitted.

Nothing tests whether the layer has spatial structure at all. On white noise (n = 300 on a 1000 m square, 30 draws) the estimate was a finite, spurious range (57–533 m) in 8 draws and a refusal in the rest, mostly as past the fitted lags or not converged, and only once as no model fitted; with detrend = "reml" it was finite in 13 of 30 (21–453 m), and 16 of the refusals were ranges of 0.18–12.5 m, too short for 30 pairs of points to lie inside them. A spurious range errs towards larger blocks, so the harm is mostly lost training data, but a caller who needs to know whether there is any structure should look at the variogram (plot() on the result) rather than at whether the answer is NA.

When the all-pairs fit is singular or did not converge and two or more directions reached a sill, the directional maximum is returned in its place (anisotropy_used = TRUE), and a log warning names the directional ranges when their ratio exceeds 1.5. When the all-pairs estimate is used and the directional ranges vary by more than 1.5, a log note (INFO) names them instead.

The returned range is in the coordinate units of the (projected) data and can be passed directly to make_folds(block_size = ...) so that CV blocks are at least as wide as the autocorrelation range.

Value

A single number, of class sac_range in the first two of the three shapes below and an unclassed NA in the third; all three behave as an ordinary number. The shapes carry different attributes:

Success

A positive effective range in projected coordinate units, with the fit attached as attributes directional (the 0°, 45°, 90° and 135° ranges, named by azimuth; NA where that direction's fit was unusable), anisotropy (largest over smallest), anisotropy_used (logical: TRUE only when the all-pairs fit was singular or did not converge and the directional maximum stands in for it), directional_status (per azimuth, why a direction is NA in directional: "ok", "over_cutoff" (its range ran past the largest lag fitted), "not_converged" or "no_fit"), directional_fitted (the range each direction's fit reported whether or not it was usable, so a refused directional range stays recoverable) and directional_fits (a list by azimuth of each direction's empirical variogram and fitted model, NULL where there is none, and NULL altogether unless keep_directional_fits = TRUE), detrended (logical: whether the variogram is of the residuals on predictor_vars or of the raw response. A missing predictor is an error, and a failed detrending fit warns and falls back to the raw response with this set to FALSE), detrend_method ("ols" or "reml" when detrended, NA otherwise), reml (with detrend = "reml": a list with n_used, subsampled, nugget_prop and sigma2 from the REML fit; NULL otherwise), crs (the projected CRS the variogram was fitted in: the unit of the range), max_dist, cutoff_dist, variogram (the empirical variogram), variogram_model (the fitted gstat model, or with detrend = "reml" a gstat model built from the REML parameters) and nugget (that model's nugget variance; see sac_nugget), so the fit can be inspected and need not be taken on trust.

Rejected range

NA_real_ when a range was fitted but is not identified: it exceeds range_frac * cutoff * max_dist (see range_frac), which applies to the all-pairs fit even when some directions reached a sill; or too few pairs of points lie inside it to identify it, because it is shorter than the shortest lag the empirical variogram resolves (the mean separation in its first bin) or, with detrend = "reml", than the distance within which 30 pairs of the points the REML fit used lie, when that is shorter (the REML range is fitted to the point pairs, not to the bins). A structure that short cannot be told from a nugget, and the bound is about identification, not a test for spatial structure (see above). Or the model did not converge; or the empirical variogram decreases with distance over its shorter lags (a net fall of more than 15 percent of the mean semivariance there, weighted by pairs), which is the shape of a periodic, hole-effect structure or of a variance that differs between a dense cluster and the rest of the layer. Sampling noise in the short-lag bins of a small sample can make that fall too: on exponential fields it refused 7–9 of 60 draws at n = 30, 3–6 at n = 50 and 0–1 at n = 100 (effective range 300 on a 1000 m square), and 15–16 of 60 at n = 30 with a range of 150. An unremoved trend makes the variogram rise instead; when it rises past the fitted lags the first test catches it, but a milder trend only lengthens the fitted range and passes, which is what predictor_vars is for. Last, the fitted range can be non-positive. It is classed sac_range as well, so it prints as NA without dumping its attributes, and it carries max_dist, cutoff_dist, variogram, variogram_model and nugget (the evidence for the rejection), plus rejected_range (the value that was refused), rejected_reason (one of "fitted range exceeds the largest lag fitted", "fitted range is below the shortest lag fitted", "variogram model did not converge", "empirical variogram decreases with distance", "fitted range is non-positive or non-finite", "no variogram model could be fitted (singular fits)"), crs (so the units the rejected number was in stay recoverable, which is what plot() labels its axis from), detrend_method, reml (as on success) and, for "fitted range is below the shortest lag fitted", range_floor (the distance the refused range fell short of). The reason names one ground for the refusal: a falling variogram first, then a range past the largest lag, and the fit not having converged only when neither holds. A fit that did not converge is therefore most often recorded as past the largest lag, and the count of "variogram model did not converge" understates how many did not. The model itself says: attr(attr(range, "variogram_model"), "converged") is FALSE when no start of either model family converged (gstat stopped at its iteration limit), whatever the reason, and its nugget and sill are then where the optimiser halted, as its range is. A model built from REML parameters carries no such attribute. resolution_profile() reads it. The refused range carries directional, anisotropy, anisotropy_used, directional_status, directional_fitted and, with keep_directional_fits = TRUE, directional_fits as well: the directional sweep runs whatever becomes of the all-pairs fit, and its per-azimuth outcome is what says whether any direction reached a sill the pooled variogram did not, or whether every direction ran past the fitted lags alike. The same shape, with rejected_range = NA, variogram_model = NULL and nugget = NA, is returned when no variogram model could be fitted at all (both the exponential and the spherical fit singular, which a flat, nugget-only variogram can produce, though on white noise it was the outcome in only 1 of 30 draws: see above); rejected_reason says so and the empirical variogram is still attached.

No fit

An unclassed NA_real_ when estimation could not be attempted at all, whose one attribute, rejected_reason, says why: "package 'gstat', which the variogram needs, is not installed", "<n> points, fewer than the 30 a variogram range is estimated from", "<n> point(s) with a finite value to model, fewer than the 30 a variogram range is estimated from", "the response is constant", "the residuals on predictor_vars are constant: the predictors explain the response exactly" or "the points have no extent (the largest distance between them is zero or could not be computed)". Nothing is fitted, so none of the other attributes above exist: no variogram, variogram_model or rejected_range, which is what tells it from a range that was fitted and refused.

Attributes and the class do not affect is.na() or is.finite(), so every downstream guard treats all three the same way it always did.

Detrending and the residual-variogram bias

Fitting a variogram to the residuals of a least-squares trend underestimates both the sill and the range, because the trend fit absorbs part of the long-wavelength spatial variation (Lark, Cullis and Welham 2006). Blocks sized from that range are then too small and a blocked validation is less conservative than it claims. How large the effect is depends on how smooth the trend terms are in space, and on how many there are. Measured for this estimator on simulated exponential fields (n = 300, true effective range 300, nugget 0.2, 40–60 draws), as the median ratio of the estimate from the trend-removed data to the estimate from the true field:

• a white-noise covariate: OLS 1.00, REML 0.99 (no bias to speak of);

• a spatially smooth covariate (a random field with range 300 or 1000): OLS 0.97, REML 0.95–0.98;

• a linear trend in the coordinates: OLS 0.92, REML 1.04;

• a quadratic trend in the coordinates (five terms): OLS 0.75, REML 1.06.

So for ordinary covariates the OLS bias is a few percent, and for trend surfaces in the coordinates it is large. Iterating between a GLS trend fit and a variogram refit (Neuman and Jacobson 1984) recovers only part of it (0.80 in the quadratic case), because the variogram of GLS residuals is biased too; REML does not fit a variogram to residuals at all, which is why detrend = "reml" returns its own range estimate. Its price is a single family (exponential with nugget), a cubic cost in n, and a somewhat wider sampling spread. The default stays "ols" so that existing scripts return what they did; a script that detrends on smooth or coordinate-based terms should pass detrend = "reml", or size its blocks with a margin.

Count and other non-Gaussian responses

The empirical variogram assumes second-order stationarity: a variance that is the same everywhere, so that semivariance depends on separation alone. A count response breaks that assumption by construction, because its variance tracks its mean, and so does any response with a mean-variance relationship (rates, proportions, skewed amounts). On such data the variogram mixes distance-dependence with mean-dependence, and the fitted range is not the thing it claims to be: where the mean is high the semivariance is inflated whatever the distance, which flattens the curve and can lengthen or shorten the apparent range depending on where the high-mean regions sit.

Supplying predictor_vars to detrend helps, because it removes the part of the mean the covariates explain, but it does not fix the variance structure: the residuals of an OLS fit to a count still have a variance that tracks the fitted mean. The principled remedy is a variogram of Pearson residuals from a model in the right family, which this function does not compute. Until it does, treat the range from a count response as no more than an order of magnitude, size blocks conservatively from it, and prefer make_folds(method = "nndm"), which does not depend on a fitted range at all.

References

Lark, R. M., Cullis, B. R. and Welham, S. J. (2006). On spatial prediction of soil properties in the presence of a spatial trend: the empirical best linear unbiased predictor (E-BLUP) with REML. European Journal of Soil Science, 57(6), 787–799. doi:10.1111/j.1365-2389.2005.00768.x

Neuman, S. P. and Jacobson, E. A. (1984). Analysis of nonintrinsic spatial variability by residual kriging with application to regional groundwater levels. Mathematical Geology, 16(5), 499–521. doi:10.1007/BF01886329

See Also

sac_nugget for the nugget behind the estimate, plot.sac_range to see the variogram the estimate rests on.

Other cross-validation: area_of_applicability(), cv_bayes(), cv_block_size_sweep(), cv_gwr(), cv_rf(), cv_spatial(), fold_separation(), gwr_model_selection(), make_folds(), sac_nugget(), select_features_forward()

Examples

if (requireNamespace("gstat", quietly = TRUE)) {
  library(sf)
  # A Gaussian random field with an exponential covariance: range
  # parameter 100, so the true effective range is 3 x 100 = 300, plus a
  # small nugget.
  set.seed(9)
  n <- 150
  xy <- data.frame(x = 5e5 + runif(n, 0, 1000), y = 5e6 + runif(n, 0, 1000))
  D  <- as.matrix(dist(xy))
  xy$z <- as.numeric(t(chol(exp(-D / 100) + diag(0.1, n))) %*% rnorm(n))
  pts <- st_as_sf(xy, coords = c("x", "y"), crs = 32632)
  r <- estimate_sac_range(pts, response_var = "z")
  # print() throughout, because only the last value of a braced block is
  # shown on its own, and every one of these is worth reading.
  print(r)                              # the effective range, in metres
  print(attr(r, "directional"))         # the four directional ranges
  print(attr(r, "variogram_model"))     # the fitted gstat model behind it

  # A field whose range the data cannot pin down: the variogram never
  # reaches a sill within the lags fitted, so the answer is NA with the
  # refused value attached rather than a long range asserted.
  xy$trend <- sin(xy$x / 400) + rnorm(n, sd = 0.2)
  r2 <- estimate_sac_range(st_as_sf(xy, coords = c("x", "y"), crs = 32632),
                           response_var = "trend")
  print(r2)
  attr(r2, "rejected_range")
}

Description

Accepts a single spatial_fit object or a named list of them. Does NOT refit. Uses fitted() for in-sample and predict() for new data.

Usage

evaluate_insample(fits, newdata = NULL, ...)

Arguments

Value

A data.frame with one row per model and columns for model name, all regression metrics, and metric_basis: what the row's metrics were computed on, "in-sample" (fitted values), "out-of-bag" (an rf_fit's fitted values, see "What the metrics are computed on") or "newdata". Rows with different bases do not compare like for like. An element that is not a spatial_fit is skipped, with a logged warning, and has no row; a list in which no element is a spatial_fit is an error.

What the metrics are computed on

With newdata = NULL the metrics come from fitted(object). That is in-sample for a gwr_fit or a bayesian_fit, but out-of-bag for an rf_fit, whose fitted() method returns out-of-bag predictions (see fit_rf_model). The data.frame model_metrics() returns carries no label distinguishing the two, so check object$info$fitted_are_oob before comparing numbers across backends; evaluate_insample() and compare_models() record it per model in a metric_basis column. compare_models_cv scores every backend the same way.

$R^2$ is $1 - RSS/TSS$ with the total sum of squares taken about the mean of the response the model was fitted to. In sample that is the ordinary $R^2$. With newdata it is out-of-sample $R^2$, the convention every cv_*() function uses: the model is measured against the prediction it had to beat, the training mean, not against the new rows' own mean, which it could not have known. It is below 0 when the model predicts the new rows worse than the training mean does, and it is NA when the response does not vary about that baseline by more than rounding error (100 machine epsilons of its magnitude, whatever its units).

Percentage errors on responses with zeros

MAPE divides by the observed value and SMAPE by $|y| + |\hat{y}|$, so neither is defined where its denominator is zero. Neither returns Inf or NaN. Both are averaged over the rows whose denominator is non-zero, and are NA when no row qualifies. Non-zero is judged at the scale of the data: a denominator no larger than 100 machine epsilons times the largest one counts as zero, so the rule does not depend on the units of the response. The n_MAPE and n_SMAPE columns record how many rows that was; the n column counts finite observation/prediction pairs. Read a percentage error next to its count: when n_MAPE < n, MAPE is an average over a subset of the data, whatever its value.

This bites on any response taking exact zeros: counts, rainfall, abundance, claim amounts. On a zero-inflated response with 62 zeros out of 120, MAPE is an average over the 58 non-zero rows, which n_MAPE = 58 now says. SMAPE fails differently and more subtly: it drops the rows where observation and prediction are both near zero (which on a well-fitted zero-inflated model are the rows it got right), so it averages the harder rows only and reads worse than the fit deserves; n_SMAPE shows how many rows it kept, and the count is only a label, not a repair.

RMSE, MAE and $R^2$ use every finite row and are unaffected; prefer them whenever the response can be zero. For a Bayesian fit, cv_bayes() additionally reports CRPS and interval coverage, which are proper scoring rules and have no such failure mode.

See Also

Other model evaluation: compare_models(), compare_models_cv(), model_metrics(), residual_morans_i()

Examples

if (requireNamespace("ranger", quietly = TRUE)) {
  library(sf)
  set.seed(1)
  pts <- st_as_sf(
    data.frame(x = 5e5 + runif(60, 0, 1000), y = 5e6 + runif(60, 0, 1000),
               a = rnorm(60)),
    coords = c("x", "y"), crs = 32632
  )
  pts$z <- 2 * pts$a + rnorm(60, 0, 0.3)
  # Fit on 40 rows and keep 20 back, so the second call really is out of
  # sample; scoring the training rows again would only re-read the fit.
  fit <- fit_rf_model(pts[1:40, ], "z", "a", num_trees = 50, seed = 1)
  print(evaluate_insample(fit))                    # in-sample (out-of-bag for RF)
  evaluate_insample(fit, newdata = pts[41:60, ])   # on the 20 held-out rows
}

Description

Fits a regression whose residual spatial structure is modelled explicitly, as a Gaussian process over the coordinates, and so kept out of the errors. Two things follow, and they are the reasons to reach for this backend. First, every quantity comes with a posterior, so predictions carry calibrated intervals instead of point estimates. Score them with cv_bayes(), which reports held-out interval coverage and CRPS. Second, the fitted length-scale is itself an estimate of how far the spatial dependence reaches, a number you can read and report.

Usage

fit_bayesian_spatial_model(
  data_sf,
  response_var,
  predictor_vars,
  family = NULL,
  gp_k = NULL,
  gp_c = NULL,
  gp_iso = FALSE,
  prior = NULL,
  chains = 4,
  iter = 2000,
  warmup = floor(iter/2),
  cores = getOption("mc.cores", 1L),
  seed = 123,
  backend = c("auto", "cmdstanr", "rstan"),
  control = list(),
  compute_loo = TRUE,
  standardize_predictors = FALSE,
  check_convergence = TRUE,
  pointize = "auto",
  boundary = NULL,
  .already_prepped = FALSE
)

Arguments

Details

Choose it over fit_gwr_model() when you want one global relationship plus an explicit spatial random field, and uncertainty you can defend; choose GWR instead when the question is how a coefficient varies across the map. Choose fit_rf_model() when predictive accuracy matters more than an interpretable model and the response is non-linear in the predictors. The cost here is time: this is full MCMC via 'brms' and Stan, so it is minutes rather than seconds, and the GP is fitted through a reduced-rank basis approximation whose size (gp_k) trades fidelity against runtime.

GP basis count and boundary factor. brms::gp() builds a full tensor grid over its covariates, so a term gp(..x, ..y, k = gp_k) carries gp_k^2 basis functions: the gp_k argument is the count per dimension, not the total rank. Both gp_k and gp_c are therefore chosen from the ratio of the estimated length-scale to the domain extent, following Riutort-Mayol et al. (2023): gp_c is set large enough to contain the upper length-scale bound, and gp_k large enough to resolve the lower one. The derived value is typically 21-25 per dimension and is largely independent of n.

The domain extent used is the one brms::gp(c = ) itself multiplies: the full pooled range of the column-centred coordinates (brms:::choose_L(), taken over the unique coordinate rows, because brms:::.data_gp() reduces the covariates to unique rows first under the default gr = TRUE, so repeat visits to one location do not widen the domain), not the per-axis half-range in which Riutort-Mayol et al. state their inequalities. Both constraints are really constraints on the boundary $L = c \times S$, so expressing them in brms's units is what keeps gp_c, gp_k and $info$gp_ell_min describing the basis brms actually builds. A gp_c derived on the half-range convention and handed to brms::gp() produces a boundary twice as wide as intended, against which gp_k under-resolves by a factor of two.

The GP term is built with scale = FALSE. brms::gp() otherwise rescales its covariates so the maximum Euclidean distance between two points is 1, and reports lscale in that space; since this function already standardises the coordinates, and the length-scale prior, gp_c and the adequacy check below are all expressed in those standardised units, a second normalisation would leave every length-scale quantity in the wrong units.

After fitting, the posterior length-scale is compared against the smallest scale the chosen basis can resolve (1.75 * gp_c * S / gp_k, stored as $info$gp_ell_min); when more than 10% of the posterior mass falls below it a WARN line is logged and the share is recorded as $info$convergence_diagnostics$gp_lscale_below_resolution, which is the signal that gp_k should be raised. This runs with the other checks, so only under check_convergence = TRUE.

Coordinate scaling and anisotropy. Before fitting the GP, X and Y coordinates are each centred and divided by their own standard deviation. This is a conditioning step: easting and northing frequently span very different ranges in a projected CRS, and handing Stan raw metres samples poorly.

Because the axes are scaled independently, a single shared length-scale in the scaled space corresponds to an anisotropic kernel in the original CRS, stretched by whatever ratio sd(X)/sd(Y) happens to take. That ratio is a property of how the sampling locations are laid out, not of the process being modelled, so it is not a defensible source of anisotropy.

gp_iso = FALSE (the default) therefore fits one length-scale per axis, letting the model estimate directional structure from the data instead of inheriting it from the standardisation. Set gp_iso = TRUE to recover the previous single-length-scale behaviour.

gp_iso does not affect cost: brms::gp() builds a tensor grid either way, so the model carries gp_k^2 basis functions regardless. The stored $info$coord_scaling list records the scaling strategy, and $info$gp_iso records which kernel was used.

Value

A bayesian_fit object (inherits from spatial_fit). Supports predict(), fitted(), residuals(), coef(), summary(), and model_metrics(). Model-specific metadata lives in $info (coords: the names of the scaled coordinate columns handed to brms::gp(); coord_scaling, predictor_scaling, gp_k, gp_c, gp_iso, gp_n_basis, gp_ell_min, gp_S: the pooled centred range brms::gp(c = ) multiplies; gp_cmeans: the column means brms centred the scaled coordinates on; gp_xy_range: the training extrema of the scaled coordinates, which predict() uses, with gp_S and gp_cmeans, to hold the GP boundary at its fitted value; gp_lengthscale_bounds: the c(lower, upper) the length-scale prior was calibrated over; gp_lscale_prior: the length-scale prior brms::validate_prior() reports the model will actually use, which is not necessarily the one this function requested (several entries, semicolon-separated, if brms resolved the axes differently); loo, looic, convergence_ok (TRUE or FALSE, and NA when nothing was checked, as under check_convergence = FALSE), convergence_diagnostics: n_divergent, max_rhat, min_neff_ratio, and rhat_failed / neff_failed, the parameters that failed each check by name with their values (empty when none failed), which is what makes a failed check actionable; and n_dropped: the rows prep_model_data() removed for missing or non-finite values or a bad geometry, so $n can be read against nrow(data_sf)). The raw brmsfit is in $engine.