spatialkit

"spatialkit" namespace, which none of those touch.

Breaking changes

Statistical corrections that change results (third audit pass)

Each item below was measured, and the measurement is quoted so you can judge

whether it affects an analysis you have already run.

• **residual_morans_i() no longer puts weight on a point's own residual.** FNN::get.knn() reports a point's OWN index among its neighbours whenever exact duplicate coordinates are present, which put 1/k on the diagonal of a matrix Moran's I is only defined for with a zero diagonal. On 40 sites x 4 repeats with a response carrying **no** spatial structure, 120 of 160 rows gained a self-weight, mean I came out at +0.086 against E[I] = -0.0063, and **77% of samples were "significant" at p < 0.05** against a nominal 5%. Repeat observations at one site are exactly what make_folds(method = "leave_location_out") is for, so this was a mainstream input. The dense fallback never had the fault, so the statistic also depended silently on whether **FNN** happened to be installed; the two paths now share one neighbour lookup and agree exactly. Requesting k + 1 neighbours and dropping self is **not** sufficient on its own — the slot self occupied displaced a genuine co-located neighbour and left a farther point standing in for it (75 of 400 retained pairs sat at distance 121 where a neighbour at distance 0 existed). Duplicate coordinates are now grouped and answered exactly.
• **residual_morans_i() gains a null argument, defaulting to "auto".** Model residuals are not exchangeable — they are orthogonal to the design matrix — so the classical randomisation moments are wrong for them. At n = 120 with six smooth covariates and independent errors, OLS residuals had mean I = -0.031 against the exchangeable E[I] = -0.008, and the z-score averaged -0.54 with sd 0.90 instead of 0 and 1. The Cliff & Ord (1981, sec. 8.3) regression-residual moments restore mean z = -0.09, sd 1.03 and a 4.3% rejection rate against a nominal 5%, and **agree with spdep::lm.morantest() to machine precision** (verified at 1e-16 through the public function). "auto" applies them only when the fit's residuals really are the OLS residuals on the rebuilt design, which a forest's and a working GWR's are not; the null actually used is reported in the return value.
• **summarize_by_cell() standard errors under a design effect were too small.** s / sqrt(n / deff) corrects the mean's variance for clustering but leaves s^2 biased low by the same clustering: for exchangeable correlation rho, E[s^2] = sigma^2 (n - deff)/(n - 1). The two errors compound. Measured 95% CI coverage at n = 20: **0.905 at rho = 0.3, 0.796 at rho = 0.6, 0.632 at rho = 0.8**; after rescaling by sqrt((n-1)/(n-deff)), 0.952 / 0.952 / 0.953. Applies to deff = "kish", deff = "variogram" and a fixed numeric deff. **The default deff = 1 path is bit-identical to before.**
• **determine_optimal_levels() ranks model-aware candidates on the standardised deviate, not on |Moran's I|.** E[I] and Var(I) both depend on the cell count, so |I| shrinks as k grows whether or not the finer tessellation captures anything. 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 — an |I| ranking prefers the largest candidate for arithmetic reasons alone. Candidates are now ordered by |z| using the Cliff & Ord residual moments (exact here, since the cell-level residuals are OLS residuals by construction); over the same runs z had mean ~0, sd ~1 and a 5% rejection rate of 0.040-0.057 at every k. The "diagnostics" attribute now carries moran_z alongside moran_i.
• **estimate_sac_range() sweeps four azimuths, not two.** A +/-22.5 degree window around 0 and 90 covers exactly **90 of the 180 distinct azimuths** — every direction between 23 and 67 degrees, and between 113 and 157, fell into neither. On simulated fields with 3:1 anisotropy and a true major-axis range of 300, the estimate came back at 255 and 249 for major axes at 0 and 90 degrees but **151 and 147 at 45 and 135**. Since make_folds(auto_range = TRUE) sizes blocks from this number, a diagonally oriented field silently got blocks half as wide as the correlation they were meant to separate. c(0, 45, 90, 135) tiles all 180 azimuths; the same fields now return 255 / 245 / 249 / 228. A direction whose variogram never reaches a sill is excluded rather than taken as a long range, and the directional attribute now has four named entries.
• **fit_bayesian_spatial_model()'s calibrated length-scale prior never reached Stan.** brms::set_prior(spec, class = "lscale") with no coef is a *global* prior, and brms applies a global prior only to coefficients with no individual prior of their own — every lscale coefficient always has one. brms dropped it with a note and Stan received brms's defaults, which made gp_lengthscale_bounds(), the tail calibration and $info$gp_lscale_prior dead weight. Confirmed with brms::make_stancode(): the requested prior is absent under the global form and present under the coefficient-level form, which is now used. $info$gp_lscale_prior is read back from brms::validate_prior(), so it records what brms will actually use.
• **The GP basis was sized against the wrong domain measure.** brms builds the boundary as choose_L(x, c) = c * max(1, max(x) - min(x)) over the pooled, column-centred covariates — the **full range**, not the per-axis half-range in which Riutort-Mayol et al. state their inequalities. Recovering the boundary from make_standata()'s eigenvalues confirms L = c * full range exactly at every c, so the old convention built a boundary **twice as wide** as gp_k was sized for: the GP was under-resolved, and $info$gp_ell_min — the diagnostic meant to catch exactly that — was twice too lenient to fire. The c floor is now brms's own default 1.25 rather than 1.2.
• **fitted() on a gwr_fit could return a coefficient surface.** The search for GWmodel's fitted-value column matched the whole candidate name vector with %in% and took the first hit in the *SDF's* column order — and the local coefficients come first. A predictor named fit, pred, prediction, fitted or yhat therefore returned its own coefficient column, silently: executed in-sample R^2 was **-1.18** against a true 0.981, and residuals(), summary(), model_metrics(), compare_models() and every cv_gwr() fold consumed it without a warning. The search now runs in preference order and excludes any candidate that is also a model term; all five colliding names now give R^2 = 0.981, identical to the renamed control.
• **coef.gwr_fit() returned GWmodel's whole SDF data slot** — 15 columns for a two-predictor fit, of which 3 are coefficients and the rest are standard errors, t-values, the response, the fitted values, residuals and Local_R2. It now returns the model terms only; reach for object$engine$SDF for the rest.
• **estimate_sac_range() is reproducible, and no longer disturbs the caller's RNG.** seed now defaults to 123L rather than NULL. The n_max subsample is an internal approximation, not part of the answer, and leaving it unseeded made the returned range differ between runs on identical input (19531 / 19589 / 19605 on three calls) while silently advancing the caller's stream — and make_folds(auto_range = TRUE) sizes its blocks from that number. Pass seed = NULL for the old behaviour.
• **estimate_sac_range() rejects a non-numeric response.** as.numeric() on a factor returns its level codes, so a factor response produced a variogram of an arbitrary integer relabelling of the categories and the estimated range changed when the levels were reordered (3700 against 2497 on the same data). Factors and character columns are now an error naming the column; logicals are read as 0/1.
• **Fold sets built from a different dataset are refused.** Fold splits are lists of ..row_id values, and row IDs are seq_len(nrow()) unless supplied, so passing cv_gwr() a folds object built from another dataset of the same size applied cleanly — every ID matched, every fold was populated, and the model was scored on splits describing other observations. make_folds() now records a small projection-invariant row fingerprint in params$row_probe, and cv_gwr(), cv_bayes(), cv_spatial() and cv_rf() error rather than proceed. Fold objects from earlier versions carry no fingerprint and are passed through unchecked.
• **evaluate_insample() rejects duplicated names in fits.** model is the key compare_models() joins its metric and Moran's I tables on, so two fits called "GWR" produced a 2x2 cross-join: four rows, every one carrying the first fit's numbers, with the second fit never scored at all.
• **fit_gwr_model() rejects a non-numeric predictor.** gwr.basic() expands contrasts via model.matrix() and fits, but gwr.predict() does not and fails, so the model appeared to fit and then silently predicted all NA.
• **fit_gwr_model() no longer rejects a two-valued continuous response.** The "binary" error is now gated on the response being integer-like. A left-censored or saturated measurement (every observation at a detection limit or a ceiling) has two distinct values and is perfectly continuous; it now warns instead. The guard also runs once per fold inside cv_gwr(), where a small training fold can legitimately hold only two distinct values.
• fitted() returning the wrong length, or nothing, is now an error in summary() and model_metrics() rather than a plausible row count over an all-NA comparison. new_spatial_fit() is the documented extension point, so a subclass with a missing or mis-sized fitted() method is user-reachable.
• The cached fitted() on a bayesian_fit is stamped with the n and a digest of the data it was computed from. The cache environment has reference semantics — which is what makes it survive copy-on-modify — so fit2 <- fit gave both objects the *same* cache, and assigning different data to the copy returned the original's values at the original's length.
• make_folds() drops rows with empty or non-finite coordinates, with a logged warning naming the count, rather than letting an EMPTY POINT reach block_kfold's nearest-block rescue and die with "replacement has length zero".
• When every fold fails, the warning now names the first underlying error. Previously "all 5 folds failed" was the whole diagnosis even when the cause was simply that **brms** or **GWmodel** was not installed.
• make_folds() records the CRS the folds were built in as params$crs. Geographic input is projected by ensure_projected() to a CRS the caller never chose, and block_size, sac_range and buffer are lengths in *that* CRS.
• .morans_i_for_k() returns NA at or below nine cells, where every cell neighbours every other and Moran's I collapses to exactly -1/(k-1) for any residual vector — a function of the cell count alone.
• residual_morans_i(fit, k = 1) works on machines without **FNN**. apply() simplified the length-1 result to a vector, making the neighbour index a 1 x n matrix and every row after the first out of bounds.

