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)
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
mtry globallyBefore 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.
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
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)
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)")
plot of chunk unnamed-chunk-8
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
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")
bw.search$tested.bandwidths to inspect
the trade-off.forests = FALSE. If you only need diagnostics and not prediction
on new points, set forests = FALSE. The output object is then much
smaller (only one ranger object is kept, the global one).geo.weighted = FALSE. Disables the bi-square case weights but
keeps the local sub-sampling. Useful for ablation studies.ranger. Call
set.seed() before any of the package functions (the package itself
does not call set.seed()).For tutorials, related publications and contact details visit the maintainer’s website at https://stamatisgeoai.eu/.
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.