Evaluation and thresholds

A model gives every place a suitability; before anyone uses the map, two questions need answers. How well does it separate the places where the species is from the others? And where, if a yes-or-no map is needed, should the line between suitable and unsuitable be drawn? In this part you will compute the usual scores for the sloth's model by hand in a graph (AUC, TSS, the Boyce index and omission rates), see what each measures and what it hides, and turn the suitability into presence and absence with eight different thresholds, which give suitable areas from 6 % to 96 % of the accessible area. Last, you will score the model on its out-of-fold predictions, the predictions of each record by a model that did not see it, and use them for calibration, thresholds and a map of residuals.

We use the model part 4 arrived at: MaxNet with linear, quadratic and hinge features (classes = "lqh") and a regularisation multiplier of 2, fitted to the 83 presences and 2 000 background points of the sloth inside its accessible area of 500 km.

The idea#

Two kinds of error#

A threshold \(t\) splits the suitability into "predicted present" (\(\geq t\)) and "predicted absent". Each record then falls in one of four boxes:

The sensitivity is the share of presences predicted present, \(a / (a + c)\); the specificity the share of background points predicted absent, \(d / (b + d)\). Every threshold trades one for the other: a low threshold catches every presence and calls most of the map suitable; a high one calls little suitable and misses presences.

With presence-only data, a "false positive" is a background point predicted suitable, and the species may well live there: the background is a sample of what is available, not of absences. Specificity is then the share of the available area predicted unsuitable. This is why the scores of presence-background models depend on the background (part 3), and why none of them can reach its maximum.

AUC: ranking presences above the background#

The receiver operating characteristic (ROC) curve plots the sensitivity against 1 − specificity for every threshold. The area under the curve (AUC) is the probability that a random presence gets a higher suitability than a random background point (Hanley & McNeil, 1982): 1 for a perfect ranking, 0.5 for a random one. It needs no threshold, which made it the most reported score in species distribution modelling (Fielding & Bell, 1997).

Its limits are well known (Lobo et al., 2008; Jiménez-Valverde, 2012). It grows with the extent of the background (the farther the background reaches from the species, the easier the ranking); it weighs parts of the curve no one uses; and with background points the best achievable AUC is below 1 and unknown. Compare AUCs only between models fitted and tested on the same records.

TSS: the best threshold's skill#

The true skill statistic (Allouche, Tsoar & Kadmon, 2006) is sensitivity + specificity − 1 at a threshold: 0 for a model no better than chance, 1 for a perfect one. Unlike Cohen's kappa, it does not depend on prevalence (the share of presences among the records), which in a presence-background data set is set by how many background points you drew. sdm.fit and sdm.evaluate report the largest TSS over all thresholds.

The Boyce index: are high values where the presences are?#

The continuous Boyce index (Boyce et al., 2002; Hirzel et al., 2006) was designed for presence-only data. It asks whether the presences are more frequent than expected in places of high suitability, and less in places of low.

The Boyce index uses only the presences and the available area, so it suits presence-only data, and the shape of the \(F\) curve tells more than the number: a curve that rises only at the top says the model separates the best places and not the rest.

Omission rates: did the model memorise its records?#

Set a threshold on the training presences and count the test presences below it (Radosavljevic & Anderson, 2014; Muscarella et al., 2014). With the lowest training suitability (the minimum training presence, MTP), a model that generalises omits almost no test presences; with the 10^th percentile of the training presences, about 10 %. Much more means the model fits its training records more tightly than new ones. sdm.fit reports both per fold (or_mtp, or_10p), and sdm.omission_rate computes them from any two sets of predictions.

Thresholds#

Many uses need a map of presence and absence: a range size, an overlap with protected areas, a habitat map for a connectivity analysis. A threshold draws the line. Liu et al. (2005) compared twelve ways of choosing one and found those that balance sensitivity and specificity among the best; for presence-only data, maximising sensitivity plus specificity (the TSS) remains a sound default (Liu, White & Newell, 2013). sdm.threshold computes them all:

When to use it#

• Every model needs an evaluation, on records it was not fitted to: the folds of sdm.fit (part 4), or independent data (a later survey, another region) through the nodes of this part.
• Compare models fitted to the same records: algorithms, settings, predictors. Use several scores at once.
• Thresholds only when the use needs a binary map. Keep the continuous suitability wherever the analysis can take it: prioritization (Zonation and linear programming take amounts), connectivity (resistance from suitability), and any comparison of maps.

What the scores are not: proof that the map is right. A model with a fine AUC can be biased in the same way as its records (records along roads, background along roads); only independent, well-designed data test that. And a score of a projected model (another time, another region) says nothing about the projection (part 6).

