Geographically Weighted Random Forest with SpatialML

Stamatis Kalogirou and Stefanos Georganos

Ιούλιος 06, 2026

1. Introduction

Geographically Weighted Random Forest (GRF) is a spatial extension of the Random Forest algorithm. It was introduced by Georganos et al. (2019) and refined in Georganos and Kalogirou (2022). For each observation i the algorithm fits a local Random Forest using only the observations that fall in the neighbourhood of i (defined by the k nearest neighbours when the kernel is adaptive, or by all the observations within a distance bw when the kernel is fixed). When geo.weighted = TRUE, observations within the neighbourhood are weighted by the bi-square kernel $$ w_{ij} ;=; \left(1 - (d_{ij}/h)^{2}\right)^{2}, $$ where \(d_{ij}\) is the Euclidean distance between observations i and j and \(h\) is the largest distance retained in the neighbourhood. The final model is the collection of all local random forests plus a global random forest fitted on the whole sample.

The package exports four user-facing functions:

Function Purpose
grf() Fit a Geographically Weighted Random Forest
grf.bw() Search the optimal bandwidth
predict.grf() Predict at new spatial locations (via S3 predict())
rf.mtry.optim() Tune the global mtry parameter (OOB or k-fold CV)
random.test.data() Generate small synthetic spatial data for testing

The package uses ranger as its random-forest back-end. Undefined local out-of-bag predictions are handled with a quiet leave-one-out fallback.

library(SpatialML)

2. Quick start with synthetic data

We start with a small synthetic dataset created on a regular 8 x 8 grid. random.test.data() returns one dependent variable (dep), two random predictors (X1, X2) and the grid coordinates (X, Y).

set.seed(42)
df <- random.test.data(nrows = 8, ncols = 8, vars.no = 3)
head(df)
#>          dep        X1         X2 X Y
#> 1  1.3709584 0.2335235 0.12887216 1 1
#> 2 -0.5646982 0.7244976 0.12908928 1 2
#> 3  0.3631284 0.9036345 0.07225311 1 3
#> 4  0.6328626 0.6034741 0.05312948 1 4
#> 5  0.4042683 0.6315073 0.53187444 1 5
#> 6 -0.1061245 0.9373858 0.11230824 1 6

2.1 Tuning mtry globally

Before fitting the GRF we tune the mtry parameter on the global data. Out-of-bag error (the default) is fast and statistically valid for Random Forests.

set.seed(1)
mtry.opt <- rf.mtry.optim(dep ~ X1 + X2, dataset = df,
                          cv.method = "oob", plot.it = FALSE,
                          verbose = FALSE)
mtry.opt$best.mtry
#> [1] 1
mtry.opt$results
#>   mtry     RMSE    Rsquared SDRMSE SDRsq
#> 1    1 1.176771 -0.08096882     NA    NA
#> 2    2 1.214827 -0.15201429     NA    NA

Use cv.method = "repeatedcv" for a more rigorous evaluation when the sample is small.

2.2 Optimal bandwidth

grf.bw() evaluates a grid of candidate bandwidths and returns the one that maximises the local OOB R-squared.

set.seed(1)
bw.search <- grf.bw(dep ~ X1 + X2, dataset = df, kernel = "adaptive",
                    coords = df[, c("X", "Y")],
                    bw.min = 8, bw.max = 18, step = 2,
                    trees = 200, mtry = mtry.opt$best.mtry,
                    verbose = FALSE)
bw.search$tested.bandwidths
#>   Bandwidth      Local       Mixed   Low.Local
#> 1         8 -0.4114707 0.002040161 0.009830132
#> 2        10 -0.4996148 0.001290432 0.011996755
#> 3        12 -0.2739078 0.024000771 0.023387175
#> 4        14 -0.2203878 0.028195146 0.027563047
#> 5        16 -0.2343679 0.015383208 0.020700446
#> 6        18 -0.1874398 0.025608808 0.032172885
bw.search$Best.BW
#> [1] 18

2.3 Fitting the GRF

With both mtry and the bandwidth chosen we fit the final model. The forests = TRUE argument is required if you want to call predict() later on new data.

set.seed(1)
m <- grf(dep ~ X1 + X2, dframe = df, bw = bw.search$Best.BW,
         kernel = "adaptive", coords = df[, c("X", "Y")],
         ntree = 200, mtry = mtry.opt$best.mtry,
         forests = TRUE, print.results = FALSE, progress = FALSE)

