Spatial Projection and Clamping

Bill Peterman

library(multiScaleR)
library(terra)
library(MASS)
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:terra':
## 
##     area

1 Introduction to Spatial Projection

Once distance-weighted scales of effect have been estimated with multiScale_optim(), the next step is often to project the fitted model across the landscape to generate spatial predictions of abundance, occurrence, or related ecological responses.

In multiScaleR, this workflow is implemented through kernel_scale.raster(), which applies the fitted smoothing kernels to raster covariates and, when requested, scales and centers projected values so they match the covariate space used to fit the model. This makes spatial projection straightforward, but it also introduces a practical problem: projected covariate values may extend beyond the range represented at the sampled points.

This matters because many ecological models use nonlinear link functions. Under a log-link or logit-link, even moderate extrapolation in covariate space can produce biologically implausible predictions. For this reason, kernel_scale.raster() includes optional clamping, which constrains projected covariate values to user-defined bounds derived from the training data.

This vignette has two goals:

  1. demonstrate how to project a fitted multiScaleR model across a raster landscape,
  2. show how clamping and pct_mx alter projected covariate values and resulting predictions, and

2 Demonstrating Clamping Under Controlled Extrapolation

2.1 Conceptual setup

To make the role of clamping explicit, we use a deliberately biased sampling design. We first simulate a landscape with a known response surface, then fit the model using points drawn from only a restricted portion of the available environmental gradient. The fitted model is therefore trained on a limited range of covariate values, but is asked to make predictions across the full landscape.

This creates two domains:

Clamping acts by restricting projected covariate values in the unsupported domain back toward the sampled range. The example below is intentionally exaggerated. The objective is not to reproduce a realistic field design, but to create a transparent demonstration of what clamping does and how pct_mx changes projection behavior.

2.2 Simulating a known landscape and response surface

We begin by simulating a raster landscape with two covariates. We then smooth those covariates using known kernels to define the true scales of effect in the data-generating process.

set.seed(321)

# Simulate a raster landscape with two binary covariates
rs <- sim_rast(user_seed = 999, dim = 500, resolution = 30)
rs <- terra::subset(rs, c("bin1", "bin2"))

# Apply the TRUE smoothing used to generate the ecological signal
true_sigmas <- c(400, 200)

true_smoothed <- kernel_scale.raster(
  raster_stack = rs,
  sigma = true_sigmas,
  kernel = "gaussian",
  scale_center = FALSE,
  verbose = FALSE
)

For simulation purposes, we next standardize each smoothed covariate across the full landscape. This defines the true environmental gradient available on the map. Later, the fitted model will scale projected values using only the sampled data. That mismatch is intentional, because it helps clarify why extrapolation occurs.

z_stack <- scale(true_smoothed)

We now define a true response surface. The linear predictor includes a strong positive effect of bin1 and a negative effect of bin2. To avoid trivial demonstrations driven by explosive exponential growth, we apply a saturating transformation before converting the linear predictor to the expected count surface.

# Define true model coefficients
alpha <- 0.5
beta1 <- 1.25
beta2 <- -0.75

# Calculate the linear predictor
linpred_true <- alpha + (beta1 * z_stack$bin1) + (beta2 * z_stack$bin2)

# Saturate the linear predictor so extreme gradients remain interpretable
linpred_true <- 4 * tanh(linpred_true / 4)

# Convert to the expected count surface (Poisson mean)
true_mu <- exp(linpred_true)

2.3 Restricting sampling to a biased domain

We next impose a sampling bias by retaining only a narrow subset of the available environmental gradient. In this example, we sample from locations with moderately high values of bin1 and intermediate values of bin2. This forces the fitted model to learn from a truncated portion of predictor space.

# Define a restricted environmental envelope for sampling
q1 <- quantile(values(z_stack$bin1), probs = c(0.50, 0.90), na.rm = TRUE)
q2 <- quantile(values(z_stack$bin2), probs = c(0.25, 0.75), na.rm = TRUE)

# Create a mask that isolates this specific domain
sample_mask <- z_stack$bin1
sample_mask[] <- ifelse(
  z_stack$bin1[] >= q1[1] & z_stack$bin1[] <= q1[2] &
    z_stack$bin2[] >= q2[1] & z_stack$bin2[] <= q2[2],
  1, NA
)

# Sample only within this restricted domain
pts <- spatSample(
  sample_mask,
  size = 150,
  method = "random",
  na.rm = TRUE,
  as.points = TRUE
)

# Extract the true mean surface at sampled points and simulate counts
mu_pts <- terra::extract(true_mu, pts)[, 2]
counts <- rpois(length(mu_pts), lambda = mu_pts)

