##
## Attaching package: 'MASS'
## The following object is masked from 'package:terra':
##
## area
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:
multiScaleR model
across a raster landscape,pct_mx alter projected covariate
values and resulting predictions, andTo 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.
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.
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)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)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
##
## 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
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_muTo 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")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.
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.
pct_mx argumentThe pct_mx argument controls how strictly those clamping
bounds are enforced.
pct_mx = 0 clamps to the exact observed minimum and
maximum.pct_mx > 0 expands the admissible range and
therefore permits limited extrapolation.pct_mx < 0 contracts the admissible range and yields
more conservative predictions.Accordingly, pct_mx should be viewed as a tuning
parameter controlling how cautiously projected values are allowed to
extend beyond the sampled domain.
We now compare four projection strategies:
pct_mx = -0.20),pct_mx = 0), andpct_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_muWe 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(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")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.
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.
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.