--- title: "spatialkit: Tessellations, Spatial Cross-Validation and Models" description: > The whole pipeline end to end on North Carolina: four tessellations, aggregation, folds, models, and a surface with its area of applicability. author: "Justin Chase" date: "`r format(Sys.Date(), '%d %B %Y')`" output: rmarkdown::html_vignette: toc: true fig_width: 6 fig_height: 4 vignette: > %\VignetteIndexEntry{spatialkit: Tessellations, Spatial Cross-Validation and Models} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include=FALSE} # ggplot2 is a Suggests package, and every chunk below either draws something # or feeds something that does. R CMD build runs this file, so the whole # vignette is gated on ggplot2 being installed rather than calling library() # unguarded and failing the build on a machine without it. has_ggplot <- requireNamespace("ggplot2", quietly = TRUE) has_geom <- requireNamespace("geometry", quietly = TRUE) has_ranger <- requireNamespace("ranger", quietly = TRUE) has_gstat <- requireNamespace("gstat", quietly = TRUE) has_gwmodel <- requireNamespace("GWmodel", quietly = TRUE) && requireNamespace("sp", quietly = TRUE) knitr::opts_chunk$set( echo = TRUE, message = FALSE, warning = FALSE, fig.align = "center", out.width = "100%", eval = has_ggplot ) ``` ```{r no-ggplot, echo=FALSE, eval=!has_ggplot, results='asis'} cat("**ggplot2 is not installed, so the code below is shown but not run.**", "Install it with `install.packages(\"ggplot2\")` and rebuild.\n") ``` ## Overview This vignette generates **synthetic spatial data** over North Carolina and walks through the `spatialkit` workflow end to end: 1. Load the boundary and check which optional packages are present 2. Simulate observations with known structure over it 3. Build four tessellation types (Voronoi, hex, square, Delaunay) 4. Assign points to cells, aggregate with `summarize_by_cell()`, and map the cell means 5. Estimate the autocorrelation range, build spatial CV folds, and score a model under **blocked** versus **random** folds 6. Fit the models to keep: a random forest, GWR when it is installed, and the two compared on identical folds 7. Predict onto a surface and mark where that surface is extrapolation Everything is self-contained — the boundary comes from the `nc.shp` demo shapefile bundled with `sf`, so no external files are needed. Optional backends are checked at the top and every section that needs one is guarded, so this document renders with whatever you happen to have installed. This table always runs, even when the gate below it is shut -- it is what says which backend is missing when the code stops producing output: ```{r backend-report, echo=FALSE, eval=TRUE} data.frame( package = c("ggplot2", "patchwork", "geometry", "ranger", "gstat", "GWmodel + sp"), available = c(has_ggplot, requireNamespace("patchwork", quietly = TRUE), has_geom, has_ranger, has_gstat, has_gwmodel), used_for = c("all plots", "the stacked comparison panels", "true Delaunay triangulation", "random forest backend", "variogram / autocorrelation range", "GWR backend") ) ``` --- ## 1. Packages and boundary ```{r load-packages} library(spatialkit) library(sf) library(dplyr) library(ggplot2) set.seed(42) ``` We load the North Carolina county boundaries shipped with `sf`, dissolve them into a single state outline, and project to **NAD83 / North Carolina (ftUS)** (EPSG:2264) so distances are planar rather than angular. Every distance, bandwidth and block size below is therefore in **US survey feet**, not metres — projected does not mean metric, and the unit is whatever the CRS says it is: ```{r boundary} nc_counties <- st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE) nc_boundary <- nc_counties |> st_union() |> st_transform(2264) |> st_as_sf() ``` --- ## 2. Synthetic observations 300 points inside the state boundary with two predictors and a spatially varying response: - `elevation` — gradient increasing west to east, plus noise - `pop_density` — decays with distance from two fake "cities" - `y` — driven by the two predictors **plus a spatial field neither of them explains**. That last term is deliberate: it is the unmodelled spatial structure that makes random cross-validation optimistic, and it is what section 5 measures. ```{r fake-data} n_points <- 300 pts_raw <- st_sample(nc_boundary, size = n_points, type = "random") pts_coords <- st_coordinates(pts_raw) x_coords <- pts_coords[, 1] y_coords <- pts_coords[, 2] elevation <- scale(x_coords)[, 1] * 500 + rnorm(n_points, 3000, 400) city1 <- c(1530000, 550000) # Charlotte-ish in EPSG:2264 city2 <- c(2150000, 750000) # Raleigh-ish dist_to_city <- pmin( sqrt((x_coords - city1[1])^2 + (y_coords - city1[2])^2), sqrt((x_coords - city2[1])^2 + (y_coords - city2[2])^2) ) pop_density <- pmax(exp(-dist_to_city / 400000) * 5000 + rnorm(n_points, 200, 100), 10) spatial_field <- 20 * sin(x_coords / 250000) * cos(y_coords / 250000) y_response <- 50 + 0.01 * elevation + 0.005 * pop_density + spatial_field + rnorm(n_points, 0, 5) points_sf <- st_sf( y = y_response, elevation = elevation, pop_density = pop_density, geometry = pts_raw ) ``` ```{r quick-peek, fig.height=3.6, fig.alt="Three hundred simulated observation points scattered across North Carolina, coloured by response value. Lower values cluster in the west and higher values along the eastern coast."} ggplot() + geom_sf(data = nc_boundary, fill = "grey95", colour = "black") + geom_sf(data = points_sf, aes(colour = y), size = 1.1) + scale_colour_viridis_c(name = "Response (y)") + theme_void() + ggtitle("Raw observation points, North Carolina") ``` --- ## 3. Four tessellations ### 3a. Voronoi from ~40 k-means seeds Voronoi cells are built around *seed* points, not the observations themselves. Seeding one cell per observation would give 300 cells each containing a single point — a nearest-neighbour interpolation rather than an aggregation, with no within-cell variation to compute a standard error from. `get_voronoi_seeds()` clusters the observations first, so cell size follows sampling density. ```{r tess-voronoi} seeds <- get_voronoi_seeds( boundary = nc_boundary, sample_points = points_sf, method = "kmeans", n = 40 ) tess_voronoi <- build_tessellation( seeds, boundary = nc_boundary, method = "voronoi", clip = TRUE, quiet = TRUE ) ``` ### 3b and 3c. Hex and square grids (~50 cells) ```{r tess-grids} tess_hex <- build_tessellation( points_sf, boundary = nc_boundary, method = "hex", approx_n_cells = 50, clip = TRUE, quiet = TRUE ) tess_square <- build_tessellation( points_sf, boundary = nc_boundary, method = "square", approx_n_cells = 50, clip = TRUE, quiet = TRUE ) ``` All four methods return cells carrying a `cell_id` column; the grid methods keep `poly_id` alongside it, holding the same values. ```{r cell-id-check} names(tess_hex$cells) ``` ### 3d. Delaunay triangles `method = "triangles"` uses `geometry::delaunayn()` when `geometry` is installed. Without it, it does **not** error — it falls back to `sf::st_triangulate()` (GEOS) on the point set and logs a warning. That is still the Delaunay triangulation of the input points; only the resolution of degenerate configurations can differ. Because the two paths are not identical, this section is guarded on `geometry` rather than on the call succeeding: ```{r tess-tri, eval=has_ggplot && has_geom} tess_tri <- build_tessellation( points_sf, boundary = nc_boundary, method = "triangles", clip = TRUE, quiet = TRUE ) nrow(tess_tri$cells) ``` ```{r cell-counts} cat(sprintf( "Voronoi: %d | Hex: %d | Square: %d | Triangles: %s\n", nrow(tess_voronoi$cells), nrow(tess_hex$cells), nrow(tess_square$cells), if (has_geom) nrow(tess_tri$cells) else "skipped (install 'geometry')" )) # has_geom is requireNamespace("geometry", quietly = TRUE), set in the # hidden setup chunk with the other availability flags the table above shows. ``` --- ## 4. Cell-level aggregation with `summarize_by_cell()` `assign_features_to_polygons()` joins points to cells; `summarize_by_cell()` does the aggregation — means, standard deviations and standard errors for the response and every predictor, plus `n` and `cell_weight`. There is no need to write the `group_by()`/`summarise()` by hand, and doing so loses the design-effect machinery below. ```{r summarize} assigned <- assign_features_to_polygons(points_sf, tess_voronoi$cells, polygon_id_col = "cell_id") cell_stats <- summarize_by_cell( assigned, response_var = "y", predictor_vars = c("elevation", "pop_density"), id_col = "cell_id" ) head(as.data.frame(cell_stats)[, c("cell_id", "n", "resp_mean_y", "..sd_resp_y", "..se_resp_y")]) ``` The `..se_*` columns are **IID** standard errors at the default `deff = 1`. That is the right standard error for each cell's own mean, the value a choropleth shows, when the points are spread through the cell. As a standard error for the population (grand) mean it is anticonservative when points inside a cell are spatially correlated. For that population-level use, `deff = "kish"` applies Kish's design-effect correction from an estimated intra-class correlation, and `attr(., "deff_applied")` records what was used: ```{r summarize-deff} cell_kish <- summarize_by_cell( assigned, response_var = "y", predictor_vars = c("elevation", "pop_density"), id_col = "cell_id", deff = "kish" ) deff <- attr(cell_kish, "deff_applied") cat(sprintf("method = %s | ICC(response) = %.3f | median design effect = %.2f\n", deff$method, deff$icc_resp, stats::median(deff$deff, na.rm = TRUE))) # Inflation of the response standard error, cell by cell. summary(cell_kish$`..se_resp_y` / cell_stats$`..se_resp_y`) ``` ### Choropleths ```{r choropleth-helper} #' Assign points to cells, summarise, and draw a choropleth. make_choropleth <- function(tess, boundary, points, title = NULL) { cells <- tess$cells assigned <- assign_features_to_polygons(points, cells, polygon_id_col = "cell_id") stats_df <- summarize_by_cell(assigned, response_var = "y", id_col = "cell_id") cells <- left_join(cells, as.data.frame(stats_df)[, c("cell_id", "resp_mean_y")], by = "cell_id") plot_tessellation_map( tessellation_sf = cells, boundary = boundary, fill_col = "resp_mean_y", palette = "viridis", tile_alpha = 0.9, outline_col = "white", outline_size = 0.3, boundary_col = "grey20", boundary_size = 0.8, legend_title = "Mean y", title = title, subtitle = sprintf("%d cells | %d observations", nrow(cells), nrow(points)) ) } ``` ```{r choro-voronoi, fig.cap="Voronoi choropleth", fig.alt="The same response aggregated into forty Voronoi cells over North Carolina. Cell sizes follow the density of the observations, and the west-to-east gradient from low to high values is visible."} make_choropleth(tess_voronoi, nc_boundary, points_sf, "Voronoi tessellation") ``` ```{r choro-hex, fig.cap="Hex grid choropleth", fig.alt="The same response aggregated into a regular hexagonal grid of about fifty equal-area cells, showing the same west-to-east gradient on a fixed geometry rather than one drawn from the data."} make_choropleth(tess_hex, nc_boundary, points_sf, "Hexagonal grid") ``` ```{r choro-square, fig.cap="Square grid choropleth", fig.alt="The same response aggregated into a regular square grid of about forty-five cells, showing the same gradient with the axis-aligned artefacts a square grid introduces."} make_choropleth(tess_square, nc_boundary, points_sf, "Square grid") ``` ```{r choro-tri, fig.cap="Delaunay choropleth", eval=has_ggplot && has_geom, fig.alt="The same response aggregated onto the Delaunay triangulation of the observation points, which produces many small irregular triangles rather than compact cells."} make_choropleth(tess_tri, nc_boundary, points_sf, "Delaunay triangulation") ``` The three stacked on one shared colour scale. Each map drawn on its own gets its own scale, so a colour in one means a different value in the next; here the limits are set from all three together and `patchwork` then collects the legends into one, which is what makes reading down the page legitimate. ```{r comparison-panel, fig.width=7, fig.height=5.4, fig.cap="All three at a glance", eval=has_ggplot && requireNamespace("patchwork", quietly = TRUE), fig.alt="The Voronoi, hexagonal and square maps stacked on one shared colour scale, so the cell values can be compared down the page. All three show the same west-to-east gradient; they differ in how the cell boundaries fall."} library(patchwork) # The cell means of all three, so the colour scale can span them. cell_means <- function(tess) { a <- assign_features_to_polygons(points_sf, tess$cells, polygon_id_col = "cell_id") summarize_by_cell(a, response_var = "y", id_col = "cell_id")$resp_mean_y } lims <- range(unlist(lapply(list(tess_voronoi, tess_hex, tess_square), cell_means)), na.rm = TRUE) # One line per map instead of a title and a subtitle, and the cell count taken # from the tessellation rather than typed in. bare <- function(tess, label) { make_choropleth(tess, nc_boundary, points_sf) + ggplot2::scale_fill_viridis_c(name = "Mean y", limits = lims) + ggplot2::labs(title = sprintf("%s, %d cells", label, nrow(tess$cells)), subtitle = NULL) + ggplot2::theme(plot.title = ggplot2::element_text(size = ggplot2::rel(1))) } (bare(tess_voronoi, "Voronoi") / bare(tess_hex, "Hex") / bare(tess_square, "Square")) + plot_layout(guides = "collect") + plot_annotation(title = "Tessellation comparison, cell-level mean response") ``` --- ## 5. Spatial cross-validation This is the part the rest of the package exists to support. ### 5a. How far does correlation reach? `estimate_sac_range()` fits an omnidirectional variogram and returns the effective range in CRS units — here, US survey feet. It also fits the four principal directions and reports their ranges in the `directional` attribute (with their largest-over-smallest ratio in `anisotropy`), but only as a diagnostic: each direction sees about a quarter of the point pairs, the maximum of four such fits is biased upward, and the windows are fixed to the coordinate axes, so nothing built from them is invariant to rotating the layer. The all-pairs fit is the estimate whenever it is usable; the directional maximum stands in for it only when the all-pairs fit fails, and the `anisotropy_used` attribute is `TRUE` in that case alone. If you *know* the field is anisotropic, size blocks from the longest directional range explicitly: `max(attr(sac, "directional_fitted"), na.rm = TRUE)`, after checking `attr(sac, "directional_status")`, since `directional` is `NA` for a direction that ran past the fitted lags, which on such a field is often the long one. It returns `NA` (still classed `sac_range`, so it prints as `NA`) when the empirical variogram never reaches a sill, because an unidentified range must not be used to size blocks. ```{r sac-range, eval=has_ggplot && has_gstat} sac <- estimate_sac_range(points_sf, response_var = "y", predictor_vars = c("elevation", "pop_density")) sac if (is.na(sac)) { cat("no identified range:", attr(sac, "rejected_reason"), "\n") } else { cat(sprintf("range = %.0f ft; anisotropy ratio %.2f; used %s\n", as.numeric(sac), attr(sac, "anisotropy"), attr(sac, "anisotropy_used"))) } ``` The variogram is attached either way, so `plot(sac)` draws the curve and lets you judge it rather than trust it (`plot(fit, type = "variogram")` does the same for a fitted model's residuals). The distance axis is labelled in the units of the CRS the variogram was fitted in, and when no range was identified the subtitle says why. A rejected range must not be handed to `make_folds(auto_range = TRUE)`, which is why it comes back `NA` rather than as a long range. ### 5b. Two fold schemes on the same data ```{r folds} # No range was identified in 5a, so the block size cannot come from the # variogram. The field was built as a sine wave with a quarter-wavelength of # about 390,000 ft, the scale over which it changes sign; 400,000-ft blocks # hold out whole lobes of it. cv_block_size_sweep() is the tool for sizing # blocks by measurement when nothing is known in advance. folds_random <- make_folds(points_sf, k = 5, method = "random_kfold", seed = 42) folds_blocked <- make_folds(points_sf, k = 5, method = "block_kfold", block_size = 400000, seed = 42) c(random = folds_random$k, blocked = folds_blocked$k, blocks = folds_blocked$params$blocks_used) ``` `plot_folds()` is the fastest way to see whether the blocks actually separate the data or are smaller than the autocorrelation range and therefore leaking: ```{r plot-folds, fig.width=7, fig.height=6, eval=has_ggplot && requireNamespace("patchwork", quietly = TRUE), fig.alt="Two maps of the same points, random folds above and blocked folds below, sharing one fold legend. In the random map the five colours are interleaved everywhere; in the blocked map, with the block grid drawn behind, every block holds a single fold and each fold is two to four whole blocks: four of the five folds fall in two separate parts of the state, and the fifth's two blocks meet only at a corner."} library(patchwork) (plot_folds(folds_random, points_sf, boundary = nc_boundary) / plot_folds(folds_blocked, points_sf, boundary = nc_boundary)) + plot_layout(guides = "collect") ``` ### 5c. The number the fold scheme changes `fit_rf_model()` and `cv_rf()` need `ranger`. `include_coords = TRUE` hands the forest the coordinates, which lets it memorise the training surface — the failure mode the default (`FALSE`) exists to prevent, and the one random folds cannot see. Scoring the same model on both fold schemes shows the size of the gap: ```{r cv-contrast, eval=has_ggplot && has_ranger} cv_random <- cv_rf(points_sf, "y", c("elevation", "pop_density"), folds = folds_random, include_coords = TRUE, num_trees = 300, seed = 1) cv_blocked <- cv_rf(points_sf, "y", c("elevation", "pop_density"), folds = folds_blocked, include_coords = TRUE, num_trees = 300, seed = 1) data.frame( folds = c("random_kfold", "block_kfold"), R2 = c(cv_random$overall$R2, cv_blocked$overall$R2), RMSE = c(cv_random$overall$RMSE, cv_blocked$overall$RMSE) ) ``` ```{r cv-contrast-text, echo=FALSE, results='asis', eval=has_ggplot && has_ranger} cat(sprintf( "Random folds report R2 = %.3f; blocked folds report %.3f on the same fitted model, a drop of %.0f%% of the reported skill.\n", cv_random$overall$R2, cv_blocked$overall$R2, 100 * (cv_random$overall$R2 - cv_blocked$overall$R2) / cv_random$overall$R2)) ``` The blocked estimate is the one to report. The random one describes interpolation between points you already have, which is not the task. --- ## 6. Fitting a model ```{r rf-fit, eval=has_ggplot && has_ranger} rf_fit <- fit_rf_model(points_sf, response_var = "y", predictor_vars = c("elevation", "pop_density")) rf_fit ``` Note that `summary()` on an `rf_fit` reports **out-of-bag** metrics, not in-sample ones — `fitted.rf_fit()` returns out-of-bag predictions, and the printout says so. A `gwr_fit` or `bayesian_fit` reports genuinely in-sample metrics, so the two are not comparable; `compare_models_cv()` is. ```{r rf-summary, eval=has_ggplot && has_ranger} summary(rf_fit) rf_fit$info$importance # coef() on an rf_fit errors: a forest has no coefficients ``` The two printouts disagree in the third decimal of R² — `0.4733` above, `0.4715` here — while reporting an identical RMSE. Neither is wrong, and neither is a bug. Both read the *same* out-of-bag predictions; they differ only in the denominator of the variance they compare against. `print.rf_fit()` echoes `ranger`'s own `r.squared`, which is `1 - MSE_oob / var(y)` using the unbiased (n − 1) sample variance. `summary()` recomputes `1 - SS_res / SS_tot` from the predictions, where `SS_tot = sum((y - mean(y))^2)` — an n denominator. The unexplained fraction therefore differs by exactly the factor n / (n − 1), here 300/299, and RMSE, which involves no such comparison, matches to four decimals. If you need a figure comparable across backends, use `cv_rf()` rather than either. ### Residual diagnostics `plot()` on any `spatial_fit` maps the residuals; visible structure means unmodelled spatial autocorrelation. `residual_morans_i()` puts a number on it. ```{r rf-resid, fig.height=3.6, eval=has_ggplot && has_ranger, fig.alt="Residuals of the random forest mapped over North Carolina on a red-to-blue scale. Positive residuals cluster in the west and negative ones through the centre, which is the visible spatial structure the model has not captured."} plot(rf_fit, type = "residuals") ``` ```{r rf-moran, eval=has_ggplot && has_ranger} mi <- residual_morans_i(rf_fit) cat(sprintf("Moran's I = %.4f (z = %.2f, p = %.3g)\n", mi$observed, mi$z, mi$p_value)) ``` ### GWR, if `GWmodel` is installed ```{r gwr-fit, eval=has_ggplot && has_gwmodel} gwr_fit <- fit_gwr_model( data_sf = points_sf, response_var = "y", predictor_vars = c("elevation", "pop_density"), adaptive = TRUE, kernel = "bisquare" ) gwr_met <- model_metrics(gwr_fit) # With adaptive = TRUE the bandwidth is a number of nearest neighbours, not a # distance in feet: each local regression uses the same count of points. cat(sprintf("Bandwidth: %.0f neighbours | in-sample R2: %.3f | RMSE: %.3f\n", gwr_fit$info$bandwidth, gwr_met$R2, gwr_met$RMSE)) ``` The fit raises no collinearity warning, but its survey of the local designs (`gwr_fit$info$local_collinearity`, see "Collinearity diagnostics" in `?fit_gwr_model`) is worth a look: in 25 of the 300 windows the condition index that includes the intercept (`cn`) is above 30. The cause is `elevation`: it averages about 3,000 but has a standard deviation of only about 400 inside a 42-neighbour window, so locally it is nearly collinear with the intercept. What that leaves poorly determined is the local intercepts, not the slopes: the slope index (`cn_slopes`), which centres the predictors in each window, stays below 30 everywhere, so no slope is masked and the fit does not warn. `plot(gwr_fit, type = "coefficients", term = "Intercept")` draws those 25 intercepts hollow. Centring the predictor (`elevation - mean(elevation)`) gives the same bandwidth and fitted values and determines the intercepts too. ### Comparing backends on identical folds `compare_models_cv()` cross-validates each requested backend on the `folds` you pass. Any backend whose package is missing is dropped with a message, so the call still returns whatever could run. ```{r compare, eval=has_ggplot && has_ranger} cmp <- compare_models_cv( points_sf, "y", c("elevation", "pop_density"), models = c("RF", "GWR"), folds = folds_blocked, rf_args = list(num_trees = 300) ) cmp$overall ``` The RF row is the forest *without* the coordinates, and on the same blocked folds it scores worse than the coordinate forest in 5c did. That is not a contradiction. The field is a smooth wave that the two predictors do not carry, so the coordinates are informative out of block as well as in, and the blocked score says how much. What 5c measured was how much of the *random*-fold score was interpolation between neighbours; what this table measures is what the predictors alone are worth. `Adj_R2` is `NA` in every `cv_*()` result's `overall`: the pooled predictions come from `k` separately fitted models and have no single parameter count to adjust for. --- ## 7. From a fit to a map, and where the map applies `predict_surface()` builds a regular grid over the training extent, joins covariates from the nearest observation, clips to the boundary and predicts: ```{r surface, fig.height=3.6, eval=has_ggplot && has_ranger, fig.alt="The fitted random forest predicted onto a regular grid over North Carolina and clipped to the state boundary, with predictions lowest in the far west and highest along the eastern coast."} surf <- predict_surface(rf_fit, n_cells = 3000, covariates = points_sf, boundary = nc_boundary) ggplot() + geom_sf(data = surf, aes(colour = .pred), size = 0.6) + geom_sf(data = nc_boundary, fill = NA, colour = "grey20") + scale_colour_viridis_c(name = "Predicted y") + theme_void() + ggtitle("Predicted surface") ``` A fitted model returns a number for any location you hand it, including locations whose predictor values look nothing like the training data. `area_of_applicability()` marks where the cross-validated score actually applies. Pass the folds you validated with: ```{r aoa, fig.height=3.6, eval=has_ggplot && has_ranger, fig.alt="The same predicted surface with cells outside the area of applicability blanked out. Nothing is blanked, because the grid took its covariates from the nearest observation and so holds no combination the model had not already seen."} aoa <- area_of_applicability(surf, model = rf_fit, folds = folds_blocked) cat(sprintf("inside: %d | outside: %d | undetermined: %d | DI threshold %.3f\n", aoa$n_inside, aoa$n_outside, aoa$n_na, aoa$threshold)) # AOA is NA wherever a predictor was missing or non-finite, and `!NA` is NA, # which R silently skips in a subscripted assignment. Test for TRUE and let # anything else count as outside. inside <- aoa$aoa$AOA %in% TRUE surf$.pred_masked <- ifelse(inside, surf$.pred, NA_real_) ggplot() + geom_sf(data = surf, aes(colour = .pred_masked), size = 0.6) + geom_sf(data = nc_boundary, fill = NA, colour = "grey20") + scale_colour_viridis_c(name = "Predicted y", na.value = "grey85") + theme_void() + ggtitle("Predicted surface, extrapolations blanked out") ``` Nothing is masked here, and the reason is worth understanding rather than taking as reassurance: `predict_surface(covariates = points_sf)` copies predictor values from the **nearest observation**, so every grid cell carries a predictor vector some training point already had. Its dissimilarity index is therefore near zero by construction. The index earns its keep when the covariates come from somewhere else — a raster, a different survey, a future scenario. Hand it values outside the training range and it says so: ```{r aoa-outside, eval=has_ggplot && has_ranger} extreme <- st_sf( elevation = c(mean(points_sf$elevation), max(points_sf$elevation) * 4), pop_density = c(mean(points_sf$pop_density), max(points_sf$pop_density) * 4), geometry = st_geometry(points_sf)[1:2] ) area_of_applicability(extreme, model = rf_fit, folds = folds_blocked)$aoa[, c("DI", "AOA")] ``` --- ## Summary | Tessellation | Cells | Notes | |:---|---:|:---| | Voronoi | `r if (has_ggplot) nrow(tess_voronoi$cells) else NA` | Adapts to point density via k-means seeds | | Hex grid | `r if (has_ggplot) nrow(tess_hex$cells) else NA` | Uniform hexagons, good for regular sampling | | Square grid | `r if (has_ggplot) nrow(tess_square$cells) else NA` | Simplest regular grid | | Delaunay | `r if (has_ggplot && has_geom) nrow(tess_tri$cells) else NA` | One triangle per Delaunay triplet; finest resolution | The choropleths show how each tessellation aggregates the response. Section 5 shows the thing that matters most: on the same data and the same fitted model, the fold scheme moves the reported score substantially, and only the blocked number describes the task you actually have. ### Where to go next This page is the wide view. Each step has a vignette that stops and argues about it: | for | read | |:---|:---| | getting your own data in, and what the coordinates are in | `vignette("getting-started")` | | how many cells, and the four criteria that disagree | `vignette("resolution")` | | the five fold schemes, and sizing a block by measurement | `vignette("spatial-cross-validation")` | | residual autocorrelation, aggregation standard errors, and two ways to leak | `vignette("diagnostics")` | | handing the regions to someone else, and what to report | `vignette("reporting")` | Ten numbered scripts in `system.file("scripts", package = "spatialkit")` run the same ground at the console, printing what to look for in each figure before drawing it. ```{r session-info, echo=FALSE, eval=TRUE} sessionInfo() ```