spatialsample:

A tidy approach to spatial cross-validation

Michael J Mahoney

mjmahone@esf.edu

About Me

• Mike Mahoney

• PhD candidate in environmental science

• 2022 summer intern with Posit (spatialsample, rsample)

• These slides: mm218.dev/boston_useR_2023

Before I get into that, I want to introduce myself quickly – I’m Mike and I’m currently a PhD candidate in environmental science, but the work I’m going to talk about today actually comes from this past summer when I was an intern with the tidymodels team at Posit working on cross validation in tidymodels.

Data splitting:

Cross-validation:

And for anyone who isn’t familiar, cross validation is a model assessment approach that helps you understand how well your model will predict new data that it hasn’t seen before. Normally you’ll evaluate your final models against an independent test set, which isn’t involved in the model fitting process at all. But if you want to evaluate models that you’re still iterating on, you’ll normally split your training data into a handful of what we call “folds”. You then usually train your model on all but one of those folds, and evaluate it against that final fold, in order to get a sense of how well your intermediate models will do on independent data. Then, when you’re finished iterating, you’ll train your model on all those folds together and evaluate it against the independent test set.

rsample and friends

library(tidymodels)
rsample::vfold_cv(spatialsample::boston_canopy) |> head()

#> # A tibble: 6 × 2
#>   splits           id    
#>   <list>           <chr> 
#> 1 <split [613/69]> Fold01
#> 2 <split [613/69]> Fold02
#> 3 <split [614/68]> Fold03
#> 4 <split [614/68]> Fold04
#> 5 <split [614/68]> Fold05
#> 6 <split [614/68]> Fold06

workflow() |> 
  add_model(linear_reg()) |> 
  add_formula(canopy_area_2019 ~ land_area * mean_temp) |> 
  fit_resamples(vfold_cv(spatialsample::boston_canopy)) |> 
  collect_metrics()

#> # A tibble: 2 × 6
#>   .metric .estimator       mean     n    std_err .config             
#>   <chr>   <chr>           <dbl> <int>      <dbl> <chr>               
#> 1 rmse    standard   377089.       10 20426.     Preprocessor1_Model1
#> 2 rsq     standard        0.353    10     0.0178 Preprocessor1_Model1

And there’s this tidymodels package called rsample that handles that splitting process for you. And the objects that rsample functions return can be used with functions from across the tidymodels ecosystem, which makes it easy to use cross-validation as part of any modeling workflow you might want.

What does “new data” mean?

ggplot(spatialsample::boston_canopy, aes(fill = canopy_area_2019)) + geom_sf() + 
  scale_fill_distiller(name = "Canopy area (2019)", palette = "YlGn", direction = 1)

But a challenge is that most rsample functions assume you can assign your data to folds at random, and still get independent training and test sets. That’s true for a lot of data, but starts to break down once your data is no longer entirely independent – once you start getting autocorrelation between your observations. And for spatial data, you almost always have autocorrelation. For instance, this is a map of tree cover in Boston, and as you can kinda guess tree cover isn’t randomly distributed across the city – down in Stony Brook Park there’s a ton of tree cover, up here in Seaport there’s a bit less. The tree cover for any one of these hexagons is pretty tightly linked to how much tree cover its neighbors have.

Are these folds really unrelated?

rsample::vfold_cv(spatialsample::boston_canopy, v = 5)

And so if you split this up at random, you’re going to get pretty highly related training and testing sets, which is probably going to make your model assessments way too optimistic. Your training and testing sets are so similar that you’re practically testing with training data.

So this is where spatialsample comes in. spatialsample is a newer tidymodels package that helps you do spatial cross-validation, where you assign data to folds based on its spatial location rather than just at random.

Spatial clustering

library(spatialsample)
set.seed(1234)
spatial_clustering_cv(boston_canopy, v = 5)

So to give you a sense of what that looks like, this is the output from the spatial_clustering_cv() function. You can see that rather than having our folds all mixed in together, we’re assigning folds based on where each of these hexagons are located.

Spatial clustering

library(purrr)
walk(spatial_clustering_cv(boston_canopy, v = 5)$splits, function(x) print(autoplot(x)))

