Local weighted regression across analog neighborhoods

Fits a weighted local regression of y on covariates within each focal location's analog neighborhood. Analog neighborhoods are defined by environmental similarity, geographic proximity, or both, while covariates capture additional predictors that influence outcomes within each neighborhood but are not captured by to the search dimensions. This generalizes the weighted mean — which averages over all within-neighborhood variation — by resolving variation driven by these auxiliary predictors. Supports ordinary and ridge-penalized weighted least squares. With purely geographic neighborhoods, this is equivalent to geographically weighted regression (GWR); with environment-based neighborhoods, it extends the analog impact model (AIM) framework to incorporate local covariate effects. This function is a wrapper that calls analog_search() with "regression" included in stat.

Usage

analog_regression(
  x,
  pool,
  y,
  weight = NULL,
  covariates,
  x_covariates = NULL,
  geog = NULL,
  env = kernel("gaussian"),
  select = "all",
  k = NULL,
  lambda = 0,
  stat = c("count", "ess", "regression"),
  se = c("none", "ess", "design"),
  normalize = "auto",
  x_cov = NULL,
  coord_type = "auto",
  env_res_adj = "auto",
  geog_res_adj = "auto",
  cell_area_weight = "auto",
  n_threads = NULL,
  progress = FALSE,
  ...
)

Arguments

x

Focal locations for which regressions will be fit. Should be a matrix/data.frame with columns x, y, and environmental variables, or a SpatRaster with environmental variable layers.

pool

The reference dataset to search for analogs. Either:

Matrix/data.frame with columns x, y, and environmental variables, or SpatRaster with environmental variable layers, OR
An analog_index object created by build_analog_index() (for repeated queries).

y

Response variable(s) to model via local regression. Can be a numeric vector, matrix, data.frame, or SpatRaster. Must have exactly the same number of rows/cells as pool. A separate regression is fit for each variable.

weight

Optional pool site weights for use in aggregation. Numeric vector, single-column matrix/data.frame, or single-layer SpatRaster, with one value per row/cell of pool. For aggregation stats like "weighted_mean", "regression", etc., weights multiply through the weighted aggregation alongside any kernel weighting and cell-area weighting; they do not influence which analogs are selected by knn_* modes (selection remains distance-only). They are reported in pair mode as a user_weight column. Values must be non-negative; NA is allowed and treated as 0 (the point is excluded from aggregation). Default NULL means no user-supplied weights.

If you want to exclude a static subset of pool sites entirely, masking pool (and any associated y / covariates) upfront is more efficient than passing weight = 0 for those sites, since the lattice index will not have to scan or distance-compute against them. Use weight = 0 for cases where the mask varies per query against a shared index, or where some sites have a continuous weight and others should be excluded.

covariates

Predictor variables for local regression, supplied at pool locations. Can be a numeric vector (single covariate), matrix, data.frame, or SpatRaster. Must have exactly the same number of rows/cells as pool. Column/layer names carry through to output. These variables are NOT used for the analog search itself – only for regression within each neighborhood.

x_covariates

Optional predictor values at focal (x) locations, at which to evaluate the local regressions. When supplied, the output includes one additional column (pred) or one column per response variable (pred_{yname} for multi-y), giving the fitted value at the focal. Must have exactly the same number of rows/cells as x, and the same column/layer names as covariates. Default NULL returns only coefficients.

env, geog

Per-family distance treatment, each a kernel() object (or NULL). A kernel bundles the hard distance threshold, the weighting kernel shape, and the kernel's scale for one family: environmental (env) or geography (geog). kernel(weight, theta, max, min) where:

