Fitting models: MaxNet, GLM and random forests

The records are clean, the predictors chosen and the background drawn: now a model can learn what separates the places where the sloth was recorded from the places available to it. In this part you will fit three kinds of model to the same data (MaxNet, a logistic regression and a random forest), test each one on records it has not seen, read what each learned (response curves and the importance of the predictors), and tune the setting that decides how closely a model may follow its records. At the end you will have a model you can defend, five cross-validated scores to report with it, and a sense of how much the choice of algorithm matters (less than the data, more than nothing).

We continue with the brown-throated sloth as part 3 left it: the records cleaned, one per cell of 0.5° and thinned at 50 km (83 presences), the accessible area within 500 km of them, four predictors (bio5, bio7, bio16 and bio17) and 2 000 background points.

The idea#

What a model of presences and background learns#

Every record table of this series has a row per point: a presence (1, a record of the sloth) or a background point (0, a place available to it), and the value of every predictor there. A species distribution model is a function of the predictors that gives presences higher values than the background. It does not know where the sloth is absent: the background is a sample of what was available, and some background points may well hold sloths. The output is a relative suitability, a ranking of places, not a probability of presence (see part 1).

The three algorithms of sdm.fit learn that function in different ways.

MaxNet is Maxent (Phillips et al., 2006), the most used method for presence-only data, as rewritten in open source by Phillips et al. (2017). Maxent does not use the predictors directly. It builds features from them: transformations that let a response bend, rise and fall, or change abruptly.

Enlarge

Which feature classes MaxNet may use is a setting (classes): letters among l (linear), q (quadratic), p (product), h (hinge) and t (threshold). By default, as in maxnet, they depend on the number of presences: l below 10, lq below 15, lqh below 80 and lqph from 80 up. The output turns the model into a value from 0 to 1: cloglog (the default, recommended by Phillips et al., 2017) or logistic.

A GLM (generalised linear model) here is a logistic regression: the log-odds of a presence as a straight-line function of the predictors, standardised first, with a small ridge penalty (C, the inverse of its strength). It is the classic of species distribution modelling: smooth, easy to read, and predictable when projected. As fitted by sdm.fit it has one term per predictor, so every response is a rising or a falling S, never a hump.

A random forest (Breiman, 2001) grows hundreds of classification trees, each on a bootstrap sample of the records and trying a random subset of the predictors at each split, and averages their votes. It finds interactions and thresholds without being told, and it can memorise its training records: a tree grown to the end puts every presence in a leaf of its own. With far more background than presences it also learns mostly about the background. Valavi et al. (2021) showed that forests for presence-only data work well when the classes are balanced. sdm.fit weights the records by default (weights = "balanced": every presence 1, every background point the number of presences over the number of background points) for the GLM and the forest.

Testing on records the model has not seen#

A model always fits its own records well. The question is whether it predicts other places. k-fold cross-validation splits the records into k groups (folds); each fold in turn is set aside, the model is fitted to the others and scored on it. sdm.fit does this for every fold of the fold column, then fits the model once more on every record (the model it returns).

How the folds are made matters. Records close together share their climate, so a test record next to a training record is almost a copy of it, and random folds flatter the model. Spatial blocks (Valavi et al., 2019; Roberts et al., 2017) keep the test records away from the training ones: the records fall in blocks of a size larger than the range of autocorrelation, and whole blocks go to a fold.

The scores computed on each test fold are those of sdm.fit's metrics: the AUC on the training and the test records and their difference, the largest TSS, the Boyce index and two omission rates. Part 5 explains each in detail. For now: the test AUC is the chance that a test presence scores above a background point (0.5 is chance); TSS and Boyce say how well the suitability separates and orders the records; the omission rates are the share of test presences that fall below a threshold set on the training presences, and the AUC difference (train minus test) grows when a model overfits (Muscarella et al., 2014; Radosavljevic & Anderson, 2014).

When to use it#

• MaxNet is the default for presence-only records: it was designed for them, it has been benchmarked more than any other method, and with tuned regularisation it is among the best (Valavi et al., 2022).
• A GLM when interpretation matters more than the last point of AUC, when records are few, and when the model will be projected far from its data: its straight-line responses extrapolate in a predictable (if not always realistic) way.
• A random forest when you suspect interactions or sharp thresholds, and have enough records (a few hundred) to let it find them without memorising.
• Several, and compare. Valavi et al. (2022) found large differences between species in which method did best; an algorithm that disagrees with the others is worth a look at its response curves.