What you need#

• A fitted model (sdm.fit or sdm.maxnet) and the table of records it was fitted to.
• The predictions at the presences and at the background as two tables, or as maps. sdm.evaluate, sdm.threshold, sdm.boyce and sdm.omission_rate take any predictions, so they evaluate models from any software: give them a CSV of predictions read with vector.read or two maps.

Build the model#

Predictions at the presences and at the background#

The model of part 4 gives a map; the evaluation needs the predictions at the records. sdm.prepare passes on the presences it kept (presences) and the background points (background_points); sdm.dataset reads the predictors at each, and sdm.maxnet_predict, given a table instead of a stack, predicts every row.

[[nodes]]
id = "data"
op = "sdm.prepare"
inputs = { occurrences = "records.vector", env = "climate.stack" }
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 = { regmult = 2, classes = "lqh" }

[[nodes]]
id = "at_presences"
op = "sdm.dataset"
inputs = { presences = "data.presences", env = "data.env" }

[[nodes]]
id = "at_background"
op = "sdm.dataset"
inputs = { presences = "data.background_points", env = "data.env" }

[[nodes]]
id = "p_pres"
op = "sdm.maxnet_predict"
inputs = { model = "fit.model", env = "at_presences.records" }

[[nodes]]
id = "p_back"
op = "sdm.maxnet_predict"
inputs = { model = "fit.model", env = "at_background.records" }

[[nodes]]
id = "evaluation"
op = "sdm.evaluate"
inputs = { presences = "p_pres.prediction", background = "p_back.prediction" }

For independent test data, give node 3 a table of the test records (points read with vector.read) and the background of the test region to node 4: the same nodes give the test scores.

Out-of-fold predictions#

The predictions above are the model's on the records it was fitted to. sdm.fit also keeps the predictions its cross-validation made: for every record, the suitability given by the model fitted without that record's fold, on its predictions output. They are called out-of-fold predictions. With the five folds of spatial blocks of part 4, each record's prediction comes from a model fitted on the other four folds, which never saw its block.

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 = 2, classes = "lqh" }

[[nodes]]
id = "oof_pres"
op = "table.filter"
inputs = { table = "fit.predictions" }
params = { expression = "presence == 1" }

[[nodes]]
id = "oof_back"
op = "table.filter"
inputs = { table = "fit.predictions" }
params = { expression = "presence == 0" }

[[nodes]]
id = "oof_evaluation"
op = "sdm.evaluate"
inputs = { presences = "oof_pres.table", background = "oof_back.table" }
params = { column = "prediction" }

[[nodes]]
id = "oof_threshold"
op = "sdm.threshold"
inputs = { presences = "oof_pres.table", background = "oof_back.table" }
params = { column = "prediction", method = "p10" }

The first rows of the table, three presences and two background points:

To use them elsewhere (a calibration, an ensemble learned from several models' predictions, a model that combines them with an expert's range map as Merow et al. (2017) did, an analysis in R or Python), add [outputs.predictions] with from = "fit.predictions" and run the model from the command line: the table is written as predictions.csv. In the application, click the predictions port of fit to see it in the Viewer.

Thresholds and the binary map#

[[nodes]]
id = "threshold"
op = "sdm.threshold"
inputs = { presences = "p_pres.prediction", background = "p_back.prediction" }
params = { method = "max_spec_sens" }

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

[[nodes]]
id = "binary"
op = "raster.calc"
inputs = { s = "predict.prediction", t = "threshold.threshold" }
params = { expression = "1 if s >= t else 0" }

sdm.threshold returns a number on its threshold port. raster.calc takes numbers as inputs as well as maps: connect the suitability map to s and the threshold to t, and the expression 1 if s >= t else 0 makes the binary map. sdm.fit also has a threshold output: the max_spec_sens threshold on all its records.

Read the results#

The scores#

The stats of sdm.evaluate, on all the records the model was fitted to, beside the cross-validated scores of sdm.fit's summary:

The scores on all records are training scores: the model has seen these records. They are what the thresholds are set on, but the fold scores are the ones to report as the model's performance. Two things stand out. The Boyce index is high on all the records (0.92) and noisy over the folds (0.55, with a standard deviation of 0.31): each test fold holds only 14 to 20 presences, too few for a smooth \(F\) curve. And the AUC is modest, because the background is close to the records (part 3); the omission rates near 0 and 10 % say the model is not overfitted.

The ROC curve and the trade-off of thresholds#

Enlarge