The model is fit to a narrow slice of the available environmental gradient, but projection is carried out across the full raster extent.

par(mfrow = c(1, 2))

# Plot the restricted sampling mask and points
plot(sample_mask, main = "Restricted Sampling Domain")
points(pts, pch = 16, cex = 0.45)

# Plot the true mean surface and points
plot(true_mu, main = "True Mean Surface")
points(pts, pch = 16, cex = 0.35)

par(mfrow = c(1, 1))

2.4 Fitting the multiscale model

We now prepare the inputs required by multiScaleR, fit a base Poisson model, and optimize scales of effect. All scaling and centering used during fitting are now derived from the sampled points, not the full landscape.

# Prepare multiscale data objects based on the sampled points
kernel_inputs <- kernel_prep(
  pts = pts,
  raster_stack = rs,
  max_D = 1500,
  kernel = "gaussian",
  verbose = FALSE
)

# Combine count data with kernel predictors
dat <- data.frame(counts = counts, kernel_inputs$kernel_dat)

# Fit the base GLM
mod <- glm(counts ~ bin1 + bin2, family = poisson(), data = dat)
# Optimize scales of effect
opt_mod <- multiScale_optim(
  fitted_mod = mod,
  kernel_inputs = kernel_inputs
)
## 
## 
## Optimization complete
summary(opt_mod)
## 
## Call:
## multiScale_optim(fitted_mod = mod, kernel_inputs = kernel_inputs)
## 
## 
## Kernel used:
## gaussian
## 
## ***** Optimized Scale of Effect -- Sigma *****
## 
##          Mean       SE     2.5%    97.5%
## bin1 357.5333 20.04240 317.9248 397.1418
## bin2 193.2692 38.01108 118.1504 268.3880
## 
## 
##   ==================================== 
## 
## ***** Optimized Scale of Effect -- Distance *****
## 90% Kernel Weight
## 
##        Mean   2.5%  97.5%
## bin1 588.09 522.94 653.24
## bin2 317.90 194.34 441.46
## 
## 
##   ==================================== 
## 
##  *****     Fitted Model Summary     *****
## 
## 
## Call:
## glm(formula = counts ~ bin1 + bin2, family = poisson(), data = dat)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.37798    0.04350  31.679  < 2e-16 ***
## bin1         0.47140    0.03579  13.172  < 2e-16 ***
## bin2        -0.29172    0.04991  -5.845 5.05e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 379.23  on 149  degrees of freedom
## Residual deviance: 180.21  on 147  degrees of freedom
## AIC: 644.87
## 
## Number of Fisher Scoring iterations: 5

3 The Extrapolation Problem

When kernel_scale.raster() is called with multiScaleR = opt_mod and scale_center = TRUE, it projects smoothed covariate rasters using the same centering and scaling parameters used to fit the model. This is essential for prediction, but it also reveals the extrapolation problem: raster cells can contain projected covariate values well beyond the range represented in the training data.

We first generate an unclamped projection.

# Project covariates without clamping
r_unclamped <- kernel_scale.raster(
  raster_stack = rs,
  multiScaleR = opt_mod,
  scale_center = TRUE,
  clamp = FALSE,       # Clamping disabled
  verbose = FALSE
)

# Predict the expected counts
pred_unclamped <- predict(r_unclamped, opt_mod$opt_mod, type = "response")

# Calculate error (Predicted - True)
unclamped_error <- pred_unclamped - true_mu

To make the problem explicit, we compare the range of each projected covariate in the training data to the range of the projected raster.

# Extract training values and projected raster values
train_vals <- opt_mod$opt_mod$model[, c("bin1", "bin2")]
proj_vals_unclamped <- as.data.frame(values(r_unclamped))

# Construct a table comparing the ranges
range_tab <- data.frame(
  covariate = c("bin1", "bin2"),
  train_min = c(min(train_vals$bin1, na.rm = TRUE), min(train_vals$bin2, na.rm = TRUE)),
  train_max = c(max(train_vals$bin1, na.rm = TRUE), max(train_vals$bin2, na.rm = TRUE)),
  raster_min_unclamped = c(min(proj_vals_unclamped$bin1, na.rm = TRUE),
                           min(proj_vals_unclamped$bin2, na.rm = TRUE)),
  raster_max_unclamped = c(max(proj_vals_unclamped$bin1, na.rm = TRUE),
                           max(proj_vals_unclamped$bin2, na.rm = TRUE))
)

