---
title: "spatialkit: Tessellations, Spatial Cross-Validation and Models"
author: "Justin Chase"
date: "2026-08-26"
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. Build four tessellation types (Voronoi, hex, square, Delaunay)
2. Assign points to cells and aggregate with `summarize_by_cell()`
3. Draw choropleths of the cell-level mean response
4. Estimate the autocorrelation range and build spatial CV folds
5. Fit a model, and score it under **blocked** versus **random** folds
6. 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:

```{r backend-report, echo=FALSE}
data.frame(
  package   = c("ggplot2", "geometry", "ranger", "gstat", "GWmodel + sp"),
  available = c(has_ggplot, has_geom, has_ranger, has_gstat, has_gwmodel),
  used_for  = c("all plots", "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}
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')"
))
```

---

## 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`,
which is anticonservative when points inside a cell are spatially correlated.
`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"}
make_choropleth(tess_voronoi, nc_boundary, points_sf, "Voronoi tessellation")
```

```{r choro-hex, fig.cap="Hex grid choropleth"}
make_choropleth(tess_hex, nc_boundary, points_sf, "Hexagonal grid")
```

```{r choro-square, fig.cap="Square grid choropleth"}
make_choropleth(tess_square, nc_boundary, points_sf, "Square grid")
```

```{r choro-tri, fig.cap="Delaunay choropleth", eval=has_ggplot && has_geom}
make_choropleth(tess_tri, nc_boundary, points_sf, "Delaunay triangulation")
```

```{r comparison-panel, fig.width=7, fig.height=2.8, fig.cap="All three at a glance", eval=has_ggplot && requireNamespace("patchwork", quietly = TRUE)}
library(patchwork)

(make_choropleth(tess_voronoi, nc_boundary, points_sf, "Voronoi") |
 make_choropleth(tess_hex,     nc_boundary, points_sf, "Hex") |
 make_choropleth(tess_square,  nc_boundary, points_sf, "Square")) +
  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 `max(attr(sac, "directional"))` explicitly. It returns `NA`
(still classed `sac_range`, so it prints as a bare `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 fit is attached either way, so `plot(fit, type = "variogram")` on a fitted
model can draw the curve and let you judge it rather than trust it — the
distance axis is labelled in the units of the CRS the variogram was fitted in,
and a fit that did not converge says so in the caption. 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}
folds_random  <- make_folds(points_sf, k = 5, method = "random_kfold", seed = 42)
folds_blocked <- make_folds(points_sf, k = 5, method = "block_kfold",  seed = 42)

c(random = folds_random$k, blocked = folds_blocked$k)
```

`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=3, eval=has_ggplot && requireNamespace("patchwork", quietly = TRUE)}
library(patchwork)
plot_folds(folds_random,  points_sf, boundary = nc_boundary) +
  plot_folds(folds_blocked, points_sf, boundary = nc_boundary)
```

### 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}
rf_args <- list(include_coords = TRUE, num_trees = 300, seed = 1)

cv_random  <- do.call(cv_rf, c(list(points_sf, "y", c("elevation", "pop_density"),
                                    folds = folds_random),  rf_args))
cv_blocked <- do.call(cv_rf, c(list(points_sf, "y", c("elevation", "pop_density"),
                                    folds = folds_blocked), rf_args))

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 R<sup>2</sup> = %.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}
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)
cat(sprintf("Bandwidth: %.1f | in-sample R2: %.3f | RMSE: %.3f\n",
            gwr_fit$info$bandwidth, gwr_met$R2, gwr_met$RMSE))
```

### 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
```

---

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

```{r session-info, echo=FALSE, eval=TRUE}
sessionInfo()
```