The Boyce curve#

A curve that rises steadily is what a good model of habitat suitability gives. A curve that is flat and then jumps, or that falls at the top, would say that the model separates only part of the range of suitability, or that its highest values are in places the records avoid.

Thresholds and the maps they draw#

The eight methods of sdm.threshold on the same predictions:

Enlarge

Threshold method

The suitable area changes fifteen-fold between methods, from the same model. Which one?

• For a range map (where could the species be?), keep most presences: p10 or a sensitivity of 0.9 to 0.95. A sloth recorded in a place is there, and a map that leaves out half the records misses known populations.
• For the core of the range (where conditions are best, for instance for a reintroduction), max_spec_sens. It is the best-skilled threshold, and with presence-background data it is stricter than it looks: it calls 26 % of the area suitable and keeps only 53 % of the presences, because "specificity" here rewards calling the available area unsuitable, whether or not the sloth lives there.
• max_kappa and equal_prevalence depend on the prevalence, which with background points is your choice of how many points to draw (83 presences against 2 000 points: 4 %). They make little sense for presence-background data: here one calls 6 % of the area suitable and the other 96 %.
• no_omission and mtp let the single least suitable record decide, often a doubtful one.

What the out-of-fold predictions say#

Scores on what the model did not see, with every presence at once. The scores of sdm.fit's metrics are computed fold by fold, each on 14 to 20 presences. Pooling the out-of-fold predictions of all five folds scores all 83 presences, each predicted by a model that did not see it (node 7):

The pooled AUC and TSS are a little below the training ones, as test scores should be. The pooled Boyce index (0.93) is far above the mean of the folds (0.55): with all the presences in one curve, the predicted-to-expected ratio is no longer noisy, and it says the model orders unseen places well. Keep in mind what pooling does: it puts side by side the predictions of five different models, so a fold whose model predicts higher everywhere shifts the pooled scores. Report the fold scores, and the pooled ones beside them.

Calibration. A model is well calibrated if its predictions mean what they say: for presence-only data, whether places of twice the suitability hold twice the share of presences. Phillips and Elith (2010) showed how to read this from the predicted-to-expected ratio of the Boyce index (their presence-only calibration plot): for a model whose suitability is proportional to the probability of presence, the ratio grows in proportion to the suitability. Drawn from out-of-fold predictions, the curve tests the calibration on records the model did not see:

Enlarge

The curve is not a straight line: flat through the middle of the range and steep at the top. Read the map's values as a ranking (the top of the range is clearly better than the rest) rather than as proportional to the chance of finding the sloth. Phillips and Elith (2010) used such curves to recalibrate badly calibrated models; that is a model fitted on the out-of-fold predictions, done outside Terraflux on the exported table.

Thresholds from out-of-fold values. A threshold set on the training predictions is set on records the model has seen, which it predicts higher than it would predict new ones. Set on the out-of-fold predictions, it says where new presences fall. For the sloth the difference is small: the 10^th percentile of the presences is 0.29 out of fold (node 8) and 0.33 on the training predictions; the threshold of the largest TSS is 0.65 and 0.66. Apply the threshold to the map of the model fitted on every record, as before (connect oof_threshold.threshold to t of the binary node).

Residuals on a map. A presence with a low out-of-fold suitability is one the model would have missed had it not seen its block: its residual, 1 minus its prediction, is large. Residuals that cluster in space say that something about that region is not in the model (a missing predictor, conditions found nowhere else in the records) or not right in the records:

Enlarge

The circled presences share a climate: a larger annual range of temperature (bio7 of 14.1 to 17.9 °C, against a median of 11.4 °C over all presences) and, all but the one in Venezuela, a hotter warmest month (bio5 of 31.7 to 33.8 °C, against 31.6 °C). They are at the seasonal edge of the climate the sloth was recorded in. Only 12 of the 83 presences have an annual range above 14 °C, and 5 of them are in fold 5: leaving that fold out removes almost half the records of that edge at once. All 12 get an out-of-fold suitability of 0.50 or less. If these records are right (check them first), the model's picture of that edge rests on a handful of records, and its map there is the least certain.

Reporting#

