Biosecurity Risk Mapping / estimating the establishment likelihood of pests 


Understanding where a harmful invasive species might arrive is critical for biosecurity decision support, especially given the finite resources available for surveillance. A risk map should also weigh the likely harm that would be caused if a pest was to establish itself in different areas. Analysis of potential harm is not part of this workflow, so from here on we will focus on estimating the likelihood of a pest establishing itself (“establishment likelihood”). 


Camac et al. 2020 provide a more comprehensive guide to producing scientifically robust models, and we suggest you review that report to both deepen your understanding of the methodology and the ways data and models might be improved. A more recent report Camac et al. 2021 provides important updates to this methodology. 





What is an establishment likelihood map? 


An establishment likelihood map is made of a grid of equal-sized grid cells over the area of interest, often over the whole of Australia. The establishment likelihood at each grid cell is calculated by multiplying up to three values in three grids: 1) in each grid cell the value in the abiotic suitability grid is usually the probability that conditions in that grid cell are not too hot, cold, dry or wet for the pest to survive; 2) in each grid cell the value in the biotic suitability grid is an estimate of how likely it is that the biological material or host required by the pest exists in that grid cell; and 3) in each grid cell the value in the pest arrivals grid describes the proportion of viable pest arrivals likely to occur in each grid cell. 


In more detail: 


Establishment likelihood = Abiotic suitability × Biotic Suitability × Probability of Pest Arrival  


Abiotic Suitability = the probability that a species can survive and/or reproduce given the abiotic conditions. 


Usually, abiotic suitability is defined by areas known to be too cold, hot, wet, or dry for the pest to survive. If the pest species is likely to do well in all the known climatic conditions available within Australia this step could be skipped, and you could upload a grid where all values were set to 1. If the species cannot survive in areas that reach freezing temperatures but can survive in other climatic conditions, it might be appropriate to simply mask out those areas that freeze (i.e., assign 0 probability) and give everything else a probability of 1. Similarly, if there is a high level of certainty that the suite of climatic conditions available to the species in its native range would define the locations used in Australia, an environmental envelope would capture those conditions, and the mean values could be given more weight if binary scoring was not of interest. If you the user want to prioritise capturing the variation in abiotic suitability, then you may want to try to generate your abiotic suitability using a species distribution model using machine learning or statistical methods that produce a wider breadth of probabilities defining the suitable climatic niche. 


Most people take a conservative approach to this step as they want to make sure not to rule out locations on a map where the establishment of a pest is possible. Keep in mind how many species of plants or animals are found outside their native climatic niches within botanic gardens or zoos throughout the globe. The key here is to identify those locations where a species will not survive, or is less likely to survive, and that really depends on the assumptions made based on your literature review. 


Biotic Suitability = the presence/absence or relative suitability (scaled to the range 0-1) of hosts or host material, including food or habitat. 


There are a variety of ways to calculate the value in each grid cell representing biotic suitability.  If the pest is a livestock disease, biotic suitability may simply be the number of livestock animals in each grid cell. If the disease can be transmitted by an insect, the number of livestock animals might be multiplied by the probability that the insect would be found in that grid cell based on a species distribution model.  


If a plant pest is likely to use any green vegetation as host material, a grid of vegetation greenness (NDVI) might be a good layer to use. Land use layers are also likely to be useful to identify which grid cells are likely to have suitable host material. The Australian Land Use and Management classification (ALUM), which maps everything from where different crops are found to natural area locations, is perhaps the most commonly used land-use layer for all of Australia. 


A critical assumption when computing biological suitability relates to the likelihood that a pest will jump to an Australian host or habitat. Often there is no data on how pests from other continents might interact with Australia’s unique ecology. 


Probability of Pest Arrival = the estimated probability of one or more viable pest arrivals in each grid cell. 


Most pests will have many ways in which they may enter or spread throughout Australia, i.e., tourists, mail, shipping containers, etc. Each of these different modes of entry is considered a pathway. While the biosecurity system aims to prevent entry of pests and other biosecurity hazards, these controls are not 100% effective. Pests that evade the biosecurity controls correspond to “leakage” events (see below and Camac et al. 2021). The probability of a pest arriving at a given location (i.e., a grid cell) in a given time period depends on: (1) the leakage rate (average number of leakage events per time period); (2) the probability that a leakage event will result in the pest’s viable establishment somewhere in susceptible habitat (e.g., in a location that contains host material); and (3) the spatial distribution of the probability of establishment for the relevant pathway (i.e., assuming a viable establishment occurs somewhere, what is the probability that it will occur in a particular cell, relative to other cells).  