max: hard upper distance threshold — candidates beyond it (in that family's distance) are excluded. For env, max may be a single Euclidean radius or a per-variable vector of absolute-difference thresholds (length equal to the number of environmental variables); scalar environmental thresholds are in Mahalanobis units when x_cov is supplied. For geog, max is a single radius (kilometers when coord_type = "lonlat", projected units otherwise).
min: hard lower distance threshold — candidates closer than it are excluded, so retained candidates form an annulus min <= d <= max. Supported only for geog (a single radius, same units as the geographic max); setting min on env is an error. Mainly used to impose a spatial buffer for cross-validation (see analog_cv()).
weight: kernel shape for weighted aggregations — "uniform" (no distance weighting), "gaussian" (exp(-d^2 / (2 theta^2))), or "inverse" (1 / (1 + d / theta)). The overall kernel weight is the product of the two families' weights, so shapes may be mixed (e.g. an inverse environmental kernel with a Gaussian geographic kernel).
theta: the kernel's scale (Gaussian bandwidth, or inverse half-weight distance). See kernel_params() for calibrated values.

A NULL kernel (the default for both) applies no threshold and no weighting for that family. See kernel() for details.

select

Character string specifying the analog selection strategy. One of:

"all" (default): Select all analogs that satisfy the max_env and max_geog constraints.
"knn_env": For each focal, select up to k analogs with smallest environmental distance, subject to filters.
"knn_geog": For each focal, select up to k analogs with smallest geographic distance, subject to filters.

k

Number of nearest analogs to return per focal location for kNN selection modes. Required when select is "knn_geog" or "knn_env"; must be NULL for select = "all".

lambda

Ridge penalty parameter (default: 0, giving ordinary weighted least squares). Higher values shrink high-variance coefficients toward zero, causing the intercept to approach the weighted mean as lambda -> Inf. Useful when some neighborhoods have few analogs relative to the number of covariates, or when covariates are strongly inter-correlated.

stat

Statistic(s) to compute. "regression" is always included. Additional stats like "count", "ess", and "weighted_mean" can be requested alongside regression coefficients. Default includes "count" and "ess" for diagnostics.

se

Standard-error framing to apply to SE-supporting stats ("weighted_mean" and "regression"). One of:

"none" (default): no SE columns are returned.
"ess": effective-sample-size framing. For weighted_mean, SE = sqrt(var_w(y) / n_eff), where n_eff = (Σw)² / Σw² is Kish's effective sample size and var_w(y) = Σwy²/Σw - ȳ_w². For regression, Var(β̂) = σ²_ess · (X'WX + λI)⁻¹, with residual variance corrected using n_eff - p degrees of freedom.
"design": design-based framing (no assumption that weights are precisions). For weighted_mean, SE = sqrt(Σ w²(y - ȳ_w)²) / Σw.

normalize

One of TRUE, FALSE, or "auto" (default). Only used if stat includes "sum_weights" or "tabulate" and pool is a raster. When active, results for these stats are divided by a global scalar so that they represent a fraction of a theoretically "perfect" scenario where the full search area within max_geog is occupied wall-to-wall by cells whose environment exactly matches x. See details under analog_search() for more info.

x_cov

Optional focal-specific covariance matrices for Mahalanobis distance calculations. Should be a matrix or data.frame with one row per focal location and one column per unique covariance component, or a SpatRaster with a layer for each component. For n environmental variables, there are n*(n+1)/2 unique components, ordered as: variances first (diagonals), then covariances (upper triangle by row).

coord_type

Coordinate system type:

"auto" (default): Automatically detect from coordinate ranges.
"lonlat": Unprojected lon/lat coordinates (uses great-circle distance; assumes max_geog is in km).
"projected": Projected XY coordinates (uses planar distance; assumes max_geog is in projection units).

env_res_adj, geog_res_adj

Control the lattice search-index resolution of the environmental and geographic families, each a multiplier on a data-dependent default (targeting ~50 pool points per occupied bin, split between families by effective dimensionality, so it scales with pool size). Each is either:

A non-negative number: 1 uses the default for that family, larger values are finer, smaller are coarser, and 0 deactivates it.
"auto" (the default for both): tune a single overall resolution scale by optimizing compute time on a subsample of focal points. If focal has relatively few rows, tuning is skipped. Not supported when downsample < 1 (set explicit numeric values instead).

