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)

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)
#> 
#> 
#> Optimization complete

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:
#> multiScaleR::multiScale_optim(fitted_mod = mod, kernel_inputs = kernel_inputs)
#> 
#> 
#> 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 = data)
#> 
#> 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 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.

prof <- profile_sigma(opt, n_pts = 15, verbose = FALSE)
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 <- kernel_scale.raster(raster_stack = land_rast,
                                multiScaleR  = opt,
                                scale_center = TRUE,
                                clamp        = TRUE
                                )
#> 
#> Smoothing spatRaster 1 of 2: land1 at sigma = 238
#> 
#> Smoothing spatRaster 2 of 2: land2 at sigma = 499
plot(r_scaled)

Kernel-smoothed raster layers at optimized scales.

pred <- terra::predict(r_scaled, opt$opt_mod, type = "response")
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.