Leakage and viability rates have been estimated for a range of biosecurity hazards by the federal Department of Agriculture, Fisheries and Forestry, based on expert elicitation (Hemming et al. 2018) and building on previous efforts such as the department’s Risk-Return Resource Allocation (RRRA) model. The latter might be informed by the relative volume of containers arriving in each postcode (for a container pathway) or human population density (e.g., if the pathway is mail or residents returning from international travel), or may be a function of distance from points of entry such as ports (e.g., vessels pathway) or airports (e.g., tourists pathway). 


Given the leakage and viability rates, and the spatial arrival weights, the probability of a location encountering one or more arrivals is calculated for each relevant pathway, and the resulting grids are combined by calculating their union to produce a single grid giving the overall probability of establishment. 


Assumptions & caveats


This workflow does not provide guidance on how to prioritise different pests or diseases, how to optimally deploy surveillance or how to account for climate change. This workflow does result in a map of establishment likelihood, but the accuracy of that map is dependent on the data used, and the assumptions made while constructing the model. It is easy to produce results with the tools available on Biosecurity Commons, however, the utility of those results will depend fully on you the user. Also, it is worth remembering that results from these workflows are relative to the data and modelling decisions made, so while they provide a map of relative establishment likelihood, the extent to which quantitative comparisons can be made across species will depend on the extent to which the components themselves (e.g. the underlying estimates of abiotic suitability) can be validly compared. There is no single best method, so take care, evaluate your assumptions, and validate your model estimates whenever possible. A full description of considerations to include when attempting to generate a useful establishment likelihood map is presented in recent reports (Camac et al. 2020, Camac et al. 2021). 


Remembering that:


Establishment likelihood = Abiotic suitability × Biotic Suitability × Probability of Pest Arrival


This workflow does not provide an estimate of the uncertainty in each variable being multiplied. However, we know that the overall error in the establishment likelihood estimate is a combination of the individual errors. Therefore, if you are highly uncertain about one of these estimates, that uncertainty carries through to your final estimate. While the platform requires you to enter a grid (raster) for each of the components (abiotic, biotic and pest arrival), if you are highly uncertain of one of these components it may be best to upload a raster where all the values in the study area have been set to one. The results would then need to be interpreted accordingly.  Research into optimal ways to include uncertainty in results is ongoing and should be encouraged.


Until uncertainty is captured quantitatively it is vital to capture the variety of assumptions underlying an establishment likelihood map and document the data used as well as how those data were derived. These steps are vital to communicating the uncertainty in the estimate, but these steps also make it easier to improve future estimates. Measures of uncertainty are commonly reported when generating more complex species distribution models, but those reported levels or uncertainty can only be trusted once the model is tested with independent data and any underlying assumptions are thoroughly tested. The same applies to biotic suitability and probability of pest arrivals.


It is worth reiterating that establishment likelihood values are relative.


If the abiotic suitability is simply set to 1 for areas deemed to have suitable conditions and 0 for unsuitable areas, then when this binary abiotic suitability is multiplied by biotic suitability or pest arrival probabilities, the resulting establishment likelihood values will be more extreme (i.e., they would be lower when suitability thresholding leads to a 0) than they would have been if the abiotic suitability had been scaled with a range of values between zero and one. This is just one simple example of why strict comparisons between numeric values of establishment likelihood from different models or species should be avoided.


Additional assessment of the potential financial, ecological, or societal harm of a potential establishment of an invasive species should be done separately.


Finally, while results arising from this approach may in some cases have a relatively high level of uncertainty, they provide an estimate of the relative likelihood of pest establishment which can be shared and improved upon. The tools available here will improve the transparency, reproducibility and consistency of workflows while also speeding up the modelling process. Importantly, by sharing workflows the platform promotes collaboration, enabling users to improve on existing estimates rather than reinventing them.


8 practical steps to estimating establishment likelihoods on Biosecurity Commons


