multiScale_optim
We are now ready to optimize the scales of effect with
multiScale_optim(). This function searches over σ values
for each raster layer (and simultaneously updates regression
coefficients) to find the combination that maximizes the likelihood of
the fitted model.
Parallelization: For this vignette, the function is
run without parallelization. In practice, optimization will be
substantially faster when you specify n_cores (e.g.,
n_cores = 4). The optimization below takes approximately 50
seconds on a single core.
opt1 <- multiScale_optim(fitted_mod = mod0,
kernel_inputs = kernel_inputs)
We get two helpful warning messages:
- The estimated σ for one raster layer is near
max_D,
meaning the optimization is pushing against the boundary of the search
space. This suggests that either the true scale of effect is large, or
that the variable has a weak or absent relationship with the response
and optimization is poorly constrained. The warning suggests increasing
max_D in kernel_prep().
- The standard error for one or more σ estimates is large relative to
the point estimate: a sign of imprecision, often caused by weak signal
or too small a sample size.
Let’s look at our results.
##
## Call:
## multiScale_optim(fitted_mod = mod0, kernel_inputs = kernel_inputs,
## n_cores = 5)
##
##
## Kernel used:
## gaussian
##
## ***** Optimized Scale of Effect -- Sigma *****
##
## Mean SE 2.5% 97.5%
## land1 248.1714 34.19032 180.2951 316.0478
## land2 538.4724 50.76299 437.6951 639.2496
## land3 538.9632 309.97857 20.0000 1154.3485
##
##
## ====================================
##
## ***** Optimized Scale of Effect -- Distance *****
## 90% Kernel Weight
##
## Mean 2.5% 97.5%
## land1 408.21 296.56 519.85
## land2 885.71 719.94 1051.47
## land3 886.52 32.90 1898.73
##
##
## ====================================
##
## ***** Fitted Model Summary *****
##
##
## Call:
## glm(formula = counts ~ site + land1 + land2 + land3, family = poisson(),
## data = dat)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.37648 0.09260 4.066 4.79e-05 ***
## site 0.06417 0.09891 0.649 0.5165
## land1 -0.50803 0.08013 -6.340 2.30e-10 ***
## land2 0.75707 0.09939 7.617 2.60e-14 ***
## land3 0.18076 0.09994 1.809 0.0705 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 224.24 on 99 degrees of freedom
## Residual deviance: 106.99 on 95 degrees of freedom
## AIC: 313.66
##
## Number of Fisher Scoring iterations: 5
##
##
## WARNING!!!
## The estimated scale of effect extends beyond the maximum distance specified.
## Consider increasing max_D in `kernel_prep` to ensure accurate estimation of scale.
##
##
## WARNING!!!
## The standard error of one or more `sigma` estimates is >= 50% of the estimated mean value.
## Carefully assess whether or not this variable is meaningful in your analysis and interpret with caution.
How to read the summary() output:
The top of the output shows the estimated σ (scale parameter) for
each raster layer, with standard errors and 95% confidence intervals.
Below each σ estimate is the effective distance: the distance
that captures a specified proportion (default: 90%) of the total kernel
weight. This converts σ into an ecologically interpretable distance in
map units (meters in this case). A Gaussian kernel with σ = 250 has 90%
of its cumulative weight within approximately 412 m of the focal
point.
The bottom of the output is the standard model summary from your
regression (GLM in this case), showing regression coefficients, standard
errors, and p-values for each covariate, exactly as you would see from
summary(glm(...)).
Confidence intervals on σ: delta method vs. profile
likelihood
By default, confidence intervals on σ are computed using the
delta method: a first-order Taylor approximation that
treats the log-likelihood surface as locally quadratic near the optimum
and derives the CI from the curvature (Hessian) of that surface.
Delta-method CIs are fast to compute and work well when the likelihood
surface is smooth and roughly symmetric around the estimate. They are
reported alongside the standard error in the default
summary() output.
However, the log-likelihood surface for σ is often asymmetric,
especially when the scale parameter is near max_D, when
sample sizes are small, or when the effect size is weak. In these
situations the delta-method CI can be misleading: the lower and upper
bounds may not accurately reflect the true uncertainty in each
direction. Profile-likelihood CIs address this by
tracing the likelihood surface directly: for each boundary value of σ,
the optimizer finds the configuration of all other parameters that
maximizes the likelihood, then identifies where that profile drops by
the critical chi-squared quantile (1.92 log-likelihood units for a 95%
CI). Profile CIs are therefore more reliable when the surface is
irregular, at the cost of additional computation.
To request profile-likelihood CIs, pass profile = TRUE
to summary():
summary(opt2, profile = TRUE)
Profile CIs are recommended whenever:
- a delta-method CI is wide or asymmetric relative to the point
estimate,
- an estimated σ is near
max_D (boundary
constraint),
- sample size is small (fewer than ~75 sites), or
- you intend to report formal uncertainty bounds in a manuscript.
Wide CIs on σ in either method indicate that the data do not strongly
constrain the scale of effect for that variable. This may mean the
variable has a weak relationship with the response, the sample size is
insufficient, or the variable should be excluded from the model.
For our full model: site has no apparent effect,
land1 has a negative effect, land2 has a
positive effect, and land3 shows a weak positive effect.
Crucially, the σ for land3 is large and
imprecisely estimated (wide CI). This is a red flag: when a variable has
no true relationship with the response, the optimization often fails to
find a well-defined σ and will wander toward large values.
With this simulated data, we know that land3 has no
effect on counts. Let’s fit a model without it.
## New model
mod0_2 <- glm(counts ~ site + land1 + land2,
family = poisson(),
data = df)
## Optimize
opt2 <- multiScale_optim(fitted_mod = mod0_2,
kernel_inputs = kernel_inputs)
No warnings! Let’s look at our results.
##
## Call:
## multiScale_optim(fitted_mod = mod0_2, kernel_inputs = kernel_inputs,
## n_cores = 5)
##
##
## Kernel used:
## gaussian
##
## ***** Optimized Scale of Effect -- Sigma *****
##
## Mean SE 2.5% 97.5%
## land1 238.4505 33.14563 172.6569 304.2441
## land2 499.5337 47.47022 405.3061 593.7614
##
##
## ====================================
##
## ***** Optimized Scale of Effect -- Distance *****
## 90% Kernel Weight
##
## Mean 2.5% 97.5%
## land1 392.22 284.00 500.44
## land2 821.66 666.67 976.65
##
##
## ====================================
##
## ***** Fitted Model Summary *****
##
##
## Call:
## glm(formula = counts ~ site + land1 + land2, family = poisson(),
## data = dat)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.39293 0.09098 4.319 1.57e-05 ***
## site 0.17476 0.08246 2.119 0.0341 *
## land1 -0.48620 0.07922 -6.137 8.40e-10 ***
## land2 0.62567 0.07245 8.636 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 224.24 on 99 degrees of freedom
## Residual deviance: 109.57 on 96 degrees of freedom
## AIC: 314.24
##
## Number of Fisher Scoring iterations: 5
The default summary() uses delta-method confidence
intervals. Because this model is well-behaved (no boundary warnings,
reasonable sample size), the delta-method CIs should be reliable here.
For a more conservative check, you can compare them to
profile-likelihood CIs:
summary(opt2, profile = TRUE)
The true data-generating values were (see section 3):
land1 σ |
250 m |
Check CI |
land2 σ |
500 m |
Check CI |
land1 β |
−0.5 |
Check CI |
land2 β |
0.7 |
Check CI |
Overall, we do a good job recovering both the regression coefficients
(β) and the scales of effect (σ) within reasonable uncertainty bounds.
This is reassuring but also a good reminder: with real data, you will
not have ground truth to check against. Careful model selection and
biological reasoning are your guides.
Visualizing the kernel decay. The
plot() method shows how the influence of each landscape
variable falls off with distance. The vertical line marks the distance
enclosing 90% of the cumulative kernel weight (the effective
distance).
## Kernel function
plot(opt2)
## Kernel function; 99% contribution
plot(opt2, prob = 0.99)
Visualizing marginal covariate effects.
plot_marginal_effects() shows the estimated effect of each
covariate on the response, holding other covariates at their mean
values. These are the standard partial effect plots from the regression
model, applied to the optimally-scaled covariate values.
plot_marginal_effects(opt2)
Profiling the likelihood surface.
profile_sigma() provides a diagnostic view of how the model
fit changes as σ varies across its search range for each covariate,
while all other σ values are held at their optimized values. Rather than
reporting a single point estimate with an interval, it shows the full
shape of the AICc (or log-likelihood) surface across the parameter
space.
This is useful in several situations:
- Verifying a well-defined optimum: A clear trough in
the AICc surface centered near the optimized σ suggests that the scale
of effect is precisely estimated. A flat or monotone profile indicates
the data do not strongly constrain σ for that variable.
- Asymmetry diagnostics: If the surface drops steeply
on one side and gradually on the other, the likelihood is asymmetric and
delta-method CIs may not capture the uncertainty accurately.
Profile-likelihood CIs (via
summary(..., profile = TRUE))
are more appropriate in this case.
- Boundary concerns: When the trough is near
max_D, the surface may be truncated by the search boundary
rather than identified by the data. This is a visual confirmation of the
warning message from multiScale_optim().
Intervals are log-spaced so that evaluation is denser at smaller σ
values, where the surface tends to change more rapidly, and sparser at
larger values where it is typically flatter.
prof2 <- profile_sigma(opt2, n_pts = 15, verbose = FALSE)
plot(prof2)
For opt2, the profile confirms a clear, well-defined
minimum in AICc near the optimized σ for both land1 and
land2, consistent with the absence of boundary warnings and
the narrow confidence intervals reported by summary().
Compare this to the full model (opt1), which included the
uninformative land3 covariate: that profile would show a
shallow or absent trough, reflecting the poor identifiability of σ when
a variable has no true relationship with the response.
We can also apply the optimized kernel to the full raster surfaces.
This is a key step if you want to generate spatially explicit
predictions across the landscape.
kernel_scale.raster() smooths each raster layer using
the estimated σ from the optimization. When
scale_center = TRUE, the smoothed values are standardized
using the mean and SD from the sample points, matching the scaling used
during model fitting. This ensures that terra::predict()
receives covariate values in the same space as the training data.
rast_opt <- kernel_scale.raster(raster_stack = land_rast,
multiScaleR = opt2)
##
## Smoothing spatRaster 1 of 2: land1 at sigma = 238
##
## Smoothing spatRaster 2 of 2: land2 at sigma = 499
Note: Because site is a non-spatial
covariate (a site-level variable, not a raster), a constant dummy raster
is created for it. This allows the spatRaster object to be
passed directly to terra::predict(), which requires all
model covariates to be present. See the
Spatial Projection and Clamping vignette for a detailed
treatment of projection, including how to handle cases where projected
covariate values fall outside the range of the training data.