And so rather than training and testing with random data, we’re testing on data that’s geographically separated from our training data, and so hopefully less related to the testing data.

Spatial blocking

spatial_block_cv(boston_canopy, v = 5, n = c(10, 10))

There’s a handful of other methods in spatialsample that are also useful; for instance, spatial_block_cv() lets you split your data up using a regular grid, which is a super popular method in ecology studies.

Spatial blocking

walk(spatial_block_cv(boston_canopy, v = 5, n = c(10, 10))$splits,
     function(x) print(autoplot(x)))

Spatial blocking

spatial_block_cv(boston_canopy, v = 5, n = c(10, 10), 
                 method = "continuous", relevant_only = FALSE)

Spatial LODO

folds <- spatial_buffer_vfold_cv(boston_canopy, v = Inf, radius = 1500, buffer = 1500)
walk(folds$splits, function(x) print(autoplot(x)))

And we’ve also got a method we call “leave one disc out”, or LODO, where you use all the observations within a certain distance of some point as your test set, and leave all the points within a “buffer” distance of that out entirely. And that buffer is really useful to make sure that you aren’t including correlated points in both your training and testing data.

Buffering

spatial_clustering_cv(boston_canopy, v = 5, buffer = 1500)

And so rather than training and testing with random data, we’re testing on data that’s geographically separated from our training data, and so hopefully less related to the testing data.

tidymodels integration

workflow() |> 
  add_model(linear_reg()) |> 
  add_formula(canopy_area_2019 ~ land_area * mean_temp) |> 
  fit_resamples(vfold_cv(spatialsample::boston_canopy)) |> 
  collect_metrics()

#> # A tibble: 2 × 6
#>   .metric .estimator       mean     n    std_err .config             
#>   <chr>   <chr>           <dbl> <int>      <dbl> <chr>               
#> 1 rmse    standard   378993.       10 18001.     Preprocessor1_Model1
#> 2 rsq     standard        0.354    10     0.0171 Preprocessor1_Model1

workflow() |> 
  add_model(linear_reg()) |> 
  add_formula(canopy_area_2019 ~ land_area * mean_temp) |> 
  fit_resamples(spatial_clustering_cv(spatialsample::boston_canopy)) |> 
  collect_metrics()

#> # A tibble: 2 × 6
#>   .metric .estimator       mean     n     std_err .config             
#>   <chr>   <chr>           <dbl> <int>       <dbl> <chr>               
#> 1 rmse    standard   934397.       10 549480.     Preprocessor1_Model1
#> 2 rsq     standard        0.348    10      0.0468 Preprocessor1_Model1

Mahoney, MJ, Johnson, L. K., Silge, J., Frick, H., Kuhn, M., and Beier, C. M. In Review. Assessing the performance of spatial cross-validation approaches for models of spatially structured data. https://doi.org/10.48550/arXiv.2303.07334

And these methods are all typically going to give you a more accurate model assessment when you’re working with spatial data than non-spatial cross-validation. This graph is from a preprint we put out last month – that green zone is the “target” range for RMSE, and you can see that clustered, LODO, and blocked all have more of their distribution in the green zone than normal V-fold cross-validation.

Other features:

Works with projected & geographic CRS

Arguments accept explicit units

Aware of CRS units, functions do unit conversion

Handles all geometry types\(^*\)

And last but not least, I wanted to mention that spatialsample is really designed around the idea of helping users fall into what you’d call a pit of success, and so handles a lot of the common edge cases that come up with spatial data – spatialsample can work with geographic coordinates, with mismatched coordinate reference systems, with different units, with points, polygons, and lines, and generally does what you’d hope spatial software would do. And by building on top of rsample, the objects and functions from spatialsample automatically integrate with the rest of the tidymodels ecosystem too.

Thank you!

Find me online:

mm218.dev

@mikemahoney218

@MikeMahoney218@fosstodon.org

Slides available at mm218.dev/boston_useR_2023

More spatialsample: https://spatialsample.tidymodels.org/

So, with that, I want to say thanks! If you have any questions or just want to talk, feel free to say hi after; otherwise, you can find me online at MikeMahoney218 or my website, mm218 dot dev. Thanks!