range_tab
##   covariate  train_min train_max raster_min_unclamped raster_max_unclamped
## 1      bin1 -1.6751655  3.265743           -3.5852971             5.217388
## 2      bin2 -0.7226202  4.078686           -0.7226202            10.321714
par(mfrow = c(1, 3))
plot(true_mu, main = "True Mean Surface")
plot(pred_unclamped, main = "Unclamped Prediction")
plot(unclamped_error, main = "Unclamped Error")

par(mfrow = c(1, 1))

In this example, the projected covariates extend well beyond the range represented in the training data. Under a log-link, those unsupported predictor values can generate very large predictions. The model is forced into covariate space where it lacks empirical support.

4 Understanding and Applying Clamping

Clamping acts on projected covariate values before prediction. When a raster cell falls outside the training-data range for a covariate, its value is truncated to user-defined bounds. This limits marginal extrapolation and can greatly reduce extreme artifacts in projected response surfaces.

Clamping is best viewed as a practical safeguard, not as a correction for model misspecification. It does not solve omitted variables, residual spatial structure, or novel combinations of predictors. It simply constrains individual projected covariate values so they do not extend arbitrarily far beyond the domain represented in the sampled data.

4.1 The pct_mx argument

The pct_mx argument controls how strictly those clamping bounds are enforced.

Accordingly, pct_mx should be viewed as a tuning parameter controlling how cautiously projected values are allowed to extend beyond the sampled domain.

4.2 Comparing different clamping settings

We now compare four projection strategies:

  1. no clamping,
  2. contracted bounds (pct_mx = -0.20),
  3. strict clamping (pct_mx = 0), and
  4. expanded bounds (pct_mx = 0.20).

The values used here are intentionally strong so that the effect of pct_mx is easy to see.

# Contracted bounds
r_cn20 <- kernel_scale.raster(
  raster_stack = rs,
  multiScaleR = opt_mod,
  scale_center = TRUE,
  clamp = TRUE,
  pct_mx = -0.20,
  verbose = FALSE
)

# Strict clamping
r_c0 <- kernel_scale.raster(
  raster_stack = rs,
  multiScaleR = opt_mod,
  scale_center = TRUE,
  clamp = TRUE,
  pct_mx = 0,
  verbose = FALSE
)

# Expanded bounds
r_cp20 <- kernel_scale.raster(
  raster_stack = rs,
  multiScaleR = opt_mod,
  scale_center = TRUE,
  clamp = TRUE,
  pct_mx = 0.20,
  verbose = FALSE
)

# Generate predictions for each clamped raster
pred_cn20 <- predict(r_cn20, opt_mod$opt_mod, type = "response")
pred_c0   <- predict(r_c0,   opt_mod$opt_mod, type = "response")
pred_cp20 <- predict(r_cp20, opt_mod$opt_mod, type = "response")

# Calculate errors
error_cn20 <- pred_cn20 - true_mu
error_c0   <- pred_c0   - true_mu
error_cp20 <- pred_cp20 - true_mu

We can again compare the projected covariate ranges under each setting.

proj_vals_cn20 <- as.data.frame(values(r_cn20))
proj_vals_c0   <- as.data.frame(values(r_c0))
proj_vals_cp20 <- as.data.frame(values(r_cp20))

range_tab_clamp <- data.frame(
  covariate = rep(c("bin1", "bin2"), each = 4),
  setting = rep(c("unclamped", "pct_mx = -0.20", "pct_mx = 0", "pct_mx = 0.20"), times = 2),
  min_value = c(
    min(proj_vals_unclamped$bin1, na.rm = TRUE),
    min(proj_vals_cn20$bin1, na.rm = TRUE),
    min(proj_vals_c0$bin1, na.rm = TRUE),
    min(proj_vals_cp20$bin1, na.rm = TRUE),
    min(proj_vals_unclamped$bin2, na.rm = TRUE),
    min(proj_vals_cn20$bin2, na.rm = TRUE),
    min(proj_vals_c0$bin2, na.rm = TRUE),
    min(proj_vals_cp20$bin2, na.rm = TRUE)
  ),
  max_value = c(
    max(proj_vals_unclamped$bin1, na.rm = TRUE),
    max(proj_vals_cn20$bin1, na.rm = TRUE),
    max(proj_vals_c0$bin1, na.rm = TRUE),
    max(proj_vals_cp20$bin1, na.rm = TRUE),
    max(proj_vals_unclamped$bin2, na.rm = TRUE),
    max(proj_vals_cn20$bin2, na.rm = TRUE),
    max(proj_vals_c0$bin2, na.rm = TRUE),
    max(proj_vals_cp20$bin2, na.rm = TRUE)
  )
)

