--- title: "Diagnostics" description: > Residual autocorrelation, aggregation standard errors, kriging adequacy, the area of applicability, two ways a score can leak, and the warnings about pooled cell values. 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{Diagnostics} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include=FALSE} # gstat fits the variograms behind the design effect and the kriging check, and # ggplot2 draws the area of applicability. Both are Suggests, so the vignette is # gated on them and R CMD build succeeds without either. # `eval` below gates CHUNKS only: knitr does not apply it to inline `r ...` # code, which runs regardless. Any inline reference to something a chunk # defines has to guard itself -- `r if (runs) x else NA` -- or the vignette # fails to build on a machine without these packages instead of skipping. has_gstat <- requireNamespace("gstat", quietly = TRUE) has_ggplot <- requireNamespace("ggplot2", quietly = TRUE) runs <- has_gstat && has_ggplot knitr::opts_chunk$set( echo = TRUE, message = FALSE, warning = FALSE, fig.align = "center", out.width = "100%", eval = runs ) ``` ```{r missing-deps, echo=FALSE, eval=!runs, results='asis'} cat("**gstat or ggplot2 is not installed, so the code below is shown but not", "run.** Install them with `install.packages(c(\"gstat\", \"ggplot2\"))`", "and rebuild.\n") ``` *This article needs two optional packages: **gstat**, which fits the variogram behind the autocorrelation range, and **ggplot2** for the figures. When either is missing the code is shown but not run, and a note at the top says so.* ## Four questions a score does not answer A cross-validated RMSE tells you how far the predictions fell from the held-out values. It is silent on four things that decide whether the analysis holds up: whether the model left spatial structure in its residuals, whether the standard errors on aggregated values are the width they claim, whether a cell mean is the best estimate of that cell available, and where on the map the score applies at all. A fifth section shows two ways the score itself can be wrong, and a sixth reads the warnings `summarize_by_cell()` raises about pooled cell values. ## A fixture and a model that is wrong in a known way The response is driven by one real predictor and by a spatial field the model never sees. `b` is noise. A linear model on `a` and `b` therefore leaves the field in its residuals, which is what the first diagnostic should detect. ```{r fixture} library(sf) library(spatialkit) set.seed(11) n <- 350 xy <- data.frame(x = 5e5 + runif(n, 0, 1000), y = 5e6 + runif(n, 0, 1000)) D <- as.matrix(dist(xy)) field <- as.numeric(t(chol(exp(-D / 70) + diag(1e-8, n))) %*% rnorm(n)) xy$a <- rnorm(n) xy$b <- rnorm(n) xy$z <- 1.5 * xy$a + field + rnorm(n, sd = 0.3) pts <- st_as_sf(xy, coords = c("x", "y"), crs = 32632) ``` ```{r learner} lm_fit <- function(train_sf, response_var = "z", predictor_vars = c("a", "b")) { d <- st_drop_geometry(train_sf) # With no predictors the model is intercept-only: "z ~ 1", not "z ~ ". rhs <- if (length(predictor_vars)) paste(predictor_vars, collapse = " + ") else "1" fml <- as.formula(paste(response_var, "~", rhs)) new_spatial_fit( subclass = "lmsurf_fit", engine = lm(fml, data = d), formula = fml, response_var = response_var, predictor_vars = predictor_vars, data_sf = train_sf ) } predict.lmsurf_fit <- function(object, newdata = NULL, ...) { if (is.null(newdata)) return(fitted(object$engine)) as.numeric(predict(object$engine, newdata = st_drop_geometry(newdata))) } residuals.lmsurf_fit <- function(object, ...) residuals(object$engine) fitted.lmsurf_fit <- function(object, ...) fitted(object$engine) for (g in c("predict", "residuals", "fitted")) registerS3method(g, "lmsurf_fit", get(paste0(g, ".lmsurf_fit"))) fit <- lm_fit(pts) ``` ## 1. Structure left in the residuals ```{r moran} mi <- residual_morans_i(fit) mi ``` A large positive Moran's I says nearby residuals resemble each other, so the model is missing something that varies smoothly across the map. Here that is the simulated field, by construction. Three parts of the print earn attention. The **null** says which moments the p-value came from: `residual` moments condition on the design matrix and are exact for OLS residuals, `randomisation` moments treat the residuals as exchangeable and are an approximation for any model. `null = "auto"` chooses per fit, so read what it chose. The **design rank** is the rank of the rebuilt model matrix, which is where the residual degrees of freedom come from. The **kurtosis** says how far the residuals sit from the Gaussian case the moments assume, with 3 as the reference. The weight matrix is \(n \times n\) and is described instead of returned. Ask for it when you want to reuse the same neighbours: ```{r keep-weights} mi_w <- residual_morans_i(fit, keep_weights = TRUE) dim(mi_w$weights) mi$weights_summary[c("n", "storage", "neighbours", "kept")] ``` `weights_summary` is present either way, so a `NULL` in `weights` is always explained. ## 2. Standard errors on aggregated values Aggregate the points into cells and the naive standard error of each cell mean is \(s/\sqrt{n}\), which treats the points in a cell as independent. For the cell's own mean that is the right standard error when the points are spread through the cell. For the cell mean as an estimate of the population (grand) mean it is not: on an autocorrelated field the points in a cell share the cell's departure from that mean, and the effective sample size is smaller than the count. ```{r cells} hull <- clip_target_for(pts, expand = 0.02, quiet = TRUE) tess <- build_tessellation(pts, boundary = hull, method = "square", approx_n_cells = 24, quiet = TRUE) asg <- assign_features_to_polygons(pts, tess$cells) c(cells = nrow(tess$cells), points = nrow(asg), median_per_cell = median(table(asg$poly_id))) ``` `summarize_by_cell()` takes a design effect, and three settings are worth comparing on the same cells. ```{r deff} sac <- estimate_sac_range(pts, "z") iid <- summarize_by_cell(asg, "z", cells_sf = tess$cells, deff = 1) kish <- summarize_by_cell(asg, "z", cells_sf = tess$cells, deff = "kish") vgm <- summarize_by_cell(asg, "z", cells_sf = tess$cells, deff = "variogram", sac = sac) data.frame(n = iid$n, iid = round(iid$..se_resp_z, 3), kish = round(kish$..se_resp_z, 3), variogram = round(vgm$..se_resp_z, 3))[1:5, ] ``` ```{r deff-inflation, echo=FALSE, results='asis'} ik <- median(kish$..se_resp_z / iid$..se_resp_z, na.rm = TRUE) iv <- median(vgm$..se_resp_z / iid$..se_resp_z, na.rm = TRUE) cat(sprintf(paste("The median cell's standard error is %.2f times wider under", "`kish` and %.2f times wider under `variogram` than the", "independent-sampling one. As an estimate of the population", "mean, a confidence interval built on the naive number is", "about half the width it should be on this field. For the", "cell's own mean the naive standard error is the right one,", "and section 3 (`kriging_adequacy()`) covers that", "quantity.\n"), ik, iv)) ``` `deff = "kish"` estimates an intra-class correlation from the within-cell and between-cell variance and applies Kish's factor. `deff = "variogram"` takes the correlation from the fitted variogram model, so it uses the spatial structure rather than the cell labels. When no usable variogram exists that request falls back to `"kish"` with a warning, and the `deff_method` column names the correction each row carries. The choice is recorded: ```{r deff-attrs} str(attr(vgm, "deff_applied")) attr(kish, "icc")$resp ``` A cell whose design effect reaches its point count carries no information about the variance at all, and its standard error is `NA`. That is a statement about the cell, so the fix is bigger cells or more points, not a different setting. ## 3. Is the cell mean the best estimate of the cell? The plain mean weights every point in a cell equally and ignores everything outside it. Block kriging uses the fitted variogram, borrows from neighbouring cells, and reports a variance derived from the spatial model. `kriging_adequacy()` runs both and compares them. ```{r kriging} ka <- kriging_adequacy(asg, "z", tess$cells, k = 4, seed = 1) ka ``` Four readings. `plain_var` is the error variance of the plain mean under the fitted variogram; a naive `se^2` below it says the naive standard error was optimistic. `kr_gain` is the share of that variance kriging removes in the cell. `kr_shift` is how far the kriged estimate moved from the plain mean, over the standard deviation the variogram gives that difference: about a third of cells beyond 1 and one in twenty beyond 2 is what a right variogram produces, and more says the means and the model disagree. The blocked cross-validation at the end standardises each held-out error by its kriging variance: a variance of 1 means the kriging variance is right, and above 1 means it is understated. ```{r kriging-cols} head(st_drop_geometry(ka)[, c("n", "mean", "se", "plain_var", "kr_pred", "kr_var", "kr_gain", "kr_shift")], 4) ``` An empty cell still gets a kriged estimate and variance, which is the practical reason to run this: a tessellation fine enough to be useful usually has cells with no observation in them. This grid has none, as the report's `empty cells` line says, because 25 cells over 350 points leaves every cell populated; the reporting vignette's grid over the North Carolina counties has five. ## 4. Where the score applies A cross-validated number describes the region the training data covers. For predictions made anywhere else, it describes nothing. `area_of_applicability()` marks the difference, following Meyer and Pebesma (2021). ```{r aoa} grid <- create_grid_polygons(hull, target_cells = 200, type = "square") gp <- suppressWarnings(st_centroid(grid)) set.seed(2) gp$a <- rnorm(nrow(gp)) gp$b <- rnorm(nrow(gp), mean = 3) # shifted well outside the training range aoa <- area_of_applicability(gp, model = fit) aoa ``` ```{r aoa-plot, fig.height=4.5, fig.alt="Cumulative distribution of the dissimilarity index. The prediction locations rise far more slowly than the training points, each measured to its nearest other training point since no folds were passed, and most of them sit beyond the dashed threshold, so the cross-validated score does not apply to most of the grid."} plot(aoa) ``` Shifting one predictor three standard deviations moves most of the grid outside the domain the model was trained on. The dissimilarity index `DI` is on `aoa$aoa`, alongside the logical `AOA`, and `threshold` is the cutoff derived from the training data. One trap worth stating plainly. `AOA` is `NA` wherever a predictor was missing or non-finite, `!NA` is `NA`, and R skips `NA` positions in a subscripted assignment. Test for `TRUE`: ```{r aoa-na} inside <- !is.na(aoa$aoa$AOA) & aoa$aoa$AOA c(inside = sum(inside), outside_or_unknown = sum(!inside)) ``` ## 5. Two ways to leak **Coordinates as predictors.** A model given `x` and `y` can reproduce the training surface by memorising location, and random folds will not detect it. `fit_rf_model(include_coords = TRUE)` logs a caution (once per session) for this reason. Score any such model with blocked folds, and read `vignette("spatial-cross-validation")` for what the blocking has to be wide enough to do. **Selection outside the fold.** Choosing predictors on the whole data set and then cross-validating the chosen set gives the selection a look at every test row. `select_features_forward()` does the selection inside the resampling, and it takes the fold scheme as an argument for the same reason the outer loop does. Blocked folds only keep test rows apart from training rows if the blocks are at least as wide as the autocorrelation range. `make_folds()` estimates that range for what the candidates leave unexplained, about 320 m on this field, and warns when the blocks are narrower. The default grid here has blocks about 250 m across, so `block_size = 400` is passed: ```{r selection} # select_features_forward() calls fit_fn(train_sf, vars), so the learner above # needs an adapter that puts `vars` in the predictor slot. Step 0 fits it with # no predictors at all, which is why lm_fit() builds "z ~ 1" for an empty set. sel_fn <- function(train_sf, vars) lm_fit(train_sf, "z", vars) sel <- select_features_forward(pts, "z", c("a", "b"), fit_fn = sel_fn, k = 4, method = "block_kfold", block_size = 400, seed = 1) sel$selected sel$history ``` `history` has one row per candidate scored at each step. Step 0 is the intercept-only model, so `a` had to beat it to be added, and at step 2 adding the noise variable `b` made the cross-validated RMSE slightly worse, so the selection stopped at `a`. A candidate is added when its pooled `gain` exceeds `tol` and the mean of its per-fold gains (`fold_gain`) exceeds one standard error of that mean (`fold_gain_se`), from step 1 on; `rule = "tol"` tests the pooled gain alone. Its inner folds default to the scheme you pass. Passing `random_kfold` there logs a caution, because a leaky inner loop can select a variable for being spatially close to the response and an honest outer loop will then report a respectable number for a dishonestly chosen feature set. ## 6. Warnings about pooled cell values `summarize_by_cell()` returns by default, beside each cell's mean, a pooled value: the mean pulled toward the mean of the cells around it. What it has to say about those values, and about whether the cells are the ones their count was chosen for, comes as R warnings with a class. Each is raised once per call and reason and carries its numbers as fields, so it can be caught and read without matching its text; a warning about the pooled columns also names the columns it applies to. This vignette hides warnings, so the helper below runs a call and returns the classes it raised: ```{r raised} raised <- function(expr) { seen <- character(0) withCallingHandlers(expr, warning = function(w) { seen <<- c(seen, class(w)[1]) invokeRestart("muffleWarning") }) seen } ``` The cells are built for pooled values at 4.5 observations per cell, by the calls `vignette("resolution")` and `vignette("reporting")` walk through. Those are the defaults, so nothing is named but the response; `estimate = "mean"` in the summary gives the plain means alone, and `vignette("resolution")` lists the profile's earlier defaults: ```{r pooled-cells} prof <- resolution_profile(pts, "z") size <- select_resolution(prof) seeds <- get_voronoi_seeds(method = "kmeans", n = size, sample_points = pts) vor <- build_tessellation(seeds, method = "voronoi", quiet = TRUE) asg_v <- assign_features_to_polygons(pts, vor$cells, polygon_id_col = "cell_id") raised(summarize_by_cell(asg_v, "z", id_col = "cell_id", cells_sf = vor$cells)) ``` A summary of the cell value the count was chosen for, on the cells and rows that were scored, raises nothing. Each of the five warnings below is a departure from that. **`spatialkit_estimate_mismatch`: the count was chosen for the other cell value.** The selection wrote down that the count was set for pooled values, the cells and the assigned layer carried that record, and this call asks for plain means: ```{r mismatch} w <- tryCatch( summarize_by_cell(asg_v, "z", id_col = "cell_id", cells_sf = vor$cells, estimate = "mean"), spatialkit_estimate_mismatch = function(w) w) unlist(unclass(w)[c("chosen_for", "reported", "rule", "cells")]) ``` The summary is returned as asked. Report the value the count was chosen for, or choose the count for the value reported: `resolution_profile(estimate = "mean")` with `select_resolution(criterion = "loo")`. **`spatialkit_cells_differ`: these are not the cells and rows that were scored.** The summary fits the response again on what it is given and compares the rows, the cells that hold one and the leave-one-out error with the record. Here twenty rows were dropped after the assignment: ```{r differ} w <- tryCatch( summarize_by_cell(asg_v[-(1:20), ], "z", id_col = "cell_id", cells_sf = vor$cells), spatialkit_cells_differ = function(w) w) w$what rbind(scored = unlist(w$scored), found = unlist(w$found)) ``` `what` names the first comparison that failed: `"n"` for another number of rows with a response, `"cells"` for another number of cells that hold one, `"loo"` for the same rows in the same number of cells with another error, which means rows sit in other cells, and `"record"` when the assigned layer and `cells_sf` carry different records. The usual causes are rows filtered after the assignment, a boundary that leaves points in no cell, and a layer other than the one profiled. On gridded locations a row can lie the same distance from two centres, on the boundary between their cells. The profile scores it in the cell of the lower-numbered centre, and `assign_features_to_polygons()` gives a point there the same cell when the cells carry the record of the centres scored, so a lattice passes the check. `build_tessellation()` draws those cells in the profile's CRS whatever CRS the seeds come back in. Cells built in another projected CRS, named with its `crs`, are the nearest-centre cells of that CRS, so rows near a boundary, on it or not, can sit in another cell, and the warning then speaks of them. The values returned describe the cells as given, and `attr(, "count_rule")$check` keeps the outcome. Cells the caller brings carry no record and are not checked. **`spatialkit_no_structure`: the cells predict no better than one mean.** `b` is noise by construction: ```{r no-structure} noise <- summarize_by_cell(asg_v, "b", id_col = "cell_id") raised(summarize_by_cell(asg_v, "b", id_col = "cell_id")) as.data.frame(attr(noise, "pooled")$columns$resp_b[ c("loo", "loo_one_cell", "structure_p", "own_weight", "structured")]) ``` `structured` is `FALSE` when the leave-one-out error of one mean for all rows (`loo_one_cell`) is within 2 percent of that of the pooled values (`loo`) or below it, or when the F test of the cell means against one mean has `p >= 0.05`. The values are returned, and a map of them still looks smooth, because each is mostly the mean of its neighbours. On 200 simulated layers with white noise as the response the median own weight was 0.11 and a cell's pooled value had a correlation of 0.79 with the mean of its 8 nearest cells' values, where the plain means had one of -0.02 (`dev/sim-pooled-values.R` in the package's repository; `?summarize_by_cell` gives the settings). Read the warning as: one mean describes this column, and the map shows the pooling. **`spatialkit_pool_degenerate`: the values were computed and carry no pooling.** Its `reason` says which of three cases it is. A column that does not vary inside any cell is one of them: ```{r degenerate} asg_v$flat <- 1 w <- tryCatch(summarize_by_cell(asg_v, "flat", id_col = "cell_id"), spatialkit_pool_degenerate = function(w) w) w$reason ``` With `"no spread inside any cell"` and with `"one cell"` every own weight is 1 and the pooled value is the plain mean. With `"no signal between cells"` the cell means differ from their neighbours' by no more than the spread inside the cells explains: every own weight is 0, and each pooled value is the mean of the cell's neighbours and holds none of the cell's own observations. That last case is not the test for a response without structure. In the simulation it came on 20 of the 200 white-noise layers and `spatialkit_no_structure` on all 200. **`spatialkit_pool_few_cells`: 2 to 9 cells.** The neighbours of a cell are then all the other cells, so its target is the mean of every observation outside it: ```{r few-cells} few <- create_grid_polygons(hull, target_cells = 4, type = "square") asg_f <- assign_features_to_polygons(pts, few, quiet = TRUE) raised(four <- summarize_by_cell(asg_f, "z")) as.data.frame(four)[, c("poly_id", "n", "resp_mean_z", "..pooled_resp_z", "..pooled_w_resp_z", "..pooled_target_resp_z")] ``` The overall mean lies between a cell's own mean and that target. A cell whose own weight is below its share of the observations therefore gets a pooled value on the far side of the overall mean from its plain mean, and two cells can change order. The `..pooled_w_*` column shows how close each cell is to that, and `estimate = "mean"` gives the plain means. `?summarize_by_cell` lists these classes and the others of the pooled values under "Conditions", with the fields each carries. ## Next `vignette("spatial-cross-validation")` covers fold schemes and block sizing. `vignette("resolution")` covers choosing the cell count these diagnostics are computed over. `?residual_morans_i` documents the two nulls and when each is exact, `?summarize_by_cell` derives the design-effect correction and reports its measured coverage, and `?area_of_applicability` covers the weighting and the threshold.