Fourth audit pass (adversarial): corrections that change results

Three reviewers were set the task of making the package fail or silently

misbehave on valid input, with their own reproductions. Every finding below

was reproduced here before it was touched; the figures quoted are from those

reproductions.

• **summarize_by_cell(deff = "kish") under-estimated the predictor ICC by about a factor of m.** The pooled one-way ANOVA grouped the m z-scored predictor columns under the same cell label, so independent per-column cell effects averaged away in the shared cell mean and the between-cell sum of squares shrank by ~1/m. Measured at true rho = 0.5 with m = 4: pooled ICC **0.12** against 0.495 per column, so every predictor SE was ~44% too small. The pooled group is now (variable, cell), which recovers 0.49.
• **Design effects are built from what a column actually observes.** A cell of 10 rows with 2 finite responses had its response SE formed at the 10-row design effect, then applied to a 2-observation mean: adding 8 NA-response rows moved the SE from 8.46 to **26.38**. Each column's design effect now uses its own non-missing count, with the mean pairwise correlation recomputed over the observed locations when a cell has NAs; cell_weight is the effective count of the primary variable, not of rows (n still counts rows).
• **CRS-less coordinates get ONE interpretation, wherever they enter.** prep_model_data() assumed EPSG:4326 for CRS-less data that looked like lon/lat and projected it, while every predict() method passed the fit's CRS as a target — a branch that *stamped* it onto the raw numbers. The same rows sat in two different places, and predict(fit, newdata = training rows) disagreed with fitted(fit) by up to one response SD (R² 0.98 in-sample, 0.64 via newdata). The heuristic is now a single function used by both branches; a fit records the assumption it was built under and replays it on CRS-less newdata. Two further symptoms of the same split — CRS-less LINESTRINGs aborting in coerce_to_points() ("crs not found"), and hex/square build_tessellation() refusing input voronoi accepted — are fixed with it. The assumption is now announced with a real R warning.
• **residual_morans_i() refuses a malformed weights matrix** instead of silently substituting the default k-NN(8) matrix (I = 0.874 returned for four malformed shapes against 0.805 for the weights actually supplied).
• **The fold-provenance fingerprint no longer refuses the caller's own data.** Three defects in the version introduced last pass: a character ..row_id was coerced to all-NA (and matched row 1 everywhere); coordinates were compared as "%.7g" strings and flipped on the ~1 in 5000 that a reprojection moved by 5e-9°; and polygon input was probed *after* pointization, so a different pointize in the cv call read as different data. Both sides now probe the geometry as supplied, keep IDs in their own type, and compare numerically within 1e-6°.
• **n_folds_attempted counts the folds supplied.** A fold whose test rows were all removed as incomplete vanished before fitting and was absent from both counts, so five supplied folds reported 4/4. It is announced with a real warning.
• **determine_optimal_levels()'s elbow uses the signed deviation below the chord.** abs() let a concave bump *above* the chord — a k where k-means fell into a worse local optimum — win with the same magnitude.
• **estimate_sac_range() returns NA for a constant variable** (an exactly explained response, or a constant one) instead of a fitted "range" of 168 or 673 from a variogram that is identically zero. make_folds(auto_range = TRUE) no longer re-opens the unseeded subsample by forwarding its own seed = NULL.
• **fit_bayesian_spatial_model() attaches a user-supp