Fitting a model is not a substitute for presence-absence surveys when you need a probability of occurrence, nor for a mechanistic model when you need to know why the species is limited. And no algorithm repairs biased records or a background drawn from the wrong area (part 2, part 3).

What you need#

• A records table with the columns presence (1 or 0), the predictors and fold: the records output of sdm.prepare or of sdm.dataset (with a folds node after it). The columns x, y, presence, fold, weight and cell are never used as predictors; every other column is, unless predictors lists them.
• The predictors as a stack, to map the model: the env output of sdm.prepare, the predictors it kept.
• The data of the series: bradypus.csv and bioclim.tif.

Build the model#

Prepare, fit, predict#

The preparation of parts 2 and 3 fits in one sdm.prepare node (part 7 explains all its settings): records cleaned and one per cell, thinned at 50 km, an accessible area of 500 km, the predictors below a VIF of 10, 2 000 background points and five folds of spatial blocks of 500 km.

format = "terraflux/1"

[[nodes]]
id = "records"
op = "vector.read"
params = { path = "bradypus.csv" }

[[nodes]]
id = "climate"
op = "raster.read"
params = { path = "bioclim.tif" }

[[nodes]]
id = "data"
op = "sdm.prepare"
inputs = { occurrences = "records.vector", env = "climate.raster" }
params = { thin = 50000, area = 500000, background = 2000, folds = { size = 500000, k = 5 }, seed = 1 }

[[nodes]]
id = "fit"
op = "sdm.fit"
inputs = { records = "data.records" }
params = { model = "maxnet", regmult = 1 }

[[nodes]]
id = "predict"
op = "sdm.maxnet_predict"
inputs = { model = "fit.model", env = "data.env" }

[[nodes]]
id = "curves"
op = "sdm.response_curves"
inputs = { model = "fit.model", records = "data.records" }

[[nodes]]
id = "importance"
op = "sdm.importance"
inputs = { model = "fit.model", records = "data.records" }

From the command line, add [outputs.summary] with from = "fit.summary" (and the same for fit.metrics, curves.curves and importance.importance) and run terraflux run sloth-fit.toml --out results: each table is written as a CSV file.

The other algorithms#

The same node, other settings. The GLM:

[[nodes]]
id = "glm"
op = "sdm.fit"
inputs = { records = "data.records" }
params = { model = "glm", C = 1 }

C is the inverse of the strength of the ridge penalty (1 by default; smaller, a flatter model); standardize (true by default) scales the predictors first so that the penalty treats them alike.

A random forest:

[[nodes]]
id = "forest"
op = "sdm.fit"
inputs = { records = "data.records" }
params = { model = "random_forest", n_estimators = 500, max_depth = 3 }

n_estimators is the number of trees (100 by default); min_samples_leaf, max_depth and min_samples_split limit how finely each tree may split the records; max_features how many predictors a split tries (the square root of their number by default); random_state the seed. The forest is scikit-learn's, tree for tree.

Other folds#

sdm.prepare makes spatial blocks by default. To fit the same records over other folds, put a folds node between the preparation and the fit: it returns the records with a new fold column.

[[nodes]]
id = "random"
op = "sdm.random_folds"
inputs = { records = "data.records" }
params = { k = 5, seed = 1 }

[[nodes]]
id = "fit_random"
op = "sdm.fit"
inputs = { records = "random.records" }

[[nodes]]
id = "clusters"
op = "sdm.spatial_clusters"
inputs = { records = "data.records" }
params = { k = 5, seed = 1 }

[[nodes]]
id = "fit_clusters"
op = "sdm.fit"
inputs = { records = "clusters.records" }

sdm.spatial_blocks (size in metres on longitude and latitude records, k, hexagons or squares), sdm.quadrant_folds, sdm.checkerboard_folds (which also needs the predictors on env, for their grid) and sdm.jackknife_folds (a fold per presence, for very few records) work the same way.

Read the results#

The scores over the folds#

The metrics of the MaxNet model, a row per test fold:

Two folds (3 and 5) are predicted much worse than the others: the model fitted without their blocks does not know their conditions well. This spread is information, not noise to hide: report the standard deviation with the mean. A test AUC of 0.69 looks modest beside the 0.85 of the first model of part 1, but the task is harder: the background now comes from places the sloth could reach, most of them just as tropical (part 3 showed the effect).

The omission rate at the 10^th percentile should be near 10 % for a model that generalises; 14 % is close. The threshold output (0.61 here) is the suitability that best separates presences from background on all the records: part 5 shows how to choose one.