First a reminder: In the workflows available on Biosecurity Commons when you run a final model or click the run button for any of the steps, you have the option to rerun each step if you change something e. g. parameter or input.  If you do rerun a step, the previous run will be overwritten. You can save any step by using the “Export Result” button. Exporting a Result will make it available under 'My Results' where it will be permanently available for future use.  It will also be retained in this Project view under 'My Exported Results'. 


Step 1 – Review the literature  

A literature review might reveal:

  • An existing risk map or establishment likelihood map that could be updated or improved upon.
  • Pest biology - physiological tolerances indicating climatic conditions where the species would not survive.  
  • Susceptible hosts - understanding the hosts the pest may require, where those hosts are distributed, and the potential for the pest to use uniquely Australian hosts.
  • Risk pathways: where are viable pests likely to come into Australia?
  • Existing distribution models of the pest and/or data useful for predicting distributions.
  • Clean, ready to use occurrence data. 
  • Existing data and worked examples on Biosecurity Commons.


Step 2 – Collate and clean data 


Biosecurity Commons will have a growing set of data that users will be able to access, but those data will often need to be supplemented with additional data not available on the platform and most data will likely need to be cleaned. Data will include rasters or grids which will in some cases need to be generated from other spatial data (e.g., rasterising vector datasets). 


At Biosecurity Commons we strongly believe that the best data will be uploaded by you the users.  Over time, at Biosecurity Commons we expect the available data that can inform establishment likelihoods will grow and as it is shared it will be improved. Much of the data available on the platform have not been fully validated with independent data or had underlying assumptions fully tested.  Over time if you upload some data, other users may improve those data.


The most common point data you are likely to use are from occurrence data. Large databases like GBIF provide excellent sources of occurrence data for developing abiotic distribution models. However, these data often have a wide range of records that would be inappropriate for any analysis. Cleaning such data requires the removal of records with taxonomic errors, temporal errors, irrelevant records, and records from non-established populations. If bringing in coordinate data from multiple sources, you may also need to transform the coordinate system so all latitude and longitude values use the same system (see step 3 for further explanation). These data cleaning steps are critical but sometimes time-consuming. 


Currently, you will often need to download the data, clean it and re-upload it. There are a variety of filtering options directly available on GBIF, and that filtered data is given a DOI when you download it. You can take these steps directly on Biosecurity Commons if using Atlas of Living Australia (ALA) data, but they have not been implemented for GBIF data yet.


Further, it may be best if the occurrence data used for the abiotic species distribution model are restricted to the species' native range, an area where the species is more likely to have occupied all the suitable areas. However, if species records outside the native range are found in areas with different climatic conditions, excluding these records could inappropriately remove areas from the establishment likelihood map that are areas at risk.


The importance of accounting for survey bias and choosing an extent are important considerations that are well covered elsewhere (Camac et al. 2020, EcoCommons).


Step 3 – Select the extent (Study Region)


The next step involves selecting your study region. This region will be represented by a template grid that has three characteristics that will be used for all other rasters used in workflows or when displaying the results.  These three characteristics include the extent (inclusive of the boundaries of the study area), resolution (single size of all grid cells in the study area), and coordinate reference system (CRS, the coordinate system and spatial projection used to turn the round earth into a flat map). 


The default template is a 1 km2 resolution raster with an extent inclusive of the entire of Australia, with the GDA94 / Australian Albers (EPSG:3577) coordinate system. This is the most commonly used projected coordinate system for Australia. Where necessary, inputs are automatically projected and resampled to the extent, resolution, and CRS of the chosen template layer ("Study Region"). However, note that these transformations require that each spatial dataset contains metadata defining its CRS. This is particularly relevant if uploading your own datasets. 


Raster data that you upload will need to be in GeoTiff format with the CRS, extent and resolution defined in the file you are going to write. 


Coordinates that you upload will need to be in decimal degrees (WGS 84, EPSG:4326) and the column headings will have to be lat and lon in a table that is saved in CSV format.


lat                         lon

-27.38888            152.87777

-27.38788            152.87888


Step 4 – Abiotic Pest Suitability


Results of abiotic suitability should be scaled from zero to one, or occasionally a binary result of zero or one is appropriate.


Again, if the pest is thought equally likely to survive in any of Australia’s climatic conditions, you could upload a raster where all values were set to 1.


