2 Demonstrating Clamping Under Controlled Extrapolation
2.1 Conceptual setup
To make the role of clamping explicit, we use a deliberately biased sampling design.
The idea: We simulate a landscape with a known response surface (we define the true ecological pattern ourselves). Then we pretend to “survey” it, but we only place sample points in a narrow slice of the available environmental gradient, leaving most of the landscape unsampled. We fit the model using only those biased samples, then project it across the entire landscape.
This creates two clearly distinguishable domains:
- a supported domain (the slice where we sampled): projected covariate values fall within the training range; the model is interpolating.
- an unsupported domain (the rest of the landscape): projected covariate values fall outside the training range; the model is extrapolating.
Because we simulated the data, we know the true response surface and can directly measure how well each projection approach reproduces it and how much clamping helps.
The example is intentionally exaggerated. The sampling bias is much stronger than you would encounter in most field studies. This is deliberate: we want to make the mechanism of clamping easy to see. In your own data, the effects may be subtler, but the principle is identical.
2.2 Step 1: Simulate a known landscape and response surface
We begin by simulating a raster landscape with two covariates
(bin1 and bin2: binary habitat
presence/absence surfaces). We then smooth those covariates using
known kernel parameters to define the true scales of effect.
This smoothed version represents the “ground truth” environmental
gradient that actually drives the ecological response.
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
)We standardize each smoothed covariate across the full landscape. This sets the true environmental gradient: the values a model with perfect knowledge would use. Later, when we fit the model using only biased samples, the scaling will be done using only the sampled subset. That mismatch between “full-landscape scaling” and “sample-based scaling” is the source of the extrapolation problem we want to illustrate.
We now define the true response surface using known coefficients: a
strong positive effect of bin1 (habitat presence is good)
and a negative effect of bin2 (its presence is bad). We
apply a saturating transformation (tanh) to the linear
predictor before converting to counts; this prevents the predicted
surface from becoming arbitrarily large at the extremes, which would
swamp the visualization of what clamping actually does.
# 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 Step 2: Impose a biased sampling design
We now simulate the “field survey.” Instead of sampling randomly
across the full landscape, we restrict sample points to a narrow slice
of the available environmental gradient: locations with moderately high
bin1 values and intermediate bin2 values.
Think of this as a field survey where access constraints, land ownership, or opportunistic design led you to sample in a particular portion of the landscape, one that does not represent the full range of available conditions. The model you fit from these points will work well within that sampled range, but you are asking it to predict across a much wider range.
# 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 left panel shows the restricted area where we sampled; the right
panel shows the true mean count surface. Notice that the sampled area
(left) covers only a portion of the landscape, missing both the
high-count areas driven by high bin1 and the
low-count areas driven by low bin1 or high
bin2. The model will have no direct information about these
extremes.
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)
2.4 Step 3: Fit the multiscale model
We now prepare the inputs required by multiScaleR, fit a
base Poisson GLM, and optimize scales of effect with
multiScale_optim(). A critical detail: all scaling and
centering is now derived from the sampled points only,
not the full landscape. This means the standardized covariate values on
the full raster will be centered and scaled relative to the sample mean
and SD, which may be quite different from the full-landscape mean and
SD.
# 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 <- multiScaleR::multiScale_optim(
fitted_mod = mod,
kernel_inputs = kernel_inputs
)
#>
#>
#> Optimization complete
summary(opt_mod)
#>
#> Call:
#> multiScaleR::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 = data)
#>
#> 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