The three algorithms#

The forest with default settings is the textbook case of overfitting: trees grown to the end memorise every training presence (training AUC 1.0), its test AUC is the lowest of all, and the threshold set on its training presences (their minimum, 1.0) is one no test presence reaches. Limiting the trees fixes it: with leaves of 20 records the test AUC is 0.71, with a depth of 3 it is 0.68 and the omission falls further. The GLM, the simplest model, does as well as any (test AUC 0.71, Boyce 0.62). With 83 presences and four predictors there is not much more to learn than smooth trends.

What the models learned#

Enlarge

A response curve shows the suitability as one predictor varies and the others stay at their mean (fixed = "median" holds them at the median). Read them as the model's partial view of the niche:

• bio7, the annual range of temperature, matters most: suitability is high where temperature varies by about 10 °C over the year or less (the humid tropics) and falls to nothing beyond 20 °C. This is ecologically sound: a sloth cannot regulate its body temperature well.
• bio16, the rain of the wettest quarter, raises suitability up to about 1 500 mm.
• bio5 and bio17 change little. MaxNet's upturn at the right end of bio17 rests on two records: a curve where there are no records is the model's guess (look at the ticks).

sdm.importance shuffles each predictor in turn (n_repeats times, 5 by default, with a seed) and measures how much the AUC falls: a predictor the model relies on breaks it when shuffled. Its predictors output lists them, the most important first. With collinear predictors, the importance is shared between them in arbitrary ways, one more reason for the selection of part 3.

Reporting a model#

In the methods: the algorithm and its settings (for MaxNet the feature classes and the regularisation multiplier; for a forest the trees, leaf size or depth; for a GLM its terms and penalty), how the records were weighted, the folds (design, number and size), and the output (cloglog). In the results: the mean and standard deviation of the test AUC, TSS, Boyce and the omission rates over the folds (the summary table), the response curves of the main predictors, their importance, and the map. ODMAP, a standard protocol for reporting distribution models (Zurell et al., 2020), lists the rest.

Settings that change the answer#

The regularisation multiplier#

The multiplier scales MaxNet's penalty on its coefficients. Too low, and the model follows every record (and every quirk of the sampling); too high, and it flattens into a trend. The slider shows the same records with multipliers from 0.25 to 8.

Enlarge

PlayRegularisation multiplier (regmult)

The numbers tell a consistent story. As the multiplier rises from 0.25 to 4, the gap between training and test AUC shrinks, the test AUC and the Boyce index rise and the omission rate falls towards 10 %: the model stops fitting noise. At 8 the model has only 4 features left and the Boyce index drops: it has become too simple to order the places well. A multiplier of 2 to 4 is the best compromise for these data, above Maxent's default of 1, as is common with few records (Radosavljevic & Anderson, 2014).

This is what ENMeval does (Muscarella et al., 2014): fit a grid of multipliers and feature classes over spatial folds, and choose the simplest model with low omission and a small AUC difference, or the best by AICc. In Terraflux, put several sdm.fit nodes side by side, or a control.repeat loop over regmult, and compare their summary tables. Choose on the test scores, never on the map you like best.

The feature classes#

More classes give a closer fit to the training records and slightly worse test scores: with these data, linear and quadratic features are enough. Hinges earn their place with many records and responses with sharp turns.

The folds#

Here random folds score about as well as spatial blocks. The thinning at 50 km and the coarse grid have already removed most of the near-duplicates that make random folds optimistic; on fine grids with clustered records the difference is often large (Roberts et al., 2017; Valavi et al., 2019). Clusters give the most extreme folds: each holds a whole region (Central America, the Atlantic forest), and omits a third of its presences, a test of transfer to regions the model has not seen. Choose the design for the use of the model: blocks for mapping within the region, clusters or quadrants when the model will be projected elsewhere.

Common pitfalls#