If climatic tolerances are well understood for the pest species, possibly discovered during the literature review, a simple reclassification of a temperature grid, or rainfall grid might be best. For example, a review may indicate that the pest only survives in areas with a minimum temperature of 10 °C, or a maximum temperature of 35 °C. Currently, such reclassification is not available on the platform, but it could be done in R, Python, or GIS software. Reclassification would simply reclassify all areas with temperatures above or below these thresholds to zero, and any temperatures between those thresholds as 1. This kind of layer would mask out all the zero cells from the final establishment likelihood. 


If climatic tolerances are less binary, and species are expected to be more likely to survive or thrive as ranges of bioclimatic variables change, then one of the environmental profile models might be best. Range-bagging, surface range envelopes or Bioclim models all tend to result in conservative estimates of where the pest might occur, often resulting in more geographic space included as suitable in the resulting geographic predictions. 


If pest species are expected to be more likely to survive or thrive as conditions approach a climatic optimum, then more complex or correlative species distribution models may be best.  However, such models tend to be less conservative in that they tend to predict a smaller geographic area as being potentially suitable. These models are often trained in the areas of the globe where the pest species is native. It is important to realise that the climatic conditions where a species is found on one continent, may not correspond to where they occur in Australia because pest occurrences in its native range were correlated to climatic conditions, but not actually limited by those conditions [literature review may help determine if climatic conditions are limiting].  Some species will show more tolerance to broader climatic conditions than others.  Understanding the biology, physiology, and ecology of potentially invasive species will be critical to appropriately setting model thresholds and arguments, and interpreting model results. “Accurate prediction to the region in which a model was fitted does not guarantee accurate prediction outside this range (Elith, 2017; Fourcade et al., 2017).” – (Camac et al., 2020)


When developing these models there are a variety of decisions to be made.  First, decisions around which occurrence points to include when training the model can have a large impact on results.  Second, the selection of background points (Maxent) or pseudo-absence points (most other statistical or machine learning models) can have a massive impact on the results (Phillips et al., 2009; Warton & Shepherd, 2010; Syfert et al., 2013). These background or pseudo-absence points need to be selected in places where the species could have dispersed and in some models, it helps if they exhibit the same survey bias as was present in the occurrence data (Camac et al., 2020). Selection of points where similar species were observed but the target species was not, or selecting points randomly from a bias layer are two ways to generate such data (see SDMs in R; Step 1 module).


Thresholding resulting probabilities from these kinds of models to generate binary outputs is not recommended as it degrades the resolution of information available to decision-makers.


A variety of species distribution models could be run, including:


Environmental Envelope Models (Most commonly used in risk mapping)

Bioclim

Surface Range Envelope

Range Bagging (Coming soon)

Climatch (Coming soon)


Machine Learning Models (Maxent does not require absence or pseudo-absence data, but all others do)

Maxent (widely used, and good predictive performance)

Boosted Regression Tree (also a good predictor)

Artificial Neural Network (also a good predictor, requires many data points)

Random Forest 

Classification Tree 


Statistical Models (alrequire absence or pseudo-absence data)

Flexible Discriminant Analysis

Generalized Additive Model

Generalized Linear Model

Multivariate Adaptive Regression Splines


Geographic Models (simple models that may be used as a preliminary step when modelling)

Circles

Convex Hull

Geographic Distance

Inverse-Distance Weighted Model (a relatively simple way to generate a bias layer when using all survey locations)

Voronoi Hull Model


Within the platform you can add an abiotic suitability layer in a variety of ways.  First, you can calculate an abiotic species distribution model using the "species distribution modelling" workflow within the Biosecurity Commons "Analysis Hub". Once completed you can add the result by selecting the "Choose from my results" option by selecting the ADD button under the Abiotic Suitablity "Input Parameters" section.


You can also select or explore datasets that you have uploaded / imported by clicking on “Explore My Datasets”.  You can also explore the thousands of curated datasets by clicking on “Explore Curated Datasets” or click on the “Import / upload data” to add your own data.  


If needed you can also perform a variety of functions on the data you input, including aggregating categories, or applying a distance weight layer (see below for more details).


Step 5 – Biotic Pest Suitability