2.4 Inspecting the fit

The fitted object is an S3 object of class "grf". Useful slots:

class(m)
#> [1] "grf"
m$LocalModelSummary$l.r.OOB              # local OOB R-squared
#> [1] -0.3795875
summary(m$LGofFit$LM_ResOOB)             # local OOB residuals
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#> -4.3458 -0.5309  0.1373  0.1017  1.0240  3.1554
head(m$Local.Variable.Importance)        # local importance per observation
#>          X1        X2
#> 1 12.175976 15.507384
#> 2 13.318965 13.472106
#> 3 10.623393 10.898932
#> 4 10.081469  8.279829
#> 5  8.818826  7.961262
#> 6  6.776464  6.550084

The columns of m$Local.Variable.Importance are the predictors and the rows correspond, in order, to the observations of dframe. To explore the spatial pattern you can map them with any plotting framework you like:

imp <- m$Local.Variable.Importance$X1
plot(df$X, df$Y, pch = 19, cex = 2,
     col = grDevices::grey(1 - imp / max(imp)),
     xlab = "X", ylab = "Y",
     main = "Local importance of X1 (darker = more important)")
Local importance of predictor X1 across the synthetic 8 x 8 grid.

plot of chunk unnamed-chunk-8

2.5 Predicting at new locations

predict() dispatches to predict.grf() because m has class "grf". For each new observation the local random forest fitted at the geographically nearest training location is used.

new.df <- random.test.data(2, 2, vars.no = 3)
predict(m, new.df, x.var.name = "X", y.var.name = "Y")
#> [1] -0.2474032 -0.0510753  1.0622908  0.5590549

By default the global random forest receives weight zero (global.w = 0). Setting global.w and local.w to non-zero values returns a linear blend of the two predictions and is a useful sensitivity test when the local model is noisy.

predict(m, new.df, x.var.name = "X", y.var.name = "Y",
        local.w = 0.5, global.w = 0.5)
#> [1] -0.2635139 -0.4262487  0.7139588  0.5709697

3. Real-world example: Greek municipal income

The Income dataset (shipped with the package) contains 325 municipalities of Greece with their centroid coordinates and four socio-economic variables. Fitting a full GRF on this data is heavier than the toy example above; the chunk below is therefore not evaluated inside the vignette but copy-paste it into an interactive session.

data(Income)
Coords <- Income[, 1:2]

# 1. Search the optimal bandwidth (be patient)
bw <- grf.bw(Income01 ~ UnemrT01 + PrSect01,
             dataset = Income, kernel = "adaptive",
             coords = Coords, bw.min = 30, bw.max = 80, step = 5)

# 2. Fit the final model
m <- grf(Income01 ~ UnemrT01 + PrSect01, dframe = Income,
         bw = bw$Best.BW, kernel = "adaptive", coords = Coords)

# 3. Local R-squared
m$LocalModelSummary$l.r.OOB

# 4. Map the residuals
plot(Coords[, 1], Coords[, 2], pch = 19,
     col = ifelse(m$LGofFit$LM_ResOOB > 0, "red", "blue"),
     xlab = "X", ylab = "Y", main = "GRF OOB residuals")

4. Practical tips

For tutorials, related publications and contact details visit the maintainer’s website at https://stamatisgeoai.eu/.

References

Georganos, S., Grippa, T., Niang Gadiaga, A., Linard, C., Lennert, M., Vanhuysse, S., Mboga, N., Wolff, E. and Kalogirou, S. (2019) Geographical Random Forests: A Spatial Extension of the Random Forest Algorithm to Address Spatial Heterogeneity in Remote Sensing and Population Modelling. Geocarto International. DOI: 10.1080/10106049.2019.1595177.

Georganos, S. and Kalogirou, S. (2022) A Forest of Forests: A Spatially Weighted and Computationally Efficient Formulation of Geographical Random Forests. ISPRS International Journal of Geo-Information, 11(9), 471. DOI: 10.3390/ijgi11090471.

Wright, M. N. and Ziegler, A. (2017) ranger: A Fast Implementation of Random Forests for High Dimensional Data in C++ and R. Journal of Statistical Software, 77(1), 1-17. DOI: 10.18637/jss.v077.i01.