range_tab_clamp
##   covariate        setting  min_value max_value
## 1      bin1      unclamped -3.5852971  5.217388
## 2      bin1 pct_mx = -0.20 -1.3401324  2.612595
## 3      bin1     pct_mx = 0 -1.6751655  3.265743
## 4      bin1  pct_mx = 0.20 -2.0101986  3.918892
## 5      bin2      unclamped -0.7226202 10.321714
## 6      bin2 pct_mx = -0.20 -0.5780962  3.262949
## 7      bin2     pct_mx = 0 -0.7226202  4.078686
## 8      bin2  pct_mx = 0.20 -0.7226202  4.894423

Notice how pct_mx smoothly regulates the hard boundaries applied to the projected rasters. We can visually inspect the predictions and errors below:

par(mfrow = c(2, 2))
plot(pred_unclamped, main = "Unclamped")
plot(pred_cn20, main = "Clamped: pct_mx = -0.20")
plot(pred_c0, main = "Clamped: pct_mx = 0")
plot(pred_cp20, main = "Clamped: pct_mx = 0.20")

par(mfrow = c(1, 1))
par(mfrow = c(2, 2))
plot(unclamped_error, main = "Error: Unclamped")
plot(error_cn20, main = "Error: pct_mx = -0.20")
plot(error_c0, main = "Error: pct_mx = 0")
plot(error_cp20, main = "Error: pct_mx = 0.20")

par(mfrow = c(1, 1))

These comparisons make two points clear. First, unclamped projection can generate extreme predictions where projected covariate values extend well beyond the sampled range. Second, pct_mx governs how aggressively those extremes are suppressed. Negative values produce the most conservative maps, pct_mx = 0 anchors projection to the observed range, and positive values allow limited re-expansion of the predictor bounds.

4.3 Quantifying prediction error inside and outside the sampled domain

Visual comparisons are helpful, but the effect of clamping is easier to interpret if we separate prediction error in the supported and unsupported domains.

# Create masks representing areas inside and outside the sampled domain
inside_mask <- sample_mask
outside_mask <- sample_mask
inside_mask[] <- ifelse(!is.na(sample_mask[]), 1, NA)
outside_mask[] <- ifelse(is.na(sample_mask[]), 1, NA)

# Helper function to compute RMSE
rmse <- function(pred, truth, mask) {
  p <- values(pred, mat = FALSE)
  t <- values(truth, mat = FALSE)
  m <- values(mask, mat = FALSE)
  
  idx <- !is.na(m) & is.finite(p) & is.finite(t)
  sqrt(mean((p[idx] - t[idx])^2, na.rm = TRUE))
}

# Compile RMSE scores into a table
rmse_tab <- data.frame(
  model = c("unclamped", "pct_mx = -0.20", "pct_mx = 0", "pct_mx = 0.20"),
  RMSE_inside = c(
    rmse(pred_unclamped, true_mu, inside_mask),
    rmse(pred_cn20, true_mu, inside_mask),
    rmse(pred_c0, true_mu, inside_mask),
    rmse(pred_cp20, true_mu, inside_mask)
  ),
  RMSE_outside = c(
    rmse(pred_unclamped, true_mu, outside_mask),
    rmse(pred_cn20, true_mu, outside_mask),
    rmse(pred_c0, true_mu, outside_mask),
    rmse(pred_cp20, true_mu, outside_mask)
  )
)

rmse_tab
##            model RMSE_inside RMSE_outside
## 1      unclamped   0.5951653     3.025820
## 2 pct_mx = -0.20   0.6011070     1.972351
## 3     pct_mx = 0   0.5951541     1.276909
## 4  pct_mx = 0.20   0.5951653     1.482737

In most applications, differences among projection methods are modest in the supported domain because prediction is interpolation. The larger contrasts occur outside the sampled domain, where clamping prevents projected values from drifting arbitrarily far from the covariate space on which the model was fit.

A final caution is worth emphasizing. Clamping can reduce extreme marginal extrapolation, but it does not guarantee ecological realism. If the model is being projected into genuinely novel environments, the resulting map is still an unsupported inference problem. Clamping only reduces how strongly that problem expresses itself.

5 Key Takeaways

Clamping should be viewed as a practical safeguard against extreme marginal extrapolation during spatial projection.

The pct_mx argument provides a simple way to tune how conservative that safeguard should be. Values near zero are often a sensible starting point for applied analyses. In this vignette, more extreme sampling bias and more extreme pct_mx values were used deliberately so that the mechanism of clamping is easy to see.