There are many ways to generate a biotic pest suitability map (Camac et al. 2020). A simple example would be to consider areas with land use classes relating to agriculture in the ALUM land use layers as potential areas for a widespread agricultural pest (scoring cells involving agriculture as 1, and others as 0). If your pest uses orchards as hosts you may use ALUM land use to identify where orchards occur, and aggregate a finer scale NDVI layer to give a total sum of vegetation greenness within that cell, which could be used to weight each cell with orchards. You may also include urban areas to account for citrus trees growing in the odd backyard. If trying to estimate the potential distribution of a fungus like Myrtle Rust, you might develop a conservative species distribution model using occurrence points from selected species in the Myrtaceae family.  For a livestock disease you may aim to get a map of where livestock are kept, and try to estimate livestock abundance within each grid cell. Often the final layer may include combining (e.g. adding) multiple layers. ALUM, NVIS and NDVI layers are held as options you can select within the biotic pest suitability's add a new input option.


Biotic suitability results should be scaled from zero to one.


We will be adding some more example biotic suitability layers, but the biggest improvements in available layers will come from you the experts in biosecurity.


Here you can also, select previous workflow results to use here or explore datasets that have been uploaded or imported by others.  You can also explore the thousands of curated datasets or import / upload your own data.  


Step 6 – Pest Suitability


Pest suitability default option is to multiply (intersect) the Abiotic Suitability by the Biotic Suitability. Probabilities should be multiplied if the two probabilities are independent, in other words the probability of Abiotic Suitability (climate) is not related to the probability of Biotic Suitability (host material). However, it is also possible to add layers together instead of multiplying them. You would only add layers if the probability at a location was based on either Abiotic Suitability or Biotic Suitability, but not both.  For example, if in the areas Abiotic Suitability was zero the Biotic Suitability was greater than zero and vice versa. A union of probabilities adds the two probabilities, but then subtracts the locations where the probabilities intersect which is appropriate if the two probabilities are not independent:


( 1- (1 – Abiotic Suitability probability)*(1 – Biotic Suitability probability) )


Step 7 – Pest Arrivals


A variety of grids have been generated based on the examples given by (Camac et al. 2020) of the proportion of units (tourists, mail parcels, containers etc) that might reach any grid cell. For the number of tourists entering the country, the total number of tourists are divided among the airports within Australia based on the percentage of international arrivals at each airport. Then data indicates that the majority of tourists stay relatively close to the airport or tourist accommodation, and a distance function is applied to indicate that very few tourists go far from the airport or tourist accommodation, with the number dropping off as you go further away.  For each pathway, and each priority pest the number of units in any grid cell is used in a function that requires two additional inputs (parameter estimates).  The functions that generate the final arrivals probability that can be made up of multiple pathways can include up to eleven formulae (see Camac et al. 2021, pages 11 –12).


The first additional input is the leakage rate. The leakage rate corresponds to the upper and lower bounds of the number of consignments expected to arrive in Australia each year that could constitute a leakage event. In other words, how many events (containers, mail packages etc.) got past existing biosecurity protections and made it into Australia. This estimate is formed using a combination of interception rates and expert opinion (Hemming et al. 2018), and the estimates are usually made by the federal Department of Agriculture Fisheries and Forestry. The estimate includes the 95% confidence interval of the estimate (lower 2.5th and 97.5th percentiles). Examples of a leakage event includes the arrival of an international passenger with a plant disease on an undeclared plat, or molluscs on the bottom of an arriving boat.


The second additional input required from users includes the estimated 95% confidence interval (lower 2.5th and upper 97.5th percentile) of the expected viability of leakage events. Estimates of viability are often based on a combination of considerations related to how likely would it be for a pest to survive a leakage event.  This can include how likely is it that the pest would survive the trip along the pathway, will there be enough pests to permit establishment, will the end point of the pathway be a place the pest can survive.  As with the leakage rates, these estimates are made based on expert elicitation (Hemming et al. 2018) by the federal Department of Agriculture, Fisheries and Forestry.


Most pests can enter Australia from multiple pathways, so the final pest arrival probability is usually the union of each pathway probability of entry grid.  The union of individual pathway probabilities is usually what is used because arrival pathway probabilities are usually not independent. Other methods may also be used, and better estimates of arrivals can be uploaded by you the user.


The first step to generating a pest arrivals pathway is to click on the “Add New Input” button on the Input Parameters tab.


There are the building blocks for you to build nine different pathways including: Residents Pathway, Tourists Pathway, Torres Straight Arrivals Pathway, Vessels Pathway, Containers Pathway, Mail Pathway, Agriculture Pathway and Machinery Pathway.  You can also choose to build your own with a “Generic Pathway”, or choose layers from your results, datasets, or the curated datasets. Select one of these options at a time by choosing "Add New Pathway" from the input parameters of Pest Arrivals. Again, most probabilities of pest arrival are comprised of multiple pathways.