• Judging a model by its training AUC. The model always fits its own records; only the test scores mean anything, and a large difference between the two is the sign of overfitting.
• A random forest with default settings. It memorises presence-only data. Limit the trees (min_samples_leaf, max_depth) and keep the records balanced (weights = "balanced", the default).
• Tuning on the map. Choose settings by their cross-validated scores, then look at the map, not the other way round.
• Comparing scores across different backgrounds or folds. Fit the models you compare on the same records table; part 3 showed how much the background moves the AUC.
• Response curves where there are no records. Read them together with the distribution of the records (the ticks of the figure); the ends of a curve are often extrapolations.
• Reading cloglog suitability as probability. It is a relative index of suitability; only presence-absence data can calibrate a probability.
• Extra columns in a records table. sdm.fit, and the lower-level sdm.maxnet given a records table on its data, leave out the columns that are not predictors (x, y, fold, weight, cell...) and use every other column. A column of your own in a table you bring (a site's code, a year) would become a predictor: keep only the predictors and those columns. sdm.maxnet fits once, on every record; sdm.fit adds the cross-validation.

Going further#

• Part 5: Evaluation and thresholds: AUC, TSS, Boyce and omission in detail, and how to turn suitability into presence and absence.
• Part 6: Projection and extrapolation: the model in another climate, and where not to trust it.
• Part 7: The whole preparation in one node.
• Suitability maps feed conservation plans: Systematic conservation planning.
• The operators: sdm.fit, sdm.maxnet, sdm.maxnet_predict, sdm.response_curves, sdm.importance.

References#

• Breiman, L. (2001). Random forests. Machine Learning, 45(1), 5–32. https://doi.org/10.1023/A:1010933404324
• Elith, J., Phillips, S. J., Hastie, T., Dudík, M., Chee, Y. E., & Yates, C. J. (2011). A statistical explanation of MaxEnt for ecologists. Diversity and Distributions, 17(1), 43–57. https://doi.org/10.1111/j.1472-4642.2010.00725.x
• Fithian, W., & Hastie, T. (2013). Finite-sample equivalence in statistical models for presence-only data. The Annals of Applied Statistics, 7(4), 1917–1939. https://doi.org/10.1214/13-AOAS667
• Muscarella, R., Galante, P. J., Soley-Guardia, M., Boria, R. A., Kass, J. M., Uriarte, M., & Anderson, R. P. (2014). ENMeval: an R package for conducting spatially independent evaluations and estimating optimal model complexity for Maxent ecological niche models. Methods in Ecology and Evolution, 5(11), 1198–1205. https://doi.org/10.1111/2041-210X.12261
• Phillips, S. J., Anderson, R. P., Dudík, M., Schapire, R. E., & Blair, M. E. (2017). Opening the black box: an open-source release of Maxent. Ecography, 40(7), 887–893. https://doi.org/10.1111/ecog.03049
• Phillips, S. J., Anderson, R. P., & Schapire, R. E. (2006). Maximum entropy modeling of species geographic distributions. Ecological Modelling, 190(3–4), 231–259. https://doi.org/10.1016/j.ecolmodel.2005.03.026
• Radosavljevic, A., & Anderson, R. P. (2014). Making better Maxent models of species distributions: complexity, overfitting and evaluation. Journal of Biogeography, 41(4), 629–643. https://doi.org/10.1111/jbi.12227
• Roberts, D. R., Bahn, V., Ciuti, S., Boyce, M. S., Elith, J., Guillera-Arroita, G., et al. (2017). Cross-validation strategies for data with temporal, spatial, hierarchical, or phylogenetic structure. Ecography, 40(8), 913–929. https://doi.org/10.1111/ecog.02881
• Valavi, R., Elith, J., Lahoz-Monfort, J. J., & Guillera-Arroita, G. (2019). blockCV: an R package for generating spatially or environmentally separated folds for k-fold cross-validation of species distribution models. Methods in Ecology and Evolution, 10(2), 225–232. https://doi.org/10.1111/2041-210X.13107
• Valavi, R., Elith, J., Lahoz-Monfort, J. J., & Guillera-Arroita, G. (2021). Modelling species presence-only data with random forests. Ecography, 44(12), 1731–1742. https://doi.org/10.1111/ecog.05615
• Valavi, R., Guillera-Arroita, G., Lahoz-Monfort, J. J., & Elith, J. (2022). Predictive performance of presence-only species distribution models: a benchmark study with reproducible code. Ecological Monographs, 92(1), e01486. https://doi.org/10.1002/ecm.1486
• Warton, D. I., & Shepherd, L. C. (2010). Poisson point process models solve the "pseudo-absence problem" for presence-only data in ecology. The Annals of Applied Statistics, 4(3), 1383–1402. https://doi.org/10.1214/10-AOAS331
• Zurell, D., Franklin, J., König, C., Bouchet, P. J., Dormann, C. F., Elith, J., et al. (2020). A standard protocol for reporting species distribution models. Ecography, 43(9), 1261–1277. https://doi.org/10.1111/ecog.04960