Tropical forests are increasingly degraded by industrial logging, urbanization, agriculture, and infrastructure, with only 20% of the remaining area considered intact. However, this figure does not include other, more cryptic but pervasive forms of degradation, such as overhunting. Here, we quantified and mapped the spatial patterns of mammal defaunation in the tropics using a database of 3,281 mammal abundance declines from local hunting studies. We simultaneously accounted for population abundance declines and the probability of local extirpation of a population as a function of several predictors related to human accessibility to remote areas and species’ vulnerability to hunting. We estimated an average abundance decline of 13% across all tropical mammal species, with medium-sized species being reduced by >27% and large mammals by >40%. Mammal populations are predicted to be partially defaunated (i.e., declines of 10%–100%) in ca. 50% of the pantropical forest area (14 million km2), with large declines (>70%) in West Africa. According to our projections, 52% of the intact forests (IFs) and 62% of the wilderness areas (WAs) are partially devoid of large mammals, and hunting may affect mammal populations in 20% of protected areas (PAs) in the tropics, particularly in West and Central Africa and Southeast Asia. The pervasive effects of overhunting on tropical mammal populations may have profound ramifications for ecosystem functioning and the livelihoods of wild-meat-dependent communities, and underscore that forest coverage alone is not necessarily indicative of ecosystem intactness. We call for a systematic consideration of hunting effects in (large-scale) biodiversity assessments for more representative estimates of human-induced biodiversity loss.
Citation: Benítez-López A, Santini L, Schipper AM, Busana M, Huijbregts MAJ (2019) Intact but empty forests? Patterns of hunting-induced mammal defaunation in the tropics. PLoS Biol 17(5): e3000247. https://doi.org/10.1371/journal.pbio.3000247
Academic Editor: Andy P. Dobson, Princeton University, UNITED STATES
Received: October 23, 2018; Accepted: April 10, 2019; Published: May 14, 2019
Copyright: © 2019 Benítez-López et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: Data are available at https://figshare.com/projects/Intact_but_emtpy_forests_Patterns_of_hunting-induced_mammal_defaunation_in_the_tropics/31118.
Funding: The authors received no specific funding for this work.
Competing interests: The authors have declared that no competing interests exist.
Abbreviations: BIC, Bayesian Information Criterion; BII, Biodiversity Intactness Index; CAR, Central African Republic; DHS, Demographic and Health Surveys; DI, defaunation index; DRC, Democratic Republic of Congo; HDE, Humanitarian Data Exchange; HF, human footprint; HPD, human population density; IF, intact forest; IFL, intact forest landscape; IUCN, International Union for Conservation of Nature; MESS, multivariate environmental similarity surface; ML, maximum likelihood; OSM, OpenStreetMap; PA, protected area; REML, restricted maximum likelihood; RR, response ratio; WA, wilderness area
Tropical forests are increasingly degraded by industrial logging, urbanization, agriculture, and infrastructure [1, 2], with only 20% of the remaining area considered intact. Yet, this figure does not include other, more cryptic but pervasive forms of degradation, such as losses of wildlife due to overhunting . Although humans in tropical areas have hunted for millennia to secure food and income, the current hunting rates are unsustainably high across the tropics due to the demand from growing human populations, an increasing commercialization of wild meat, and higher human accessibility to otherwise remote areas [4–8]. As a consequence, many tropical forests in the developing world are becoming “empty” (sensu Redford, 1992 ), with wildlife populations reduced or locally extirpated (i.e., defaunated), resulting in concomitant catastrophic consequences for ecosystems and the services and livelihoods that they provide . Despite the increasing evidence of overhunting in the tropics, there is virtually no information about the spatial variation of hunting-induced defaunation and the areas where impacts might be most severe. Mapping “hotspots” of defaunation is crucial as a first step towards identifying and quantifying possible consequences for ecosystem functioning, as well as designing more targeted conservation measures.
Overhunting, as opposed to deforestation, is undetectable by remote-sensing techniques , and results from local studies have limited applicability to other areas [12–15]. However, upscaling local data with models based on quantitative relationships between impacts on wildlife populations and the main drivers of hunting pressure represents an unexplored approach to predict large-scale defaunation, particularly in understudied areas. Some of the main drivers of hunting include hunters’ accessibility to wildlife resources via road development and settlement establishment , hunters’ preferences for certain species [3, 10], and proximity to urban markets . Additional factors are human population growth and subsequent increases in wild meat demand, socioeconomic status, food security, and governmental controls on hunting via law enforcement in protected areas (PAs) . While the main drivers of such a multifaceted phenomenon were recently identified , the spatial pattern of hunting pressure on wildlife populations at the pantropical scale remains elusive.
Here, we model and project the spatial patterns of hunting-induced mammal defaunation in the tropics and identify areas where hunting impacts on mammal communities could be high. To this end, we developed a modelling framework based on a suite of important socioeconomic drivers of hunting pressure and taking into account the vulnerability of species to hunting . Our models were based on the most extensive database to date on hunting impacts on mammal populations, consisting of 3,281 mammal abundance estimates in hunted and non-hunted areas extracted from 163 studies (S1 Table, https://figshare.com/projects/Intact_but_emtpy_forests_Patterns_of_hunting-induced_mammal_defaunation_in_the_tropics/31118). We related the abundance declines to socioeconomic drivers of hunting, including the distance to hunters’ access points, accessibility of urban markets, protection status (whether hunting occurred inside or outside PAs), human population densities, poverty levels, and access to domestic meat. Species body mass and diet were included as proxies of vulnerability to hunting  (S1 and S2 Figs). We modelled abundance declines with hurdle models to simultaneously incorporate the probability of being locally extirpated as well as abundance reductions (Methods). Subsequently, we used our models to map defaunation gradients across the tropics and quantify the magnitude and spatial extent of the population declines of 3,923 mammal species. We averaged the declines across species into a defaunation index (DI) ranging from 0 (intact mammal assemblage) to 1 (fully defaunated mammal assemblage). We conservatively consider areas with a DI >0.1 (more than 10% average reduction in mammal abundance across all species) to be partially defaunated (hereafter defaunated), and areas with DI >0.7 to be severely defaunated. We identified defaunation hotspots in areas where at least one third of the species had declines >70%. We overlaid our defaunation maps with intact forest (IF)  and human footprint (HF)  maps to assess the extent to which these pristine landscapes could be defaunated. Both initiatives indicate the location and extent of IFs, but do not consider the effect of hunting . Finally, we assessed potential hunting-induced mammal abundance declines in PAs (PA, International Union for Conservation of Nature [IUCN] I–IV categories).
Results and discussion
Distance to hunters’ access points, species body mass, and human population density (HPD) were the most important predictors of defaunation, followed by stunting and the protection status of the area (S1 Text, S5 Table, S3, S4, S5 Figs). More than half of the pantropical forest area (60%) is located within 10 km of the nearest human settlement, and >80% is within 20 km (S2 Fig). These are typical distances that hunters travel in the tropics , thereby suggesting that most of the tropical forest area is relatively accessible for hunters. Despite the complexity of the phenomena being modelled, the predictive performance of our models was satisfactory, with an overall sensitivity and specificity of 0.5 and 0.7, respectively, and an overall balanced accuracy of 0.6 for the different defaunation categories (S6B and S6C Fig). Pseudo-R2 values were 0.32 and 0.24 for the full and cross-validated models, respectively. Model performance was highest for the high defaunation level (DI > 0.7, S6B and S6C Fig); hence, the models are particularly useful to identify potential defaunation hotspots.
Across all mammal species and the entire pantropical area, we estimated an average DI of 0.13 ± 0.1 (mean ± SD, median: 0.09, IQR: 0.17, N = 30,004,854 grid cells) (Fig 1A). Results were highly similar if we removed areas outside the socioeconomic domain covered by our data (DI: 0.12 ± 0.1, median: 0.08, IQR: 0.16, N = 29,030,794 grid cells). For large mammals (>20 kg), we predicted an average decline of more than 40% across the tropics (DI: 0.42 ± 0.3, median: 0.43, IQR: 0.53, Fig 1D). An average reduction of 27% was predicted for medium-sized mammals (1–20 kg) (DI: 0.27 ± 0.2, median: 0.21, IQR: 0.42, Fig 1C), whereas impacts on small mammals (<1 kg) were negligible (DI: 0.05 ± 0.2, median: 0, IQR: 0.04; Fig 1B), reflecting that these are usually not hunted . Our estimates are lower than the previously reported average decline of 83% , as here we are predicting hunting pressure in areas that are less hunted or intact, and for all mammals in the tropics, including a larger proportion of small species than in the original data set based on empirical studies (70% versus 20%).
The insets represent the total area (y-axis) under different levels of defaunation (x-axis, from DI = 0 [blue] to DI = 1 [red]). Note that the y-axes in the four insets have different scales. Available at https://figshare.com/projects/Intact_but_emtpy_forests_Patterns_of_hunting-induced_mammal_defaunation_in_the_tropics/31118. DI, defaunation index.
We predict that approximately 47% of the pantropical forest area (ca. 14 million km2) is defaunated (DI > 0.1, Fig 1A), with mammal population declines of at least 50% in 3.5% (540 thousand km2) of the tropical forests. The average DI was highest in countries from West and Central Africa, particularly in The Gambia, Ghana, Togo, and Cameroon (DI range: 0.3–0.5, Fig 1A and S7 Fig), followed by some of the Asian countries (Thailand and Bangladesh). Nigeria, Burundi, Rwanda, Sri Lanka, and Java (Indonesia) also had high DIs, but these estimates were lower when we removed areas where local human population densities exceeded those covered in our data set (S8 Fig). For the rest of the tropics, our predictions were well within the range of observed values of the socioeconomic predictors (S8 Fig).
We identified hotspots of hunting-induced defaunation in West and Central Africa (Cameroon, Guinea, and Cote D’Ivoire), Central America (Panama, Mexico, Costa Rica, Guatemala, and Honduras), Northwest South America (Colombia, Venezuela), and some areas in Southeast Asia (Thailand, Malaysia, and Southwest China, Fig 2). In the West African countries, particularly in Cameroon, more than 50% of the total number of species are predicted to have their populations reduced by 70%–100% because of hunting activities. Non-defaunated refugia (DI ≤ 0.1) were identified in the Guiana shield (Suriname, Guyana, and French Guiana) and the Brazilian Amazon, which represent inaccessible regions and sparsely populated areas (Fig 1A and S7 Fig).
Colors range from green (low relative number of species with DI >0.7) to red (high relative number of species with DI >0.7). Orange to red areas showcase hotspots of hunting-induced mammal defaunation (with at least one third of the species with DI >0.7). Available at https://figshare.com/projects/Intact_but_emtpy_forests_Patterns_of_hunting-induced_mammal_defaunation_in_the_tropics/31118. DI, defaunation index.
Our maps with hunting-induced declines of medium- and large-sized mammals resemble spatial patterns of biomass harvest reported for Central Africa [21, 22], with high hunting pressure predicted in west, central, and north Cameroon, the Albertine Rift, south and north Democratic Republic of Congo (DRC), south Congo, and Gabon. Similarly to these earlier studies, we also identify hunting-free areas in central Congo, north Gabon, central DRC, and east Central African Republic (CAR). Furthermore, our defaunation maps for medium- and large-sized species suggest that high hunting pressure is expected along the Amazon River network and, particularly, in the east Amazon basin, which matches previous results on the spatial extent of hunting pressure for Ateles spp. (approximately 8.5 kg) .
Our results are consistent with the idea that hunting is downsizing tropical mammal communities [3, 10]. The hunting-induced alteration of body size distributions across the tropics could trigger shifts in ecological functioning by impairing key ecological processes such as seed dispersal, predation, and herbivory, for which large-bodied species play a significantly more important role than smaller species [23–26]. We further estimated that declines were more severe for carnivores and herbivores (DI: 0.24 ± 0.2, median: 0.19, IQR: 0.37 and 0.22 ± 0.2, median: 0.17, IQR: 0.28, respectively) than for frugivores and insectivores (DI: 0.09 ± 0.1, median: 0.03, IQR: 0.1 and 0.06 ± 0.1, median: 0.02, IQR: 0.07, Fig 3). These results reflect differences in body mass distributions between feeding guilds (frugivores and insectivores are, on average, smaller than carnivores and herbivores), as well as differences in the spatial distribution and extent of the geographic ranges of the different guilds, with each species range encompassing diverse values of the socioeconomic drivers of hunting. Overall, the loss of large carnivores and herbivores may diminish top-down and bottom-up regulation, which can in turn trigger trophic cascades, resulting in the destabilization of tropical ecosystems and eventually leading to a net loss of diversity [27, 28].
The dashed gray line indicates the mean DI across the pantropical forest zone. The y-axes have different scales. Available at https://figshare.com/projects/Intact_but_emtpy_forests_Patterns_of_hunting-induced_mammal_defaunation_in_the_tropics/31118. DI, defaunation index; sq., square.
We estimated that approximately 9% of IFs and 11% of wilderness area (WAs) are defaunated when full mammal assemblages are considered (Fig 4A and 4C). Large-bodied mammals could be, however, defaunated in more than half of the remaining IFs (52%, 2.8 million km2) and WAs (62%, 4.3 million km2), with barely any remaining intact areas in Central Africa (Fig 4B and 4D). Compared with degraded forests, IFs support more imperilled biodiversity and exceptional environmental values, including carbon sequestration and storage, water provision, indigenous cultures, and the maintenance of human health . We now show here that even the last of the wild areas and the IF landscapes could be partly devoid of large mammal populations, potentially resulting in downsized mammal assemblages. As mammal complexity and diversity is reduced due to hunting, the net carbon storage of IFs could be compromised. For example, large frugivores are effective seed-dispersal agents, and their contribution to carbon storage in tropical forests is 2-fold: via discarded fruits that contribute to biomass accumulation in soil  and by being the main agent of dispersal of high-carbon large-seeded tree species . Additionally, the loss of large herbivores results in lower herbivory rates in defaunated forests, which allows fast-growing herbivore-sensitive wind-dispersed plants (i.e., lianas) to outcompete slower-growing animal-dispersed trees. As lianas become more abundant, the net aboveground carbon uptake can be substantially reduced , thereby jeopardizing the role of IFs as carbon sinks .
The insets represent the total area of forest (y-axis) predicted to be defaunated (DI > 0.1, red) and intact (DI ≤ 0.1, green) in the respective graphs. Available at https://figshare.com/projects/Intact_but_emtpy_forests_Patterns_of_hunting-induced_mammal_defaunation_in_the_tropics/31118. DI, defaunation index; IF, intact forest; WA, wilderness area.
Finally, we found that mammal populations could be defaunated in 20% of the IUCN PA when all mammal species are considered (S9 Fig), and that this pattern is more conspicuous in the case of large mammals (57% of the total PA). Most of the PAs predicted to be at risk of defaunation are located in Benin, Burundi, Bangladesh, Thailand, and India, with DI >0.3 (S10 Fig). While area-based conservation seems to be effective for conserving forest habitats, the evidence remains inconclusive regarding the effectiveness for maintaining wildlife populations [32–35], especially in the absence of anti-poaching regulations.
Our models conservatively assess the spatial extent of hunting, as potential time-lagged changes in mammal population abundances were not accounted for in the analysis. That is, we provide a snapshot of the spatial extent of defaunation by upscaling local studies performed between 1980 and 2017 to pantropical estimates. Hunting pressure, however, is sustained in time, and hunters will move towards more remote grounds once the larger and medium-sized species have disappeared in the proximity of their villages . Thus, the mammal densities reported by local studies may be far from equilibrium and could have continued to decline years after the study was completed. In this sense, our estimates are relatively optimistic, yet they suggest a more dramatic picture of biodiversity loss than depicted in previous studies [2, 18], with possibly profound ramifications for ecosystem functioning and the livelihoods of wild-meat-dependent communities. Furthermore, we did not explicitly account for cultural beliefs and taboos regarding the consumption of some species in our analyses, as this information is extremely sparse in the literature. However, this is unlikely to have affected our conclusions because cultural beliefs and taboos do not prevent commercial hunting , and, in addition, their importance is declining due to rapid socioeconomic changes (i.e., increasing wealth and human mobility ). Finally, we did not consider that hunting impacts could be exacerbated by the synergies derived from habitat degradation, deforestation, and fragmentation . Because hunter accessibility increases with increasing fragmentation via road development [40, 41], future projections of hunting-induced defaunation should rely on forecasts of HPD and include, when possible, the establishment of new roads that would provide access to the last remaining remote areas in the tropics.
We advocate for the inclusion of hunting-induced defaunation in large-scale biodiversity assessments, in which it has been routinely ignored due to data paucity. Other potential applications of our defaunation maps include species extinction risk assessments, conservation planning , the quantification of the effects of hunting on body size and animal biomass distributions , and progress evaluations to achieve global biodiversity targets . For example, by combining our defaunation maps with spatially explicit mammal density estimates , we could gauge the remaining standing mammalian biomass in tropical forests. Special attention should be given to assess the integrity of those IFs and WAs that are at risk of defaunation, as their status could shift from “intact” to “half-empty” or “empty” forests. Retaining the integrity of intact tropical forests will not be possible if global and national environmental strategies do not address ongoing hunting practices.
The relationship between mammal species abundance and hunting pressure was quantified using data from peer-reviewed and non-peer reviewed literature selected through a systematic literature search (S1 Table). We expanded upon the data set in Benítez-López and colleagues  that included studies that assessed the impact of hunting on wildlife abundance. Specifically, we searched for studies in which species abundance was reported in at least one hunted area and one unhunted control area, and at increasing distance from access points. Studies with potential confounding effects due to other disturbances (e.g., hunted and logged area versus unhunted unlogged area) were discarded (see Benítez-López and colleagues  for details on the search strategy and the study inclusion criteria). For updating the database, we specifically searched for studies performed in countries not included in the original database. In total, we included 163 studies covering 296 mammal species and 3,281 mammal abundance ratios, which corresponds to a 70% increase in the number of records in comparison with the original database of 1,938 ratios. Each study was georeferenced, and predictors were either recorded from the studies or extracted from spatially explicit raster maps (S1 Fig). Changes in abundance due to hunting pressure were expressed as the response ratio (RR) between the abundance of each species in hunted (Xh) and unhunted (Xc) sites within each study (RR = log(Xh/Xc)) . Some ratios were zero for species that were completely extirpated in areas close to hunters’ access points (mean abundance equals zero), precluding log-transformation. Therefore, we converted our response variable into a binary variable (zero and nonzero abundance ratios) and a continuous variable (RRs calculated for abundance ratios >0) and modelled using a two-stage or hurdle model approach (see Modelling section). RRs are negative (RR < 0) or positive (RR > 0) if the abundance estimates are lower or higher, respectively, because of hunting pressure.
As predictors of hunting-induced defaunation, we used proximity to hunters’ access points (km), accessibility of urban markets (travel time to major towns, min), human population density (ind/km2), availability of domestic meat as an alternative food source (kg/km2), prevalence of stunting in children <5 years old (which indicates insufficient growth as a consequence of malnutrition and is an indicator of poverty ), literacy rate (as an indicator of access to qualified jobs), protection status (protected, nonprotected) , and species traits (body mass, kg, and diet—frugivores, herbivores, insectivores, omnivores, and carnivores) (S2 Table, S1 and S2 Figs). All spatially explicit predictors were calculated within the extent of present-day (sub-) tropical forest ecosystems (i.e., “forest zones”) based on the global year 2000 tree canopy cover data set [1, 48].
Distance to hunters’ access points.
The distances to hunters’ access points (settlements, roads) were extracted from each study . When the distance was not explicitly stated (e.g., comparison of close hunted area versus remote control area), we georeferenced the study locations and calculated the distance to the nearest settlement. Most of our studies (88%) referred to settlements as the main access point, and thus we generated a distance to the nearest settlement raster map for our model projections. We downloaded settlement locations for all countries in the tropical forest zone. The settlement data were collated from different sources, including national databases, the Humanitarian Data Exchange (HDE, https://data.humdata.org/), and OpenStreetMap (OSM, https://download.geofabrik.de) (S3 Table). Settlement data were visually inspected to assess country coverage and alignment with satellite imagery. When there were two or more data sets per country, we used the data set with the most extensive coverage. We merged data sets when they had similar coverages but differed at a few settlement points. Usually, the people that engage in hunting activities come from rural (and relatively) remote areas; i.e., decaying hunting pressure is usually related to the distance to small towns and villages, mostly because the livelihoods of people in urban areas do not depend on wild meat acquisition. Therefore, we excluded as access points settlements that referred to urban areas (those overlapping the built-up category in ESA Land Cover Maps for the year 2015 [https://www.esa-landcover-cci.org/] and those that were categorized as cities in OSM). We also filtered out location points that referred to the name of the country, county, region, or island (OSM data set). Subsequently, we calculated a raster map of the distance to the nearest settlement across the whole pantropical forest zone at a resolution of 30” (ca. 1 km; S2 Fig). Caribbean Islands, Oceanic Islands, and Papua New Guinea were discarded from the analyses due to the paucity of data on settlement locations.
Human population density and travel time to major towns.
HPD is an indicator of wild meat demand and hunting pressure. We obtained raster maps of HPD (1-km resolution) for the period 1990–2015 from SEDAC [49, 50]. Per study location, estimates of HPD were extracted that matched the study year. Because urban areas may act as hubs where wild meat is commercialized in urban markets or transported elsewhere, we included travel time to major towns as a proxy for accessibility of urban markets and potential urban demand. Estimated travel times to the nearest town with more than 50,000 inhabitants were used as a proxy for accessibility of urban markets (1-km resolution, see ). Travel times were extracted from the accessibility map for the year 2000  for studies performed before that year. For studies performed between 2000 and 2015, travel times were interpolated from accessibility maps for the years 2000  and 2015 . For more recent studies, the accessibility map from 2015 was used.
Stunting and literacy rate.
As a poverty indicator, we used the prevalence of stunting, which represents insufficient child growth due to persistent dietary deficiencies and/or illness susceptibility. Stunting is considered a better indicator of economic and social deprivation than estimates of per capita income, as it indicates chronic failure to alleviate poverty . To determine stunting, we used a spatial database on the prevalence of stunting in children <5 years old . This database was outdated for some subnational areas, and thus we updated the spatial information with the WHO Global Database on Child Growth and Malnutrition (https://www.who.int/nutgrowthdb/about/en/), Demographic and Health Surveys (DHS), UNICEF MICS, and national surveys (e.g., ENSIN in Colombia), producing a new global map of stunting (S2 Fig). We recorded estimates for different time periods, and then we extracted stunting values per location by matching the values with the year of study (e.g., if the study year was 2007, we used stunting estimates corresponding to or close to 2007).
Education has been shown to correlate with the potential to access the labor market and, thus, alternative livelihoods that are less dependent on wild meat. We used literacy rate per country from the World Bank database (https://data.worldbank.org/) as a proxy for educational attainment.
We included livestock biomass as a proxy of accessibility of domestic meat as an alternative protein source to wild meat. We used global livestock density maps  to estimate the amount of domestic meat available per grid cell. The densities were transformed into livestock biomass per unit area (kg/km2) based on the average weights of cattle, sheep, pigs, and chicken extracted from the literature (S4 Table). The average weights were calculated based on a large number of animals (range: 1,531–4,951 individuals, depending on livestock type) across multiple studies (range: 24–37 studies).
Prior to modelling, we assessed the collinearity among the explanatory variables, which was low overall (S6 Fig). The highest correlation (Pearson’s rho = −0.59) was between stunting and literacy rate, and thus we kept the former, which is more intimately linked to poverty, for further analyses.
We used a hurdle or "two-stage" mixed model to accommodate the distribution of our response variable, which included local extirpations in ca. 14% of the cases (N = 408 out of 3,281) and abundance declines (or increases) compared with control areas in the rest of the data set (N = 2,873). Hurdle models use a binomial distribution to specify the probability of getting a 0 or a positive value, and then fit a zero-truncated probability density function to the nonzero data . Hence, we used a binomial model to predict whether mammal populations were locally extinct or not (i.e., whether the relative abundance values were 0 or >0) and then fitted a Gaussian model through the nonzero RRs. We specified as random effects Country, Study, and Species to account for between-country variation in hunting laws and policies, culture, taboos, and traditions [56, 57] and to control for nonindependence in the data from the same study or species. All continuous variables included as fixed effects (see Predictors section) were standardized before model selection, and predictions were done with the models refit with unstandardized predictors. We fit models using maximum likelihood (ML) for model selection and restricted maximum likelihood (REML) for coefficient estimation . Model selection was conducted for the binomial and continuous model based on the Bayesian Information Criterion (BIC), with models with a ΔBIC ≤2 considered supported and used for inference and for spatial predictions (S5 Table). We assessed the explained variance of the best models using the marginal R2 (fixed effects) and the conditional R2 (fixed and random effects) (see S1 Text). We also assessed variance partitioning among random effects and calculated semi-partial R2 values for each fixed effect included in the best binomial and continuous models (S5 Fig). We evaluated the predictive accuracy of the hurdle model with 5-fold cross-validation with an 80%/20% training/testing set. We split our predictions into three defaunation intensity categories of high (DI > 0.7), moderate (DI = 0.1–0.7), and low (DI ≤ 0.1), which roughly resemble the categories used in the Planetary Boundaries framework (Biodiversity Intactness Index [BII] = 90%, equivalent to DI = 0.1, and 90%–30%, equivalent to DI = 0.1–0.7) . We assessed the accuracy of our model for predicting these categories of defaunation using sensitivity, specificity, and balanced accuracy. We also squared the correlation coefficient between the “Predicted” and the “Observed” data obtained from all repetitions of the 5-fold cross-validation to calculate a pseudo-R2, which was used as the predictive performance metric for the continuous range of observed abundance declines.
Maps of hunting-induced defaunation
We extracted IUCN species ranges for all mammal species with distributions that overlapped the tropical forest area (N = 3,923 mammal species). This area was based on the global “forest zone” following Potapov and colleagues (2017) , who classified each grid cell from Landsat imagery (30-m resolution) with tree canopy cover greater than 20% in the year 2000  as forest (S2A Fig). IUCN species ranges were gridded to match the spatial resolution of the predictors (1 km with Mollweide equal area projection), and, subsequently, we used our model to project the hunting-induced decline in abundance for each species. Our projections were based on the taxonomic identity of the species (captured by the random-effect intercept “Species”), the country where it was located (random-effect intercept “Country”), and its body mass (species vulnerability to hunting pressure), combined with the distribution of context-dependent drivers of hunting pressure (distance to settlements, population density, network of PAs, etc.) within the species range. Our empirical data cover 7.5% of all tropical mammal species and 30% of medium-sized and large-sized mammal species (i.e., those that are generally hunted ) included in the model projections. The number of species included in our database was proportional to the number of mammalian species included in the model projections (S11 Fig), and the extrapolations covered the range of body mass included in our data (range: 0.018–3,940 kg). The best represented taxa were elephants, armadillos, anteaters, and sloths, followed by ungulates, tapirs, primates, and carnivores. We back-transformed our predicted RR into a defaunation intensity index per species (DIs) as the reverse of the exp(RR), i.e., DIs = 1 − exp(RR). Species-specific defaunation maps were then aggregated to create a composite map of hunting-induced defaunation by averaging the DIs values across all species per grid cell (, with S being the number of species in a grid cell). We present our results in the form of defaunation gradients that range from 0 (not defaunated) to 1 (fully defaunated) , and consider areas with average DI > 0.1 as partially defaunated (hereafter, defaunated). We also calculated the proportion of species with DIs > 0.7 to identify hotspots of defaunation caused by hunting. Because hunting is known to be a size-differential pressure [3, 10], we generated defaunation maps for small (<1 kg), medium (1–20 kg), and large (>20 kg) mammal species separately, in addition to an overall defaunation map. Additionally, we quantified defaunation for specific trophic groups (carnivores, herbivores, frugivores, insectivores) that play key roles in ecosystem functioning via seed dispersal, top-down, or bottom-up regulation [25, 28]. Finally, to identify geographic areas outside the range of socioeconomic predictor variables covered by our data, we calculated and mapped the multivariate environmental similarity surface (MESS) . This analysis indicates areas where our DI estimates should be interpreted with caution, as they are based on extrapolation beyond the socioeconomic values used to fit our models (e.g., areas where HPD is higher or lower than the range of values included in the model fitting).
We then estimated the degree to which intact forest landscapes (IFLs) and WAs are defaunated (DI > 0.1) by overlapping our defaunation maps with the IFLs map, as defined by Potapov and colleagues (2017) , and the HF map, in which HF ≤ 2 is considered low and corresponds to WAs . We calculated the total area that is defaunated (DI > 0.1) and intact (DI ≤ 0.1) within the WA and IF for all species and large mammal species. We used similar procedures to assess the risk of defaunation of IUCN PAs (cat. I–IV).
All analyses were conducted in R 3.4.1 . The package “lme4”  was used to run the mixed models, “MuMIn”  was used for model selection and to calculate the marginal and conditional R2 of the models, “data.table”  was used to manipulate the large databases, “caret”  was used to calculate the accuracy metrics, “ModEvA”  was used to calculate the MESS, and “raster”  and “rgdal”  were used for GIS operations. Spatial analyses were conducted in Mollweide equal-area projection at a resolution of 1km using R and ArcGIS .
S1 Fig. Model scheme.
S2 Fig. Spatial distribution of hunting studies and socioeconomic drivers of hunting pressure.
(A) Location of 163 studies (in blue) with 3,281 abundance estimates for mammals in areas under hunting pressure. (B) Distance to the nearest rural settlement (km), (C) livestock biomass (kg/km2), (D) HPD (ind/km2), (E) travel time to major cities, (F) prevalence of stunting among children under five by the lowest available subnational administrative unit, varying years. Based primarily on the WHO Global Database on Child Growth and Malnutrition (https://www.who.int/nutgrowthdb/about/en/). Available at https://figshare.com/projects/Intact_but_emtpy_forests_Patterns_of_hunting-induced_mammal_defaunation_in_the_tropics/31118. HPD, human population density.
S3 Fig. Partial plots of the relationships between the probability of a species/population being locally extirpated (0) or not (1) due to hunting, and the retained predictors in the best model.
(A) Distance to hunters’ access points, (B) HPD, (C) PA status (yes, no), (D) body mass, and (E) prevalence of stunting. CIs (95%) are shown in gray. The scale of the y-axis has been adjusted to enhance visualization of the fitted lines. Available at https://figshare.com/projects/Intact_but_emtpy_forests_Patterns_of_hunting-induced_mammal_defaunation_in_the_tropics/31118. HPD, human population density; PA, protected area.
S4 Fig. Partials plots of the relationships between the RR (change in species abundance) and the retained predictors in the best model.
The dashed gray line indicates that hunting pressure has no effect on species abundance (RR = 0). Positive values indicate an increase in species abundance, whereas negative values indicate a negative effect on species abundance. (A) Distance to hunters’ access points, (B) body mass, (C) interaction between body mass and distance, and (D) HPD. CIs (95%) are shown in gray. In (C), dark blue: 0.1 kg, e.g., Oryzomys spp.; light blue: 1 kg, e.g., Sylvilagus brasiliensis; yellow: 10 kg, e.g., Alouatta spp.; orange: 100 kg, e.g., Panthera onca; red: 4,000 kg, e.g., Loxodonta africana. Available at https://figshare.com/projects/Intact_but_emtpy_forests_Patterns_of_hunting-induced_mammal_defaunation_in_the_tropics/31118. HPD, human population density; RR, response ratio.
S5 Fig. Effects and relative importance of predictors on mammal abundance declines due to hunting pressure.
Standardized coefficient estimates of the variables retained in the best (A) binomial (extinct/no extinct) and (B) Gaussian models (RR). Explained variance by (C) the random effects and the (D) fixed effects of the binomial and Gaussian models. Available at https://figshare.com/projects/Intact_but_emtpy_forests_Patterns_of_hunting-induced_mammal_defaunation_in_the_tropics/31118. BM, body mass; Dist, distance to hunters’ access points; HPD, human population density; PA, protected area; RR, response ratio; Stunt, stunting.
S6 Fig. Collinearity test between explanatory variables and predictive performance of the models.
(A) Correlation plot between explanatory variables, (B) predictive performance metrics (mean ± SD) for three categories of defaunation (low, DI < 0.1; intermediate, DI = 0.1–0.7; high, DI = 0.7–1.0). (C) Predicted versus observed categories of defaunation intensity obtained with the best hurdle model for the cross-validated data set. Size of the squares relative to the size of the grid indicates the proportion of the observed data of a given DI category (columns) to match with the prediction of a particular DI category (rows). Available at https://figshare.com/projects/Intact_but_emtpy_forests_Patterns_of_hunting-induced_mammal_defaunation_in_the_tropics/31118. BM, body mass, DI, defaunation index; Dist, distance to hunters’ access points; HPD, human population density; Literacy, literacy rate; LivestockBio, biomass of domestic livestock; Stunt, stunting; TravTime, travel time to major towns.
S7 Fig. Mean DI per country and 95% CI (black lines) for (A) all pantropical area and (B) after excluding areas outside the socioeconomic domain covered by our data.
Colors denote different regions. Available at https://figshare.com/projects/Intact_but_emtpy_forests_Patterns_of_hunting-induced_mammal_defaunation_in_the_tropics/31118. CAR, Central African Republic; DI, defaunation index; DRC, Democratic Republic of Congo.
(A) Geographic areas inside and outside the socioeconomic domain covered by our data, as estimated by the MESS. The values represent the similarity between each grid cell in pantropical range and those in the reference data set used to fit the models. Values range from positive (green) to negative (red). Positive values represent interpolation areas with similar socioeconomic factors (distance to hunters’ access points, HPD, and prevalence of stunting) than those used to fit the models that are covered by our data set. Negative values indicate localities where at least one socioeconomic variable is outside the range of socioeconomic variables in our data set. (B) Main variable that is dissimilar in each grid cell compared with the socioeconomic domain in our data set. Orange, distance to the nearest rural settlement; green, HPD; blue, prevalence of stunting. Available at https://figshare.com/projects/Intact_but_emtpy_forests_Patterns_of_hunting-induced_mammal_defaunation_in_the_tropics/31118. HPD, human population density; MESS, multivariate environmental similarity surface.
S9 Fig. Geographic variation and spatial extent of hunting-induced defaunation in IUCN PAs (I–IV) for (A) all species and (B) large-bodied species.
Available at https://figshare.com/projects/Intact_but_emtpy_forests_Patterns_of_hunting-induced_mammal_defaunation_in_the_tropics/31118. IUCN, International Union for Conservation of Nature; PA, protected area.
S10 Fig. Mean DI and 95% CI (black lines) within IUCN PAs (I–IV) per country.
Colors denote different regions. Available at https://figshare.com/projects/Intact_but_emtpy_forests_Patterns_of_hunting-induced_mammal_defaunation_in_the_tropics/31118. CAR, Central African Republic; DI, defaunation index; DRC, Democratic Republic of Congo; IUCN, International Union for Conservation of Nature; PA, protected area.
S11 Fig. The relationship between the number of species represented in our database (N = 296) and the number of tropical species for which we extrapolated our models (N = 3,293) for 14 orders.
Lines show 10% (dotted), 50% (solid), and 90% (dashed) representations of the predicted species in our data set. Available at https://figshare.com/projects/Intact_but_emtpy_forests_Patterns_of_hunting-induced_mammal_defaunation_in_the_tropics/31118.
S1 Table. List of data sources included in our analyses and the associated metadata: Author and year, type of source (SP, MT, DT, TR, BC), location, habitat, order, type of access point, type of hunting, legality status, number of studies, and methods used in each source.
BC, book chapter; DT, doctoral thesis; MT, master thesis; SP, scientific publication; TR, technical report.
S2 Table. Overview of explanatory variables included in the hurdle models.
S3 Table. Sources of settlement location data.
S4 Table. Number of animals and studies used to estimate average body weights for cattle, sheep, goats, pigs, and chickens.
S5 Table. Model selection results for (a) the binomial model (0/1, extirpated versus not extirpated) and (b) the continuous model.
Models were ranked according to BIC. We only show models with a BIC weight >0.01. The best model (ΔBIC < 2, in bold) was used in the cross-validation analyses and for spatial predictions. BM, body mass; BIC, Bayesian Information Criterion; Dist, distance to hunters’ access points, PA, protected area; PopDens, human population density; Stunt, stunting; TravTime, travel time to major towns.
This study contributes to the GLOBIO 4.0 project (www.globio.info). We are grateful to Eddy Scheper for calculating the distances to nearest settlement raster maps and to Leonie Lautz for sharing the data on body weights of different livestock. This study would have not been possible without the tremendous effort by a large number of authors who have assessed local-scale hunting impacts on mammal populations during the past four decades.
- 1. Potapov P, Hansen MC, Laestadius L, Turubanova S, Yaroshenko A, Thies C, et al. The last frontiers of wilderness: Tracking loss of intact forest landscapes from 2000 to 2013. Science Advances. 2017;3(1). pmid:28097216
- 2. Venter O, Sanderson EW, Magrach A, Allan JR, Beher J, Jones KR, et al. Sixteen years of change in the global terrestrial human footprint and implications for biodiversity conservation. Nature Communications. 2016;7.
- 3. Benítez-López A, Alkemade R, Schipper AM, Ingram DJ, Verweij PA, Eikelboom JAJ, et al. The impact of hunting on tropical mammal and bird populations. Science. 2017;356(6334):180–3. pmid:28408600
- 4. Bennett EL. Is there a link between wild meat and food security? Conservation Biology. 2002;16(3):590–2.
- 5. Peres CA. Effects of subsistence hunting on vertebrate community structure in Amazonian forests. Conservation Biology. 2000;14(1):240–53.
- 6. Robinson JG, Redford KH. Neotropical Wildlife Use and Conservation. Chicago: Chicago University Press; 1991.
- 7. Wilkie DS, Carpenter JF. Bushmeat hunting in the Congo Basin: an assessment of impacts and options for mitigation. Biodiversity and Conservation. 1999;8(7):927–55.
- 8. Fa JE, Peres CA, Meeuwig J. Bushmeat exploitation in tropical forests: an intercontinental comparison. Conservation Biology. 2002;16(1):232–7.
- 9. Redford KH. The empty forest. BioScience. 1992;42:412–22.
- 10. Ripple WJ, Abernethy K, Betts MG, Chapron G, Dirzo R, Galetti M, et al. Bushmeat hunting and extinction risk to the world’s mammals. Royal Society Open Science. 2016;3(10):160498. pmid:27853564
- 11. Peres CA, Barlow J, Laurance WF. Detecting anthropogenic disturbance in tropical forests. Trends in Ecology & Evolution. 2006;21(5):227–9.
- 12. Peres CA. Effects of hunting on western Amazonian primate communities. Biological Conservation. 1990;54(1):47–59. https://dx.doi.org/10.1016/0006-3207(90)90041-M.
- 13. Peres CA, Nascimento HS. Impact of game hunting by the Kayapo of south-eastern Amazonia: implications for wildlife conservation in tropical forest indigenous reserves. Biodiversity and Conservation. 2006;15(8):2627–53.
- 14. Koerner SE, Poulsen JR, Blanchard EJ, Okouyi J, Clark CJ. Vertebrate community composition and diversity declines along a defaunation gradient radiating from rural villages in Gabon. Journal of Applied Ecology. 2016.
- 15. Kuehl HS, Nzeingui C, Yeno SLD, Huijbregts B, Boesch C, Walsh PD. Discriminating between village and commercial hunting of apes. Biological Conservation. 2009;142(7):1500–6.
- 16. Milner-Gulland EJ, Bennett EL. Wild meat: the bigger picture. Trends in Ecology & Evolution. 2003;18(7):351–7. https://dx.doi.org/10.1016/S0169-5347(03)00123-X.
- 17. Allan JR, Venter O, Watson JE. Temporally inter-comparable maps of terrestrial wilderness and the Last of the Wild. Scientific Data. 2017;4:170187. pmid:29231923
- 18. Watson JE, Evans T, Venter O, Williams B, Tulloch A, Stewart C, et al. The exceptional value of intact forest ecosystems. Nature ecology & evolution. 2018:1.
- 19. Abernethy K, Coad L, Taylor G, Lee M, Maisels F. Extent and ecological consequences of hunting in Central African rainforests in the twenty-first century. Philosophical Transactions of the Royal Society of London B: Biological Sciences. 2013;368(1625):20120303. pmid:23878333
- 20. Nasi R, Taber A, Van Vliet N. Empty forests, empty stomachs? Bushmeat and livelihoods in the Congo and Amazon Basins. International Forestry Review. 2011;13(3):355–68.
- 21. Fa JE, Olivero J, Farfán MA, Lewis J, Yasuoka H, Noss A, et al. Differences between Pygmy and Non-Pygmy Hunting in Congo Basin Forests. PLoS ONE. 2016;11(9):e0161703. pmid:27589384
- 22. Ziegler S, Fa JE, Wohlfart C, Streit B, Jacob S, Wegmann M. Mapping Bushmeat Hunting Pressure in Central Africa. Biotropica. 2016:n/a-n/a.
- 23. Peres CA, Emilio T, Schietti J, Desmoulière SJ, Levi T. Dispersal limitation induces long-term biomass collapse in overhunted Amazonian forests. Proceedings of the National Academy of Sciences. 2016;113(4):892–7.
- 24. Camargo-Sanabria AA, Mendoza E, Guevara R, Martínez-Ramos M, Dirzo R. Experimental defaunation of terrestrial mammalian herbivores alters tropical rainforest understorey diversity. Proc R Soc B. 2015;282(1800):20142580. pmid:25540281
- 25. Ripple WJ, Estes JA, Beschta RL, Wilmers CC, Ritchie EG, Hebblewhite M, et al. Status and ecological effects of the world’s largest carnivores. Science. 2014;343(6167):1241484 pmid:24408439
- 26. Estes JA, Terborgh J, Brashares JS, Power ME, Berger J, Bond WJ, et al. Trophic downgrading of planet Earth. Science. 2011;333(6040):301–6. pmid:21764740
- 27. Terborgh JW. Toward a trophic theory of species diversity. Proceedings of the National Academy of Sciences. 2015;112(37):11415–22.
- 28. Ripple WJ, Newsome TM, Wolf C, Dirzo R, Everatt KT, Galetti M, et al. Collapse of the world’s largest herbivores. Science Advances. 2015;1(4):e1400103. pmid:26601172
- 29. Sobral M, Silvius KM, Overman H, Oliveira LF, Raab TK, Fragoso JM. Mammal diversity influences the carbon cycle through trophic interactions in the Amazon. Nature ecology & evolution. 2017;1(11):1670.
- 30. Van Der Heijden GM, Powers JS, Schnitzer SA. Lianas reduce carbon accumulation and storage in tropical forests. Proceedings of the National Academy of Sciences. 2015;112(43):13267–71.
- 31. Lewis SL, Lopez-Gonzalez G, Sonké B, Affum-Baffoe K, Baker TR, Ojo LO, et al. Increasing carbon storage in intact African tropical forests. Nature. 2009;457(7232):1003. pmid:19225523
- 32. Craigie ID, Baillie JE, Balmford A, Carbone C, Collen B, Green RE, et al. Large mammal population declines in Africa’s protected areas. Biological Conservation. 2010;143(9):2221–8.
- 33. Geldmann J, Barnes M, Coad L, Craigie ID, Hockings M, Burgess ND. Effectiveness of terrestrial protected areas in reducing habitat loss and population declines. Biological Conservation. 2013;161:230–8.
- 34. Gray CL, Hill SL, Newbold T, Hudson LN, Börger L, Contu S, et al. Local biodiversity is higher inside than outside terrestrial protected areas worldwide. Nature Communications. 2016;7:12306. pmid:27465407
- 35. Coetzee BW, Gaston KJ, Chown SL. Local scale comparisons of biodiversity as a test for global protected area ecological performance: a meta-analysis. PLoS ONE. 2014;9(8):e105824. pmid:25162620
- 36. Coad L, Schleicher J, Milner-Gulland EJ, Marthews TR, Starkey M, Manica A, et al. Social and Ecological Change over a Decade in a Village Hunting System, Central Gabon. Conservation Biology. 2013;27(2):270–80. pmid:23369059
- 37. Kideghesho JR. Co-existence between the traditional societies and wildlife in western Serengeti, Tanzania: its relevancy in contemporary wildlife conservation efforts. Biodiversity and Conservation. 2008;17(8):1861–81.
- 38. Jenkins RK, Keane A, Rakotoarivelo AR, Rakotomboavonjy V, Randrianandrianina FH, Razafimanahaka HJ, et al. Analysis of patterns of bushmeat consumption reveals extensive exploitation of protected species in eastern Madagascar. PLoS ONE. 2011;6(12):e27570. pmid:22194787
- 39. Brook BW, Sodhi NS, Bradshaw CJ. Synergies among extinction drivers under global change. Trends in ecology & evolution. 2008;23(8):453–60.
- 40. Peres CA. Synergistic effects of subsistence hunting and habitat fragmentation on Amazonian forest vertebrates. Conservation biology. 2001;15(6):1490–505.
- 41. Laurance WF, Peletier-Jellema A, Geenen B, Koster H, Verweij P, Van Dijck P, et al. Reducing the global environmental impacts of rapid infrastructure expansion. Current Biology. 2015;25(7):R259–R62. pmid:25754645
- 42. Tulloch VJ, Tulloch AI, Visconti P, Halpern BS, Watson JE, Evans MC, et al. Why do we map threats? Linking threat mapping with actions to make better conservation decisions. Frontiers in Ecology and the Environment. 2015;13(2):91–9.
- 43. Santini L, González-Suárez M, Rondinini C, Di Marco M. Shifting baseline in macroecology? Unravelling the influence of human impact on mammalian body mass. Diversity and Distributions. 2017;23(6):640–9.
- 44. Tittensor DP, Walpole M, Hill SL, Boyce DG, Britten GL, Burgess ND, et al. A mid-term analysis of progress toward international biodiversity targets. Science. 2014;346(6206):241–4. pmid:25278504
- 45. Santini L, Isaac NJ, Maiorano L, Ficetola GF, Huijbregts MA, Carbone C, et al. Global drivers of population density in terrestrial vertebrates. Global Ecology and Biogeography. 2018;27(8):968–79.
- 46. FAO. Chronic undernutrition among children: an indicator of poverty. Rome: FAO/EDRN & ESNA; 2003.
- 47. The World Database on Protected Areas (WDPA) [Internet]. 2016 [cited 2016 Feb 10]. www.protectedplanet.net.
- 48. Hansen MC, Potapov PV, Moore R, Hancher M, Turubanova SA, Tyukavina A, et al. High-Resolution Global Maps of 21st-Century Forest Cover Change. Science. 2013;342(6160):850–3. pmid:24233722
- 49. Center for International Earth Science Information Network—CIESIN—Columbia University, International Food Policy Research Institute—IFPRI, The World Bank, Centro Internacional de Agricultura Tropical—CIAT. Global Rural-Urban Mapping Project, Version 1 (GRUMPv1): Population Density Grid. Palisades, NY: NASA Socioeconomic Data and Applications Center (SEDAC); 2011.
- 50. Center for International Earth Science Information Network—CIESIN—Columbia University. Gridded Population of the World, Version 4 (GPWv4): Population Density. Palisades, NY: NASA Socioeconomic Data and Applications Center (SEDAC); 2016.
- 51. Nelson A. Estimated travel time to the nearest city of 50,000 or more people in year 2000. Ispra, Italy: Global Environment Monitoring Unit—Joint Research Centre of the European Commission.; 2008.
- 52. Weiss D, Nelson A, Gibson H, Temperley W, Peedell S, Lieber A, et al. A global map of travel time to cities to assess inequalities in accessibility in 2015. Nature. 2018.
- 53. UNICEF. Progress for Children, No.2 –A report card on gender parity and primary education. New York: 2005 April 2005. Report No.
- 54. Robinson TP, Wint GRW, Conchedda G, Van Boeckel TP, Ercoli V, Palamara E, et al. Mapping the Global Distribution of Livestock. PLoS ONE. 2014;9(5):e96084. pmid:24875496
- 55. Zuur AF, Ieno EN, Walker NJ, Saveliev AA, Smith GM. Mixed Effects Models and Extensions in Ecology with R. New York: Springer; 2009.
- 56. Ngoufo R, Yongyeh N, Obioha E, Bobo K, Jimoh S, Waltert M. Social norms and cultural services-community belief system and use of wildlife products in the Northern periphery of the Korup National Park, South-West Cameroon. Change and Adaptation in Socio-Ecological Systems. 2014;1(1).
- 57. Bobo KS, Aghomo FFM, Ntumwel BC. Wildlife use and the role of taboos in the conservation of wildlife around the Nkwende Hills Forest Reserve; South-west Cameroon. Journal of ethnobiology and ethnomedicine. 2015;11(1):2.
- 58. Steffen W, Richardson K, Rockström J, Cornell SE, Fetzer I, Bennett EM, et al. Planetary boundaries: Guiding human development on a changing planet. Science. 2015;347(6223):1259855. pmid:25592418
- 59. Dirzo R, Mendoza E, Ortíz P. Size-related differential seed predation in a heavily defaunated neotropical rain forest. Biotropica. 2007;39(3):355–62.
- 60. Elith J, Kearney M, Phillips S. The art of modelling range-shifting species. Methods in ecology and evolution. 2010;1(4):330–42.
- 61. R-Core-Team. R: A language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing; 2017.
- 62. Bates D, Maechler M, Bolker B, Walker S. lme4: Linear mixed-effects models using Eigen and S4. R package version. 2014;1(7):1–23.
- 63. Bartoń K. MuMIn: multi-model inference. R package version 1.10. 0. See https://CRAN.R-project.org/package=MuMIn. 2014.
- 64. Dowle M, Srinivasan A. data.table: Extension of `data.frame. R package version 1.10.4–3. https://CRAN.R-project.org/package=data.table. 2017.
- 65. Kuhn M, Wing J, Weston S, Williams A, Keefer C, Engelhardt A, et al. caret: Classification and regression training. R package version 6.0–21. CRAN, Vienna, Austria. 2015.
- 66. Barbosa A, Brown J, Jimenez-Valverde A, Real R. modEvA: Model evaluation and analysis. R package version 1.3. 2. 2016.
- 67. Hijmans RJ, Van Etten J. raster: Geographic data analysis and modeling. R package version 2.2–31. Google Scholar. 2014.
- 68. Keitt TH. rgdal: Bindings for the Geospatial Data Abstraction Library, R package version 0.6–28. https://cran.r-project.org/package=rgdal. 2010.
- 69. ESRI. ArcGIS Desktop: Release 10.3.1. Redlands, CA: Environmental Systems Research Institute; 2015.