Remember when probabilities are not independent, use of a union is appropriate. 


( 1- (1 – Tourists Pathway probability)*(1 – Mail Pathway probability) )


Step 8 – Pest Establishment Likelihood


The default for this final step is to multiply (intersect) the Pest Suitability probability (remember this was the Abiotic Suitability multiplied by the Biotic Suitability) by the Pest Arrivals probability grid. Probabilities should be multiplied if the two probabilities are independent, in other words the probability of Pest Suitability is not related to the probability of pest arrivals. However, it is also possible to add layers together instead of multiplying them. You would only add layers if the probability at a location was based on either Pest Suitability or Pest Arrivals, but not both. A union of probabilities would be appropriate if the Pest Suitability and Pest Arrivals were not independent. Union of probabilities adds the two probabilities, but then subtracts the locations where the probabilities intersect:


( 1- (1 - Pest Suitability probability)*(1 – Pest Arrivals probability) )


It is also possible to leave out abiotic suitability, biotic suitability or pest arrivals if there is no helpful information for any of those three estimates. Do this by uploading a grid where all the values in your Study Region are set to 1. The variety of ways the user can compute an overall Pest Establishment Likelihood is one reason the values of results run on different pests cannot be quantitatively compared.  In other words, results are relative or dependent on the decisions made during the construction of the model. Still the relative probability maps are useful for understanding where pest establishment is most or least likely.


Functions described further


Distribute Features: Distributes feature polygon values across a template or mask raster spatial layer. The total value assigned to each polygon is divided across all corresponding raster cells when a template (single value plus NAs) raster is provided, or across non-zero cells when a mask (zero and non-zero values plus NAs) raster is provided, which is used as a binary mask (regardless of the variety of non-zero values).


Aggregate Categories: Aggregate the categories within a spatial layer to a courser resolution, that of another (template) layer, based on the presence of each of the user-selected categories in each cell. The resulting aggregated layer cells may be binarized to one or zero, or indicate the proportion of selected categories present in the cell.


Distance weight layer: Calculates a distance-weighted spatial layer based on the proximity of each cell to a series of points (features), or via a pre-calculated distance layer. Each cell is weighted via a negative exponential function applied to the distance of the nearest point or pre-calculated distances.”  (also one way to generate a bias layer)



Literature cited:


Camac, J. S., Baumgartner, J., Robinson, A., & Elith, J. (2020). Developing pragmatic maps of establishment likelihood for plant pests. Technical Report for CEBRA project 170607.


Camac, J. S., Baumgartner, J. B., Hester, S., Subasinghe, R., & Collins, S. (2021). Using edmaps & Zonation to inform multi-pest early-detection surveillance designs. Technical Report for CEBRA project 20121001.


Elith, J. (2017). Predicting distributions of invasive species. In: Invasive species risk assessment and management (eds. Robinson AP, Walshe T, Burgman MA, Nunn M), pp. 93–129. Cambridge University Press.


Fourcade, Y., Besnard, A. G., & Secondi, J. (2017). Paintings predict the distribution of species, or the challenge of selecting environmental predictors and evaluation statistics. Global Ecology and Biogeography, 27, 245–256.


Hemming, V., Burgman, M. A., Hanea, A. M., McBride, M. F., & Wintle, B. C. (2018). A practical guide to structured expert elicitation using the IDEA protocol. Methods in Ecology and Evolution, 9(1), 169-180.


Phillips, S. J., Dudik, M., Elith, J., Graham, C. H., Lehmann, A., Leathwick, J., & Ferrier, S. (2009). Sample selection bias and presence-only distribution models: implications for background and pseudo-absence data. Ecological Applications, 19, 181–197.


Syfert, M. M., Smith, M. J., & Coomes, D. A. (2013). The Effects of Sampling Bias and Model Complexity on the Predictive Performance of MaxEnt Species Distribution Models. PLoS ONE, 8, e55158–10.


Warton, D. I., & Shepherd, L. C. (2010). Poisson point process models solve the “pseudoabsence problem” for presence-only data in ecology. The Annals of Applied Statistics, 4, 1383–1402.