How should ICE slopes be computed for local marginal effects in spatial analysis?

I am replicating a methodology that uses Individual Conditional Expectation (ICE) plots to explore heterogeneity in predictor effects across neighbourhoods (UK LSOA). The paper states

> In this study, ICE was applied to each demographic indicator to assess whether its association with low Mapillary coverage was consistent in direction but varied in intensity, or whether it exhibited directional reversals across neighbourhoods. To summarise local effects, the slope maps of ICE curves at the LSOA centroids were computed, which highlight the spatial distribution of marginal effects.

Using the `iml` package, I currently fit a linear regression across each ICE curve (sampled at 5 grid points, for simplicity) and take the slope coefficient. This gives me an average slope over the feature’s domain.

Does slopes at the LSOA centroids mean I should instead compute the local derivative of the ICE curve at the observed feature value for each neighborhood, rather than the average slope across the whole curve?

Here is a simplified version of my code:

```
library(randomForest)
library(iml)
library(dplyr)

# Dummy dataset
data(mtcars)

# Fit a random forest
rf_fit2 <- randomForest(mpg ~ ., data = mtcars)

# Wrap in iml Predictor
predictor2 <- Predictor$new(
model = rf_fit2,
data = mtcars[, -1], # predictors only
y = mtcars$mpg
)

# Function to compute average slope across ICE curve
compute_ice_slope2 <- function(feature_name, predictor){
ice_obj <- FeatureEffect$new(
predictor,
feature = feature_name,
method = "ice",
grid.size = 5 # for simplicity
)

ice_obj$results |>
group_by(.id) |>
summarise(
slope = coef(lm(.value ~ .data[[feature_name]]))[2]
)
}

# Example: slope for 'hp'
slopes_hp2 <- compute_ice_slope2("hp", predictor2)
head(slopes_hp2)

# Plot ICE curves
plot(ice_obj)
```
Session Info:

```
> sessionInfo()
R version 4.5.2 (2025-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 26200)

Matrix products: default
LAPACK version 3.12.1

locale:
[1] LC_COLLATE=English_United States.utf8 LC_CTYPE=English_United States.utf8 LC_MONETARY=English_United States.utf8
[4] LC_NUMERIC=C LC_TIME=English_United States.utf8

time zone: Europe/Budapest
tzcode source: internal

attached base packages:
[1] parallel stats graphics grDevices utils datasets methods base

other attached packages:
[1] randomForest_4.7-1.2 iml_0.11.4 GPfit_1.0-9 janitor_2.2.1 lubridate_1.9.5 forcats_1.0.1
[7] stringr_1.6.0 readr_2.2.0 tidyverse_2.0.0 patchwork_1.3.2 reshape2_1.4.5 treeshap_0.4.0
[13] future_1.69.0 fastshap_0.1.1 shapviz_0.10.3 kernelshap_0.9.1 tibble_3.3.1 doParallel_1.0.17
[19] iterators_1.0.14 foreach_1.5.2 ranger_0.18.0 yardstick_1.3.2 workflowsets_1.1.1 workflows_1.3.0
[25] tune_2.0.1 tidyr_1.3.2 tailor_0.1.0 rsample_1.3.2 recipes_1.3.1 purrr_1.2.1
[31] parsnip_1.4.1 modeldata_1.5.1 infer_1.1.0 ggplot2_4.0.2 dplyr_1.2.0 dials_1.4.2
[37] scales_1.4.0 broom_1.0.12 tidymodels_1.4.1

loaded via a namespace (and not attached):
[1] RColorBrewer_1.1-3 rstudioapi_0.18.0 jsonlite_2.0.0 magrittr_2.0.4 farver_2.1.2 fs_1.6.6
[7] ragg_1.5.0 vctrs_0.7.1 memoise_2.0.1 sparsevctrs_0.3.6 usethis_3.2.1 curl_7.0.0
[13] xgboost_3.2.0.1 parallelly_1.46.1 KernSmooth_2.23-26 desc_1.4.3 plyr_1.8.9 cachem_1.1.0
[19] lifecycle_1.0.5 pkgconfig_2.0.3 Matrix_1.7-4 R6_2.6.1 fastmap_1.2.0 snakecase_0.11.1
[25] digest_0.6.39 furrr_0.3.1 ps_1.9.1 pkgload_1.5.0 textshaping_1.0.4 labeling_0.4.3
[31] timechange_0.4.0 compiler_4.5.2 proxy_0.4-29 remotes_2.5.0 withr_3.0.2 S7_0.2.1
[37] backports_1.5.0 DBI_1.2.3 pkgbuild_1.4.8 MASS_7.3-65 lava_1.8.2 sessioninfo_1.2.3
[43] classInt_0.4-11 tools_4.5.2 units_1.0-0 future.apply_1.20.2 nnet_7.3-20 Metrics_0.1.4
[49] doFuture_1.2.1 glue_1.8.0 callr_3.7.6 grid_4.5.2 sf_1.0-24 checkmate_2.3.4
[55] generics_0.1.4 gtable_0.3.6 tzdb_0.5.0 class_7.3-23 data.table_1.18.2.1 hms_1.1.4
[61] utf8_1.2.6 pillar_1.11.1 splines_4.5.2 lhs_1.2.0 lattice_0.22-9 sfd_0.1.0
[67] survival_3.8-6 tidyselect_1.2.1 hardhat_1.4.2 devtools_2.4.6 timeDate_4052.112 stringi_1.8.7
[73] DiceDesign_1.10 pacman_0.5.1 codetools_0.2-20 cli_3.6.5 rpart_4.1.24 systemfonts_1.3.1
[79] processx_3.8.6 dichromat_2.0-0.1 Rcpp_1.1.1 globals_0.19.0 ellipsis_0.3.2 gower_1.0.2
[85] listenv_0.10.0 viridisLite_0.4.3 ipred_0.9-15 prodlim_2025.04.28 e1071_1.7-17
```

I am not sure it's correct because “at the LSOA centroids” is almost certainly spatial language: I have one ICE curve per LSOA, and I map a scalar summary (here: a slope / marginal effect) to the centroid location of that LSOA.

Because each ICE curve is f(x_j, x_{-j,i}) as a function of feature x_j while holding the other features fixed at LSOA i’s values, a “slope map” usually aims to show a local marginal effect for each LSOA, that is the derivative (or a local finite-difference approximation) at the observed feature value for that LSOA, then plotted at that LSOA’s centroid.

I think my current approach (a linear regression over the whole ICE curve) is an average slope over the feature’s domain, which is a different estimand and can hide sign changes or curvature (exactly the thing ICE is meant to reveal).

How can I create the "slope map"?

I apologize if the quesiton is not suitable for GitHub. Feel free to close it.