4.1 Load the package and data

multiScaleR ships with a simulated data set of Poisson counts at 100 survey locations and three landscape raster layers.

data("landscape_counts")
dat <- landscape_counts

data("surv_pts")
pts <- vect(surv_pts)

land_rast <- rast(pkg_extdata("landscape.tif"))

The landscape_counts data frame contains counts and a site-level covariate (site). The raster stack has three layers: land1 (binary habitat), land2 (continuous), and land3 (continuous). Counts were simulated with land1 (σ = 250 m, β = −0.5) and land2 (σ = 500 m, β = 0.7) as true predictors; land3 has no effect.

plot(land_rast)

Landscape rasters with survey locations.

plot(land_rast$land1, main = "land1 with survey points")
points(pts, pch = 19, cex = 0.6, col = 'red')

Landscape rasters with survey locations.

4.2 Prepare inputs with kernel_prep()

kernel_prep() extracts raster values within max_D meters of each survey point and computes initial distance-weighted covariate values. Set max_D to a value that comfortably exceeds the scale of effect you expect. If optimization later warns that σ is near max_D, re-run kernel_prep() with a larger value.

kernel_inputs <- kernel_prep(
  pts          = pts,
  raster_stack = land_rast,
  max_D        = 1700,
  kernel       = "gaussian",
  verbose      = FALSE
)

Combine the kernel-weighted covariate values with the response data:

df <- data.frame(dat, kernel_inputs$kernel_dat)

4.3 Fit an initial model

Fit the regression model you intend to optimize. The formula and family you specify here are inherited by multiScale_optim(), so get the model structure right before optimizing. The coefficient values at this stage are only starting values.

mod <- glm(counts ~ site + land1 + land2,
           family = poisson(),
           data   = df)

4.4 Optimize scales of effect

multiScale_optim() searches over σ values for each raster layer while simultaneously updating regression coefficients to maximize the model likelihood. Use n_cores to parallelize and speed up the search.

opt <- multiScaleR::multiScale_optim(fitted_mod    = mod,
                                     kernel_inputs = kernel_inputs)

If optimization is sensitive to the initial σ values, use the screened start strategy. This first scouts a small set of σ candidates, then launches one full optimization from the best candidate. It is most helpful when models have multiple landscape covariates, weak scale signal, or evidence that a single starting value is settling on a local optimum. The marginal σ scans are serial, but the short screening attempts also use n_cores when you request parallel optimization. On Windows, that means screened starts use PSOCK workers in the same way as the full optimization, so the same refit_fn serialization advice applies.

set.seed(1)
opt_screen <- multiScaleR::multiScale_optim(
  fitted_mod = mod,
  kernel_inputs = kernel_inputs,
  start_strategy = "screen",
  screen_n_sigma = 4,
  screen_n_jitter = 2,
  screen_maxit = 5
)

4.5 Examine results

summary() reports the estimated σ and its standard error and 95% confidence interval for each spatial covariate, followed by the standard regression summary. By default, confidence intervals are computed using the delta method. For more reliable bounds when the likelihood surface is asymmetric (e.g., when σ is near max_D or sample size is small), pass profile = TRUE.

summary(opt)
#> 
#> Call:
#> multiScale_optim(fitted_mod = mod0_2, kernel_inputs = kernel_inputs, 
#>     verbose = FALSE)
#> 
#> 
#> Kernel used:
#> gaussian
#> 
#> ***** Optimized Scale of Effect -- Sigma *****
#> 
#>           Mean       SE     2.5%    97.5%
#> land1 238.4880 33.12452 172.7364 304.2397
#> land2 499.5767 47.47055 405.3484 593.8050
#> 
#> 
#>   ==================================== 
#> 
#> ***** Optimized Scale of Effect -- Distance *****
#> 90% Kernel Weight
#> 
#>         Mean   2.5%  97.5%
#> land1 392.28 284.13 500.43
#> land2 821.73 666.74 976.72
#> 
#> 
#>   ==================================== 
#> 
#>  *****     Fitted Model Summary     *****
#> 
#> 
#> Call:
#> glm(formula = counts ~ site + land1 + land2, family = poisson(), 
#>     data = data)
#> 
#> Coefficients:
#>             Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)  0.39291    0.09098   4.319 1.57e-05 ***
#> site         0.17482    0.08246   2.120    0.034 *  
#> land1       -0.48628    0.07923  -6.138 8.38e-10 ***
#> land2        0.62571    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.25
#> 
#> Number of Fisher Scoring iterations: 5
#> 
#> 
#> 
#> ***** Follow-up recommendation *****
#> Recommended `max_D`: 1700
#> Recommended starting `sigma` values:
#>    land1    land2 
#> 238.4880 499.5767 
#> 
#> Action:
#> Use the optimized parameters as efficient starting values if refitting.

The true simulated values (σ = 250 m for land1, σ = 500 m for land2) should fall within the reported confidence intervals.

4.6 Visualize kernel decay and covariate effects

plot() shows how the kernel weight assigned to each landscape cell decreases with distance. The vertical line marks the distance enclosing 90% of the total kernel weight (the effective scale of effect).

plot(opt)

plot_marginal_effects() shows the estimated effect of each covariate on the response, holding all other covariates at their mean values.

plot_marginal_effects(opt)

4.7 Diagnose the likelihood surface

profile_sigma() evaluates AICc (or log-likelihood) across the sigma search range for each covariate, holding all other sigmas at their optimized values. A clear trough centered near the optimized sigma confirms a well-identified scale of effect; a flat or monotone profile indicates the data do not strongly constrain sigma for that variable.

You can parallelize this step with n_cores. A practical rule is to use up to the number of scaled surfaces (spatial covariates), which is often the most efficient setting when computer resources permit.

n_profile_cores <- min(length(opt$scale_est$Mean),
                       max(1L, parallel::detectCores(logical = FALSE)))
prof <- vignette_cache$quickstart$prof
plot(prof)

For a deeper explanation of how to interpret the profile and when to prefer profile-likelihood CIs over delta-method CIs, see vignette("multiScaleR_Guide", package = "multiScaleR").

4.8 Project the fitted model to the landscape

kernel_scale.raster() applies the optimized kernel to the full raster extent and, with scale_center = TRUE, standardizes the result using the same centering and scaling used during model fitting. The output can be passed directly to terra::predict().

r_scaled <- terra::unwrap(vignette_cache$quickstart$r_scaled)
plot(r_scaled)

Kernel-smoothed raster layers at optimized scales.

pred <- terra::unwrap(vignette_cache$quickstart$pred)
plot(pred, main = "Predicted counts")
plot(surv_pts, add = T, col = 'red')

Predicted counts across the landscape.

The nature of this toy example and simulated data is such that predicted counts go to zero at low levels of land2. In this example that is most of the border region. This results is simply an artifact of the small landscape and large scale of effect.