A family that the query does not constrain (no corresponding max_* and not the knn sort key) is automatically deactivated, overriding any explicit value (with a message), since binning an unconstrained family only costs time. Ignored if pool is an analog_index (uses the index's resolution).

cell_area_weight

Controls cell-area weighting when pool is a raster. One of "auto" (default; on for raster pools, off otherwise), TRUE (force on; errors if pool is not a SpatRaster), or FALSE (force off). Cell-area weights correct aggregation statistics for non-uniform cell areas (e.g. lonlat grids near the poles, or projected grids on non-equal-area projections); they are computed via terra::cellSize() and normalized to mean 1. When pool is a pre-built analog_index, this argument must agree with the index's stored configuration: cell_area_weight = FALSE errors if the index was built with cell-area weighting on (rebuild the index instead).

n_threads

Optional integer number of threads to use for the computation. If NULL (default), the global RcppParallel setting is used (see RcppParallel::setThreadOptions).

progress

Logical; if TRUE, display a progress bar during computation. Progress tracking works by splitting the focal dataset into chunks and processing them sequentially. Useful for large datasets. Default is FALSE.

...

additional arguments passed to analog_search() (usually unneeded)

Value

A data.frame (or SpatRaster if x is a SpatRaster) with one row per focal location containing:

index, x, y: Focal location identifiers
Columns for any additional stats requested (e.g., count, ess)
coef_intercept: Regression intercept
coef_{name}: Regression slope for each covariate
When se != "none": se_intercept and se_{name} giving standard errors of each coefficient
With multiple y variables: columns are named coef_{coeff}_{varname} (and se_{coeff}_{varname} when SEs are returned)
When x_covariates is supplied: pred (single-y) or pred_{varname} (multi-y) giving the fitted value at each focal.

See analog_search() for column-naming conventions across stats and metadata() for attached metadata attributes.

Details

Method

For each focal location, the function:

Selects analog pool locations based on select, env, geog, and k
Computes distance-based kernel weights for each analog (via the env / geog kernels)
Fits a weighted least squares regression of y on covariates using these weights, with optional ridge penalty lambda
Returns the regression coefficients (intercept + slopes), and optionally their standard errors (see se in analog_search()).

The math: beta = (X'WX + lambda * I_p)^{-1} X'Wy, where W is diagonal weights, X is the design matrix (intercept + covariates), and I_p penalizes covariate coefficients only (not the intercept).

Relationship to Weighted Mean

With lambda -> Inf, covariate coefficients shrink to zero and the intercept converges to the weighted mean. With lambda = 0 (the default), the full local regression is used. If covariates are centered (weighted mean zero within each neighborhood), the intercept equals the weighted mean at any lambda. The lambda parameter thus provides smooth interpolation between a simple weighted average and a full local regression.

Common Configurations

Geographically weighted regression (select = "all" or "knn_geog", with max_geog and geographic weights): Local spatial regression using geographic proximity to define and weight neighborhoods, equivalent to GWR.
Environmental-nearest regression (select = "knn_env", with max_geog and k): Fixed-size neighborhoods of the k most similar environments within geographic range.

Prediction

To evaluate fitted values at the focal locations, pass covariate values at those locations via x_covariates. The output will then include a pred column (single-y case) or pred_{yname} columns (multi-y case) alongside the usual coefficient columns.

If you only need coefficients, leave x_covariates = NULL.

Examples

if (FALSE) { # \dontrun{
# GWR-style spatial regression
gwr_result <- analog_regression(
  x = sites,
  pool = sites,
  y = sites$income,
  covariates = data.frame(education = sites$edu, access = sites$access),
  select = "knn_geog",
  k = 50,
  env = NULL,
  geog = kernel("gaussian", theta = 20),
  se = "ess"
)

# With focal-side covariates, to obtain fitted predictions
pred_result <- analog_regression(
  x = future_sites,
  pool = sites,
  y = sites$income,
  covariates   = data.frame(education = sites$edu),       # at pool
  x_covariates = data.frame(education = future_sites$edu) # at focals
)
head(pred_result[, c("coef_intercept", "coef_education", "pred")])
} # }