Report the scores on test data (the folds' means and standard deviations, or independent data), with the background and the folds they depend on: AUC, TSS, Boyce and the omission rates together, never AUC alone. For a binary map, give the method and the value of the threshold, the sensitivity it keeps, and the suitable area, and say why that method suits the use.

Settings that change the answer#

• The data you evaluate on. Training scores flatter the model; random folds may too (part 4). The same model scores Boyce 0.92 on its training records and 0.55 over spatial folds.
• The background. AUC, TSS and specificity are all measured against the background: a wider background raises them all (part 3).
• The Boyce window. sdm.boyce (and the boyce of sdm.evaluate) uses a moving window a tenth of the range wide by default (window); nclass sets classes instead. Few presences make the index noisy whatever the window.
• The threshold method, as the slider shows: a factor of 15 in the area called suitable.
• Training or out-of-fold predictions. The same scores and thresholds from the out-of-fold predictions are test scores; the pooled Boyce index (0.93) and the fold mean (0.55) differ because of the number of presences per curve, not the model.

Common pitfalls#

• Reporting the AUC alone, or comparing AUCs of models with different backgrounds or study areas.
• Evaluating on the training records and calling it validation. Use the folds, or independent data.
• Choosing the threshold that gives the map you expected. Choose the method for the use, before looking, and report it.
• Kappa with background points. It depends on the prevalence, which the number of background points sets.
• A binary map where a continuous one would do. Thresholds throw away the gradient of suitability; keep it for prioritization, connectivity and comparisons.
• Out-of-fold predictions from random folds. They are only as independent as the folds: with random folds, a record's neighbours are in the training set, and the out-of-fold predictions flatter the model as the fold scores do (part 4).
• Pooling without the fold scores. The pooled out-of-fold predictions come from five models; report them with the scores of each fold.
• Reading specificity as absence. A background point predicted suitable is not an error in the field: it may hold the species.

Going further#

• Part 6: Projection and extrapolation: the model in another climate, and where its suitability cannot be trusted (south Chile is an example).
• Part 4: Fitting models: the cross-validated scores of sdm.fit.
• Binary maps as habitat for patterns and networks: MSPA.
• The operators: sdm.evaluate, sdm.threshold, sdm.boyce, sdm.omission_rate.

References#

• Allouche, O., Tsoar, A., & Kadmon, R. (2006). Assessing the accuracy of species distribution models: prevalence, kappa and the true skill statistic (TSS). Journal of Applied Ecology, 43(6), 1223–1232. https://doi.org/10.1111/j.1365-2664.2006.01214.x
• Boyce, M. S., Vernier, P. R., Nielsen, S. E., & Schmiegelow, F. K. A. (2002). Evaluating resource selection functions. Ecological Modelling, 157(2–3), 281–300. https://doi.org/10.1016/S0304-3800(02)00200-4
• Fielding, A. H., & Bell, J. F. (1997). A review of methods for the assessment of prediction errors in conservation presence/absence models. Environmental Conservation, 24(1), 38–49. https://doi.org/10.1017/S0376892997000088
• Hanley, J. A., & McNeil, B. J. (1982). The meaning and use of the area under a receiver operating characteristic (ROC) curve. Radiology, 143(1), 29–36. https://doi.org/10.1148/radiology.143.1.7063747
• Hirzel, A. H., Le Lay, G., Helfer, V., Randin, C., & Guisan, A. (2006). Evaluating the ability of habitat suitability models to predict species presences. Ecological Modelling, 199(2), 142–152. https://doi.org/10.1016/j.ecolmodel.2006.05.017
• Jiménez-Valverde, A. (2012). Insights into the area under the receiver operating characteristic curve (AUC) as a discrimination measure in species distribution modelling. Global Ecology and Biogeography, 21(4), 498–507. https://doi.org/10.1111/j.1466-8238.2011.00683.x
• Liu, C., Berry, P. M., Dawson, T. P., & Pearson, R. G. (2005). Selecting thresholds of occurrence in the prediction of species distributions. Ecography, 28(3), 385–393. https://doi.org/10.1111/j.0906-7590.2005.03957.x
• Liu, C., White, M., & Newell, G. (2013). Selecting thresholds for the prediction of species occurrence with presence-only data. Journal of Biogeography, 40(4), 778–789. https://doi.org/10.1111/jbi.12058
• Lobo, J. M., Jiménez-Valverde, A., & Real, R. (2008). AUC: a misleading measure of the performance of predictive distribution models. Global Ecology and Biogeography, 17(2), 145–151. https://doi.org/10.1111/j.1466-8238.2007.00358.x
• Merow, C., Wilson, A. M., & Jetz, W. (2017). Integrating occurrence data and expert maps for improved species range predictions. Global Ecology and Biogeography, 26(2), 243–258. https://doi.org/10.1111/geb.12539
• 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., & Elith, J. (2010). POC plots: calibrating species distribution models with presence-only data. Ecology, 91(8), 2476–2484. https://doi.org/10.1890/09-0760.1
• 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